CN108229082A - 一种基于数据空间快速计算的联合反演方法 - Google Patents

Abstract

本发明公开了一种基于数据空间快速计算的联合反演方法，包括如下步骤：S1：建立重磁电震联合反演目标函数；S2：建立电阻率、速度、密度和磁化率的初始模型；S3：正演计算和雅克比矩阵的求取；S4：计算交叉梯度偏导数和拉格朗日算子；S5：计算获得模型空间迭代模型；S6：计算获得数据空间迭代模型；S7：对获得的数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则返回S3。该基于数据空间快速计算的联合反演方法计算速度快、占用内存小。

Description

一种基于数据空间快速计算的联合反演方法

技术领域

本发明涉及多物性联合反演的处理技术领域，特别提供了一种基于数据空间快速计算的联合反演方法。

背景技术

综合地球物理方法已经广泛应用于深部地质探测和矿产、水、石油、地热等资源勘查。而联合反演是其综合地球物理方法数据处理的主流方法，国外科学家已经提出了二种或多种物性参数之间结构耦合的联合反演方法，国内也实现了二种或多种物性参数之间的结构耦合联合反演方法。

目前，对开展大数据量的联合反演数据处理已经成为重要的发展方向。现有的联合反演多采用在模型空间下进行反演计算，但是当反演物性参数种类从二个增加到多个或者反演区域网格剖分数据量过大时，对计算机的计算速度和内存占用都将成为考验，将会严重制约了综合地球物理方法的发展。

发明内容

鉴于此，本发明的目的在于提供一种基于数据空间快速计算的联合反演方法，以解决现有的联合反演方法计算速度慢、占用内存大等问题。

本发明提供的技术方案是：一种基于数据空间快速计算的联合反演方法，包括如下步骤：

S1：建立重磁电震联合反演目标函数，所述联合反演目标函数如式(1)：

约束条件：τ(m)＝0

其中：

C_d＝diag[C_d1,C_d2,C_d3,C_d4]，C_m＝diag[C_m1,C_m2,C_m3,C_m4]，

α＝[α₁,α₂,α₃,α₄]，

τ_ij＝▽m_i×▽m_j，i,j＝1,2,3,4 (2)

其中，m为模型参数，m₁,m₂,m₃,m₄分别为电阻率、速度、密度和磁化率模型参数；m⁰为先验模型参数，m⁰¹,m⁰²,m⁰³,m⁰⁴分别为电阻率、速度、密度和磁化率先验模型参数；d为观测数据，d₁,d₂,d₃,d₄分别为视电阻率、地震走时、重力异常和磁异常；C_d为观测数据d的数据协方差矩阵， C_d1,C_d2,C_d3,C_d4分别为视电阻率、地震走时、重力异常和磁异常的数据协方差矩阵；C_m为模型参数m的模型协方差矩阵，C_m1,C_m2,C_m3,C_m4分别为电阻率、速度、密度和磁化率的模型协方差矩阵；α为阻尼因子，α₁,α₂,α₃,α₄分别为大地电磁、地震初至波走时、重力和磁法的阻尼因子，f(m)表示正演响应， f₁(m),f₂(m),f₃(m),f₄(m)分别为大地电磁、地震初至波走时、重力和磁法的正演响应，▽为梯度，τ为交叉梯度，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度；为数据拟合项、为模型平滑约束项；T和-1分别为矩阵转置和矩阵求逆；

S2：建立电阻率、速度、密度和磁化率的初始模型：

在笛卡尔坐标系下沿x，z二个坐标轴分别将初始模型空间划分成为 N_x，N_z个小的矩形网格单元，间距为Δx(i)(i＝1,...,N_x)，Δz(i)(i＝1,...,N_z)，初始模型均采用均匀半空间模型；

S3：正演计算和雅克比矩阵的求取：

采用Wanamaker(1987)提出的三角形剖分有限元算法进行大地电磁正演计算；采用了Singh(2002)提出的对任意多边形棱柱进行重磁异常正演计算；采用Colin Zelt(1998)改进的有限差分方法求解Eikonal方程进行地震初至波走时正演计算；

通过公式(3)求取雅克比矩阵：

其中，A为正演响应f(m)的雅克比矩阵；

S4：计算交叉梯度偏导数和拉格朗日算子：

通过公式(4)计算交叉梯度偏导数

其中，B为交叉梯度函数τ的导数，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度，i,j＝1,2,3,4；

拉格朗日算子计算公式如下：

首先将目标函数中的交叉梯度约束项通过拉格朗日算子法加入到目标函数中，如式(5)：

其中，Ψ为加入约束条件的目标函数，Γ为拉格朗日算子，

然后，对正演响应f(m)和交叉梯度约束条件τ(m)进行泰勒级数展开，如式(6)：

将式(6)代入到式(5)中，得式(7)：

其中，

对公式(7)求极值推导出模型改变量的表达式如式(8)：

Δm＝N^-1·n-N^-1·B^T·Γ (8)

其中：

将公式(8)带入到公式(1)中的交叉梯度约束项中，求取拉格朗日算子Γ，如式(11)：

S5：计算获得模型空间迭代模型：

将式(11)得到的拉格朗日算子反代入到式(8)中，得到模型空间模型改变量式(12)：

Δm_m＝N^-1·n-N^-1·B^T·Γ (12)

之后，获得模型空间迭代模型如式(13)：

m_m＝m₀+Δm_m (13)

S6：计算获得数据空间迭代模型：

对模型空间模型改变量式(12)进行恒等式变换

将式(14)和式(15)带入到式(8)中，得到数据空间模型改变量

Δm_s＝N^-1·n-N^-1·B^T·Γ (16)

之后，获得数据空间迭代模型如式(17)：

m_s＝m₀+Δm_s (17)

S7：对获得的数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3，

其中，拟合差通过式(18)计算获得

其中，N为观测数据个数。

本发明提供的基于数据空间快速计算的联合反演方法能够将野外采集到的观测数据通过反演计算转换成地下二维电阻率、密度、磁化率和速度模型的结构信息，从四种不同物性方面反映更为真实的地下结构信息。该反演方法克服了传统模型空间联合反演方法中，计算速度慢、占用内存大等问题，将模型空间存储和计算灵敏度矩阵通过数学推导转换到数据空间进行存储和计算，由于观测数据的个数远远小于模型参数的个数，因此在内存存储量方面，数据空间方法要明显优于模型空间方法，从而提高了计算速度和效率。该反演方法在保证计算精度的条件下，计算效率很高，适用于目前国内外十分主流的的windows操作系统，并且可以在普通个人电脑上推广使用。该反演方法能够为探索成矿区带不同深度尺度上地球物理特征与成矿之间的关系，为解决深部资源勘查提供有力的依据，指明了找矿、水、石油、地热等资源的方向。

附图说明

下面结合附图及实施方式对本发明作进一步详细的说明：

图1为本发明提供的基于数据空间快速计算的联合反演方法的流程图；

图2为验证本发明提供的方法而设计的理论模型图；

图3为(a)列为单独反演结果、(b)列为模型空间联合反演结果，(c) 列为数据空间联合反演结果，其中，每列由上至下分别表示电阻率、密度、磁化率和速度反演结果。

具体实施方式

下面将结合具体的实施方案对本发明进行进一步的解释，但并不局限本发明。

如图1所示，本发明提供了一种基于数据空间快速计算的联合反演方法，包括如下步骤：

S1：建立重磁电震联合反演目标函数，所述联合反演目标函数如式(1)：

约束条件：τ(m)＝0

其中：

C_d＝diag[C_d1,C_d2,C_d3,C_d4]，C_m＝diag[C_m1,C_m2,C_m3,C_m4]，

α＝[α₁,α₂,α₃,α₄]，

τ_ij＝▽m_i×▽m_j，i,j＝1,2,3,4 (2)

其中，m为模型参数，m¹,m²,m³,m₄分别为电阻率、速度、密度和磁化率模型参数；m₀为先验模型参数，m₀₁,m₀₂,m₀₃,m₀₄分别为电阻率、速度、密度和磁化率先验模型参数；d为观测数据，d₁,d₂,d₃,d₄分别为视电阻率、地震走时、重力异常和磁异常；C_d为观测数据d的数据协方差矩阵， C_d1,C_d2,C_d3,C_d4分别为视电阻率、地震走时、重力异常和磁异常的数据协方差矩阵；C_m为模型参数m的模型协方差矩阵，C_m1,C_m2,C_m3,C_m4分别为电阻率、速度、密度和磁化率的模型协方差矩阵；α为阻尼因子，α₁,α₂,α₃,α₄分别为大地电磁、地震初至波走时、重力和磁法的阻尼因子，f(m)表示正演响应， f₁(m),f₂(m),f₃(m),f₄(m)分别为大地电磁、地震初至波走时、重力和磁法的正演响应，▽为梯度，τ为交叉梯度，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度；为数据拟合项、为模型平滑约束项；T和-1分别为矩阵转置和矩阵求逆；

S2：建立电阻率、速度、密度和磁化率的初始模型：

在笛卡尔坐标系下沿x，z二个坐标轴分别将初始模型空间划分成为 N_x，N_z个小的矩形网格单元，间距为Δx(i)(i＝1,...,N_x)，Δz(i)(i＝1,...,N_z)，初始模型均采用均匀半空间模型；

S3：正演计算和雅克比矩阵的求取：

采用Wanamaker(1987)提出的三角形剖分有限元算法进行大地电磁正演计算；采用了Singh(2002)提出的对任意多边形棱柱进行重磁异常正演计算；采用Colin Zelt(1998)改进的有限差分方法求解Eikonal方程进行地震初至波走时正演计算；

通过公式(3)求取雅克比矩阵：

其中，A为正演响应f(m)的雅克比矩阵；

S4：计算交叉梯度偏导数和拉格朗日算子：

通过公式(4)计算交叉梯度偏导数

其中，B为交叉梯度函数τ的导数，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度，i,j＝1,2,3,4；

拉格朗日算子计算公式如下：

首先将目标函数中的交叉梯度约束项通过拉格朗日算子法加入到目标函数中，如式(5)：

其中，Ψ为加入约束条件的目标函数，Γ为拉格朗日算子，

然后，对正演响应f(m)和交叉梯度约束条件τ(m)进行泰勒级数展开，如式(6)：

将式(6)代入到式(5)中，得式(7)：

其中，

对公式(7)求极值推导出模型改变量的表达式如式(8)：

Δm＝N^-1·n-N^-1·B^T·Γ (8)

其中：

将公式(8)带入到公式(1)中的交叉梯度约束项中，求取拉格朗日算子Γ，如式(11)：

S5：计算获得模型空间迭代模型：

将式(11)得到的拉格朗日算子反代入到式(8)中，得到模型空间模型改变量式(12)：

Δm_m＝N^-1·n-N^-1·B^T·Γ (12)

之后，获得模型空间迭代模型如式(13)：

m_m＝m₀+Δm_m (13)

S6：计算获得数据空间迭代模型：

对模型空间模型改变量式(12)进行恒等式变换

将式(14)和式(15)带入到式(8)中，得到数据空间模型改变量

Δm_s＝N^-1·n-N^-1·B^T·Γ (16)

之后，获得数据空间迭代模型如式(17)：

m_s＝m₀+Δm_s (17)

S7：对获得的数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3，

其中，拟合差通过式(18)计算获得

其中，N为观测数据个数。

该基于数据空间快速计算的联合反演方法能够将野外采集到的观测数据通过反演计算转换成地下二维电阻率、密度、磁化率和速度模型的结构信息，从四种不同物性方面反映更为真实的地下结构信息。该反演方法克服了传统模型空间联合反演方法中，计算速度慢、占用内存大等问题，将模型空间存储和计算灵敏度矩阵通过数学推导转换到数据空间进行存储和计算，由于观测数据的个数远远小于模型参数的个数，因此在内存存储量方面，数据空间方法要明显优于模型空间方法，从而提高了计算速度和效率。该反演方法在保证计算精度的条件下，计算效率很高，适用于目前国内外十分主流的的windows操作系统，并且可以在普通个人电脑上推广使用。该反演方法能够为探索成矿区带不同深度尺度上地球物理特征与成矿之间的关系，为解决深部资源勘查提供有力的依据，指明了找矿、水、石油、地热等资源的方向。

下面对本发明提供的基于数据空间快速计算的联合反演方法进行验证。

如图2所示，设计一个物性结构一致的双棱柱模型，左侧异常体电阻率为10^1.5Ω·m，剩余密度为-2.0g/cm³，磁化率为-0.2SI，速度为2000m/s，异常体大小600×400m；右侧异常体电阻率为10000Ω·m，剩余密度为2.0g/cm³，磁化率为0.2SI，速度为5000m/s，异常体大小为600×900m；两个异常体顶面埋深均为500m，埋藏于剩余密度为0g/cm³，磁化率为0.0001SI，电阻率为100Ω·m，速度为4000m/s的均匀半空间中.大地电磁数据包含了TM、TE两种极化模式下的视相位和视电阻率，频点个数为10个，频率范围在1～1000Hz之间，9个观测点等距离分布在测线0～6km之间；重力、磁法观测点各有30个；地震震源位于地下50m处，共17个，接收点位于两口井中，水平位置分别为1.5km和4.5km，每口井内有10个接收点，间隔为0.2km.重力、磁法、地震走时地下二维等间隔网格剖分都为140×60， MT地下二维网格剖分为174×74，需要在重磁震等间距网格剖分的基础上向外延展，延展网格以不等间隔剖分，重磁震网格剖分区域即为联合反演区域。

图3为反演结果的对比图，通过图3可以看出：无论是模型空间联合反演(图3b)，还是数据空间联合反演(图3c)相比于单独反演(图3a)，异常体的空间几何形态均有更好的恢复，重、磁联合反演结果中异常体下部高密度和高磁化率发散现象得到明显改善，异常体上下界面清晰可见，这是由于交叉梯度结构相似性的约束，重磁模型结构受到了地震和大地电磁模型结构的制约，同时电阻率和地震速度模型也由于重磁方法的结构相似性的影响，改善了水平分辨能力。同时，数据空间联合反演结果相比于模型空间联合反演结果恢复效果更好，可以看出，模型空间反演结果较为光滑，异常体边缘相对模糊，而数据空间反演结果聚焦程度更高，异常体边界清晰，异常体几何形态和真实模型基本吻合。在反演计算时间方面，两种反演方法都迭代了8次，模型空间联合反演耗时约2164s，数据空间联合反演耗时约1883s；在内存消耗方面，模型空间联合反演计算时，占用最大内存约为675MB，而数据空间大约只有177MB。可以看出，地下不同物性模型结构不一致时，交叉梯度联合反演也可以得到准确的反演结果，同时无论在计算时间还是在内存消耗方面数据空间联合反演更具有优势。

Claims (1)

1.一种基于数据空间快速计算的联合反演方法，其特征在于，包括如下步骤：

S1：建立重磁电震联合反演目标函数，所述联合反演目标函数如式(1)：

约束条件：τ(m)＝0

其中：

C_d＝diag[C_d1,C_d2,C_d3,C_d4]，C_m＝diag[C_m1,C_m2,C_m3,C_m4]，

α＝[α₁,α₂,α₃,α₄]，

其中，m为模型参数，m₁,m₂,m₃,m₄分别为电阻率、速度、密度和磁化率模型参数；m₀为先验模型参数，m₀₁,m₀₂,m₀₃,m₀₄分别为电阻率、速度、密度和磁化率先验模型参数；d为观测数据，d₁,d₂,d₃,d₄分别为视电阻率、地震走时、重力异常和磁异常；C_d为观测数据d的数据协方差矩阵，C_d1,C_d2,C_d3,C_d4分别为视电阻率、地震走时、重力异常和磁异常的数据协方差矩阵；C_m为模型参数m的模型协方差矩阵，C_m1,C_m2,C_m3,C_m4分别为电阻率、速度、密度和磁化率的模型协方差矩阵；α为阻尼因子，α₁,α₂,α₃,α₄分别为大地电磁、地震初至波走时、重力和磁法的阻尼因子，f(m)表示正演响应，f₁(m),f₂(m),f₃(m),f₄(m)分别为大地电磁、地震初至波走时、重力和磁法的正演响应，为梯度，τ为交叉梯度，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度；为数据拟合项、为模型平滑约束项；T和-1分别为矩阵转置和矩阵求逆；

S2：建立电阻率、速度、密度和磁化率的初始模型：

在笛卡尔坐标系下沿x，z二个坐标轴分别将初始模型空间划分成为N_x，N_z个小的矩形网格单元，间距为Δx(i)(i＝1,...,N_x)，Δz(i)(i＝1,...,N_z)，初始模型均采用均匀半空间模型；

S3：正演计算和雅克比矩阵的求取：

采用Wanamaker(1987)提出的三角形剖分有限元算法进行大地电磁正演计算；采用了Singh(2002)提出的对任意多边形棱柱进行重磁异常正演计算；采用Colin Zelt(1998)改进的有限差分方法求解Eikonal方程进行地震初至波走时正演计算；

通过公式(3)求取雅克比矩阵：

其中，A为正演响应f(m)的雅克比矩阵；

S4：计算交叉梯度偏导数和拉格朗日算子：

通过公式(4)计算交叉梯度偏导数

其中，B为交叉梯度函数τ的导数，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度，i,j＝1,2,3,4；

拉格朗日算子计算公式如下：

首先将目标函数中的交叉梯度约束项通过拉格朗日算子法加入到目标函数中，如式(5)：

其中，Ψ为加入约束条件的目标函数，Γ为拉格朗日算子，

然后，对正演响应f(m)和交叉梯度约束条件τ(m)进行泰勒级数展开，如式(6)：

将式(6)代入到式(5)中，得式(7)：

其中，

对公式(7)求极值推导出模型改变量的表达式如式(8)：

Δm＝N^-1·n-N^-1·B^T·Γ (8)

其中：

将公式(8)带入到公式(1)中的交叉梯度约束项中，求取拉格朗日算子Γ，如式(11)：

S5：计算获得模型空间迭代模型：

将式(11)得到的拉格朗日算子反代入到式(8)中，得到模型空间模型改变量式(12)：

Δm_m＝N^-1·n-N^-1·B^T·Γ (12)

之后，获得模型空间迭代模型如式(13)：

m_m＝m₀+Δm_m (13)

S6：计算获得数据空间迭代模型：

对模型空间模型改变量式(12)进行恒等式变换

将式(14)和式(15)带入到式(8)中，得到数据空间模型改变量

Δm_s＝N^-1·n-N^-1·B^T·Γ (16)

之后，获得数据空间迭代模型如式(17)：

m_s＝m₀+Δm_s (17)

S7：对获得的数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3，

其中，拟合差通过式(18)计算获得

其中，N为观测数据个数。