张聪聪 赵 杰 张亚敏

(1.山西省地震局临汾中心地震台，山西 临汾 041000； 2.山西省地震局应急中心，山西 太原 030012)

1 概述

模拟退火法是一种启发式的蒙特卡洛非线性反演方法[1]，它模拟退火的物理过程：物质先加热融化，然后慢慢冷却，最后达到能量最小的晶体状态。模拟退火反演算法与传统的线性反演方法相比该方法具有:不依赖初始模型的选择、能寻找全局最小点而不陷入局部极小、在反演过程中不用计算雅克比偏导数矩阵等优点, 因而在地球物理资料非线性反演中受到广泛的应用。

2 问题描述

某水库发生浅源地震，其坐标为(x，y，h)，在其附近有9个观测台站，其坐标为(xk，yk，hk)观测到了直达波P的走时tpk，地面水平，地球介质均匀，发震时刻为t，P波速度vp，则存在下面的方程：

现在我们需要求出该地震发震位置与发震时刻(x，y，h，t)，其中9个台站坐标与地震到时数据[2]见表1，假定vp=2 000 m/s。

表1 台站位置坐标及地震到时

3 模拟退火法计算(MSA)

3.1 模拟退火法步骤[3]

1)给定模型每一个参数变化范围，在这个范围内随机选择一个初始模型m0，并计算相应的目标函数值E(m0)；

2)对当前的模型进行扰动产生一个新模型m，计算相应的目标函数值E(m),得到ΔE=E(m)-E(m0)；

3)若ΔE<0,则新模型被接受，若ΔE>0，则新模型m按概率P=exp(-ΔE/kT)进行接受，其中,k为一个常数;exp为自然指数;T为温度。因为P(ΔE)的函数取值范围是(0,1)，在0～1之间随机产生一个数R，当P(ΔE)>R时，接受修改，即m0=m；

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

5)缓慢降低温度T，Tnew=(0.99)k×T；

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

3.2 重要参数

3.2.1温度及下降速率

温度T在模拟退火反演中是最重要的问题，由接受概率公式可知P=exp(-ΔE/kT)，温度实际上是一种非线性加权[3]。初始温度T太高会增加计算时间；T也不能取太小，否则不能遍历模型空间。收敛阈值Tm不能取太大，这样会达不到收敛状态；Tm也不能太小，这样会增加计算量。下降速率如果太快，虽然运算速度快，但容易陷入局部最优解；下降速率太慢，则运行效率较低。

本文选择双曲线下降行作为退火机制，且初始选择较高初始温度100，选择较低收敛温度0.000 01，下降速率0.99。通过matlab选择较好的参数。

3.2.2扰动函数的选择[4]

本文主要分析三种模型扰动方法：

1)全局扰动:

X′=xmin+(xmax-xmin) ×R。

其中，X′为模型新值;R为0～1之间的随机数。

2)Ingher提出的VFSA方法:

X′=x+(xmax-xmin) ×yi；

yi=T×sign(R-0.5)×[(1+1/Tk)|2R-1|-1]。

其中，yi为扰动因子;T为温度;X′为模型新值;x为模型旧值;R为0～1之间的随机数。

3)和温度相关的局部收敛加强型:

u=T×tan[pi×(R-0.5)];

X′=u×(xmax-xmin)+x。

其中，u为扰动因子;T为温度;X′为模型新值;x为模型旧值;R为0～1之间随机数。

4 重要参数选择

4.1 温度

取不同的初始温度，收敛阈值Tm=0.000 001，下降速率k=0.99，得到的地震定位见表2，图1为不同起始温度到时残差。

表2 不同初始温度对反演结果的影响

由表2，图1可知不同起始温度对数据观测值影响不大，会影响到迭代次数，根据接受概率公式P=exp(-ΔE/T)可知初始温度T的选择应该与ΔE同一数量级或者高一两个数量级为宜，这样才能更好的实现全局搜索，且有较高工作效率。T初始值选取，可以通过查看T取较高值时ΔE的曲线图的纵坐标可得。如图2所示，迭代次数取1 000左右较为合适，即T取1或者10较适宜。为了提高搜索效率，本次采用T=1作为初始温度。

表3 不同的收敛温度对地震反演结果的影响

T=1k=0.99迭代次数发震位置与发震时刻(x,y,h,t)xyhtTm=0.01460Tm=0.001689Tm=0.000 00191817.7734.92-200.370.474 69.2727.28-208.520.470 825.3615.67-210.240.473 541.0339.12-200.540.474 537.2437.52-203.790.472 939.2138.56-200.890.473 039.6241.62-199.980.473 739.6641.44-200.000.473 739.5641.54-199.990.473 7

由表3，图3可知，当初始温度选择T=100，下降速率为0.99时，随着温度阈值的减小，迭代次数的增加，反演结果越来越稳定，且到时残差越来越小，精度越来越高。此次实验为了较高的反演精度选择Tm=0.000 000 1。为了更好的全局收敛，选择较低的降温系数。本次实验选择T=1，k=0.99，Tm=0.000 000 1 。

4.2 三种扰动方式分析

表4 三种数据扰动方式反演结果

由表4可知，全局扰动方法收敛性较差，其余两种方法收敛性较好。图4,图5为具有较强的局部收敛性扰动函数的扰动因子收敛图，温度为T=100到Tm=0.000 000 1(下降速率0.99)，从图4可以看出该扰动因子可以较好的收敛，且扰动范围大，迭代次数少，效率较高。从图5可以看出该扰动方法可以收敛，但是需要较多的迭代次数，且扰动范围较小。虽然第二种方式扰动较大，但根据公式可知很多扰动为无效扰动，第三种方式基本全部为有效扰动。综上所述：

全局性：method 1>method 3>method 2;

收敛速度：method 2>method 3>method 1;

综合考虑，选择收敛性较强的method 2。

根据试验，采用T=1，k=0.99，Tm=0.000 000 1，收敛函数method 2对地震进行定位，然后正演获取台站地震到时，与实际地震到时观测值作对比(见图6)，符合性较好，结果可以接受。

5 结论与建议

1)模拟退火法的使用需要了解先验模型信息和一定的实际工作经验，在经验缺乏的情况下，设定高的初始温度、大的降温系数(0.99)、小的收敛阈值。

2)由接受函数公式P=exp(-ΔE/T)可知，模拟退火算法初始温度太高时虽然可以在比较大的模型空间中搜索，较大概率接受较差的解，但是会增加迭代次数，增加计算时间。通过计算与试验可知，初始温度选取与ΔE同量级或者高一量级较为适宜，ΔE数量级可以通过在较高初始温度下画ΔE收敛图，看初始幅值获得。

3)温度阈值越小，精度越高，但是会增加迭代次数，增加计算时间。

4)扰动函数有多种方式，根据实际情况选择合适的扰动方法，必要时画出扰动因子收敛曲线图，判断全局性和收敛性的情况。