模拟退火算法原理及应用

模拟退火算法原理及应用研究主讲:陈华根同济大学海洋与地球科学学院一模拟退火算法及VFSA算法模拟退火算法在反演中的应用：•非线性组合优化算法：模型扰动，模拟退火，全局寻优。•能量函数—目标函数•模拟退火过程—反演迭代传统模拟退火流程图YesNo随机选择初始模型m0计算能量函数E(m0)模型扰动产生新模型m1=m0+△m0计算能量函数E(m1)△E=E(m1)-E(m0)△E≤0?m0=m1新模型按Metropolis准则接受缓慢降低温度满足收敛条件为止Metropolis接受准则：•⊿E0，新模型无条件被接受—接受能量值较小状态•否则，)Eexp(Tkrb•产生随机数ξ∈[0，1]•若rξ，新模型被接受，否则被舍弃。—接受能量值较大状态，从而在模拟退火温度控制下全局寻优。VFSA算法分析：•模型扰动：831/115.0sgn7312'－－uiiiiiiTuTyABymm•接收概率：93/111/1－hTEhP•退火计划：103/10－NckExpTkTVFSA的温度衰减曲线：VFSA的降温速度是比较快的1020304050607080901001020304050607080901000.980.80.50.2衰减率系数递增方向初始温度：100迭代次数：100参数个数N:2温度迭代次数高温下VFSA算法模型状态分布图：高温下VFSA算法的状态空间遍历能力逊于随机数发生器的遍历能力-3-2.5-2-1.5-1-0.500.511.522.53-3-2.5-2-1.5-1-0.500.511.522.53xyVFSA方式扰动全局随机扰动VFSA算法迭代次数k与系数yi的关系示意图：低温下模型扰动的空间过大，扰动后模型被接受的机率必然降低，势必影响寻优效率，最终影响算法完成后最终解的精度-1-0.500.51yi高温带过渡带低温带010203040506070809010005101520迭代次数退火温度VFSA二改进的VFSA算法—MVFSA算法MVFSA有以下改进：•过程一：较高的初始温度，VFSA算法的退火计划，模型作全局随机扰动—搜索并锁定最优解区间；•过程二：较低的初始温度，适当回火的退火计划，模型作局部随机扰动--扰动在当前模型周围进行—在锁定最优解空间后，由于其搜索空间变得较小，以此提高模型接受效率。过程一：模拟退火，全局搜索iiiiABuBm'T=T0*EXP(-α*(j-1)1/2)过程二：回火升温，局部搜索)(/)5.0('jLABummiiiiT=T0*EXP(-α*(j-k0/β)1/2)图3-5VFSA与MVFSA算法的退火温度曲线比较0迭代次数退火温度VFSA改进算法过程一过程二MVFSA算法迭代次数k与系数yi的关系示意图-1-0.500.51yi高温带过渡带低温带010203040506070809010005101520迭代次数退火温度MVFSA模型试验—Z=f(x,y)型：模拟退火计划表表3-1算法初试温度温度衰减率叠代次数扰动次数初始位置VFSA2000.9985003x0=2.5,y0=2.5MVFSA2000.9985003x0=2.5,y0=2.5VFSA算法扰动状态分布和寻优轨迹图-3-2.5-2-1.5-1-0.500.511.522.53-3-2.5-2-1.5-1-0.500.511.522.53起点最优解xy寻优轨迹接受状态MVFSA算法扰动状态分布和寻优轨迹图-3-2.5-2-1.5-1-0.500.511.522.53-3-2.5-2-1.5-1-0.500.511.522.53起点最优解xy寻优轨迹接受状态VFSA算法目标函数之差与迭代次数关系图501001502002503003504004505000510目标函数之差迭代次数扰动状态接受状态寻优轨迹204060800510局部放大MVFSA算法目标函数之差与迭代次数关系图501001502002503003504004505000510目标函数之差迭代次数扰动状态接受状态寻优轨迹204060800510局部放大•VFSA及MVFSA算法在退火计划十分完备的情况下,表现相当完美：算法起点相同，寻优路径不同，最终找到的都是同一最优解•VFSA与MVFSA算法的模型状态均分布这个状态空间，但VFSA模型状态在最优解点出现一个十字型状态，MVFSA算法在整个最优解区域形成一个矩形，这与它们的模型扰动方式有关。•在相同的退火计划下两种算法的时间，VFSA算法约为103秒，而MVFSA算法只用时约75秒，多次试验表明：MVFSA算法计算时间约比VFSA算法少20-30%。算法稳健性试验：模拟退火计划表表3-3算法初试温度温度衰减率叠代次数扰动次数初始位置VFSA20.98301x0=2.5,y0=2.5MVFSA20.98301x0=2.5,y0=2.5VFSA算法扰动状态分布和寻优轨迹图-3-2.5-2-1.5-1-0.500.511.522.53-3-2.5-2-1.5-1-0.500.511.522.53起点xy寻优轨迹接受状态扰动状态最终解最优解位置MVFSA算法扰动状态分布和寻优轨迹图-3-2.5-2-1.5-1-0.500.511.522.53-3-2.5-2-1.5-1-0.500.511.522.53起点最优解xy寻优轨迹接受状态扰动状态稳健性试验结论：•多次试验表明：在同等退火计划下，VFSA算法较易落入了局部极值区，而MVSFA算法则比较稳健。应用•一从简单入手-重力模拟退火反演研究MT-重力联合反演研究重力正演计算：计算单元：2.5度体的多边形截面棱柱体Δg(x,y,z)=G·б·∑cosψi[F1(y2-y,i)+F1(y1-y,i)]重力目标函数：MiobsicaligMgxg12/1目标函数的含义：正演值与实测值的相对均方误差。优点：无量纲，并与测点数无关，便于与MT方法共同开展联合反演重力反演的等值效应现象km-5-4-3-2-10012345678910161820222426△g(mGal)km实际数据反演结果实际界面反演界面等效密度界面2.70(g/cm)32.30(g/cm)3消除锯齿状界面的方法：•人机联作方式•修正在反演程序执行过程中进行，不需暂停反演程序9/)232(2112'jjjjjiTTTTTT)2/()(11'nTnTTTjjji垂直侧边梯形的MVFSA反演结果ABCD-2-1.5-1-0.500102030405060708090100012345678MVFSAkmkm△g(mGal)实际模型反演结果模型计算反演结果密度差0.28(g/cm)3垂直侧边梯形组合的MVFSA反演结果ABCDEFGHIJ-2-10010203040506070809010002468模型计算反演结果kmkm△g(mGal)实际模型反演结果MVFSA密度差0.28(g/cm)3背斜-向斜模型的模拟退火反演结果-7-6-5-4-3-2-1002468101214161820246模型计算反演结果模型界面反演结果△g(mGal)(km)(km)2.70(g/cm)32.30(g/cm)32.52(g/cm)3MVFSA算法的实际资料处理：-8-6-4-202.342.542.702.542.602.6005101520253035510g(mGal)(km)(km)Q-J2.60D-SAnST-C2.60Q-J实测点拟合曲线N12新87反演界面MT解释综合解释目标函数与迭代关系曲线图200400600800100012001400123迭代次数目标函数50100150200123模型修正阶段放大二MT-重力联合反演研究联合反演研究现状：•线性反演算法居多、非线性算法少，用模拟退火算法进行联合反演研究更少。•尽管目前开展的联合反演研究已有多种，但研究内容主要集中在地震与重力、地震与MT联合反演的研究。•有关电磁测深与重力的联合反演研究只查阅到一篇相关的论文，因此对这方面的研究基本上还是空白。1.MT与重力联合可以使两方法相互弥补•MT纵向分辨率与重力横向的分辨率的互补•实际工作中MT方法的测点点距一般较稀，而野外重力数据的采集点较密。MT-重力联合反演必要性:2.充分利用野外资料：在生产实际中，非地震方法一般同时开展。MT-重力联合反演可能性:•电性与密度同源界面是两种方法联合的前提•地下结构有可能以同源或部分同源的形式出现示例:海安-盐城地区密度与电阻率统计表表5-1地层密度(×103kg/m3)电阻率(Ωm)Q16N2.138.7E2.30-2.313.9K2t4.8K2p2.4611.9T1s150T1x84P1Q2.60275C3c3000C2n1617C1h100C1g18C1j2.58600D3w2.5387S3m58S2f66S1g2.5427.5O3d118O2d5218O1I2.6710051ε3g2.70-2.732785电性与密度界面有几个是一致的：•上第三系-下第三系的物性界面•下第三系-白垩系的物性界面•侏罗系-三叠系物性界面等。开展MT与重力两者的联合反演是可能的。MT与重力联合反演技术难点：•如何处理电性与密度界面的不一致情况是首先要遇到的技术难点•如何构筑一个共同的目标函数—因为地下界面的变化使MT与重力两者的场值变化幅度是不同的MT的正演可以写为以下变分问题：经比较，本文选用有限元法作为二维MT正演方法，既保持较高的计算精度，又适应于复杂结构的地电模型。max2222211zzGdyVdydzVzVyVVI电性与密度界面的GIS编码技术•若是共同的物性界面，MT与重力均需要作正演计算•若是单独的电性界面，只参加MT正演计算•若是单独的密度界面，只参加重力正演计算对考虑三种情况：由此产生的三个问题：•如何记录由MT与重力合成的模型，这个模型中有两者一致和不一致的物性界面•如何扰动由两种不同地球物理方法合成的模型•如何让扰动后模型分别参加MT与重力正演计算以求得有关的目标函数解决办法：借鉴GIS中的属性编码思想•记录模型数据的MTG(i,j),按照MVFSA算法的具体要求进行全局或在当前模型周围扰动产生新模型•编码数组Code(k)（k=1,2…N)，对记录电性与密度界面情况的数组MTG(i,j)进行界面属性编码•各自的正演计算具体的界面属性编码方法：•电性与密度界面一致：Code(k)＝0，MT与重力均作正演计算；•只是电性界面：Code(k)＝1，只有MT作正演计算；•只是密度界面：Code(k)＝－1，只有重力作正演计算；目标函数及其加权联合：•WMT+Wg=1，体现了联合后的目标函数与Φmtφg的尺度相同，以防止目标函数的发散。•WMT、Wg的大小决定了何者为主反演方法，相应的权系数越大，则该方法为主反演方法。•落实到WMT、Wg的具体取值则是一个非常复杂的问题。ggMTMTWW背斜模型的反演结果-5-4-3-2-1020024681012141618205101520模型计算反演结果500△g(mGal)kmkm模型计算反演结果(Ω·m)(Ω·m)2.70(g/cm)32.30(g/cm)3二层模型的反演结果-6-4-2030105000246810121416182013579模型计算反演结果模型计算反演结果g(mGal)电性界面电性-密度界面2.70(g/cm)32.30(g/cm)3(Ω·m)(Ω·m)(Ω·m)(km)(km)三层模型反演结果-6-4-203002468101214161820246g(mGal)(km)(km)电性界面10500模型计算反演结果电性-密度界面密度界面(Ω·m)(Ω·m)(Ω·m)2.70(g/cm)32.50(g/cm)32.30(g/cm