直接反演多层密度界面的方法研究*

王 笋 申重阳

1)中国地震局地震研究所(地震大地测量重点实验室)，武汉 430071)

2)厦门地震勘测研究中心，厦门 361021

3)中国地震局地壳应力研究所武汉创新基地，武汉430071

1 引言

利用重力异常资料反演地下介质密度差异明显的分界面起伏形态称为界面反演，对盆地基底研究、区域构造和深部构造起伏研究都具有重要作用［1］。实际应用中往往需要对地下的多个密度界面进行反演，由于重力异常叠加在一起，导致反演的难度很大。目前反演多层密度界面主要用场分离方法，对重力异常的空间频谱特征进行分析，使特定频率成分对应于单个界面或深度范围，把多层密度界面的重力场逐层分离，将剥离后的剩余场按单层界面的方法逐层反演。但场的分离很难做到与各个界面完全对应。

将地下场源区域规则划分成很多多边形单元，通过反演确定这些单元的密度变化来约束异常形态，其反演条件自由度高，易操作，模拟的地质模型范围更广［2］。但如果想精细地描述场源的物性分布，就要增加单元数量，这样必然扩大解的空间范围［3］，而且使计算量急剧增加和出现无意义的解。虽然计算机技术的发展解决了前一个问题，但后一个问题只能通过改进反演方法来解决。

本文选用经典的二度体模型模拟各层密度界面，用地震等资料给出各层界面的深度控制，运用参数范围约束和光滑约束压缩解空间，并使用初始模型要求不高的非线性全局优化方法——遗传算法进行多层密度界面反演计算。

2 算法原理

2.1 模型的建立和参数化

考虑二维情形，若地下空间被呈层状分布的数种介质充填，不同介质具有不同的密度(图1)，其间的分界面的起伏将产生附近空间重力位场的变化。根据位场的叠加原理，地面观测到的重力扰动等于各层分界面产生的重力异常之和;要求取单个界面产生的重力异常，可以其固定的上下界(图1 上红色虚线)将附近剖面划分出来，计算此部分剖面场源产生的重力效应;为简化计算，记界面两侧密度差值为σ，将界面一侧视为密度为0 的空间，另一侧用密度σ 充填。

图1 地下密度界面的等效场源Fig.1 Equivalent mass distribution of density interface

如图2，把界面起伏视作相对深度D 的一系列水平长方柱体高度Δh 的变化，用一维数组dp0［N］来存放离散化的柱体高度值，以其作为初始模型模型参数。确定这些柱体的剩余密度值σ 后，再根据二度矩形截面体公式计算地面测点处的重力效应。

初始模型参数发生扰动，生成演进模型，随之发生的重力场变化等于以初始界面为底边、演进界面为顶边的所有长方柱体产生的引力之和(图3)。

用一维数组dpa［N］存放离散化的柱体高度值即演进界面偏离初始界面的程度，以其作为反演对象参数，dpa［N］的取值范围即界面深度约束范围，该约束值可以根据钻孔、地震、电法及地质等已知资料获得。

2.2 目标函数和光滑约束

用于反演的目标函数包含由模型参数计算得到的重力场值与观测数据的误差:考虑M 个不规则分布的重力台站的观测数据，对于第i 个重力台站Pi(xi，zi)，设重力台站点Pi的平面坐标、重力异常值和重力计算值分别为Pi(xi，zi)、gi，obs(i=1，…，M)和gi，cal(i=1，…，M)，则其均方差为:

在没有约束的情况下，重力反演是病态的。为避免出现一些明显不合理的结果，引入紧邻参数估值相似，使被反演的众多参数中空间相邻参数之间的紧密联系、数值上尽可能地接近。由于这种约束倾向于产生物性光滑过渡的场源，所以这种约束又称为光滑约束［4］。其形式为:

为使反演问题稳定，采用Tikhonov 正则化方法［5］，使包括两项的、关于模型拟合差和光滑约束项的目标函数最小，这里选择的目标函数为:

式中b 为正则化参数。

2.3 遗传算法

遗传算法是一种以概率论为基础的求解多极值非线性问题的最优化方法。

应用遗传算法进行界面反演，关键在于如何对模型参数进行编码。以往使用的二进制编码表示一个参数就需要几位，要求的分辨率越高，编码越长，当参数众多时，使用起来很不方便，故不能适应较多变量的反演计算［6］。为适应物性模型场源精细划分的需要，采用十进制整数编码方式［7］。

参数编码采用十进制整数编码，即把每个参数的编码用一个整数来表示。N 个参数的L 组模型共有N×L 个编码。参数群体中的一个编码代表一个基因，表示一个模型的N 个基因构成一条染色体，群体中共有L 条染色体。对每条染色体译码，然后计算出这个模型的目标函数。首先计算初始群体的L 个目标函数，以目标函数作为判据进行搜索。

3 理论模型反演

给定一个试验模型:三层密度界面，密度差分别为0.08、0.11 和0.43，每个界面各有一个斜坡或凹槽(图4)。应用非均匀分布的223 个观测点，其观测精度相同。不添加任何区域特征的正演计算得到的重力异常(图5)其特征为鞍形隆起。为验证反演算法对误差的敏感度，对异常加入8%随机误差。

初始模型设置为3 层水平界面，深度取平均深度。应用本文所用方法搜索解空间，求得拟合差达到设计要求的解，以确定密度界面的形态(图6)。

由图6 可见，从只有界面平均深度信息的初始模型出发，通过在光滑约束下对重力场特征的拟合，反演结果具有了和验证模型相同的特征，结构基本一致，说明本反演程序分离重力异常的效果是可靠的，添加光滑度约束对多层密度界面进行同时反演的方法是可行的。

4 实际模型反演

类乌齐-玉树-玛多剖面位于青藏高原东部，剖面所在区域存在4 个主要密度界面:沉积层底界，上地壳底界，中地壳底界，下地壳底界(Moho 面)①王夫运，等.青海玉树7.1级地震科学考察项目研究报告［R］.2010.。利用钻井和地震资料可探明结晶基底面的形态，而上、中、下地壳的底界在地震资料上不够清晰，需要利用重力资料［8］进行计算。若用频率域场分离的方法反演结晶基底面下方的3 个密度界面，最浅的上地壳底界的重力异常只保留高频成分，长波长变化信息丧失，界面只在平均深度附近跳动，与地震资料给出的此剖面上地壳底界埋深变化较平缓，总体呈南高北低的趋势的结果偏离较大。本文根据P波速度-深度剖面，以及P 波速度和密度的经验公式，建立岩石圈初始界面-密度模型，确定界面深度变化范围(图7)。先将沉积层的异常剥离，再运用前述方法同时反演这3 个密度界面的形态，重力异常拟合效果如图8，均方差约为2 ×10－5ms－2。输出模型如图9 所示，Moho 面形态与文献［7，8］得出的结果基本相同，并展现了上、中、下地壳分界面的起伏变化。

图9 类乌齐-玉树-玛多密度界面模型Fig.9 Final density model of Leiwuqi-Yushu-Maduo profile

5 结论

位场异常的分离不可简单视为某个频率段对应于特定深度。一个可行的办法是对模型使用尽量多的已知资料作控制约束，用钻井、地震及地质资料等给出深度控制，添加光滑约束消除包含虚假高频成分的解，达到压缩解空间的目的，再运用全局优化方法直接反演多层密度界面。

以往用场分离的方法反演多层密度界面，是对重力异常进行滤波，把通过高通滤波器的重力异常归为浅处的界面的，忽略了较浅界面可能的长波长变化(如整体向一侧倾的)，这与很多区域的地下情况不符，无法得到理想的反演结果;本文提出应用光滑约束，其效果与对各层密度界面进行低通滤波类似，同样使密度界面与一定频率范围的重力异常对应，但保留了浅部界面的低频异常，与实际地下情形符合得更好。另外，观测到的重力场值不可避免地混有误差和干扰，混染重力异常的频谱，严重影响场分离的效果;本文所述方法从界面的形态特征出发，有对重力异常的干扰成分不敏感的优点。

1 陈军，等.应用改进的遗传算法反演多层密度界面［J］.地球科学(中国地质大学学报)，2000，25(6):651－655.(Chen Jun，et al.Application of improved genetic algorithm to inversion of multi-layer density interface［J］.Earth Science-Journal of China University of Geosciences，2000，25(6):651－655)

2 姚长利，等.重磁遗传算法三维反演中高速计算及有效存储方法技术［J］.地球物理学报，2003，46(2):252－258.(Yao Changli，et al.High-speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms［J］.Chinese Journal of Geophysics，2003，46(2):252－258)

3 Chen Shi，Zhang Jian and Shi Yaolin.Gravity inversion using the frequency characteristics of the density distribution［J］.Applied Geophysics，2008，5(2):99－106.

4 姚长利，郝天珧，管志宁.重磁反演约束条件及三维物性反演技术策略［J］.物探与化探，2002，26(4):253－257.(Yao Changli，Hao Tianyao and Guan Zhining.Restrictions in gravity and magnetic inversions and technical strategy of 3D properties inversion［J］.Geophysical ＆ Geochemical Exploration，2002，26(4):253－257)

5 金其年.Tikhonov 正则化的饱和性与逆结果及后验参数选取［J］.中国科学(A 辑)，1999，29(8):715－723.(Jin Qinian.Saturation and inverse-results of Tikhonov regularization and selection of posteriori parameter［J］.Science In China(Series A)，1999，29(8):715－723)

6 杨光亮，等.类乌齐-玉树-玛多剖面重力异常研究［J］.大地测量与地球动力学，2011，(5):1－4.(Yang Guangliang，et al.Study on gravity anomaly of profile riwoqe-yushumaduo［J］.Journal of Geodesy and Geodynamics，2011，(5):1－4)

7 Yang Guangliang，et al.Joint inversion of gravity and seimic data along a profile across the seimogenic fault of 2010 Yushu Ms7.1 earthquake［J］.Geodesy and Geodynamics，2011，2(4):21－27.Doi:10.3724/SP.J.1246.2011.00021.

8 王谦身，安玉林.青藏高原东部玛多-沙马地区的重力场与深部构造［J］.地球物理学进展，2001，16(4):4－10.(Wang Qianshen and An Yulin.Gravity field and deep structure of Maduo-Shama region in eastern Qinghai-Xizang(Tibetan)Plateau［J］.Progress of Geophysics，2001，16(4):4－10.)