重力高阶方向导数的可视化拟合及其应用

第 29卷 � 第 3期 物探化探计算技术 � � 2007年 5月 � � � � 基金项目: 国家 863项目 ( 2002AA132012) 收稿日期: 2006- 02- 28 文章编号: 1001� 1749( 2007) 03� 0205� 04 重力高阶方向导数的可视化拟合及其应用 田 � 舍 1, 宋宇辰 2 ( 1.北京市地质勘察技术院, 北京 � 102218; 2.中国地质大学, 北京 � 100083 ) 摘 � 要: 针对目前人们所使用的交互式重力剖面拟合软件多以正演布格重力异常为主, 在应用 效果方面存有欠缺的情况, 提出了在作重力剖面布格异常拟合的同时, 增加高阶方向导数的拟 合。其目的是发挥布格重力异常的高阶方向导数在突出浅部地质体的特征,压制区域性深部地 质因素,区分相邻地质体引起的叠加异常等方面的优势。这里采用最小二乘法演算了布格重力 异常高阶方向导数的计算过程, 对计算程序编制和使用中涉及的相关问题进行了阐述,并例举了 在地球物理勘探中,利用高阶方向导数提高重力解释效果的实例。并在最后对重力剖面进 行可视化拟合解释时应注意的要点做了分析, 强调了兼顾布格重力异常及其高阶方向导数二者 拟合曲线均基本相似,而不应追求其中某一方面精确拟合的解释原则。 关键词: 重力资料; 高阶导数; 可视化解释 中图分类号: P 631�1� � 文献标识码: A 0� 前言 随着计算机技术的不断提高,地球物理勘探资 料的分析解释逐渐向可视化方向发展,人们愈来愈 多地感受到在计算机可视化技术的支持下,利用已 知的地质特征正演拟合物探曲线所带来的便利。 目前人们所使用的重力剖面可视化实时拟合软件 多以正演布格重力异常 �g为主 [ 1] ,由于实测布格 重力异常在反映地质体细节特征方面所现出的 比较 �迟钝�,且存在区域场校正不精确的因素,实 际应用效果有所欠缺。布格重力异常的各阶导数 在不同形状地质体上有着不同的特征,它可以突出 浅而小的地质体特征, 而压制区域性深部地质因 素;也可将几个互相靠近、埋深相差不大的相邻地 质体引起的叠加异常划分开来, 其导数的阶次越 高,这种分辨能力亦越强。在实际工作中, 通常采 用差商的办法来近似地代替求导,若想提高计算精 度,可用最小二乘法 [ 2, 3]对实测重力布格值进行函 数拟合,然后对 x求导, 即可得到各阶方向导数。 作者通过在可视化条件下,拟合布格重力异常 的同时,增加重力剖面高阶水平方向导数的拟合, 提高重力场在反映地质体细节特征方面的 �灵敏 度�, 以改善可视化技术条件下交互式重力剖面拟 合软件的功能, 提高地球物理勘探的地质效果。 1� 方法原理 假设在实测重力剖面上, 由 n个测点组成的区 间中,布格重力异常 �g可以展成四阶多项式,即 �g (x )� f (x ) = a0 + a1x+ a2x2 + a3x3 + a4x4 ( 1) 则 �g ( 0 ) � a0; ��g (0)�x � a1 � 2�g ( 0) �x2 � 2a2; �3�g ( 0) �x3 � 6a3; �4�g ( 0) �x4 � 24a4 ( 2) 用区间上 n个测点值进行最小二乘拟合, 使 � = �n i= 1 (f ( xi ) - �g ( xi ) ) 2 = m in的待定系数 a0、 a1、�、a4必满足方程组: ���ak = 0 ( k = 0, 1, �, 4),展开此式,得到法方程组 a0 �n i= 1 1 + a1�n i= 1 xi + a2�n i= 1 x 2 i + a3�n i= 1 x 3 i + � � a4 �n i= 1 x 4 i = �n i= 1 �g i a0 �n i= 1 xi + a1 �n i= 1 x 2 i + a2 �n i= 1 x 3 i + a3 �n i= 1 x 4 i + � � a4 �n i= 1 x 5 i = �n i= 1 xi �g i a0 �n i= 1 x 2 i + a1 �n i= 1 x 3 i + a2 �n i= 1 x 4 i + a3 �n i= 1 x 5 i + � � a4 �n i= 1 x 6 i = �n i= 1 x 2 i �g i a0 �n i= 1 x 3 i + a1 �n i= 1 x 4 i + a2 �n i= 1 x 5 i + a3 �n i= 1 x 6 i + � � a4 �n i= 1 x 7 i = �n i= 1 x 3 i �g i a0 �n i= 1 x 4 i + a1 �n i= 1 x 5 i + a2 �n i= 1 x 6 i + a3 �n i= 1 x 7 i + � � a4 �n i= 1 x 8 i = �n i= 1 x 4 i �g i ( 3) 解此方程组即可得到 a0、a1、�、a4。当重力 剖面为等点距时,为使计算简便将 n取为奇数, 计 算区间中点置于坐标原点, 并以点距为单位,则 xi = i- 1-m ( i= 1、2、�、n, 且 n = 2m + 1)。即 xi 为区间 [ -m, m ]中的连续数,由于其对称性,故有 �n i= 1 xi = �n i= 1 x 3 i = �n i= 1 x 5 i = �n i= 1 x 7 i � 0, 于是与上述法 方程组对应的增广矩阵成为 n 0 �n i= 1 x 2 i 0 �n i= 1 x 4 i �n i= 1 g i 0 �n i= 1 x 2 i 0 �n i= 1 x 4 i 0 �n i= 1 xi�g i �n i= 1 x 2 i 0 �n i= 1 x 4 i 0 �n i= 1 x 6 i �n i= 1 x 2 i �g i 0 �n i= 1 x 4 i 0 �n i= 1 x 6 i 0 �n i= 1 x 3 i �g i �n i= 1 x 4 i 0 �n i= 1 x 6 i 0 �n i= 1 x 8 i �n i= 1 x 4 i �g i � 可将其化为二个方程组,即式 ( 4)和式 ( 5)。 na0 + a2 �n i= 1 x 2 i + a 4 �n i= 1 x 4 i = �n i= 1 �g a0 �n i= 1 x 2 i + a 2 �n i= 1 x 4 i + a4 �n i= 1 x 6 i = �n i= 1 x 2 i �g i a0 �n i= 1 x 4 i + a2 �n i= 1 x 6 i + a4 �n i= 1 x 8 i = �n i= 1 x 4 i �g i ( 4) a1�n i= 1 x 2 i + a3�n i= 1 x 4 i = �n i= 1 x i�g i a1�n i= 1 x 4 i + a3�n i= 1 x 6 i = �n i= 1 x 3 i �g i ( 5) 这种方程组可以按 �克莱姆法则 �用行列式直 接求解。 a0 = �11 �n i= 1 �g i - �21 �n i= 1 x 2 i �g i + �31�n i= 1 x 4 i �g i /� a2 = - �12 �n i= 1 �g i + �22 �n i= 1 x 2 i �g i - �32 �n i= 1 x 4 i �g i /� a4 = �13 �n i= 1 �g i - �23 �n i= 1 x 2 i �g i + �33�n i= 1 x 4 i �g i /� 亦即 a0 = �n i= 1 �11 - x2i �21 + x4i �31 � �g (x i ) = �n i= 1 p0i�g (x i ) p0i = ( �11 - x2i �21 + x4i �31 ) /� � ( i = 1, 2, �, n ) 为计算 a0时, 第 i个节点上的系数。式中 �为行 列式, 其中 �11、�21、�31则为对应行列的二阶子行 列式,当 n确定后它们都是常数。 � = n �n i= 1 x 2 i �n i= 1 x 4 i �n i= 1 x 2 i �n i= 1 x 4 i �n i= 1 x 6 i �n i= 1 x 4 i �n i= 1 x 6 i �n i= 1 x 8 i �11 = � n i= 1 x 4 i �n i= 1 x 6 i �n i= 1 x 6 i �n i= 1 x 8 i � �21 = � n i= 1 x 2 i �n i= 1 x 4 i �n i= 1 x 6 i �n i= 1 x 8 i �31 = �n i= 1 x 2 i �n i= 1 x 4 i �n i= 1 x 4 i �n i= 1 x 6 i 求 a2、a4的系数 p2i、p4i ( i= 1�n )方法同上。 而 a1和 a3的求法则要简单些,可直接写出。 a1 = �n i= 1 x 6 i � �n i= 1 xi �g i - �n i= 1 x 4 i � �n i= 1 x 3 i �g i �n i= 1 x 2 i � �n i= 1 x 6 i - �n i= 1 x 4 i � �n i= 1 x 4 i = �n i= 1 p1i�g i 206 物探化探计算技术 29卷 p1 i= xi �n i= 1 x 6 i - x 3 i �n i= 1 x 4 i �n i= 1 x 2 i � �n i= 1 x 6 i - �n i= 1 x 4 i � �n i= 1 x 4 i a3= �n i= 1 x 2 i � �n i= 1 x 3 i �g i - �n i= 1 x 4 i � �n i= 1 xi�gi �n i= 1 x 2 i � �n i= 1 x 6 i - �n i= 1 x 4 i � �n i= 1 x 4 i = �n i= 1 p3i�g i p3 i= x 3 i �n i= 1 x 2 i - x i�n i= 1 x 4 i �n i= 1 x 2 i � �n i= 1 x 6 i - �n i= 1 x 4 i � �n i= 1 x 4 i 按上诸式求出 a1、�、a4后,再按式 ( 2)求得重 力各阶方向导数值。据此可编制计算程序,计算时 在实测剖面上取小范围的奇数点作最小二乘拟合, 仅计算所取区间中点的导数值,然后同步移到下一 点。在实际应用中,我们要特别注意重力实测曲线 经常存在随机干扰成份, 要想取得好的导数效果, 就要求原始曲线比较光滑, 精度高且干扰小,否则 误差会很大。事实上我们在作最小二乘拟合时,得 到的 a0既为圆滑重力值, 亦可在作重力导数计算 前先对实测重力曲线进行圆滑,再求其导数。 关于计算点数 n的选取,要根据实测剖面的异 常特征和精度进行试验。一般情况下, n值较大 时对干扰有一定的消除或压制,但 �边部损失 �也 大, n值较小高阶导数异常有可能奇变。图 1是 n 值分别为 9、11、13时的三阶方向导数曲线, 随着 n 值的增大,导数异常曲线亦越光滑且异常形态变化 不大。 图 1� 不同 n值三阶方向导数比较 (单位 : 10- 9 /m2 s2 ) F ig. 1� The com parison o f the th ird d irectiona l der iva tive w ith var ious n 2� 应用实例 以津京地区地热勘探为例。津京地区蕴藏着 丰富的地热资源,地热储层主要为古生界、中上元 古界碳酸盐岩地层 (蓟县系白云岩和寒武 - 奧陶 系灰岩 ) ; 中生界、新生界的巨厚沉积对下部热储 层起到隔热保温作用,成为热储盖层;切入深度大、 属张性特征的深大断裂具有导热导水作用,是地热 地质构造的重要特征。因此, 在地热勘探中分析研 究断裂的形态、走向、产状等特征至关重要。 图 2为津京某地前期地热勘查布格重力异常 剖面。由图 2可见, 剖面二侧布格异常值相差较 大、落差较陡、曲线顺滑,反映出左侧抬升、右侧下 降的断裂构造形态。图 3为剖面的一阶方向导数, 它明显表现出在剖面 5 km处存在一条断裂。曲线 左侧的微弱上弯,使我们误认为是断裂左倾引起的 异常反映。按照常规的分析判断, 我们将按照一条 断裂模型去拟合地层的特征。 图 2� 布格重力异常剖面 (单位: 10 g� u) F ig. 2� The bouguer grav ity abnom a l section 图 3� 重力一阶方向导数异常 (单位: 10 E ) F ig. 3� The first d irectiona l der ivative anom a ly 但与重力剖面位置相同的直流电测深剖面却 反映出了不尽相同的结果 (图 4), 电测深 �s断面 图上清晰地勾画出了剖面左侧第四系厚度薄而右 侧厚,这与布格重力异常的特征是一致的,但 �s断 面反映的断裂形态是一组断阶,位置分别在剖面的 5 km和 8 km 处附近, 二断裂间明显存在一个台 阶。图 5(见下页 )是布格重力异常剖面的三阶方 向导数,它所反映的断裂特征与 �s断面完全一致, 显示出高阶导数在区分微弱叠加异常方面的明显 优势。 图 4� 直流电测深剖面 (单位: � � m ) F ig. 4� The vertica l electr ic sound ing section by DC me thod 图 6(见下页 )是重力正演拟合的地质断面,拟 2073期 田舍等: 重力高阶方向导数的可视化拟合及其应用 图 5� 重力三阶方向导数异常 (单位: 10- 9 /m2� s2 ) F ig. 5� The th ird d irec tiona l der iva tive anoma ly 合前对布格重力异常进行了必要的区域场校正,根 据电测深剖面和布格重力三阶方向导数反映的特 征,建立初始地质模型。根据前述计算过程,我们 改进了目前常用的重力剖面可视化正演拟合软件。 经反复调整,最终达到地质模型的正演值及其三阶 方向导数与实测剩余布格重力异常及三阶方向导 数形态相一致。图 7为前例重力三阶方向导数拟 合情况。经后期钻探证实,该地前期的物探分析和 推断解释是完全正确的。 图 6� 拟合地质剖面 F ig. 6� The fitting geo log ical section 图 7� 重力三阶方向导数拟合形态 F ig. 7� The fitting shape o f the th ird d irec tiona l der ivative 3� 结论 重力剖面正演拟合软件一般采用二度半体重 力正演公式计算地质模型的场值, 通过不断调整其 形态与实测布格重力异常进行拟合。软件程序正 演计算的是近地表有限范围内,且具有二度体特征 的,地质模型体所产生的重力场值, 而实测布格重 力异常还包含了地下深部,以及测区周边地质因素 的影响。因此, 在做正演拟合时需先进行区域场校 正,但在实际工作中很难准确做到。在这种情况 下,如果我们一味追求重力异常的精确拟合, 可能 会导致将区域地质因素添加到模型地质体中而造 成错误的结果。 重力勘探的目的是了解基岩 (或目标 )地质体 的特征,布格重力异常的高阶导数在突出浅部地质 体的特征, 压制区域性深部地质的因素, 区分相邻 地质体引起的叠加异常等方面具有独特的优势,能 为研究者进行正演拟合提供很多有用信息。但如 果一味追求高阶导数的精确拟合, 同样会导致错误 的结果,因为实际工作中带来的随机干扰亦有可能 在消除不净的情况下被放大。 因此在利用可视化重力剖面分析解释软件建 立和调整地质模型时,应兼顾布格重力异常及其各 阶导数的特点, 不可过份追求某一方面的精确拟 合,而应以二者拟合曲线的基本形态相似为原则, 尽可能地从其它物探方法和地质资料中,寻找出相 关信息加以佐证,这样才有可能取得预期的效果。 参考文献: [ 1] � 姚长利, 黎益仕 ,管志宁. 重磁异常正反演可视化实 时方法技术改进 [ J]. 现代地质, 1998, 12( 1): 115. [ 2] � 王家林, 王一新, 万明浩, 等. 石油重磁解释 [M ]. 北 京:石油工业出版社, 1991. [ 3] � 李祺. 物探数值方法导论 [ M ]. 北京: 地质出版社, 1991. [ 4] � 刘云生,译. 三维地质及位场数据的正演模拟 [ J]. 物 化探译丛, 1998, (增刊 ): 98. [ 5] � 田黔宁, 吴文鹂 ,管志宁. 任意形状重磁异常三度体 人机联作反演 [ J]. 物探化探计算技术, 2001, 23( 2): 125. [ 6] � 吴文鹂,管志宁. 可视化技术及混合优化算法在重磁 三维反演中的应用 [M ] .北京: 中国地质大学, 1997. [ 7] � 何昌礼.复杂形体的高精度重力异常正演方法 [ J]. 物化探计算技术, 1988, 10( 2): 121. [ 8] � BARNETT C T. Theoretica l modeling of the m agne tic and g rav itationa l fie lds o f an arb itrarily shaped three�di� m ens iona l body [ J] . Geophy sics, 1976, 41( 6): 1353. [ 9] � POHANKA V, OPT IMUM. Expression fo r compu tation o f the g rav ity fie ld of a hom ogeneous po lyhedra l body [ J]. Geophysical Prospecting, 1988, 36( 7) : 733. [ 10] � THOMAS ENMARK. A versatile ine ractive com puter program for com putation and autom atic optim iza tion o f grav ity m ode ls[ J]. Geoexploration, 1981, 19( 1): 47. 作者简介: 田舍 ( 1957- ) ,男, 高级工程师,长期从 事石油地质、地热地质、区域地质等领域地球物理 勘探工作。 208 物探化探计算技术 29卷 of se ism ic data, and determ ine the litho logy variation of coa l- bearing strata, and the lithology o f coa l beds and roo f/bo ttom o f coal beds. Key words: impedance inversion; litho log ical se is� m ic survey ing; simulation log data THE CALCULATION OF FORMATION STR� ESS BY CONVENTIONAL WELL LOGGING AND FM I LOGG ING DATA LIU Gao�bo1, FENG Wen�guang1, TAO X iao�hong2, et a.l ( 1. Co llege of Energy Resources, Chengdu Un i� versity o f Techno logy, Chengdu 610059, Ch ina ; 2. Cen tra l Sichuan M in ing D istric,t Southw est O il and Gas F ield Branch Company, Petroch ina, S ichuan Suining 629001, China ). COMPUTING TECH � N IQUES FOR GEOPH YSICAL AND GEOCHEM ICAL EXPLORATION, 2007, 29( 3) : 200 The ca lculation of formation stress by conven tion� alw e ll logg ing data is of the advantages of low cos,t easy to access data and acquire successive formation section of fo rmation stress. The paper introduces the methods o f interval transit tim e fitt ing of S�wave w ith P�wave and Po isson rat io calcu lation. A new model for fitt ing compressive resistance by integrat ing clay content and P- w ave impedance is put forw ard. The paper a lso calcu lates tectonic stress by FM I logging and hydro fracturing data. F ina lly, the form ation stress sections o f tw enty five w ells are ca lculated w ith the appropriate ca lculat ing mode.l Key words: fo rmation stress; conven tiona l w ell log� ging; FM I; rock mechan ics parameters VISUALIZATION FITTING AND APPLICA� TION ON HIGH - ORDER DERIVATIVES IN GRAVITY TIAN She 1 , SONG Yu�chen2 ( 1. Beijing Institute o f Geo log ica l Explorat ion and Techno logy, Be ijing 102218, Ch ina; 2. Ch inaUn iversity ofG eosciences, Be jing 100083, Ch ina ) . COMPUTING TECH � N IQUES FOR GEOPH YSICAL AND GEOCHEM ICAL EXPLORATION, 2007, 29( 3) : 205 A rm ed at the defects in the app licat ion o f the conventional interact ing g rav itationa l profilem ainly u� sing forw ard of Bouguer grav ity anomaly, F itting w ith high- order derivat ives is adopted w ith the advantages to feature low er geo log ical bodies, suppress reg ional deeper geo log ica l facto rs and recogn ize the deference of anomalous caused by borders on geo log ical bod ies. An examp le is g iven to show the improvem ent o f grav� ity data interpretation using high- order deriva tives in geophysica l prospect ing. F inally, it is ana lyzed for fitt ing exp lanation on v isualization in grav itational pro� files, and at the same t ime it is stressed that it is sim� ilarly in fitt ing curves both Bouguer gravity anoma ly and high - order derivative, wh ile it is no t pursued the exp lanat ive principle in a sing le fitt ing precision. Key words: grav ity data; h igher derivatives; expla� nat ion on v isualization STUDY ON DATA COMBINATION AND THE METHOD OF CORRECTING ERRORS OF MULTIPLE ARRAYS IN RESISTIVITY IMAG� ING METHOD YAO Jian, ZENG Zhao�fa, LI Fang ( Co llege of G eo� Exploration Science and Techno logy, Jilin Un iversity , Changchun 130026, Ch ina) . COMPUTING TECH � N IQUES FOR GEOPH YSICAL AND GEOCHEM ICAL EXPLORATION, 2007, 29( 3): 209 As the resistiv ity im ag ing method data o f two ad� jacent arrays have the overlaps, w e present a method to combine the arrays and remove the error. To con� tro l the effect o f error, w e apply a method to so lve the data combinat ion o f arrays system ic error correct ion based on the least squares me thods. By applcating, it is proved to be pract ical and e ffect ive in so lv ing the error in the measuremen.t Key words: resistivity imag ing; comb ination of ar� rays; system error- correction; least squares methods SIMULATION OF ADAPTIVE ANTENNA FOR GPR APPLICATION ON THE LAYERED SO IL YUAN Guo�sheng, LIYong�m ing( TheK ey Labo rato� ry of H igh Vo ltage Eng ineering & E lectrical New Techno logy, M in istry of Educat ion, E lectrical Eng i� neering Co llege of Chongqing U niversity, Chongq ing 400044, Ch ina ). COMPUTING TECHNIQUES FOR GEOPH YSICAL AND GEOCHEM ICAL EXPLORA� TION, 2007, 29( 3): 213 A new type of adaptive antenna for g round pene� trat ing radar is described in th is paper. The adaptive antenna is based on a strip bow - tie structure since this structure make it possible to sw itch on /o ff its strip elements in order to form a bow - t ie an tenna w ith certa in flare ang le. By adjusting the flare ang le in this w ay, the shape and size o f the subsurface foot� print and the inpu t impedance of the antenna can be varied. The result o f the simu lation w ith the fin ite in� tegration techn ique indicates that the strip bow - tie antenna can achieve the above design ob jec.t Key words: ground penetrating radar; fin ite integra� �2� ABSTRACTS� � � � � � � � � � � � � � � � � � � Vo .l 29� No. 3