（固体地球物理学专业论文）基于重力异常的密度界面反演研究.pdf

摘要 摘要 重力异常界面反演是重力勘探中的一个重要研究课题，它是一个求解第一 类f r e d h o l m 非线性积分的问题。目前情况下直接求解这类问题比较困难，般 采用拟合法或逐步线性化办法实现反演计算。而在计算过程中，由于反演的不 适定性，往往不能很好地获得反演解。 本文先从二度体重力异常正演公式出发，研究了重力异常密度界面反演问 题，导出地下界面迭代反演的计算公式；然后利用广义似然函数对模型参数进 行优化选择，从而获得模型参数修正量迭代计算方法，实现二度体单一密度界 面的稳健迭代( 最小二乘加权迭代) 反演。模型计算表明方法可以很好反演出 地下界面。由于实际资料不能完全分离出正常场，为了消除正常场的残余量对 密度界面反演的影响，我们又用重力异常水平梯度束反演密度界面。先取重力 异常的水平梯度，然后利用广义似然函数对模型参数进行优化选择，从而得到 地下模型参数修正量的迭代计算公式，将l s q r 方法引入到该方法中来，用它 求解大型方程组，得到模型参数的修正量，从而实现了二度体单一密度界面的 反演，模型计算表明此方法是行之有效的。 最后讨论了正常场对重力异常密度界面反演和重力异常水平梯度密度界面 反演的影响。 关键词：界面反演重力异常稳健迭代水平梯度l s q r 方法 敏前皮 k 8 7 ，甲 a b s t r a c t a b s t r a c t i n v e r s i o no fd e n s i t yi n t e r f a c eb yu s i n gg r a v i t ya n o m a l yi sa l li m p o r t a n tr e s e a r c h s u b j e c ti ng r a v i t yp r o s p e c t i n g ，w h i c hi sap r o b l e mt os o l v et h ef i r s tf r e d h o l m n o n l i n e a ri n t e g r a le q u a t i o n i ti sd i f f i c u l tt os o l v et h ep r o b l e md i r e c t l y , a n dg e n e r a l l y f i t t i n gp r o c e s sa n dg r a d u a ll i n e a r i z i n gp r o c e s s e sa r eu t i l i z e dt oi n v e r s et h ei n t e r f a c e i nc o m p u t a t i o n a lp r o c e s s ，i ti sh a r dt oo b t a i ng o o di n v e r s i o nr e s u l t sd u et oi 1 1 p o s e d p r o b l e m i nt h i st h e s i sw ef i r s te d u c et h ei t e r a t i v ei n v e r s i o nf o r m u l af r o mg r a v i t ya n o m a l y f o r m u l ao f2 一dd e n s i t yi n t e r f a c e ，a n dt h e no p t i m i z et h em o d e lp a r a m e t e rb yu s i n g g e n e r a l i z e dl i k e l i h o o df u n c t i o n t h u s ，t h ei t e r a t i v em e t h o df o rc o m p u t i n gc o r r e c t i o n v a l u eo fm o d e lp a r a m e t e ri so b t a i n e d ，a n dr o b u s ti t e r a t i v ei n v e r s i o n ( 1 e a s ts q u a r e w e i g h t i n g i t e r a t i v ei n v e r s i o n ) o f2 - dd e n s i t yi n t e r f a c ei sa c h i e v e d m o d e l c o m p u t a t i o ns h o w st h a tt h ei n t e r f a c ec o u l db ei n v e r s e dw e l l 。a sr e a ld a t ac a nn o tb e c o m p l e t e l ys e p a r a t e df r o mn o r m a lg r a v i t yf i e l d ，i no r d e rt oe l i m i n a t et h ei m p a c to f t h er e m n a n t so ft h en o r m a lg r a v i t yf i e l do nd e n s i t yi n t e r f a c e ，w eg e tt h eh o r i z o n t a i g r a d i e n to fg r a v i t yt o i n v e r s i o nw i t ht h em e t h o d t h e nw eo p t i m i z et h em o d e l p a r a m e t e r sb yu s i n gg e n e r a l i z e dl i k e l i h o o df u n c t i o n a f t e rt h el a r g e s c a l e de q u a t i o n i so b t a i n e d ，w ec h o o s et h el s q rm e t h o dt os o l v ei t q u i c k l ya n dr o b u s ti t e r a t i v e i n v e r s i o n ( 1 e a s ts q u a r ew e i g h t i n gi t e r a t i v ei n v e r s i o n ) o f2 dd e n s i t yi n t e r f a c ei s a c h i e v e d m o d e lc o m p u t a t i o ns h o w st h a tt h em e t h o di se f f i c i e n t i nt h ee n d ，w ed i s c u s st h ee f f e c to ft h en o r m a lg r a v i t yo nt h ei n v e r s i o no f2 - d d e n s i t yi n t e r f a c eb a s e do ng r a v i t ya n o m a l ya n dt h eh o r i z o n t a lg r a d i e n to fg r a v i t y a n o m a l y k e y w o r d s ：i n v e r s i o no f2 一dd e n s i t yi n t e r f a c e ，g r a v i t ya n o m a l y , r o b u s ti t e r a t i o n ， t h eh o r i z o n t a lg r a d i e n to fg r a v i t y , l s q rm e t h o d 1 1 中国科学技术大学学位论文原创性和授权使用声明 本人声明所呈交的学位论文，是本人在导师指导下进行研究工作 所取得的成果。除已特别加以标注和致谢的地方外，论文中不包含任 何他人已经发表或撰写过的研究成果。与我一同工作的同志对本研究 所做的贡献均已在论文中作了明确的说明。 本人授权中国科学技术大学拥有学位论文的部分使用权，即：学 校有权按有关规定向国家有关部门或机构送交论文的复印件和电子 版，允许论文被查阅和借阅，可以将学位论文编入有关数据库进行检 索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。 保密的学位论文在解密后也遵守此规定。 作者签名：迅 o a 年7 月幺日 前言 1 厶j l 日| j 舌 重力异常是地下各物性界面起伏与岩性不均匀变化的综合反映。在一般情 况下，地下岩石密度的不均匀往往和某些地质构造的变异或某些矿产分布有关， 因此利用重力资料研究地壳密度界面分层，对于了解地壳内部结构、研究区域 地质及探测金属和非金属矿产，都有十分重要的实际意义。 关于地球内部的信息，绝大多数来源于对地表地球物理观测资料的解释。 利用勘探来源于对地表地球物理观测资料的处理和解释。利用勘探地球物理资 料推算地下结构、构造和物性变化，一直是勘探地球物理资料处理和解释的重 要研究课题，它是一个十分复杂的理论和应用问题。为了正确地进行解释，必 须具备两个方面的知识：一、根据已知地质体的形状、产状，研究它所引起异 常的特征、数值大小、分布和变化规律等，即地球物理学的j 下演问题。二、根 据异常特征、数值大小、分布和变化规律，确定地质体的形状和空间位置，即 地球物理学的反演问题。任何一种地球物理观测，最终都要求刻画地球内部介 质的结构和构造，根据各种地球物理观测资料，来模拟推算地球内部的结构形 态和物质成分，定量计算地球物理参数，解释地球内部的奥秘。 重力异常密度界面反演是已知重力异常在地面某区间上的观测值，柬确定 异常源的某个密度界面。密度界面的传统反演方法都要求界面平均深度比界面 起伏大得多，6 0 7 0 年代盛行的最优化法是一种非线性解法，但由于多解性使得 该法的应用受到一定的限制。此后，约束条件下的最优化求解界面深度方法出 现，在一定成都上解决了上述问题( 陈钟琦等，1 9 9 1 ) 然而这种约束解也存在 初值和特征点的选定问题。 p a r k e r ( 1 9 7 3 ) 在位场界面正演计算中引入快速傅立叶变换，首次提出频率 域的重力异常正演公式，为密度界面深度分布的重力非线性频率域反演提供了 理论依据，其方法具有速度快、功能强的特点，成为重磁异常定量解释中的工 具，适合于含油气构造区的地球物理资料综合解释及联合反演；缺点是为避免 迭代计算发散而进行异常及反演深度值滤波，所以深度分辨率受到损失。 其后o l d e n b u r g ( 1 9 7 4 ) 在此基础上，提出了确定起伏密度界面的迭代反演 法。由于该方法中包含了一个不稳定因子向下延拓算子，导致了本方法的刁i 前缶 稳定性，通常的做法在反演过程中加入一个低通滤波器，以限制高频振荡 ( g o m e z o r t i z ，2 0 0 5 ；张作新等，1 9 9 8 ) ，但这种做法在消除高频假波的同时必然 会对有效信号进行一定的压制，造成反演精度的降低。 为了提高反演的精度，一些学者对f f t 进行各种改进：在向下延拓因子上 附加正则化因子以提高精度( 张风旭等，2 0 0 5 ) ，但向下延拓的精度和稳定性， 没有得到实质性的提高。 随着新的观测手段特别是卫星技术的应用，地球重力场的观测信息越来越 丰富。除传统的地面重力外，还有卫星测高、g p s 水准、卫星跟踪卫星( s s t ) 、 卫星重力梯度( s g g ) 和航空重力等观测值。综合利用这些观测资料，将有助 于提高重力异常的计算精度( 王海瑛等，2 0 0 0 ) ，上述任何一种完全覆盖的观测 值都可唯一确定地球重力场。采用多种数据的求解模型将是一个超定模型，求 解重力异常的超定模型将具有实用价值。据此给出了基于第一类f r e d h o l m 型积 分方程离散型解法求解重力异常超定模型，给出求解病念方程的正则化算法。 采用正则算法能有效提高参数估计的精度，其中，应用最广泛的是t i k h o n o v 正 则化方法( g r o e t s c h ，1 9 8 4 ；李浩等，1 9 9 2 ；李功胜。2 0 0 3 ；刘继军，2 0 0 5 ； 张 瑞等，2 0 0 6 ) 。 还有作者提出样条函数法( 汪炳柱等，1 9 9 5 ；王硕儒等，1 9 9 6 、2 0 0 0 ) ，不 像频率域的p a r k e r 公式那样只能考虑有限项的影响，因而不存在忽略高次项所 带来的误差，样条函数法比较灵活。李会元等( 1 9 9 9 ) 将确定异常源的密度界 面问题归结为求非线性泛函极小化的变分问题，并在h i l b e r t 空蒯中用 g a u s s n e w t o n 迭代方法进行重力异常界面反演，并采用隐式迭代，增强了解的 稳定性。g u s p i ( 1 9 9 3 ) 在o l d e n b u r g 公式的基础上，做了改进，对幂级数系数 的表达式，只列出了前4 项，在界面反演时太少了，并且没有详细讨论幂级数 的收敛和滤波问题。赵冠华等( 1 9 9 8 ) 又在g u s p i 的基础上做了很大改进，较好 地进行了密度界面的菲线性非迭代反演计算方法研究。 一般而言，反演总会因实际资料的观测误差及其数量的有限，而表现为不 适定性。重力异常界面的反演问题，常归结为求解个关于未知界面函数的非 线性泛函的最优化问题或者变分问题。求解此类问题目前多采用迭代方法，其 中，非线性模型的线性化或局部线性化是整个求解过程的关键。线性化的结果 是，把重力异常的小扰动量近似地表示为界面的小扰动量的f r e d h o l m 第一类积 分方程关系：而在考虑它的离散实现时，又多化为有关迭代修正量的一个线性 前言 代数方程组来求解。 目前，多数的线性化迭代方法及其离散求解过程都基于若干离散数掘之上， 有的明显地缺乏严格的理论依据，有的虽然求得简单的模型，却因过多地舍弃 了有用信息，从而导致了不收敛过程。 高尔根等( 2 0 0 0 ) 提出了二维声波方程稳健迭代速度反演方法，我们从中 得到很大启发，将稳健迭代的反演方法推广到了密度界面反演中去，提出了稳 健迭代的密度界面反演方法。我们从二度体重力异常正演公式出发，引入界面 起伏，得到地下界面迭代反演的计算公式。为避免反演的不适定性，定义广义 似然函数 pp f ( c ，p ) = l p - ( n 7 芦 。 2 3 i ( i 、 对模型参数进行优化选择，即误差分布满足模型参数的优化选择时，广义似然 函数取极大值。然后对积分方程离散化，迭代求解参数的修正量，这实际上是 一个迭代加权最小二乘问题。b y r d 和p y n e ( 1 9 7 9 ) 研究表明：迭代加权最小二 乘问题。当范数p ( 1 尸 2 ) 在满足相当弱的情况下是完全收敛的；都可保证迭 代最小二乘收敛，即可以实现二度体单一密度界面的稳健迭代反演。由于实际 资料不能完全分离出正常场，为了消除正常场的残余量对密度界面反演的影响， 我们又用重力异常水平梯度来反演密度界面。先取重力异常的水平梯度，然后 利用广义似然函数对模型参数进行优化选择，从而得到地下模型参数修f 量的 迭代计算公式，用l s q r 方法求解大型方程组，得到模型参数的修形量迭代计 算公式，再由界面起伏进行迭代计算，便可实现了二度体单一密度界面的反演。 经模型试算，表明此稳健迭代的反演方法在用重力异常反演密度界丽和用重力 异常水平梯度反演密度界面时都适用，并取得良好的结果。 第l 章重力场理论方法 第1 章重力场理论方法 重力学作为- j 3 古老、经典的地球科学一直是地球物理研究方法中不可或 缺的重要手段。自从牛顿发现力有引力以来，到今天重力学理论已经发展了数 百年，许多著名的数学物理学家都对地球重力学理论的发展做出了重要贡献。 拉普拉斯l a p l a c e 引进了位的概念并提出了位场理论的基本方程一拉普拉斯方程 和球谐分析方法，这些理论都是用重力方法研究地球形状的基础。欧拉( e u l e r ) 提出固体旋转理论，研究了当地球转动轴偏离主转动惯量时产生的晃动或自山 章动。斯托克( s t o k e s ) 研究大地水准面的形状，确定重力异常与大地水准面起 伏的关系，得到了斯托克斯定理，为大地测量学的发展提供了坚实的理论基础。 重力学的基础理论研究到现在已经发展的非常成熟了，作为今后的发展方向主 要是在方法技术和应用方面的研究( 王谦身等，2 0 0 3 ) 。 重力场是反应地球介质密度变化和各种环境( 固体地球潮汐、内部热流、 固体和液体之间的质量交换及迁移、表面负荷和地震构造运动等) 下动力学特 征的最基本和最直接的物理量。重力异常反演可以定性和定量地推断地下客观 存在的密度界面、异常地质结构和地质体的形状、产状及剩余密度分布，是人 们认识地球的重要手段之一。 1 1 地球物理反演方法发展概述 地球是人类赖以生存和发展的家园，以美国b a c k u s 和g i l b e r t 为首的一批地 球物理学家和应用数学家研究了地球物理数据和地球模型的共性，揭示了地球 物理反问题内在的非唯一性，并提出用算子谱分析构成反问题解估计的方法， 以及评价解估计的各项准则。 6 0 年代末至7 0 年代初，广义线性反演方法被广泛应用于地球振荡、地震波 与地球波速结构大地电磁和重磁反问题中。 7 0 年代发展起来的“以生存和发展的家园，为了不断丌发和利用地球，实 现人与自然的和谐发展就必须更加深刻的认识他。开始，人们用地震、地磁、 第1 章莺力场理论方法 地热、重力等方法去推测地球的内部结构和变化，推动了地球物理反演的不断 发展。 8 0 年代以来偏微分方程反问题研究进一步深化了反演的内涵，根据运动方 程和本构方程，大多数地球物理现象都可以用二阶双曲线型、椭圆型或抛物线 型偏微分方程的边值问题描述，因此地球物理反问题本质上是相应类型的偏微 分方程系数项或源项的反演。同时期以走时反演为基础的地震层析成像取得了 巨大的成功，由层析图揭示了波速的横向变化，由波速的横向变化了解地幔密 度和温度的横向变化，即低速异常与低密度、高温的地幔物质对应。由层析图 上地幔波速环状分布模式可推测地幔内热流的为位置，从而问接地显示出地幔 的对流图形。近年来，地球物理界已不满足于波速成像，丌始向密度各向异性 及粘度等更多参数成像发展，促进人们对深部岩性及流变性质的了解和相关地 球动力学问题的解决。 到上世纪末，非线性反演方法发展突飞猛进( k g a l l a g h e r & m s a m b r i d g e ，1 9 9 4 ) 除了一些传统的非线性方法( 也称非启发式) 方法，如梯度 法( g r a d i e n tm e t h o d ) ，牛顿法( n e w t o nm e t h o d ) 和蒙特卡洛法( m o n t eg a r l o m e t h o d ) 外，一些称为启发式的反演方法，如模拟退火( s i m u l a t e da n n e a l i n g ， s a ) 、遗传算法( g e n t i ca l g o r i t h m ，g a ) 、人工神经网络法( a n i 一c i a ln e u r a l n e t w o r k s ，a n n ) 小波分析法( w a v e l e t a n a l y s i s ，w a ) 等应运而尘，随着并行 计算机大出现，需要大量计算时间的非线性反演方法如鱼得水，有了迅速发展 的前提条件。近年来非线性反演方法初步显示了强大的生命力( s e n 等，1 9 9 2 ) 。 2 0 世纪8 0 年代后，联合反演方法得到了迅速的发展，尤其基于重力反演结 果存在多解性( 不唯一性) 的固有问题，应用重力和地震数据的联合反演已经 非常广泛了。杨辉等( 2 0 0 4 ) 利用相对准确的地震勘探结果作为分离重力场的 先验信息，然后用分离后的剩余场来反演地震反射不清楚的界面( 剥离法反演) 。 关小平( 1 9 9 4 ) 等研究了传统重震联合反演存在的问题，建议充分利用地震资 料作为形体参数进行场分离，对分离出来的目的曾为的重力效应再利用p a r k e r 公式进行反演( o l d e n b u r g ，1 9 7 4 ) ，以求出较深的或没有可靠地震资料的界而，在 此基础上利用速度和密度参数之间的关系，进行联合反演，取得了较好的效果。 a n d e r s o n ( 1 9 9 8 ) 等用顺序法对地震和重力资料联合反演速度和密度，为深度偏 移成像提供了准确的速度模型。 第l 章重力场理论方法 1 2 重力位场反演方法 重力位场反演对于确定地下地质体的物性参数及几何形念的研究起着重要 的作用，其方法得到广泛深入的研究( j i a n g h a ix i a 等，1 9 9 5 ；杨长福，2 0 0 4 ) 。 1 2 1 位场反演问题的定义 重力异常反演问题，简单地说就是，由实测重力异常及其导数的数值大小、 空间分布和变化规律，定性和定量推断地下客观存在的异常地质体结构、构造 和地质体的形状、产状及剩余密度分布等( 陈善，1 9 8 5 ) 。 任意形态物体g ( x ，j ，) 表达式可写为 g c x ，y ，= g a oj 0 己三萎 器s d 刀d f = g ，i i i ( q ) p ( q ) k ( p ，q v e d r l d f 其中：丌( s ，r ，f ) 为物体的形态函数，由物体边界确定，也可称物体边界位 置函数，该函数在物体内值为l ，在物体外为0 ，呈阶梯状。p ( e ，7 7 ，f ) 为物体密 度分布函数，它们共同确定了物体状态。 由上式，反演问题的数学定义为( 1 ) 由观测面上重力异常g ( x ，y ) 分布， 在给定物体边界位置函数的条件下，求解物体的密度分布函数：( 2 ) 由观测而 上重力异常g ( x ，) ，) 分布，在给定物体密度函数的条件下，求解物体的边界位置 函数；( 3 ) 由观测面上重力异常g ( x ，y ) 分布，在给定特殊约束( 如设物体密度 均匀、形态规则) 条件下，求解物体密度参数和几何参数。在这里，给定的函 数和特殊约束成为反演问题的定解条件。 1 2 2 重力异常反演方法概述 重力异常反演理论方法分为线性反演和非线性反演两大类。在给定定解条 件，使反问题具有唯一解的条件下，就反演定义和反演算法来说，密度分御反 演问题属于线性反演：而几何边界反演问题属于非线性反演；如果采用特殊约 束，把反演问题简化为有限个参数反演，此时的反演一般为非线性反演。 6 第l 章重力场理论方法 由上述分析可见，重力异常反演以非线性反演为主。非线性反演比较困难， 因此，研究者们经常由重力异常分布特征出发，采用各种简化算法，包括线性 化方法，把非线性反演变为简单反演，甚至变为线性反演。于是就出现了下面 各类别的反演方法。 依据参考文献的分类原则，把具体的重力异常反演方法分为如下四类。 ( 1 ) 解析反演方法： 此类方法用于求解均匀规则物体总质量( 或密度) 和几何参数反演问题。 其中包括特征点法、任意点法、切线法、参量图法4 种方法。 ( 2 ) 积分反演方法： 此类方法，既可用于求解均匀规则物体反演问题，又可用于求解任意密度 与横截面有限二度体和任意密度与形态三度体的总质量、质量重心位置。其中 包括空间域反演方法和波数域反演方法。这旱把波数域反演方法也列入积分法， 是因为正反演傅立叶变换是积分变换。当然，波数域反演方法也单独列为一类。 此类反演方法也只能用于异常数据位于直线或平面上这一简单情况。 选择法和最优化反演方法： 此类反演方法，可以用于异常数据位于起伏地形上和物体为任意形念与非 均匀密度等复杂情况。其中包括量板选择法、人机联作选择法、经典最优化选 择法和遗传算法等最优寻查方法。 全方位反演方法： 此类反演方法，可以在异常数据位于起伏地形上、物体为任意形态与非均 匀密度等复杂条件下应用。由于此方法采用“全空间解析延拓”关键技术，因 此它包含了传统的“场的空间分布法”。 重力异常反演理论和方法研究一直是重力异常解释的核心问题，贯穿重力 学发展的全过程，因此大量的反演理论和方法不断地涌现出来，使反演理论和 方法丰富多样。 1 2 3 反演问题的多解性和不稳定性 重力异常反演问题的多解性是地球物理反演中普遍存在的问题。产生多解 性的主要原因是不同场源体的等效性，即多个显然不同的地质体可以产z 卜完令 相同的重力场，如一个均质实心球体的在球体外部空间产生的重力场可以等效 于一个在球心处相同质量的质点产生的重力场。另外，观测和计算误差及观测 7 第1 章重力场理论方法 数据的有限性，也是造成反演结果多解性的原因。消除多解性的主要途径是综 合利用各种先验信息以增强约束条件，以及采用新的更合理的反演算法。 重力异常反演的另一个问题是解的不稳定性，即观测或计算误差在反演巾 被急剧方法，使得反演结果不收敛。解的不稳定性主要是由高频噪声和浅部= = 一- 扰引起的，所以在反演前要对重力异常进行位场分解。 8 第2 章重力异常稳健迭代密度界面反演研， 第2 章重力异常稳健迭代密度界面反演研究 2 。1 方法原理 利用重力异常进行密度界面反演，考虑到地一f 物质的分御极其复杂，为讨 论问题的方便，对地下界面模型做如下简化。如图2 i 所示，地下一二维密度界 面，界面深度f ( f ) 是水平坐标ft ：j o - 值函数，界面的平均深度为，上、下层 密度分别是q 和0 2 。由界面起伏所引起的重力异常为： a g ( x ，) = g ( q - c r 2 ) e l n ( x 手) 2 + ( z f ) 2 螈d 乡 ，- g ( q 训e l n 铲告 ( 2 1 ) f ( 舌) 图2 。l 二维层状界面模型 晃面深度f ( 步) 是孝的单值函数，界面平均深度为 ：、下层密度足qo 2 。 在z = o 时，同时设需要反演的界面段以外的两端水乎界面的深度为，从而可 得地表的重力异常： 第2 章重力异常稳健迭代密度界面反演f j f 究 g ( x ，。) = g ( q 一吒) 上：l n 蔓三挚孝 = g ( q 一吒) l int,+00 ( z 一孝) 2 + 9 2 ( 孝) ( x 一孝) 2 + | i z d 2 ( x 一孝) 2 + j f l d 2( x 一孝) 2 d = 2 万g ( 仃，一仃：) + g ( 仃。一吼) f ：l n 量专辫孝( 2 2 ) 议夼回刖_ 起伏农坯瓦刀 ( 亏) 3 。( 亏，+ a g 。( 亏) ，且 u ，力疗仄彳父f j 后炯夼 面起伏函数，可得： 。r n 、，r。、r l 。( x 一善) 2 + 2 ( 孝) + 2 ( 孝) ( 孝) + 2 ( 亭) 芦 蜮五o ) = 2 n g 缸一o - , ) h 0 + g ( c r t 吨) 1 n 竖业等篙搿4 严, 黜 m ，y 一，i 。 。 = 2 万g ( ( 7 i - - 0 2 ) + g ( q c r 2 ) e l n 鱼甚静孝+ +g(盯一)工：，n(t+三篁盖笔莹竺鬻)d孝 将右端第三项展开，舍去二次小量和高次项得： 蚺，o ) = 2 碳q 训嘶刮e i n 黹告+ 2 g ( m ) 器孝 对于实际重力反演问题，剖面长度取为有限段，设总长度为2 l ，则得： 删) _ 2 损q 刮懈训c i n 帮即g ( q 刮c 器f 由于求得幺( 孝) 后，积分c l n 等砉糌孝为已知值，可令： k f x 芒r1 ：2 g ( o - 一仃，、篁圭l ( 2 ) k ( x ，善，色) = - , 一仃z i 了二 粝 2 3 g。(x，。)=2万g(仃，一仃：)+g(仃。一口：)cln垒篆辫毒 ( 2 4 ) 1 0 第2 章重力异常稳健迭代密度界血反演t i j d 4 9 ( x , 0 ) = 每伍o ) + e k 专) 戋( 告嵫 ( 2 5 ) 在测得地下界面产生的地表重力异常五( x ) t 寿n t ，通过求解上式，可以求得 幺( 善) ，利用幺+ ，( 孝) = 厶( 孝) + 幺( 孝) 关系，可以实现地下界面的迭代反演计算。 2 2 模型参数的优化选择 ( 2 5 ) 式为第一类f r e d h o l m 积分方程，一般情况下，由于测量数据，( 石) 足 不精确的，并且选取的计算区间【- i j ，+ l 有限，其积分核可能是病态的。为此， 我们受到高尔根等( 2 0 0 0 ) 声波方程稳健迭代方法的启发，把误差函数表示为 6 ( x ) = 疋( x ) 一g ( x ，0 ) ，并且定义一个广义似然函数，使误差分布满足模型参数 的优化选择时，广义似然函数取极大值。这个误差函数可以通过广义幂指数在 概率分布簇中适当地表示为： m ，艄2 两p 丽p 州伽) 这里r 是g a m m a 函数，p 和p 都是正参数。上式中，当p = l 时，它简化为二项指 数分布：而当p = 2 时，它变为具有2 = 2 0 - 2 ( 仃2 是方差) 的普通高斯分钿，这 时r ( 1 2 ) = 玎；对于p 接近无穷时( 尸_ ) ，则变为均匀分布。它的优点在于： 它所描述的误差分布不仅仅是高斯分布，而是一个常见的重力异常干扰的多频 分量分布，它同时也提供了一个误差概率分布情况下，获得优化模型参数的准 则。在本方法中，这个广义似然函数具体为： p 芝l 正( t ) 一g ( o ) l 产缸 f ( a g ，p , p ) 2 而e x p 一生厂 上式中，i x 为重力异常沿地表x 轴离散采样总数，缸为离散点间距。优选反演 参数( 善) 过程，就是让上述似然函数极大化，它相当于指数取极小值。因此， 似然函数极大化等价于： 第2 章重力异常稳健迭代密度界面反演研究 兰l ( ) 一g ( z ，o ) 1 p 缸：m i n ( 2 6 ) 为了从上式中推得模型参数的变换量修正值幺( 孝) ，我们对( 2 5 ) 式采用梯形 法进行离散化，可得： a g ( x , ，o ) = ( 五，o ) 十芝k ( 薯，彭，幺) 幺( 彭) f + 去k ( 一，乒心，幺) ( 善一m ) 孝+ + 7 k ( x i ，芎扎，毛0 厶专，嗨n l 、) 专 = a g 。( 薯，o ) + 艺足( 葺，彭，色) ( 白) 孝 ( 2 7 这里= 一l ，- 2 ，- 1 ，0 ，1 ，2 ，n l 是地下界面的离散采样点编号，善为每 段的长度。当= 一n l 和n l 时，霞( 薯，乞，幺) = k ( 一，白，幺) 2 ；- 为其它值时， 霞( 葺，六，) = k ( ，乞，毛) 。将上式代入( 2 6 ) 式，并对( 2 6 ) 式作关于模型参量 厶六蟛) 的导数，再令其为零，同时令d ( 一) = 正( 薯) 一瓴( ，o ) ，可得： 挚墨g ：益：幺2 皇亟! ；挚竺! 二! ：兰：! ! ! 差。兰! 二：坌：量兰! ! 苎兰三 ( 2 8 ) 白r 2 - p 厶 p2 一” 式中k = 一m ，- 2 ，- 1 ，0 ，1 ，2 ，m 。并且由下式给定 巧= d ( 五) 一芝足( ，乞，) ( 乞) 厶f ( 2 9 我们取下列矩阵形式： 霞= 露( ，鼠，六) 孝) i , , i , k 肌- m 是核矩阵： w = d 垤扯i 卜2 ：是残差对角阵； 。= ( 彭) 2 一，是模型参数修正列向量； d = d ( i ) ) 竺是地表测量重力异常与参考界面起伏计算异常差。 1 2 第2 章重力异常稳健迭代解度界面反演研究 这样一来，( 2 7 ) 式可进一步写成下列矩阵形式： 露t w g , a ；。= 露t w d ( 2 1 0 ) 在上述方程中，权矩阵w 是一个依赖于误差分布的残差矩阵。对于p = 2 时， 亦就是当误差满足普通高斯分布时，w = i ， 变为一个具有等权的单位矩阵。 在这种情况下，上述方程( 2 1 0 ) 变为关于修正列向量的最小二乘解问题： 霞t 霞 ；。= 霞t d 。对于p = 1 情况，误差分布被假设为具有二项指数分布。这时， 按照残差绝对值的倒数，对权系数矩阵的主对角元素分配权值。残差绝对值越 大，分配的权值越小；反之，分配更大的权值。 通过求解( 2 1 0 ) 式，获得模型参数的修正量。实现迭代过程，实质上是 一个迭代加权最, i x - - 乘问题。b y r d 和p y n e ( 1 9 7 9 ) 研究表明：迭代加权最小二乘问 题，当范数尸( 1 p f l ，2 ，i r ( k l 占 g 是某个给定的小正数，k 是迭代次数。通过这样修改，这个算法完全收敛到某 个固定点i ；。原则上，p 可以取l 和2 之间的任意值，都可保证迭代加权最小 二乘法收敛，即可以实现稳健迭代反演。 2 3 数值模拟 为了验证本算法的有效性和计算精度，我们对模型数据进行了试算。以半 圆模型、正弦模型、梯形等地质地球物理模型对该方法作试验，考虑到反演计算 过程的稳定性，结合高尔根等( 2 0 0 0 ) 的结论，计算中选取p = 1 3 5 进行反演计 算。 第2 章重力异常稳健迭代密度界萄反演研究 模型一 中间有半圆型突起的水平密度界面模型( 如图2 2 所示) 。 地下有一密度界面，其二层介质的剩余密度) 勺- 5 0 0 k g m 3 j 界面的深度均为 l o 砌，水平方向总长度为1 0 0 砌，地下界面的离散采样间距为l 拥， 沿地表 x 轴离散采样间距为0 5 砌。( c ) 、( d ) 图分别是以起算深度为1 0 拥、5 砌进 行反演得到的结果两次计算深度与真值得均方误差分别为o 4 2 七m 和 0 3 5 砌。 图2 2 半圆突起界面的反演结果及重力异常对比 ( c ) 、( d ) 图中实线代表理论密度界面模型，虚线代表反演结果；( a ) 、( b ) 图巾殳线和虚线分 别代表理论密度界面和反演界面引起的重力异常。( a ) 、( c ) 图以起算深度为l o k m 进行反演：( b ) 、 ( d ) 图以起算深度) q s k m 进行反演 1 4 第2 章重力异常稳健迭代密度界面反演缈究 模型二弦波形起伏的密度分界面模型( 如图2 3 所示) 地下有一密度界面，其二层介质的剩余密度为一5 0 0 豫m 3 ，界面的深度均为 1 0 k i n ，水平方向总长度为l o o k m ，地下界面的离散采样间隔为l 砌，沿地表 x 轴离散采样间距为o 5 砌。( c ) 、( d ) 图分别是以起算深度为1 0 锄、5 拥进 行反演得到的结果两次计算深度与真值得均方误差分别为o 1 0 砌和 0 1 4 o n ，计算深度已很接近真值。 图2 3 正弦起伏界面的反演结果及重力异常对比 ( c ) ( d ) 图中实线代表理论密度界面模型，虚线代表反演结果；( a ) ( b ) 图巾实线和虚线分别 代表理论密度界面和反演界面引起的重力异常。( a ) ( c ) 图以起算深度为1 0 k r a 进行反演：( b ) ( d ) 图以起算深度为5 k m 进行反演。 1 5 第2 章重力异常稳健迭代密度界面反演研究 模型三中间有一梯形下凹的水平密度分界面模型( 如图2 4 所示) 地下有一密度界面，其二层介质的剩余密度为一5 0 0 k g i m 3 ，界面的深度均为 1 0 k r a ，水平方向总长度为1 0 0 k m ，地下界面的离散采样间隔为l 砌，沿地表 x 轴离散采样间距为0 5 k i n 。( c ) 、( d ) 图分别是以起算深度为l o k m 、5 砌进 行反演得到的结果。两次计算深度与真值得均方误差分别为o 4 0 砌和 0 4 2k m 。 图2 4 梯形凹陷界面的反演结果及重力异常对比 ( c ) ( d ) 图中实线代表理论密度界面模型。虚线代表反演结果：( a ) ( b ) 图巾实线和虚线分别代 表理论密度界面和反演界面引起的重力异常。( a ) ( c ) 图以起算深度为i o k m 进行反演：( b ) ( d ) 图以 起算深度为5 k i n 进行反演。 1 6 第2 章重力异常稳健迭代街度界面反演研究 2 4 关于正常场影晌的探讨 前面所推导的反演方法，对于较好分离出局部场情况下，可以直接用来进 行迭代反演计算，但是在实际的重力勘探资料处理中，由于野外测量的莺力异 常大都采用相对重力异常，对于一个异常体来说，由于采用的f 常场背景的影 响，其数值上经常会相差一个常数；另方面，我们在利用重力资料进行m o h o 面计算时，经常采用地表重力异常资料的向上延拓某个高度的结果作为m o h o 界 面产生的异常进： t m o h o 界面的计算。这样处理的结果，也会造成实际资料实际 上是m o h o 界面的异常和一个常数背景场的和。对于有地震测深资料约束情况卜， 采用压缩质面法( 刘元龙等，1 9 7 7 ) 等方法在计算m o h o 界面起伏形状时，常数场 可以在计算中被相互抵消，不影响计算结果。然而，对于勘探来说，并不能保 证每个计算区域都能得到很好的地震测深约束资料，为了克服因为讵常场的选 取所造成的台阶误差，进而造成反演结果的不精确问题，我们考虑把萨常场的 补偿修正也作为反演计算中的一个参数，进行综合反演计算，避免不必要的误 差影响。 首先，我们从( 2 5 ) 式出发，把( 2 5 ) 式写成下列形式： 姆( 为0 ) = 喙( o ) + l ，k ( 墨六乞) ( 倒孝+ c ( 2 11 ) 式中c 为常数，它是正常场在消除中的残余量。 我们按照推导( 2 1 0 ) 式的办法，最终可以推得下列矩阵方程： k r 形霞f n1 ：k d ( 2 1 2 ) 一c 夕 上式中i ：f fe 1 ，e 为元素为1 的。阶列向量，其它量的涵义与( 2 1o ) 式柏 同。通过求解( 2 1 2 ) 式，我们可以获得地下赛面的起伏形状和f 常背景场的残 余量的迭代反演结果。利用幺+ ( f ) = 六( 孝) + ( f ) 关系，可以实现地下界面的 迭代反演计算。 模型计算： 地下有一阶梯模型，其二层介质的剩余密度为5 0 0 k g 朋，界面的深度均为 1 0 枷，水平方向总长度为l o o k m ，地下界面的离散采样间隔为l 枷， 沿地表z 轴离散采样间距为0 5 砌，正常场的残余量c = o 4g a l ，来进行反演计算。 j 7 第2 章重力异常稳健迭代密度界面反演研究 下两图为加入c ，起算深度为5 k m 反演计算时，界面模型对比图和重力异常 对比图。 重力异常对比田( s k i n ) 图2 ，5 n 入c 时，阶梯模型重力异常反演对比图 ( 图为起算深度为5 k m 进行反演，理论密度界面( 加c 和米加c ) 和反演密度界面0 l 起的章力异常对比) 界面模型对比田( s k i n ) 图2 6 加入c 时，阶梯模型界面反演对比图 ( 图为起算深度为5 k m 进行反演，理论密度界面和反演密度界囱的对比) 18 第2 章重力异常稳健迭代密度界面反演研究 下两图为加入c 时，起算深度为l o k m 反演计算时，界面模型对比图和重力异 常对比图。 图2 7 加入c 时，阶梯模型重力异常反演对比图 ( 图为起算深度为l o k m 进行反演，理论带度界面( 加c 和朱力c ) 和反演街度界面，l 起的系力异常对比) 图2 8 加入c 时， ( 图为起算深度为l o k m 进行反演。 阶梯模型界面反演对比图 理论桁度群岛和反演密度界血的对比) - 1 9 第2 章重力异常稳健迭代密度界面反浈研究 下两图为不加入c 时，分别以l o k m 、5 k i n 为起算深度反演计算时，界面模型 对比图和重力异常对比图。 面 p 暑 x k m 图2 9 不加c 时，阶梯模型界面模型及重力异常对比图 ( 图为不加c 时。起算深度为i o k r a 进行反演，j ：圈为理论桁度界血和反演密度界面对比，下幽为理论密 度界面私反演密度界面弓l 起的重力异常对比) e 嵩 x k m 图2 1 0 不加c 时，界面模型及重力异常对比图 ( 图为不加c 时，起算深度为5 k i n 进行反演。上图为理论镪度界面和反演带度界卣对c = l 二，下鳖为理论俄度 界面和反演密度界面引起的重力异常对比) 2 0 第3 章重力异常水f 梯度密度界面反演彬f 究 第3 章重力异常水平梯度密度界面反演研究 第二章我们进行了重力异常密度界面反演，在考虑到实际重力场影响时， 我们加入了正常场残余量c 进行反演计算，可以消除正常场的影响。如果在重 力异常表达式中加入残余量c ，对重力异常求导，也能消除萨常场残余量对计算 结果的影响。因此，本章我们对重力异常取水平梯度，用上一章重力异常稳健 迭代密度界面的反演方法，实现了重力异常水平梯度密度界面反演，具体内容 如下。 3 1 方法原理 重力异常水平梯度密度界面反演与第二章的重力异常密度界面反演方法在 重力异常的推导部分是相同的，为简便起见，本章我们从公式( 2 5 ) 出发，对 重力异常取水平梯