CN108132486A - 一种重磁梯度与地震数据联合界面反演的优化模拟退火法 - Google Patents

Abstract

本发明公开了一种重磁梯度与地震数据联合界面反演的优化模拟退火法，包括以下步骤：模型确定，目标函数f(m)计算，随机扰动产生新的模型，目标函数f(m_i)计算，Δf值运算，i值判断，收敛条件判断，输出模型。油气汇聚大多为多个断裂的交叉位置，因此断裂划分的准确性直接决定钻井的位置，断裂位置是分区板块构造单元、地球动力学的主要组成部分，其准确性直接决定地球模型分布；本发明能够准确的计算出结果，可有效地降低油气勘探的成本，避免不必要的损失。

Description

一种重磁梯度与地震数据联合界面反演的优化模拟退火法

技术领域

本发明涉及地球科学技术领域，具体为一种重磁梯度与地震数据联合界 面反演的优化模拟退火法。

背景技术

重磁勘探具有快速、经济、范围大等优点，是地球物理勘查的常用手段 之一。利用观测得到的重磁异常可反演得到地下界面(层位)的具体分布， 进而可为下一步的详细勘探计划提供依据。

现今进行界面反演的常用方法是Parker-Oldenburg法，该方法具有计 算速度快和精度高的优点，在实际数据解释中应用十分广泛。重磁梯度测量 是近几年兴起的一种高精度地球物理勘探手段，其具有精度高、参数多、抗 噪性强等优点，测量结果能更好地描述地下异常体的特征，具有更高的分辨 率，能更好地描述界面的细节特征，对于界面的突变特征具有更高的分辨率。 从重磁基本原理出发，依据重磁异常界面正反演公式，推导出利用重磁梯度 数据进行界面反演的基本公式，并对各个梯度的反演结果进行加权，其中权因子的分配是依据各个分量的测量精度来设定。

重磁梯度与地震数据反演时所针对的参数具有线性对应性，因此这两类 数据可以有效地结合进行界面(地层)的反演，提高反演结果的准确性。模 拟退火算法(SA)是近几年来发展起来的一种全局最优化算法[1]，其主要优 点是，不用求目标函数的偏导数，不用解大型矩阵方程组，即能找到一个全局 最优解，且易于加入约束条件，方法易于移植，克服了局部反演方法依赖于初 始模型的选取以及解容易落入局部极值的缺陷。重磁梯度数据具有高分辨率 和高精度的优点，因此将重磁梯度与地震数据通过模拟退火法实现组合，进而通过寻找满足两种数据特征的全局最优解，有效地提高了对地层突变特征 的分辨率，获得更加准确的界面分布。

界面特征对于分析油气分布、区域构造特征均具有重要的意义，其准确 性直接制约资源潜力评价。因此将重磁梯度与地震数据联合约束界面反演可 有效地提高反演结果的准确性，准确的界面分布对于油气资源潜力评价具有 重要的价值，对于下一步精细勘探计划具有重要的意义。

现有界面反演技术是依据重磁原始异常，但对于界面的突变特征反映不 明显。地震勘探在火成岩覆盖区和构造破碎区难以获得良好的反射结果，从 而无法准确地获得层位分布。

依据地震测量结果难以完整地呈现地层的真实形态，因此利用重磁梯度 数据来补充地震未明区的反演工作。重磁梯度是重磁异常的变化率，能更明 显地反映界面的变化特征。重磁梯度与地震数据良好地互补，可准确地获得 界面的形态，对于区域油气资源潜力评价具有重要的价值。

本技术是将重磁梯度与地震数据通过模拟退火法实现组合，进而通过寻 找满足两种数据特征的全局最优解，有效地提高了对地层突变特征的分辨率， 获得更加准确的界面分布。

发明内容

本发明的目的在于提供一种重磁梯度与地震数据联合界面反演的优化模 拟退火法，以解决上述背景技术中提出的现有界面反演技术是依据重磁原始 异常，但对于界面的突变特征反映不明显。地震勘探在火成岩覆盖区和构造 破碎区难以获得良好的反射结果，从而无法准确地获得层位分布的问题。

为实现上述目的，本发明提供如下技术方案：

一种重磁梯度与地震数据联合界面反演的优化模拟退火法，包括以下步 骤：

S1、模型确定，对重磁梯度与地震数据进行收集数据，基于已知信息确 定模型空间及初始模型，设置模型参数；

S2、目标函数f(m)计算，重磁正演及梯度异常、地震旅行时和目标函数f(m)，f(m)＝W_gf_g(m)+W_gzf_gz(m)+W_sf_s(m)；

S3、随机扰动产生新的模型m_i；

S4、目标函数f(m_i)计算，计算重磁及其梯度异常、地震旅行时和目标函 数f(m_i)，f(m_i)＝W_gf_g(m_i)+W_gzf_gz(m_i)+W_sf_s(m_i)；

S5、Δf值运算，Δf＝f(m_i)-f(m)，当Δf≤0时完全接受模型，当Δf＞0 时按照概率P接受模型；

S6、i值判断，当i≥M时，按照降温原则降低温度，当i＜M时，重新返 回S3步骤进行重新测算；

S7、收敛条件判断，计算结果是否满足收敛条件判断，当判断为不满足 时，以i＝1重新代入S3步骤中进行运算；

S8、输出模型，当S7中满足收敛条件判断为满足时，输出模型并结束运 算。

优选的，密度界面重磁梯度异常的理论计算公式如下：

单一密度界面在水平面(xoy面)内的引力位为：

引力位波谱为：

根据Erdelyi(1954)积分表，上式可变为：

ξ和η的积分限是-∞～∞，而ζ的积分限是(h₀+Δh)～h₀，其中Δh为界面 相对平均深度h₀的起伏，界面在h₀以上Δh为负，在h₀以下Δh为正，而Δh又是 ξ和η的函数，故：

已知重磁异常垂直导数波谱与引力位波谱之间的关系为

则在地面引起的重磁异常垂直导数的波谱为：

因为且[-Δh(w)]⁰/0！＝1，所以上式可以改写离散形式：

通过以上的过程获得了界面重磁梯度的计算方法。

优选的，目标函数：采用重磁梯度异常单独进行反演的目标函数为：

M为测点数；L为模型界面的层数；与为第i点上模型计算得到的 重磁梯度异常与观测值；λ＝[v₁，v₂，...，v_L-1，σ_L1，σ_L2，...，σ_LM，h_L1，h_L2，...h_(L-1)M]^T为模型参数 矢量；W_gz为重磁梯度误差φ_gz的权系数；

模型扰动：模拟退火算法中新模型的产生是对当前模型进行扰动得到的； 采用Cauchy分布方法来产生新模型，即

x_i′＝x_i+y_i(B_i-A_i)

y_i＝T sgn(r-0.5)[(1+1/T)^|2u-1|-1]

式中x_i为当前模型的第i个参数，x_i′为随机扰动后的模型参数，[A_i，B_i]是 x_i的取值范围，r为(0，1)之间的随机数，T为当前温度，u是确定非均匀性 程度的常数，也称之为形状因子，sgn为符号函数，y_i称为扰动因子；

采用依赖于温度的似Cauchy分布产生新模型的优点是在高温情况下进行 大范围的搜索，在低温时仅在当前模型附近进行搜索，由于似Cauchy分布有 一平坦的“尾巴”，使其更易于快速地跳出局部极值，加快了模拟退火方法 的收敛速度；

接收概率：多维分形中的一种新的熵，即

式中，K为常数；q为任意实数；p_i为内循环的接受概率；由上 式推导出的广义Gibbs分布为

β为常数。当q→1时上式即为

p_i＝exp(-βE_i)/Z₁

E_i为***的能量，式(14)、(15)为模拟退火算法的理论基础；

新的接收概率的计算公式为：

P＝[1-(1-h)ΔE/T]^1/(1-h)

式中ΔE为扰动得到的新模型的目标函数E(m)与当前模型的目标函数E(m₀)之差，即ΔE＝E(m)-E(m₀)；T为温度(T＝1/β)；h为实数；新模型按式 (16)计算得到的概率进行接受；

显然，当h→1时式(16)为

P＝exp(-ΔE/T)

降温方式：模拟退火方法的降温方式

T(K)＝T₀exp(-CK^1/N)

式中，T₀为初始温度，K为迭代次数，C为给定的常数，N为待反演的 参数的个数，上式可改写成其中0.7≤α≤1，在应用中，1/N采用 0.5或1来代替。

与现有技术相比，本发明的有益效果是：油气汇聚大多为多个断裂的交 叉位置，因此断裂划分的准确性直接决定钻井的位置，断裂位置是分区板块 构造单元、地球动力学的主要组成部分，其准确性直接决定地球模型分布； 本发明能够准确的计算出结果，可有效地降低油气勘探的成本，避免不必要 的损失。

附图说明

图1为本发明模拟退火算法反演重磁梯度异常的基本流程结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行 清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而 不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做 出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1，本发明提供一种技术方案：一种重磁梯度与地震数据联合界 面反演的优化模拟退火法，包括以下步骤：

S1、模型确定，对重磁梯度与地震数据进行收集数据，基于已知信息确 定模型空间及初始模型，设置模型参数；

S2、目标函数f(m)计算，重磁正演及梯度异常、地震旅行时和目标函数 f(m)，f(m)＝W_gf_g(m)+W_gzf_gz(m)+W_sf_s(m)；

S3、随机扰动产生新的模型m_i；

S4、目标函数f(m_i)计算，计算重磁及其梯度异常、地震旅行时和目标函 数f(m_i)，f(m_i)＝W_gf_g(m_i)+W_gzf_gz(m_i)+W_sf_s(m_i)；

S5、Δf值运算，Δf＝f(m_i)-f(m)，当Δf≤0时完全接受模型，当Δf＞0 时按照概率P接受模型；

S6、i值判断，当i≥M时，按照降温原则降低温度，当i＜M时，重新返 回S3步骤进行重新测算；

S7、收敛条件判断，计算结果是否满足收敛条件判断，当判断为不满足 时，以i＝1重新代入S3步骤中进行运算；

S8、输出模型，当S7中满足收敛条件判断为满足时，输出模型并结束运 算。

具体的，密度界面重磁梯度异常的理论计算公式如下：

单一密度界面在水平面(xoy面)内的引力位为：

引力位波谱为：

根据Erdelyi(1954)积分表，上式可变为：

ξ和η的积分限是-∞～∞，而ζ的积分限是(h₀+Δh)～h₀，其中Δh为界面 相对平均深度h₀的起伏，界面在h₀以上Δh为负，在h₀以下Δh为正，而Δh又是 ξ和η的函数，故：

已知重磁异常垂直导数波谱与引力位波谱之间的关系为

则在地面引起的重磁异常垂直导数的波谱为：

因为且[-Δh(w)]⁰/0！＝1，所以上式可以改写离散形式：

通过以上的过程获得了界面重磁梯度的计算方法。

具体的，目标函数：采用重磁梯度异常单独进行反演的目标函数为：

M为测点数；L为模型界面的层数；与为第i点上模型计算得到的 重磁梯度异常与观测值；λ＝[v₁，v₂，...，v_L-1，σ_L1σ_L2，...，σ_LM，h_L1，h_L2，...h_(L-1)M]^T为模型参数 矢量；W_gz为重磁梯度误差φ_gz的权系数；

模型扰动：模拟退火算法中新模型的产生是对当前模型进行扰动得到的； 采用Cauchy分布方法来产生新模型，即

x_i′＝x_i+y_i(B_i-A_i)

y_i＝T sgn(r-0.5)[(1+1/T)^|2u-1|-1]

式中x_i为当前模型的第i个参数，x_i′为随机扰动后的模型参数，[A_i，B_i]是 x_i的取值范围，r为(0，1)之间的随机数，T为当前温度，u是确定非均匀性 程度的常数，也称之为形状因子，sgn为符号函数，y_i称为扰动因子；

采用依赖于温度的似Cauchy分布产生新模型的优点是在高温情况下进行 大范围的搜索，在低温时仅在当前模型附近进行搜索，由于似Cauchy分布有 一平坦的“尾巴”，使其更易于快速地跳出局部极值，加快了模拟退火方法 的收敛速度；

接收概率：多维分形中的一种新的熵，即

式中，K为常数；q为任意实数；p_i为内循环的接受概率；由上 式推导出的广义Gibbs分布为

β为常数。当q→1时上式即为

p_i＝exp(-βE₁)/Z₁

E_i为***的能量，式(14)、(15)为模拟退火算法的理论基础；

新的接收概率的计算公式为：

P＝[1-(1-h)ΔE/T]^-1(1-h)

式中ΔE为扰动得到的新模型的目标函数E(m)与当前模型的目标函数 E(m₀)之差，即ΔE＝E(m)-E(m₀)；T为温度(T＝1/β)；h为实数；新模型按式 (16)计算得到的概率进行接受；

显然，当h→1时式(16)为

P＝exp(-ΔE/T)

降温方式：模拟退火方法的降温方式

T(K)＝T₀exp(-CK^1/N)

式中，T₀为初始温度，K为迭代次数，C为给定的常数，N为待反演的 参数的个数，上式可改写成其中0.7≤α≤1，在应用中，1/N采用 0.5或1来代替。

界面特征对于分析油气分布、区域构造特征均具有重要的意义，其准确 性直接制约资源潜力评价。重磁梯度与地震数据反演时所针对地参数具有线 性对应性，因此这两类数据可以有效地结合对界面(地层)进行反演，提高 反演结果的准确性。重力梯度数据具有高分辨率和高精度的优点，因此将重 力梯度与地震数据联合反演可有效地提高对地层突变特征的分辨率，可获得 更加准确的界面分布。模拟退火算法(SA)是近几年来发展起来的一种全局最 优化算法，其主要优点是，不用求目标函数的偏导数，不用解大型矩阵方程组，即能找到一个全局最优解，且易于加入约束条件，方法易于移植，克服了局部 反演方法依赖于初始模型的选取以及解容易落入局部极值的缺陷。基于模拟 退火法的优势，将综合考虑重力梯度与地震数据，获得最优的反演结果。此 外，在进行界面反演的过程中针对反演迭代过程中的算法进行了优化，避免 了以往该类方法所存在的大量重复计算，能更有效地完成大数据量的反演工 作，从而为大面积石油资源勘查以及区域构造普查提供有力的技术支撑。

实施例

根据中国南海重磁垂直梯度计算海底地形，为了得到中国南海海域重磁 垂直梯度，首先利用Geosat，ERS-1/2，T/P，Jason-1，EnviSat-1卫星测高数据 采用垂线偏差计算大地水准面的方法计算出高分辨率的海洋大地水准面，然 后采用下式计算出重磁垂直梯度：

g_z＝γ₀N_zz＝-γ₀(N_xx+N_yy).............................................(5.13)

其中，γ₀为参考椭球面正常重磁值，N_xx、N_yy和N＝分别为海洋大地水准 面在3个坐标轴方向上的二阶导数。

计算得到的中国南海海域重磁垂直梯度异常，采用快速模拟退火法对该 异常进行反演来获得中国南海海底地形，其初始模型深度为5647m，最终迭代 均方误差为1.558×10-4mGal/m，迭代次数为511次。反演得到的海底地形具 有较高的分辨率，能很好地反应细节地形变化。为了验证其结果的正确性， 将其与LDEO船测海深进行比较，本文方法的反演结果与实际结果的最大差距 为240m，而在大部分地区差距均较小，结果可信度较高，因此利用模拟退火 法进行梯度异常的界面反演在实际数据解释工作中也有良好的应用效果。根 据海底地形的特征可油气沉积有利位置。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而 言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行 多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限 定。

Claims (3)

1.一种重磁梯度与地震数据联合界面反演的优化模拟退火法，其特征在于：包括以下步骤：

S1、模型确定，对重磁梯度与地震数据进行收集数据，基于已知信息确定模型空间及初始模型，设置模型参数；

S2、目标函数f(m)计算，重磁正演及梯度异常、地震旅行时和目标函数f(m)，f(m)＝W_gf_g(m)+W_gzf_gz(m)+W_sf_s(m)；

S3、随机扰动产生新的模型m_i；

S4、目标函数f(m_i)计算，计算重磁及其梯度异常、地震旅行时和目标函数f(m_i)，f(m_i)＝W_gf_g(m_i)+W_gzf_gz(m_i)+W_sf_s(m_i)；

S5、Δf值运算，Δf＝f(m_i)-f(m)，当Δf≤0时完全接受模型，当Δf＞0时按照概率P接受模型；

S6、i值判断，当i≥M时，按照降温原则降低温度，当i＜M时，重新返回S3步骤进行重新测算；

S7、收敛条件判断，计算结果是否满足收敛条件判断，当判断为不满足时，以i＝1重新代入S3步骤中进行运算；

S8、输出模型，当S7中满足收敛条件判断为满足时，输出模型并结束运算。

2.根据权利要求1所述的一种重磁梯度与地震数据联合界面反演的优化模拟退火法，其特征在于：密度界面重磁梯度异常的理论计算公式如下：

单一密度界面在水平面(xoy面)内的引力位为：

引力位波谱为：

根据Erdelyi(1954)积分表，上式可变为：

ξ和η的积分限是-∞～∞，而ζ的积分限是(h₀+Δh)～h₀，其中Δh为界面相对平均深度h₀的起伏，界面在h₀以上Δh为负，在h₀以下Δh为正，而Δh又是ξ和η的函数，故：

已知重磁异常垂直导数波谱与引力位波谱之间的关系为

则在地面引起的重磁异常垂直导数的波谱为：

因为且[-Δh(w)]⁰/0！＝1，所以上式可以改写离散形式：

通过以上的过程获得了界面重磁梯度的计算方法。

3.根据权利要求1所述的一种重磁梯度与地震数据联合界面反演的优化模拟退火法，其特征在于：目标函数：采用重磁梯度异常单独进行反演的目标函数为：

M为测点数；L为模型界面的层数；与为第i点上模型计算得到的重磁梯度异常与观测值；λ＝[v₁，v₂，...，v_L-1，σ_L1，σ_L2，...，σ_LM，h_L1，h_L2，...h_(L-1)M]^T为模型参数矢量；W_gz为重磁梯度误差φ_gz的权系数；

模型扰动：模拟退火算法中新模型的产生是对当前模型进行扰动得到的；采用Cauchy分布方法来产生新模型，即

x_i′＝x_i+y_i(B_i-A_i)

y_i＝T sgn(r-0.5)[(1+1/T)^|2u-1|-1]

式中x_i为当前模型的第i个参数，x_i′为随机扰动后的模型参数，[A_i，B_i]是x_i的取值范围，r为(0，1)之间的随机数，T为当前温度，u是确定非均匀性程度的常数，也称之为形状因子，sgn为符号函数，y_i称为扰动因子；

采用依赖于温度的似Cauchy分布产生新模型的优点是在高温情况下进行大范围的搜索，在低温时仅在当前模型附近进行搜索，由于似Cauchy分布有一平坦的“尾巴”，使其更易于快速地跳出局部极值，加快了模拟退火方法的收敛速度；

接收概率：多维分形中的一种新的熵，即

式中，K为常数；q为任意实数；p_i为内循环的接受概率；由上式推导出的广义Gibbs分布为

β为常数。当q→1时上式即为

p_i＝exp(-βE_i)/Z₁

E_i为***的能量，式(14)、(15)为模拟退火算法的理论基础；

新的接收概率的计算公式为：

P＝[1-(1-h)ΔE/T]^1/(1-h)

式中ΔE为扰动得到的新模型的目标函数E(m)与当前模型的目标函数E(m₀)之差，即ΔE＝E(m)-E(m₀)；T为温度(T＝1/β)；h为实数；新模型按式(16)计算得到的概率进行接受；

显然，当h→1时式(16)为

P＝exp(-ΔE/T)

降温方式：模拟退火方法的降温方式

T(K)＝T₀exp(-CK^1/N)

式中，T₀为初始温度，K为迭代次数，C为给定的常数，N为待反演的参数的个数，上式可改写成其中0.7≤α≤1，在应用中，1/N采用0.5或1来代替。