多孔隙介质参数识别的正则化模拟退火法①

王迎

(黑龙江工业学院机电工程学院，黑龙江鸡西158100)

随着科学技术的飞速发展和人们生活水平的不断提高，世界范围内面临着油气资源需耗量大幅增加的趋势，而现有的油气资源是有限的，如果不去进一步勘察探寻油气资源的分布，那么势必会经受资源急剧短缺的严峻问题.可见，资源勘探技术尤为重要，其中涉及到在介质的波动方程中识别参数的问题.研究人员在勘探资源之前需要预先知道一些参数，诸如孔隙率、饱和度等来还原对应的介质波动方程，以此来预测所研究区域的地质层面结构，而要预知这些参数，就要通过参数识别方法来完成[1].另外，所勘探的地层通常是多孔介质，更接近真实的地层介质，因此其研究更具有实际意义[2].本文中有效结合了具有全局收敛特性的模拟退火法以及求解非适定问题的正则化方法，形成了二维多孔隙介质波动方程的参数孔隙率识别问题的正则化-模拟退火算法.

1 参数识别模型

描述二维多孔隙介质波动方程的模型为

固体和流体物质的物理参数及各个参数间的关系如下所示，这些参数均为关于x和z的变量.

Kb：骨架体变模量；ρs：固体密度；Gb：骨架切变模量；m：有效流体密度；Ks：固体体变模量；ρf：液体密度；λb：固体骨架的Láme系数；λ⁀：双相介质的Láme系数；μ：孔隙介质的切变模量；Kf：液体体变模量；β：孔隙率；α⁀、M：饱和多孔隙介质压缩性系数，且系数满足

边界条件：

初始条件：

(1)～(10)式构成了二维双相介质波动方程的正问题模型，即由已知孔隙率β(x,z)可以求解出u(x,z,t).若孔隙率β(x,z)未知，给出附加条件

若对二维双相介质波动方程(1)～(4)式考虑平行波入射时，f=g(t)δ(x)，其中g(t)为震源函数，δ(x)为δ函数，且满足g(t)=0,t＜0，则(1)～(11)式就构成了二维多孔隙介质波动方程一般意义下的参数识别问题.

2 模型的离散化与参数识别算法

2.1 离散化模型设即分别为沿x,z方向的步长，τ为时间步长，记对方程(1)～(11)进行差分离散.为了简单方便，将离散结果统一记作(12)式.

其中

其中

2.2 参数识别的算法记

对β(i,j)及分别按适当顺序排列，可排列成向量形式，分别记为B,G.

令

由(12)式定义了一个非线性向量值函数A:B→G.假设已知也进行对应排列，组成向

于是反演βi,j的问题转化成求解泛函

的非线性优化问题.

为了克服反问题的不适定性，进而达到数值求解的稳定性，利用正则化方法在(13)式中引入光滑泛函，获得泛函

其中，α为正则化参数，M1、M2分别为x方向、z方向的二阶光滑矩阵[4].

通过求解满足下式的最优解B*，将要求解的二维流体饱和多孔隙介质波动方程反问题再次转化为求解非线性优化问题

Jα(B)的梯度为

其中，A′(B)是A(B)关于B的Jacobi矩阵.

模拟退火法作为随机化算法[4]的通用方法之一，主要思想是通过设置随机迭代函数来产生随机解，最终在有限时间内搜寻到符合终止条件的最优解.根据所研究的问题给出一个初始猜测值B0，通过(16)式迭代产生B1，B2，….我们期望存在一个Bk是(15)式的解或者存在某一收敛于解的点列.

其中c为产生扰动的系数，rand(q+1,p+1)为随机数矩阵.

在上述模型的基础上，结合正则化方法，得到多孔隙介质波动方程反演的正则化模拟退火算法，其算法流程为：

步骤1:给定初始温度T0，初始点B0为对应阶初值矩阵，及ε＞0，正则化参数α＞0,k⇐0.

否则，如果k+1=N，则Tn+1=σTn，B0⇐Bk+1,k⇐0，转步骤2；

否则，转步骤4.

步骤 4：计算 Δm=Jα(Bk+1)-Jα(Bk)，

若Δm＜0，则Bk=Bk+1，k⇐k+1，转步骤2；

否则，转步骤2.

上述算法中的kB=1.380 7×10-23J·K-1为Boltzmann常数，rand( )q+1,p+1为随机数矩阵，c为产生扰动的系数，Δm为目标函数之差（评价函数的增量），σ为指数降温的系数，N为B的维数.我们可以用图1简要描述此算法流程：

图1 正则化模拟退火法流程图

3 数值模拟

数值模拟使用形式为g(t)=sin(2πvt)exp(-2πvt)的雷克子波，其中入射波的频率v取为75Hz.研究区域取为水平方向距离为260m、竖直方向为260m的区域，即a=260，b=260.在进行有限差分离散时，取水平差分步长为hx=10m，竖直差分步长为hz=10m，接收器在地表的每个差分点上.正则化参数α=0.000 95，相对误差ε=1%，T0=1 050，c=1.5，σ=0.85，τ=0.002 s，参数识别时取孔隙率的初始值B0为所有元素都是0.2的相应阶矩阵.

所选模型如图1所示，它是一种异常体，其中两块具有尖点，其相应的各介质参数如表1所示.

表1 多孔隙介质模型的参数表

假设背景介质的孔隙率的精确值为0.4，位于上方的异常体内的孔隙率的精确值为0.25，位于下方的异常体内的孔隙率的精确值为0.3.对此模型使用正则化-模拟退火算法进行反演得到的结果如图2所示.

4 结论

基于全局收敛的模拟退火方法与求解不适定性问题的正则化方法的思想，将两者有机结合构造了正则化-模拟退火算法.为增加参数识别的收敛性和稳定性，利用此算法识别二维多孔隙介质波动方程的参数孔隙率.数值模拟的结果显示出此算法能克服一般的参数识别方法的局部收敛性，能有效地解决具有非适定性、非线性性质的参数识别问题.

图1 孔隙率模型

图2 正则化模拟退火算法对孔隙率模型的参数识别结果

[1]王立荣.反问题的数值解法[D].深圳:哈尔滨工业大学深圳研究生院，2011.

[2]田迎春，章梓茂，赵成刚，等.黏弹性流体饱和孔隙介质动力反应分析的显式有限元法[J].科学技术与工程，2010，10(5):1112-1112.

[3]张新明，刘克安，刘家琦.流体饱和多孔隙介质波动方程多尺度反演[J].应用力学学报，2007，24(1):87-90.

[4]黄锋振，袁艳，张修宝，等.基于模拟退火的傅里叶变换成像光谱仪干涉图相位修正[J].光谱学与光谱分析，2011，31(7):2012-2012.