第八讲 重力异常反演.docx

1、第八讲 重力异常反演应用重力学第八讲重力异常反演 正问题(Forward Calculation) 已知地质体的形状、产状和剩余密度等，通过理 论计算来求得异常的分布和规律。m d 正演 反问题(Inversion) 已知异常的分布特征和变化规律，求场源的赋存 状态（如产状、形状和剩余密度等）。? 解正问题是解反问题的基础，解反问题是目的。 仅从地质角度，解重力反演问题的目标9 矿体类问题:寻找、研究或推断金属或非金属 矿体；9 构造类问题:研究地质构造，包括控矿构造， 如含石油、天然气、煤的构造以及区域性的深 部构造等。 从地球物理角度，解重力反演问题的目标9 矿体类问题:确定地质体的几何和

2、物性参数；9 构造类问题:确定物性分界面的深度及起伏；9 密度分布问题:确定密度的分布。一、计算地质模型体的几何及物性参数(一)直接法 直接利用由反演目标引起的局部异常，通过某种积分运算和函数关系，求得与异常分布有关地质 体的某些参量。 三度体剩余质量的求法 m n m ij g .u. m g .u. M = 2.386xyg + 2 R2 g (R) i =1 j=1 三度体重心水平坐标的求法2.386 m n y x g x0 x m i m ij g .u. M t2.386 yi =1 j=1 m n y g y0 x m i m ij g .u. M ti =1 j=1 二度体横截

3、面积的求法nx g i = i =1 4Gtg 1 xD 二度体横截面重心水平坐标的求法x0 =n ( x g )i i =1 n g ii =1(二)特征点法 根据异常曲线上的一些点或特征点（如极大值 点、零值点、拐点）的异常值及相应的坐标求取 场源体的几何或物性参数; 仅适用于剩余密度为常数的几何形体。异常曲线形态分类 第一类是单峰异常，零值点在无穷远处 如球体的g曲线、台阶的Vxz曲线等； 第二类是具有极大值、极小值和一个零值点 如球体的Vxz曲线、台阶的Vzz、Vzzz曲线； 第三类是具有一个极大值、两个极小值和两个零值点 如球体、水平圆柱体的Vzz和Vzzz曲线； 第四类是台阶的g曲

4、线，一边高一边低的形态应用条件 对异常作平滑处理，尽量准确确定原点的位置； 对异常曲线作分离处理，获得单纯由研究对象引 起的异常； 对剩余（局部）异常进行分类，判明该异常的场 源体接近于何种可能的几何形体，然后选用相应 的反演公式。球体g的反演(第一类曲线)单峰异常，零值点在无穷远处GMD GMDg1 = 2 2 3 / 2 g2 = 2 2 3 / 2 ( x1 + D ) ( x2 + D ) g ( x2 + D2 )3 / 2 n = 1= 2 2 1 g ( x2 + D2 )3 / 2 1/ 2 x2 x2 n2 / 3 D = 21 n2/ 3 1 D = 12 x1/ n n2

5、/3 x1/ n 1当n=2时D = 1.305 x1 / 2= 0.6524( x1 / 2 x1/ 2 )当n=3时当n=4时D = 0.9622 x1 / 3D = 0.8111x1 / 4= 0.4811( x1 / 3= 0.4056( x1 / 4 x1/ 3 ) x1/ 4 )maxM = 14.99D2 g 1/ 3 1/ 3 R 3M0.62 M t = = m 4 ( 0 ) h = D R g / cm3 M 实 = g / cm3M tm m D= 1.305x1/ 2 M = 14.99D2 g t m max1/ 3g .u. mR =0.62 M t gc/ m3

6、 DR = 1g / cm3m m D= 1.305x1/ 2 = 1.305 153 = 199.67mM = 14.99D2 g 200mt m maxg .u. = 14.99 199.67 199.67 6.98 = 4.1714 106 t1/ 3 mR =0.62 M t gc/ m3 = 0.62 (4.1714 106 / 1)1/ 3 = 99.81m100m球体Vxz的反演(第二种类型的曲线)具有极大值、极小值和一个零值点1 1 xmax= D, xmin= D 2 2 (V ) = (V ) = 48GMxz maxxz min25 5D3D = 2 xmin= 2 xm

7、ax= xmin xmaxxzM = 0.01746D3 (V )max球体Vzz的反演(第三种类型曲线)具有一个极大值、两个极小值和两个零值点(V ) =2GMzz max D3x0 = 2Dm , x0 = 2Dm zzM = 0.00749D3 (V )max0 zz max = 0.00265x3 (V )球体Vzzz的反演(第三种类型曲线)(V ) =6GMzzzmax D4D = 1.225 x0M = 2.499D4 (V )103 = 5.622x4 (V )zzzmax 0zzzmax台阶曲线的反演(第四种类型曲线)台阶的g曲线，一边高一边低的形态gmax= g ( x) +

8、 g ( x) = 2 G2 G ( 2 1x1/ n ) = G + tg n D D = x1/ n tg (2 n)2n1 = g ( x) + g ( x) 2 Gh x1/ n 1 ( x) + g ( x) 4 G = (2 n)tg 2nH x1/ n 1 g ( x) + g ( x) 4 G = + (2 n)tg 2n(三)选择法 根据异常分布和变化特征，结合地质和其他地球 物理和物性等资料，给出初始地质体模型; 进行正演计算，将理论异常与实测异常对比;m k k k k 1 2 n = gk =1-f(x ,y ,z ,b ,b , b )2= min 若两者偏差较大，对

9、模型进行修改，重算其理论 异常计算，再次进行对比; 如此反复进行，直至两种异常的偏差达到事前要 求的误差范围为止，则这最后的理论模型就可作 为所求的解答了。修改模 型参数否初始模型计算理论重力异常与实测异常对比是否满足精度是反演结果选择法特点 异常可以是整条剖面或整个测区的数据，受个别 点误差的影响较小，抗干扰的能力较强。 所求的地质体可以是一个或几个复杂的不规则几 何形体,密度分界面,或者密度的分布。 需要重复而复杂的正演计算, 可编制相应程序由 计算机来自动完成。 解释复杂重力资料时，能够考虑研究区已知的地 质构造资料，在反演过程中利用这些资料，控制 或约束计算结果，使得到的地质模型更接近

10、实际 的地质体。(四)人机交互式反演方法(姚长利-重磁异常正反演解释系统)二、计算密度分界面的深度 密度分界面与区域构造和储油构造有密切的关系，因此计算密度分界面的起伏或深度的变化在 区域构造研究和石油勘探中具有重要的意义。(一)线性回归法 如果界面起伏平缓，可以认为重力变化与界面的 起伏近似呈线性关系。h = a + b g n i i (a, b) = (h h )2i =1= min 令 = = 0, = 0, 则 a b 2 h n a = =gi hi gi gi i b = = gi hi hi gi ,n g 2 (g )2 n g 2 ( g )2 i i i i (二)压缩质

11、面法(刘云龙,1977)条件：界面起伏较小，埋藏深度较大。基本原理 将界面从最小深度h和最大深度H处向中间挤压，使之在 界面平均深度D=(h+H)/2上压缩成一个面密度不均匀分布 的物质面， 将该物质面剖分成局部面密度均匀分布的水平物质带(二维)或物质片，面密度 j = hj 计算物质带或物质片的正演，或迭代反演物质带或物质片 的面密度，进而求出界面深度。g (i) 2G tg 1 2(j 1) i 2( j 2) 11 n = tg jj =1 2 2 (三)迭代法(Cordell,1968)gyOt Q长方体zxP测点 由无限平板重力异常公式给出t初值Mt1,q= gobs ( q ) 2

12、 G 计算模型初值的重力异常gcalc,n, p= q =1gGf (P, Q, tn,q , , D) 计算t的下一个修改值t n +1,q =t n,q( obs ,q )gM obs ,qcalc,qcalc ,n,q 目标函数 = ( g g )2 = min q =1 判断计算结果满足要求否，不合要求则转到第2步 继续计算；否则停止计算。t1,q= gobs ( q ) 2 Ggcalc,n, pM= q =1Gf (P, Q, tn,q , , D) n=n+1t n +1,q =t n,qg obs ,q( )gcalc ,n,qM obs ,qcalc ,q = ( g g )

13、2 = min? q =1否 是反演结果理论界面深度 正演的重力异常m d=Gm理论界面深度 反演界面深度(迭代50次)m m=G-1d(四)频率域反演法(Parker-Oldenburg,1973) 假定在x-z直角坐标系中，重力异常用g(x)表示，场源层的上部边界为z=0，下部边界为z=h(x)，这个边界显示界k面的起伏。 根据帕克的二维傅里叶变换公式得到重力异常傅立叶变换F g ( x) =2G e kz0 n=1n1k F h n ( x) n! 从上式的无限和式中提出n=1的项，并重新排列，得到F h(z0x) = =F =g ( x)e n 1k F h n( x) 2 Gn =2

14、 n! 假定已知地层与下部介质之间的密度差及参考面深度z0已知或给定，就可以应用上式进行下列迭代计算：(1) 给界面起伏h(x)的初值，例如h(x)=0，(2) 将h(x)的初值代入上式的右端项，计算右端项的傅立 叶变换，(3) 右端项的傅立叶反变换即改进的界面起伏h(x)(4) 判断：计算结果是否满足某个收敛标准，或是否达到 给定的最大迭代次数，如果是，即停止计算；否则转到第2 步，以本次迭代结果作为初值，继续迭代计算。 级数的收敛与界面的起伏大小有关。记H为界面相 对于参考面的最大起伏值，当H/z01时，级数均匀 收敛并与波数k无关；当H/z0=1时，即界面与观测 面相交时，级数不收敛，而

15、且，界面越接近观测面 收敛速度越慢。 上式右端第一项包含了一个指数项，这就会放大数 据中的高频干扰，故在右端项乘以一个低通滤波器 来加以限制。反演正演三、计算地层密度的分布 根据重力异常反演地下密度的不均匀分布。gXXZ 根据密度分布图，就可以判断是否存在高(低)密度体及其大小、位置，或发现密度的分界面，再 结合地质资料就可以作出地质上的结论。 这个方法思路简单、清楚、直观，但是实现起来 相当困难。 这种方法目前尚处于理论研究阶段或限于推断深 部地质构造问题，很少用于解决实际重力勘探问题。(一)二维密度成像 g = 2G x ln( x2 + z 2 ) + 2 1 x in nk +1 l

16、+1 = a ztg zk l Nin ngi = ain n , i = 1, 2, ., m n =1a111a211+ a121+ a221+ + a1N N = g1+ + a2 N N = g2 am11 + am 21 + + amN N = gm A P = G (二)三维密度成像gyO xz = ( 2 + 2 + 2 )1/2 g = G ln( + ) + ln( + ) tg 1= x2y2 z2 =g G ln( + ) +ln( + ) tg 1= x1 y1 z1 x2 y2 z2 x1 y1 z1 = Fg = F 下半空间剖分为M个长方体，观测面上有N个观测点g

17、ij = Fij jMgi =j =1Fij j ,i = 1, 2, ., N g1 F1F112 F1M 1 g2 = F2F1 22 2 g N FN 1 FNM M G = F P (姚长利，2007)(刘天佑，2007)(刘天佑，2007)重 计算地质模型体的几何及物性参数力 直接法、特征点法、选择法、异 人机交互反演法常反 计算密度分界面的深度演 线性回归法、压缩质面法、迭代法、方 频率域反演法法 计算地层密度的分布二维、三维密度成像四、反问题解的多解性 处在不同深度的、 具有不同起伏的界 面，还包括在这些 界面之间的许多界 面，都能够引起在 测量精度范围内的 相同的重力异常。(据

18、Skeels，1947)内特尔顿(1987)指出： 与一定的异常宽度对应的场源的可能的最大深度，就是引 起同样宽度异常的点源(球体)的深度。 在这个最大深度和地面之间，存在一个可能源的锥形区。 同一个异常可以由埋藏很浅的薄透镜体，或不太狭窄的较 厚的物体，或球体引起，而且在球体上方较浅的深度上， 却有无限多个可能的场源引起这个相同的异常。 引起同一异常的不同的场源，有一个唯一的共同的特性， 就是它们的剩余质量必须是相同的。引起多解性的原因 场的等效性； 观测数据离散、有限; 实测的异常包含一定误差的； 数据整理带来的误差。 反问题解的非唯一性是客观存在的。限制多解性的方法 除了应用重力资料外，

19、还引用工区内的地质、钻 井、物性和其他地球物理资料等，尽可能地增加 已知条件和约束条件，则反问题的解答数目就会 大大减少，甚至可以得到单一的解。如球体反问题的求解，如果已知剩余密度值，则 它的半径和顶部埋深也就唯一地被确定了。 提高仪器的测量精度。 改进各项校正的计算方法，使得能够更加精确地 得到重力异常的分布。(刘天佑，2007)(刘天佑，2007)五、重力归一化总梯度法 原苏联地球物理学家别列兹金于1967年提出了重力 归一化总梯度法，它是一种利用重力场中的特征点(如极大点或极小点)或者解析函数的奇点，探测重 力异常的场源体并估计其位置的方法。 在物体的内部，引力位V满足泊松方程2V =

20、4 G 在物体的外部，引力位V满足拉普拉斯方程2V2V2V2V = + + = 0x2 y 2 z 2 剩余质量的引力位及导数在场源体以外空间都是解析函数，而在场源处则失去解析性。 在解析函数中，失去解析性的点叫函数的奇点。 确定场源的问题，就是通过对异常解析延拓来确定 函数的奇点问题。以二度水平圆柱为例g( x, z) =2G zz 2 + x 2 在x=0处，把g(x,y)延拓到圆柱中心z=0处，在圆柱中心处的g(0,0)将变为无穷大。 圆柱中心是极点类型的奇点。寻找重力场的奇点有一些困难： 实测异常是由解析成分（有用信号）及随机干扰组 成，要用实测异常求奇点首先必须把干扰滤掉； 实测异常

21、为离散形式，而且其观测剖面长度是有限 的，在确定奇点时，需施加某些限制； 不知道实测异常的解析式，无法求其奇点。 二度水平圆柱重力异常gx( , z) = 2GDD2 + x2 向下延拓到任一点P(x,z)时g ( x, z) =2GD z 若沿Z轴(x=0)g ( z) =2G(D 1z) 2+ x 2D z 在z=0的附近将上式展开成台劳级数为g ( z) 1 n! 0)n = 2Gn+1z =0( z n=0n! (D z) n2 1 n 2G z = = G n+1 z n=0 D D n=0 D 当向下延拓趋于场源(zD)时,级数发散，即D为 g(z)的一个奇点。 当级数不取无限项时，情况就不同了，这时2GN z ng ( z)= D n =0 D 通过柱体中心(zD)，因式是一有限项级数，只有当z时才发散，否则只要N为有限值，级数 都是收敛的。 所以向下延拓直接求与场源质量有关的奇点是很 困难的。 作如下改进：2GnN z 2Gn z g ( z)= + D n =0 D D n = N +1 D 右边的第二项是级数的余项，也即舍去误差。 令=z/Dg ( ) =g N () + rN ( ) 若N取适当值，则1时，有rN()1时将有rN()gN()，用它可以 确定与场源有关的奇点。 以rN()对gN()进行归一化，则g