摘要
岩石的密度、磁性与岩石类型、区域构造有着密切联系,研究岩石的密度和磁性分布对位场异常的解释、研究大地构造分区及矿产资源调查具有重要的意义。如何根据位场异常计算地球内部岩石的密度、磁化率或磁化强度分布,称之为位场异常的物性反演。但由于位场的叠加和等效性等原因,物性反演问题一直是国内外研究的热点和难点问题。
目前位场物性反演的方法大致可分为最优化选择法和非线性反演方法两大类,具体的算法有广义最小二乘法(Gauss法)、最速下降法、阻尼最小二乘法(Marguardt法)、奇异值分解法、脊回法、人工神经网络法BP算法、遗传算法、模拟退火算法和基于Radon变换、多尺度边缘的算法等。这些算法的基本思路是,将地下场源区域一次性划分成若干形态已知、密度或磁性未知的小棱柱体单元,首先给定密度或磁性等物性参量一个初始值,然后通过各种反演算法和各种约束条件求得参量的改变值,最终求得合理的物性参量,使得模型的理论正演异常与实测异常之差的L1或L2范数等最小。
目前已有的这些算法,在理论上是合理、可行的,但在实际应用中却存在一些问题。其中,最优化选择法的主要问题,一是需要建立大型的线性代数方程组,计算时间很长;二是方程组病态严重,解的稳定性很差。非线性反演方法的主要问题是,反演数据量较大时,存在超常规的计算量、存储量及收敛速度慢的问题。总体来说,由于目前计算机的内存和速度的限制,在对位场资料的物性反演,尤其是大面积资料的三维物性反演时,已有的反演方法都很难奏效。
针对现有的问题,本文在系统研究了位场分离方法、位场向下延拓方法和频率域层源位场异常快速反演方法的基础上,提出了一种基于位场分离与延拓的三维物性快速反演方法。方法的基本思路为,首先运用位场分离的插值切割法对平面上的场进行不同深度层源的分离,得到各深度层源在地面引起的异常;然后应用大深度的位场迭代法向下延拓技术将各深度层源的地面异常延拓到各深度层的顶部;最后,根据各深度层顶部的异常利用基于棱柱体组合的频率域反演方法反演出各个深度层源的密度或磁性。新反演方法不需要对场源进行一次性剖分和解线性代数方程组,避开了制约三维反演实用化的超大内存需求、解稳定性差和收敛慢引起的特别冗长的计算时间的瓶颈问题。
通过新方法对新疆色尔特能地区、普光气田和东海及邻区重力异常的视密度反演和对霍邱铁矿磁异常的视磁性反演,表明该方法具有精度高、稳定性强和速度快的优点,适用于大面积位场资料的三维物性反演,具有一定的理论和实际意义。运用Fortran语言、c语言和C++Builder编译环境,结合动态链接库技术和多线程技术对新方法编制了Windows系统下的可视化软件,方便了方法的应用和推广。
The density and magnetization of rock are closely related to the lithology and tectonics. Therefore, the study on the density and magnetization of rock is very significant for the interpretation of potential anomaly data, tectonic unit division and mineral resources investigation. The calculating of rock density and magnitization according to potential anomaly is called physical properties inversion of potential anomaly. Because of potential superposition and equivalence principle, the research of physical properties inversion has been an interested and difficult problem at home and abroad.
At present, there are two main methods to inverse the physical properties of potential anomaly, such as optimization method and non-linear method. The algorithms of these two methods include Gauss algorithm, the steepest descent algorithm, Marquardt algorithm, singular value decomposition algorithm, ridge regression algorithm, the Back Propagation algorithm of artifical neural networks, genetic algorithm method and simulated annealing algorithm et al. The basic thinking of present methods is to divide the subsurround field sources region into many rectangular blocks with known shape and unknown density or magnetization. At first, initial values of density or magnetization are given. Then, the changed values of density or magnetization are needed to be gained. Finaly, the reasonable physical properties values that make the L1 norm or L2 norm of the difference between theory forward anomaly data and observed data to be least are obtained.
Even though the present methods are reasonable and feasible in theory, there are still many problems in practice. There are mainly two problems about optimization method. Firstly, the establishment of large linear algebraic equation is needed, it consumes large computer memory and long computation time. Secondly, the equation is ill-conditioned and solution is unstable. The main problem of non-linear inversion methods is the long calculating time cosuming and large memory occupying when the quantities of inverse data are large. Overall the present inversion methods are ineffective during inversion of potential data, especially 3-D physical properties inversion in the large area due to the restriction of memory and computing velocity of computer.
This paper presents a rapid method for 3D physical properties inversion based on separation and continuation of potential field. The new method is carried out in the following three steps. Firstly, potential feilds corresponding to strata of different depth on horizontal groud are gained through separating the potential fields on the plane by a cutting method. Secondly, the fields corresponding to strata of different depth on horizontal groud are reduced to a plane corresponding to the top surface of different depth strata by using downward continuation. Finally, according to the field on the top surface, inversion is conducted to get structure of apparent density or magnetization in frequency. The inversion technique characterized by faster computing speed, it does not need to divide the subsurround field sources region into many rectangular blocks in one time and does not need to slove linear algebraic equations as well. So, this technique can overcome the requirement of large computer memory, unstable solution and computation time bottleneck which hinders the application of 3Dinversion to practice.
Several cases, including apparent density inversion of gravity anomolies in Xinjiang, Puguang gas field and the East China Sea and its adjacent area as well apparent magnetization inversion of magnetic anomaly in Huoqiu iron mine area, are presented to demonstrate the effect of applying this method. This method is suitable for 3D inversion of potential field in large area and has some theorical and practical significance. A visual software in Windows system is designed with Fortran language, C language, DLL technology and multithread technology based C++ Builder 6.0, which provides convenience for the use and popularization of new method.
引文
[1]程方道,刘东甲,姚汝信.划分重力区域场与局部场的研究[J].物探化探计算技术,1987,9(1):3-11.
[2]程明.可视化编程中界面的设计[J].科技论坛,2005,23(1):97-98.
[3]陈建国,夏庆霖.利用小波分析提取深层次物化探异常信息[J].地球科学-中国地质大学学报,1999,24(5):509-512.
[4]陈生昌,肖鹏飞.位场向下延拓的波数域广义逆算法[J].地球物理学报,2007,50(6):1816-1822.
[5]段本春,徐世浙,阎汉杰,等.划分磁异常的插值切割法在研究火成岩体分布中的应用[J].石油地球物理勘探,1998,34(1):125-131.
[6]范朝阳,张良驹.多线程程序设计的概念与应用[J].小型微型计算机系统,1996,17(4):1-6.
[7]高德章,候遵泽,唐建.东海及邻区重力异常多尺度分解[J].地球物理学报,2000,43(6):842-849.
[8]高德章,唐建,薄玉玲.东海地球物理综合探测剖面及其解释[J].中国海上油气(地质),2003,17(1):38-43.
[9]高德章,赵金海,薄玉玲,等.东海重磁地震综合探测剖面研究[J].地球物理学报,2004,47(5):853-861.
[10]关小平,黄嘉正.重力、地震资料联合反演实例[J].石油地球物理勘探,1995,30(3):379-385.
[11]管志宁,张昌达,陈方道,等.磁法勘探重要问题理论分析与应用[M].北京:地质出版社,1993,42-56.
[12]管志宁,郝天姚,姚长利.21世纪重力与磁法勘探的展望[J].地球物理学进展,2002,17(2):237-224.
[13]郭积忠.重力测量区域校正的新方法-应用网函数插值法计算重力区域场[J].物探与化探,1981,5(03):137-152.
[14]郭玉贵,李延成,许东禹,等.黄东海大陆架及邻域大地构造演化史[J].海洋地质与第四纪地质,1997,17(1):1-12.
[15]郝天珧,刘伊克,段昶.中国东部及其邻域地球物理场特征与大地构造意义 [J].地球物理学报,1997,40(5):677-690.
[16]郝天珧,徐亚,胥颐,等.对黄海-东海研究区深部结构的一些新的认识[J].地球物理学报,2006,49(2):458-468.
[17]侯重初.位场频域向下延拓方法[J].物探与化探,1982,6(1):33-40.
[18]侯重初,蔡宗熹.剖面曲线上的位场转换系统[J].地球物理学报,1984,27(6):573-581.
[19]侯遵泽,杨文采.中国重力异常的小波变换与多尺度分析[J].地球物理学报,1997,40(1):85-95.
[20]黄卫祖.重力资料的三维密度反演[J].长春地质学院学报,1993,23(1):75-79.
[21]江为为,宋海斌,郝天珧,等.东海陆架盆地及其周边海域地质地球物理场特征[J].地球物理学进展,2001,16(2):18-27.
[22]梁锦文.分离重磁区域场与局部场的维纳滤波器[J].地球物理学报,1983,26(1),80-88.
[23]梁锦文.位场向下延拓的正则化方法[J].地球物理学报,1989,32(5):600-608.
[24]李九亭.划分重力区域异常与局部异常的变阶滑动趋势分析法[J].物探化探计算技术,1998,20(1):53-61.
[25]李盛汉.油气重力异常的特征与分离[J].石油物探,1997,36(4):70-76.
[26]刘军.重力异常频率域分离方法及应用[J].物探与化探,1998,22(6):446-452.
[27]刘光鼎.东海的地质与油气勘探[J].地球物理学报,1988,31(2):184-197.
[28]刘光鼎.中国海区及邻域地质地球物理系列图[M].北京:科学出版社,1992,5-6.
[29]刘青松,王宝仁.利用多次匹配滤波技术进行垂向位场分离[J].物探化探计算技术,1996,18(4):279-286.
[30]刘祥重.区域重力波数域解释及找油效果[J].石油地球物理勘探,1988,23(1):110-118.
[31]刘卫国,蔡旭辉.FORTRAN 90程序设计教程[M].北京:北京邮电大学出版社,2003,1-3.
[32]李伟生.深入C++Builder编程[M].西安:西安电子科技大学出版社,2001,2-3.
[33]李宗杰,杨林,王勤聪.二维小波变换在位场数据处理中的应用试验研究[J].石油物探,1997,36(3):121-125.
[34]栾锡武,高德章,喻普之,等.中国东海及邻近海域一条剖面的地壳速度结构研究[J].地球物理学进展,2001,16(2):28-34.
[35]马丽,张露.Windows动态链接库技术面面观[J].科技信息,2007,15(1):59-60.
[36]毛小平,吴蓉元.频率域位场下延的振荡机制及消除方法[J].石油地球物理勘探,1998,33(2):230-237.
[37]马永生.四川盆地普光大气田的发现与勘探[J].海相油气地质,2006,11(2):35-40.
[38]马永生.四川盆地普光超大型气田的形成机制[J].石油学报,2007,28(2):9-14.
[39]孟小红,王霞.利用Radon变换进行三度体重磁异常反演[J].地球科学,1995,20(5),594-598.
[40]穆石敏,申宁化,孙运生.区域地球物理数据处理方法及其应用[M].长春:吉林科学技术出版社,1990,79-86.
[41]宁津生,汪海洪,罗志才.基于多尺度边缘约束的重力场信号的向下延拓[J].地球物理,2005,48(1):63-68.
[42]潘作枢.位场资料处理及解释问题[M].北京:地质出版社,1992,1-2.
[43]宋景明.关于重力场分离的思考[C].杭州:重磁数据处理解释应用研讨会论文集,2008,128-129.
[44]汤井田,宋守根,何继善.多分辨分析和重磁异常的识别与分层次提取[J].中国有色金属学报,1994,4(3):6-14.
[45]童永春.一种提取重力局部异常的方法-九点差分格式异步迭代法[J].地质与勘探,1984,12(1):51-58.
[46]王邦华,王理.重磁位场的正则化向下延拓[J].物探化探计算技术,1998,20(1):30-35.
[47]汪炳柱,徐世浙,刘保华,等.多次插值切割法分场的一个实例[J].石油地球物理勘探,1997,32(3):431-438.
[48]汪炳柱,王朔儒.二维位场向上延拓和向下延拓的样条函数法[J].物探化探 计算技术,1998,20(2):126-129.
[49]王家林.对我国石油重磁勘探发展的几点思考[J].勘探地球物理进展,2006,29(2):82-86.
[50]王敏,邵定宏.动态链接库技术及其应用实例[J].软件时空,2006,22(9):272-274.
[51]王四龙,许孝庭.用三维趋势面分析分离位场异常[J].煤田地质与勘探,1993,21(5):60-64.
[52]王兴涛,夏哲仁,石磐,等.航空重力测量数据向下延拓方法比较[J].地球物理学报,2004,47(6):1017-1022.
[53]文百红,程方道.用于划分磁异常的新方法-插值切割法[J].中南矿冶学院学报,1990,2l(3):229-235.
[54]吴健生,陈冰,王家林.东海陆架区中北部前第三系基底综合地球物理研究[J].热带海洋学报,2005,24(2):8-15.
[55]徐宝慈,李春华.位场数据处理理论与问题[M].长春:吉林大学出版社,1995,1-2.
[56]徐菊生,刘锁旺,朱仲芬,等.中国东海和邻区重力测量结果及其构造意义[J].地震地质,1986,8(1):43-52.
[57]胥颐,刘建华,郝天珧,等.中国东部海域及邻区岩石层地幔的P波速度结构与构造分析[J].地球物理学报,2006,29(4):1053-1061.
[58]徐世浙.二度重磁力场解析延拓的实际方法[J].地球物理学报,1965,14(2):136-146.
[59]徐世浙,张研,文百红,等.切割法在陆东地区磁异常解释中的应用[J].石油物探,2006,45(03):316-318.
[60]徐世浙.位场延拓的积分-迭代法[J].地球物理学报,2006,49(4):1176-1182.
[61]徐世浙.位场大深度向下延拓[J].物探化探计算技术,2006,28.(增刊):29-31.
[62]徐世浙,曹洛华,姚敬金.重力异常三维反演-视密度成像方法技术的应用[J].物探与化探,2007,31(1):25-28.
[63]徐世浙.迭代法与FFT法位场向下延拓效果的比较[J].地球物理学报,2007,50(1):285-289.
[64]徐世浙,余海龙.位场曲化平的插值-迭代法[J].地球物理学报,2007,50(6):1811-1815.
[65]徐世浙,余海龙,姚敬金,等.新疆色尔特能地区视密度和视磁性反演[J].物探化探计算技术,2007,29(增刊):12-16.
[66]羊春华.筛选.趋势分析法分离区域异常与局部异常[J].物探与化探计算技术,2005,29(2):167-169.
[67]杨金玉,徐世浙,余海龙,等.视密度反演在东海及邻区重力异常解释中的应用[J].地球物理学报,2008,5l(6):1909-1916.
[68]杨文采,陈军中.地球物理反演的理论和方法[M].北京:地质出版社,1997,2-3.
[69]姚长利,管志宁.位场转换的抽样分组法[J].石油地球物理勘探,1997,32(5):696-702.
[70]姚长利,郑元满.重磁遗传算法三维反演中动态数组优化方法[J].物探化探计算技术,2002,24(3):240-245.
[71]姚长利,郝天姚.重磁遗传算法三维反演中高速计算及有效存储方法技术[J].地球物理学报,2003,46(3):252-258.
[72]于德武.用等效源法进行磁异常转换[J].物探化探计算技术,2004,26(2):133-135.
[73]于舟,李昭荣.重磁异常的方向维纳串联滤波器[J].物探化探计算技术,1989,11(2):149-155.
[74]曾华霖.重磁勘探反演问题研究评述[J].物探与化探,1990,14(3):182-190.
[75]曾华霖,阚筱玲.重磁勘探反演问题[M].北京:石油工业出版社,1991,12-13.
[76]张小路.重磁反演中的控制随机搜索法[J].物探化探计算技术,1991,13(4):275-281.
[77]朱自强,程方道.重磁异常反演的直接解法-随机搜索法[J].地质与勘探,1992,28(9):34-71.
[78]朱自强,黄国祥.重磁反演的逆人工神经网络BP算法及其应用[J].中南矿冶学院学报,1994,25(3):288-293.
[79]Abdelrahman E M,Bayoumi A I,El-Araby H M.A least-squares minimization approach to invert gravity data[J]. Geophysics, 1991, 56(1): 115-118.
[80] AgocsWB. Least-squares residual anomaly determination[J]. Geophysics, 1951, 16 (4): 686-696.
[81] Bott M H P. Solution of the linear inverse problem in magnetic interpretation with application to occanic magnetic anomalies[J]. Geophys, J. Roy. Astr. Soc. , 1967, 13 (2): 313-323.
[82] Cooper G. The stable downward continuation of potential field data[J]. Exploration Geophysics, 2004, 35 (2): 260-265.
[83] Dampancy C N G. The equivalent source technique[J]. Geophysics, 1969, 34 (1): 39-53.
[84] Emilia D A. Equivalent sources used as an analytic base for processing total magnetic field profiles[J]. Geophysics, 1973, 38(3): 339-348.
[85] Fajklewicz, Zbigniew. The use of cracovian coputation in estimating the regional gravity[J]. Geophysics, 1959, 24(3): 465-478.
[86] Fedi M, Quara T. Wavelet analysis for the regional-residual and local separation of potential field anomalies[J]. Geophysical Prospecting, 1998, 46(5): 507-525.
[87] Fedi M, Florio G. A stable downward continuation by using the ISVD method [J]. Geophys. J. Int. , 2002, 151 (1): 146-156.
[88] Grant F S. Review of data processing and interpretation methods in gravity and magnetics[J]. Geophysics, 1972, 37(4): 647-661.
[89] Jacobsen B H. A case for upward continuation as a standard separation filter for potential-field maps[J]. Geophysics, 1987, 52(8), 1138-1148.
[90] Lee CS, Shor G G, Bibee L D, et al. Okinawa Trough: Origin of a backarc basin[J]. Marine Geology, 1980, 35 (2): 219-241.
[91] LeydenR, EwingM, Murauchi S. Sonobuoy refraction measurement in the East China Sea[J]. Am. Assoc. Petrol Geol Bull, 1973, 57(12): 2396-2403.
[92] Ludwig W J, Muuauchi S, Den N, et al. Structure of East China Sea-western Philippine Sea margin off southern Kyushu. Japan[J]. Journal of Geophysical Research, 1968, 78(14): 2526-2536.
[93] Meulebrouck J V. Density and magnetic tomography of the upper continental crust: application to a thrust area of the French Masif Central[J]. Annales Geophysicae, 1984, 2(5): 579-592.
[94] Mottl J, Mottlova L. Solution of the inverse gravimetric problem with the aid of integer linear Programming[J]. Geoexploration, 1972, 10(1): 53-62.
[95] Nabighian M N, Grauch V J S, Hansen R O, et al. The historical development of the magentic method in exploration[J]. Geophysics, 2005, 70(6): 33-61.
[96] Nabighian M N, Ander M E, Grauch V J S, etal. The historical development of the gravity method in exploration[J]. Geophysics, 2005, 70(6): 63-89.
[97] Pawlowski R S, Hansen R O. Gravity anomaly separation by Wiener filtering [J]. Geophysics, 1990, 55(5): 539-548.
[98] Pawlowski R S . Preferential continuation for potential-field anomaly enhancement[J]. Geophysics, 1995, 60(2): 390-398.
[99] Ridsdill-Smith T A. Separation filtering of aeromagnetic data using filter-banks[J]. Exploration Geophysics, 1998, 29(5): 577-583.
[100] Sandwell D T, Smith W H F. Marine gravity anomaly from Geosat and ERS 1 satellite altimetry[J]. Journal of Geophysical Research, 1997, 102 (B5): 10039- 10054.
[101]Skeels D C. What is residual gravity? [J]. Geophysics, 1967, 32(5): 872-876.
[102] Spector A, Grant F S. Statistical models for interpreting aeromagnetic data[J]. Geophysics, 1970, 35 (2): 293-302.
[103] Syberg F J R. A Fourier method for the regional-residual problem of potential fields[J]. Geophysicsl Prospecting, 1972, 20(1): 47-75.
[104] Syberg F J R. Potential field continuation between general surfaces[J]. Geophysical Prospecting, 1972, 20(2): 267-282.
[105]Ulrych T J. Effect of wavelength filtering on the shape of the residual anomaly[J]. Geophysics, 1968, 33 (6): 1015-1018.
[106] Wessel P. An empirical method for optimal robust regional-residual separation of geophysical data[J]. Mathematical Geology, 1998, 30(4): 391-407.