张聪聪 赵 杰 张亚敏
(1.山西省地震局临汾中心地震台,山西 临汾 041000; 2.山西省地震局应急中心,山西 太原 030012)
模拟退火法是一种启发式的蒙特卡洛非线性反演方法[1],它模拟退火的物理过程:物质先加热融化,然后慢慢冷却,最后达到能量最小的晶体状态。模拟退火反演算法与传统的线性反演方法相比该方法具有:不依赖初始模型的选择、能寻找全局最小点而不陷入局部极小、在反演过程中不用计算雅克比偏导数矩阵等优点, 因而在地球物理资料非线性反演中受到广泛的应用。
某水库发生浅源地震,其坐标为(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 台站位置坐标及地震到时
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.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之间随机数。
取不同的初始温度,收敛阈值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 三种数据扰动方式反演结果
由表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),符合性较好,结果可以接受。
1)模拟退火法的使用需要了解先验模型信息和一定的实际工作经验,在经验缺乏的情况下,设定高的初始温度、大的降温系数(0.99)、小的收敛阈值。
2)由接受函数公式P=exp(-ΔE/T)可知,模拟退火算法初始温度太高时虽然可以在比较大的模型空间中搜索,较大概率接受较差的解,但是会增加迭代次数,增加计算时间。通过计算与试验可知,初始温度选取与ΔE同量级或者高一量级较为适宜,ΔE数量级可以通过在较高初始温度下画ΔE收敛图,看初始幅值获得。
3)温度阈值越小,精度越高,但是会增加迭代次数,增加计算时间。
4)扰动函数有多种方式,根据实际情况选择合适的扰动方法,必要时画出扰动因子收敛曲线图,判断全局性和收敛性的情况。