密度界面反演方法很多，主要有迭代法、直接法、统计分析法、非线性反演方法等.一直以来，由于非线性方法的计算效率问题，人们往往采用线性化方法或多次迭代的逐级线性优化方法进行反问题的求解.然而，随着计算机计算能力的提高和非线性反演方法算法的改进，非线性反演方法在实际地球物理问题的求解中逐渐应用起来.自适应模拟退火反演方法便是发展较好的方法之一，它具有不用求目标函数的偏导数及解大型矩阵方程组，即能找到一个全局最优解，而且易于加入约束条件的优点.众所周知，为了减少重、磁反演的多解性，利用带有约束信息约束界面反演结果，在减少反演的多解性、提高反演精度方面起到了积极的作用.而综合利用已有控制点资料作为约束信息，进行地球物理联合反演研究，可获得更为合理的结果，成为近年来的发展趋势.然而针对有些地区缺少甚至没有相应的约束信息的情况，本文提出了引入剩余密度模型的改进的自适应模拟退火密度界面反演方法.

地球物理学的反问题大多是一个多参数、非线性、多极值、非连续函数的优化问题.对于非线性反演问题，用局部方法反演的结果容易陷入局部极值解，而非全局最优解.所以，用全局寻优反演方法来求此类问题的最优解就显得越来越重要.在全局寻优反演方法中，最典型的算法为模拟退火方法(SA).模拟退火方法是通过模拟物质退火的物理过程，建立起来的一种非线性反演方法.其基本思想是：生成一系列参数向量模拟粒子的热运动，通过缓慢地减小一个模拟温度的控制参数，使模拟的系统最终冷却结晶达到系统能量最小值的过程.其实质是利用了地球物理反演问题求解过程与熔化固体退火过程的相似性，开辟了地球物理反演的新途径.这种方法的主要优点是：不用求目标函数的偏导数及解大型矩阵方程组，即能找到一个全局最优解，而且易于加入约束条件.本文在原有基础上引入了剩余密度模型，利用频率域变密度正演迭代的方法进行约束求解.针对因研究区范围大且地形复杂的地区的界面反演的不确定性的问题，而提出的一种改进的变剩余密度的自适应模拟退火密度界面反演方法.

重、磁数据的处理解释中，由于实测数据范围有限一般会导致处理后发生边部损失.为尽量减少这种现象，更充分地利用实测资料进行地质解释，通常需对实测数据进行扩边，以满足后续处理的需要.

重、磁异常场中区域异常一般展布范围较大、变化相对简单，而局部异常展布范围小、变化较复杂.根据重、磁场的这种特点，段本春曾提出一种区域场扩边方法，该方法的基本思想是只对区域场进行扩边计算.此方法不仅简便易行，而且在压制边界效应方面效果好，对一般特点变化的异常普适性强. 本文设有网格化后的异常数据 ZY(1：m′，1：n′).需将数据扩边为ZK(1-m′：m+ m′，1-n′：n+n′)，使得ZK(i，j)= ZY(i，j)，(1≤i≤m，1≤j≤n).

区域场四个边的扩边计算公式分别为

以本文实例而言，用于反演的格莱尼重力异常数据网度为20 km×20 km，先确定每条测线上两侧各扩166.7 km即可，应该是两侧各扩8个点，并两侧各扩8条线.利用各测线边部9个实测数据求出相邻两点间异常差的平均值，这个平均值可以作为与边部9实测数据相应的区域场的变化梯度，然后利用这个区域场变化梯度作为增减因子向外扩 166.7 km即可.针对重力异常反演的现状，采用矩形单元体平面剖分组合模型，实现重力正演.重力正演迭代的过程是本方法中关键的一步，采用方法为频率域变密度正演方法.

结合已知地质地球物理资料构建剩余密度文件后，重力异常反演的目标函数为

该建模方法适应性强，适合于实际工作中常见的密度横向变化的复杂模型.使用改进的全局寻优的快速模拟退火算法，对重力异常进行反演，结合这种灵活的密度建模方法，反演过程中只需要反演各矩形单元体的密度参数，即可同时得到地质体的界面或形态以及密度值分布.在地震等先验信息约束下，该重力反演方法提高了反演精度并减少了多解性，可有效解决界面反演不确定性的难题.

整个反演方法按以下步骤进行：

(1)给定模型每一参数变化范围，包括剩余密度初始模型ρ₀，并在这个范围内随机选择一个初始模型m₀，并计算相应的目标函数值E(m₀)；

(2)对当前模型m0、ρ0进行扰动产生一个新模型m、ρ，利用频率域变密度正演计算相应的目标函数值E(m)，得到ΔE= E(m)-E(m₀)；

(3)若ΔE＜0，则新模型m、ρ被接受，若ΔE>0，则新模型m按概率P=exp(-ΔE/T)进行接受，T为温度.当模型被接受时，置m=m0，E(m0)=E(m)；

(4)在温度T下，重复一定次数的扰动和接受过程，即重复步骤(2)，(3)；

(5)缓慢地降低温度T；

(6)重复步骤(2)，(5)，直至收敛条件满足为止.

具体流程如图 1所示.

图 1

Fig. 1

图 1 改进的自适应模拟退火密度界面反演流程图 Fig. 1 Flow chart of modified simulated annealing density interface inversion method

2.1 模型实例

为了进一步说明变密度自适应模拟退火三维界面反演方法的适用性，本文通过一个模型实例来进行说明.图 2是本文建立的深度模型实例.结合构建的剩余密度模型(图 3)，利用频率域变密度正演得到重力异常场(图 4).然后利用Parker反演方法和改进的变密度自适应模拟退火反演方法所得到的反演结果(图 5)及其与原始深度模型差值(图 6).

图 2

Fig. 2

图 1 深度模型实例 Fig. 1 The example of depth model

图 3

Fig. 3

图 3 建立初始剩余密度模型 Fig. 3 The foundation of density difference model

图 4

Fig. 4

图 4 模型重力异常正演结果 Fig. 4 Gravity anomaly forward results of model

图 5

Fig. 5

图 5 Parker界面反演方法反演结果(a)及其与原始深度模型差值(b). Fig. 5 The inversion results of Parker inversion method(a) and the difference between the inversion results and the example of depth model(b).

图 6

Fig. 6

图 6 改进的自适应模拟退火密度界面反演方法 反演结果(a)及其与原始深度模型差值(b). Fig. 6 The inversion results of modified simulated annealing density interface inversion method(a) and the difference between inversion results and the example of depth model(b).

通过模型两种反演结果与原始深度模型差值对比，体现本文所选用的变密度自适应模拟退火反演方法的优点及其在构造单元多、密度变化复杂的地区应用的可行性.

本文将其应用于印度-孟加拉湾研究区，采用的原始重力异常数据为格莱尼重力异常，即对布格重力异常数据进行了格莱尼改正.格莱尼改正值的与地势密切相关，从青藏高原到印度洋改正值降幅达140 mGal.它隐含于布格重力异常之中成为布格异常“区域背景”的一部分.本文研究区实例的跨度相对较大，为了消除格莱尼改正的影响，采用格莱尼重力异常(图 7)代替布格重力异常计算莫霍面深度.结合前人结果，采用向上延拓法、补偿圆滑滤波法及正则化滤波法进行对比得到研究区的区域异常场(图 8).

图 7

Fig. 7

图 7 印度-孟加拉湾研究区格莱尼异常 Fig. 7 The Grannie gravity anomaly in the study area of India-the bay of Bengal

图 8

Fig. 8

图 8 印度-孟加拉湾研究区区域重力异常 Fig. 8 The regional gravity anomaly in the study area of India-the bay of Bengal

在已有地质、地球物理资料比较充分的地区，尤其是有地震测深结果约束的地区，我们可以采用带控制点的三维界面反演方法，可以得到精度较高的反演结果.而对于已有约束点较少或者没有约束点的地区，采用单一密度界面反演的Parker法(图 10)则往往使结果偏差增大，从而影响分区拼接和整体反演精度.本文通过构建合理的剩余密度模型，进行变密度自适应模拟退火界面反演结果(图 11)，并得到了其与Parker法反演结果的差值(图 12).为了进一步说明变密度自适应模拟退火反演方法效果，并将自适应模拟退火反演方法与Parker法的反演结果与穿越喜马拉雅山构造带宽频地震测深剖面结果进行了对比，其所在位置见图 10~12中红色实线标注.

图 9

Fig. 9

图 9 建立初始剩余密度模型 Fig. 9 The foundation of density difference model

图 10

Fig. 10

图 10 Parker法反演结果(σ=0.48 g/cm3，平均深度38 km) Fig. 10 The inversion results of Parker inversion method(σ=0.48 g/cm3，average depth 38km)

图 11

Fig. 11

图 11 改进的自适应模拟退火密度界面反演方法反演结果 (σ1=0.56 g/cm3，σ2=0.43 g/cm3，σ3=0.4 g/cm3，平均深度38 km) Fig. 11 The inversion results of modified simulated annealing density interface inversion method (σ1=0.56 g/cm3，σ2=0.43 g/cm3，σ3=0.4 g/cm3，average depth 38 km)

图 12

Fig. 12

图 12 两种方法反演结果差值 Fig. 12 The difference between inversion results of Parker inversion method and modified simulated annealing inversion method

印度-孟加拉湾研究区包括孟加拉湾、印度次大陆和喜马拉雅造山带，这三个次分区莫霍面深度变化极大.孟加拉湾地区莫霍面深度约为20~30 km，印度次大陆地区莫霍面深度约为36~42 km，而喜马拉雅造山带地区莫霍面深度约为40~60 km.此研究区内三个次分区的地壳性质分别属于减薄型大陆地壳、正常大陆地壳和增厚型大陆地壳.莫霍面深度相差极大，如果采用单一密度界面方法进行反演，往往顾此失彼.例如，本文采用Parker法得到的结果，孟加拉湾及印度地区的偏差绝对值在3~6 km.通过对两种方法反演结果差值的对比，我们可以发现，自适应模拟退火反演方法极大的提高了反演精度，尤其是分区内剩余密度相差较大的地区.为了更清楚的对比两种方法的反演结果，我们将二者反演结果与实际地震测深结果进行了偏差分析对比，如图 13所示.不难发现变密度的自适应模拟退火反演方法具有更好的效果.本文通过模型与实际应用实例均体现了改进的变密度自适应模拟退火反演方法的优点和可行性.

图 13

Fig. 13

图 13 两种方法反演结果与地震测深剖面结果对比图 Fig. 13 The difference among the two inversion results and seismic sounding profile

3 结 论

非线性算法的反演将是今后发展的一个主流方向，本文结合已有地质地球物理基础，通过引入变剩余密度建模来改进模拟退火算法，能加速收敛速度，使得反演结果更好地适合了地球物理场.通过模型与应用实例进一步表明，在密度界面反演中，针对范围广而复杂的地球物理场，在缺乏相应的数据或无约束信息的区域内采用变剩余密度的模拟退火反演方法是可行的.