利用模拟退火法反演地下三维密度界面的数值模拟

利用模拟退火法反演地下三维密度界面的数值模拟
李宗浩;杨宏飞
【摘要】解的优化问题是重磁反演中很重要的一个方面.为了提高反演的精度和效率,采用一种线性网格化的模拟退火法来进行密度反演.将模拟区域划分为若干规则且密度均匀的长方体.为了减小边界条件的影响,可将研究区域中不同部分划分为规模不同的长方体.由引力位可推导出长方体外任意一点的重力异常公式;利用重力的可叠加性,计算每个观测点的重力异常.在给定长方体的规模参数后,根据巳知的重力异常,利用模拟退火法求取最优化解,从而确定密度参数.
【期刊名称】《科学技术与工程》
【年(卷),期】2014(014)004
【总页数】4页(P153-156)
【关键词】优化问题;网格化建模;模拟退火法;密度反演
【作者】李宗浩;杨宏飞
【作者单位】成都理工大学地球物理学院,成都610059;成都理工大学能源学院,成都610059
【正文语种】中文
【中图分类】P318.6
通过重力异常反演地球内部横向的密度结构，是为了了解地球内部的精细的地质结构。

通过求解出密度的不均匀分布，就可以判断是否存在高，低密度体。

同时可以
根据密度分布显示他们的位置，大小和形态等。

目前，对于多层密度界面的反演一般都是在空间域或者频率域中进行。

肖鹏飞［1］等在2007年采用徐世浙提出的
空域迭代法代替原Parker-Oldenburg［2］法中的向下延拓算子，极大地提高了
密度反演的精度和稳定性。

张风旭［3］等在2006年在Parker算法的基础上，用余弦变换代替FFT，从而进一步提高了重力界面反演的计算精度和速度，但向下延拓的精度和稳定性没有得到实质性的提高。

空间域由于计算速度慢，只适合于局部异常体的研究，频率域中的Parker-Oldenburg位场迭代反演方法具有计算速度
快的优势，但在利用Parker-Oldenburg方法进行界面反演时，未涉及物性随深
度的变化情况。

模拟退火法是一种有别于常规最优化方法的算法。

该方法可以克服目标函数局部极小的限制，从而使反演获得全局最优解，提高反演精度和效率。

本文将研究区域进行网格化划分，并结合模拟退火法对研究区域进行三维密度界面的数值模拟。

通过分析，取得了较好的反演结果。

1 方法原理及模型构建
1.1 方法原理
模拟退火法源于统计热力物理学，它模拟熔融状态下物体缓慢冷却达到结晶状态的物理过程.模拟退火反演算法的基本思想是:生成一系列参数向量模拟粒子的热运动，通过缓慢地减小一个模拟温度的控制参数，使模拟的系统最终冷却结晶达到系统能量最小值的过程。

模拟退火法与传统寻优方法在迭代过程中对模型的接收上有所不同。

模拟退火法不但接收误差能量减小的模型，同时也以某种概率接收误差能量增加的模型。

正是这种接收模式避免了迭代搜索过程陷入局部最优解的局面。

这一突破点使该方法能够收敛到全局最优解。

提高反演的精度和效率，使反演结果更理想。

1.1.1 模拟退火算法
据王山山等［4］，任义庆等［5］对算法的研究:设反演的目标函数为E(ml)，
m∈Ω，m为模型，Ω为反演的解空间(模型空间)。

假定有待反演的N个模型参数表示为:m=m1，m2，…，m N，其中每一个 m i 有一个对应的模型空间［m imin，m imax］。

其算法步骤为:
1)随机产生一个初始模型m0，设定系统的初始温度T0，模型反演空间，最大迭代次数，拟合误差阀值ε等参数。

令l=0;
2)对第l次迭代的模型参数值m l，计算其目标函数E(ml)，计算温度Tl，在模型空间里随机地修改模型参数值为 m l+1。

按下式计算修改的接收概率:
式(1)中:ΔE=E(m l+1)－E(m l)为目标函数的差值。

若ΔE＜0，表明系统朝着能量减小的方向移动，新模型m l+1可以概率1接收。

若ΔE＞0，计算概率 p(ΔE)。

显然，0 ＜ p(ΔE)＜ 1。

然后在［0，1］中产生一个随机数R。

若p(ΔE)＜R。

则接收新模型m l+1。

否则拒绝修改模型。

3)计算目标函数E(m l)，若E(m l+1)＜ε则退出，否则l=l+1，转到第二步。

对于退火过程，要求温度T随着迭代次数的增加而缓慢降低。

本文采用双曲线下降型。

式(2)中:Tl为第l次迭代的温度，T0为初始温度，l为迭代次数。

初始温度T0不能太高，否则增加计算时间，T0也不能太低，否则模型选取不能覆盖整个模型空间，只是在初始模型附近选取，不能进行全局寻优。

T0只能通过实验计算得到。

对于从模型参数值m l到m l+1的修改，我们才有对m l的每个元素进行随机扰动从而得到m l+1。

即
1.2 模型构建
设定的研究区域为400 km×60 km×60 km。

如图1所示，将研究区域划分为长
方体。

为了减小边界条件的影响，我们在研究区域内的最外层将单元格的大小设定为:10 km×10 km×10 km。

将次外层设置为5 km×5 km×5 km。

剩下的设置为
2 km×2 km×2 km。

相同大小的单元体的密度值相同，如图1所示。

据明圆圆［6］等对模型构造的划分:将图1中的任意一个单元体置于图2中。

据陈胜早［7］对单元块重力异常的计算:设其长宽高分别为a，b，c。

取每个单元体的几何中心(ε，η，ζ)为该单元的坐标。

ρ为该单元体的剩余密度。

则该单元体对其外一点p(x，y，z)产生的重力异常为
由重力的可加性可得所有长方体在p(x，y，z)处产生的总的重力异常为
图1 长方体水平层
图2 长方体外任一点的重力异常
2 数值模拟
设置单元体大小为10 km ×10 km×10 km的密度为:2 600 kg/m3，5 km×5
km×5 km的密度为:2 700 kg/m3，2 km×2 km×2 km的密度为:2 800 kg/m3。

为了编写程序的方便，我们将这些密度按从左上角的单元格为起点，顺时针旋转存储为一维的数组。

并且先存储10 km×10 km×10 km，继而是5 km×5 km×5 km，最后存储2 km×2 km×2 km的密度。

在平行于x方向，z=0分别布设三条测线，即:y=10，y=30，y=50。

每条测线上设置100个观测点，观测点之间的间距为4 km。

这样，把以上这些参数带入到编写的模拟退火法程序中。

可以得出以下结果。

图3 模型密度与反演密度对比图
图4 反演密度与理论密度之差
图5 反演密度与理论密度误差百分比
由于反演出来的密度数值的个数较多。

我们在三种不同规模的单元格中各选取一定比例的密度数值对比。

如图3～图5。

图3中:红色表示反演的密度，上层的绿色表示的是原始模型的密度。

从数值上分析，可知:他们的误差很小。

从图4可以更进
一步地分析，他们差值的最大值不到0.2 g/cm3。

其反演精度很高。

从图5可以
更进一步地看出，其误差百分比最大不超过6%。

以此说明，MSA方法的反演精
度较高。

从图6可以看出，此图6中包含了三条侧线的异常曲线对比。

其中，横坐标(0～100)表示第一条侧线，后面依次为第二，三条测线。

反演出的重力异常与理论异
常曲线吻合得非常好。

从图7中的反演异常与理论异常的相对误差来看，其中最
大的误差为2.4%，其他大部分误差都在0.2% ～1.4%之间，说明反演的精度很高。

从以上的对比可以看出:
1)从这次选用的迭代算法模拟退火法(MSA)可以看出。

其精度和误差精度都很高，这主要是MSA具有一定的概率跳出局部极小值的特性。

尽可能搜索全区的极小值。

2)从设置的模型可以看出，这次的模型考虑到边界效应。

特意设置了为了减小边界效应的边界模型。

其反演的精度和误差都有很大的提高。

3 结论
图6 三条反演重力异常测线与理论测线反演对比
图7 三条反演重力异常测线与理论测线反演异常相对误差
通过重力异常反演地球内部横向的密度结构，对于了解地球内部的精细的地质结构有着重要的意义。

对于进一步对地下资源的勘察也有着重要的意义。

本文通过构建一种利于减小边界效应的网格化的地质模型，并结合模拟退火法。

通过以上数值模拟可以看出该模型的构建对减少边界效应有很好的效果。

为密度反演提供了一种高效率的密度反演方法。

参考文献
【相关文献】
1 肖鹏飞，陈生昌，杨长福，等．油气藏潮汐重力的初步研究.石油物探，2007;46(2):202—206
2 柯小平，王勇，许厚泽，等．青藏东缘三维Moho界面的位场遗传算法反演.大地测量与地球动力学，2006;26(1):100—104
3 张凤旭，孟令顺，张凤琴，等．重力位谱分析及重力异常导数换算新方法—余弦变换.地球物理学报，2006;49(1):244—248
4 王山山，聂勋碧．基于波场计算的最优化速度反演.石油物探，1994;33(1):81—95
5 任义庆，徐仲达，马在田．应用模拟退火法反演横波速度.石油地球物理勘探，
1996;31(5):677—684
6 明圆圆，范美宁．鱼群算法的重力密度异常反演方法.物探化探计算技术，2012;34(6):666—671
7 陈胜早.多层变密度模型与壳幔结构研究.中国科学(B辑)，1989;(9):991—1000。