CN108680964A - 一种基于结构约束的归一化重磁电震联合反演方法 - Google Patents

Abstract

本发明公开了一种基于结构约束的归一化重磁电震联合反演方法，包括如下步骤：S1：建立归一化联合反演目标函数；S2：建立初始模型；S3：正演计算和雅克比矩阵求取；S4：计算归一化交叉梯度偏导数和拉格朗日算子；S5：计算获得迭代模型；S6：对获得的迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出归一化联合反演模型，否则执行S3。该基于结构约束的归一化重磁电震联合反演方法克服了单一方法的局限性、减小了地球物理反演多解性、不依赖于岩石物性关系，同时还可以有效的避免由于数量级相差较大的物性参数直接耦合反演所导致的错误反演结果。

Description

一种基于结构约束的归一化重磁电震联合反演方法

技术领域

本发明涉及重力、磁法、大地电磁和地震初至波走时四种地球物理方法同时处理的领域，特别提供了一种基于结构约束的归一化重磁电震联合反演方法。

背景技术

根据岩石物性参数的不同，地球物理方法可分为重力、磁法、地震和电法等，重力和磁法具有很好的横向分辨能力，对地质构造分布划分起到了重要的作用，但纵向分辨能力较差，且勘探深度浅；大地电磁测深法(MT)的勘探深度与频率有关，由于频率范围丰富，其勘探深度达到几十或者上百公里，同时具有不受高阻层屏蔽，对良导体反映灵敏等优点，因此在研究地球深部构造等方面具有特殊的优势，但是横向分辨能力较差；对于地震初至波走时方法，具有较高的纵向分辨能力，在划分地下结构时分辨能力敏感，但是在盐丘、高角度断裂存在时容易形成盲区，并且初至波传播深度较浅，导致勘探深度较浅。以上每种单一方法都存在一定的缺陷，每种方法都只能从单一角度来评价地下介质的岩石物理特性，单凭一种地球物理方法很难准确勘探地下结构，因此，为了得到更准确的地下信息，利用多种地球物理参数对同一地下地质体进行综合解释已经成为当今发展的趋势。

联合反演是综合地球物理方法的重要手段和工具，传统的联合反演方法是基于岩石物性关系的耦合的方法，但该方法很难建立准确的岩石物性关系，并不适用于复杂的地质情况，具有一定的局限性，基于空间分布结构性耦合的联合反演方法并不依赖于岩石物性关系，而是强调结构上的相似性，具有更普遍的适用性。

然而，目前，国内针对重力、磁法、大地电磁法和地震初至波走时方法之间较多的进行了二种方法之间的结构耦合研究，较少的考虑更多物性参数之间的耦合问题，同时没有考虑到，由于不同物性参数的量级和单位相差较大，如果直接耦合不同物性参数，会引起的错误反演结果。

发明内容

鉴于此，本发明的目的在于提供一种基于结构约束的归一化重磁电震联合反演方法，以解决现有的联合反演方法未考虑归一化数量级相差较大的不同物性参数耦合问题和无法处理类似于实际情况的复杂模型问题。

本发明提供的技术方案是：一种基于结构约束的归一化重磁电震联合反演方法，包括如下步骤：

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为交叉梯度函数τ的导数；

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

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

其中，Γ为拉格朗日算子，

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

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

其中，

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

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

S5：计算获得迭代模型：

式(9)得到的拉格朗日算子反代入到式(8)中，得到最终的迭代模型改变量式(10)：

最后获得迭代模型如式(11)：

m＝m₀+Δm (11)

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

拟合差如式(12)

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

本发明提供的基于结构约束的归一化重磁电震联合反演方法不依赖于不同物性参数之间的经验关系，在复杂不确定物性参数关系的地区也可以实用。另外，本发明可以同时对重力、磁法、大地电磁和地震初至波走时四种不同物性方法进行同步联合反演，在多交叉梯度约束下，获得多参数多约束的联合反演结果，相比于单独反演和传统联合反演，可以有效地降低反演多解性，提高反演结果的分辨率，同时还可以有效的避免由于数量级相差较大的物性参数直接耦合反演，所导致的错误反演结果，为深部资源勘查提供了有力的技术支持。

附图说明

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

图1为本发明提供的基于结构约束的归一化重磁电震联合反演方法的流程图；

图2为验证本发明提供的方法而设计的地下结构的复杂理论模型图；

图3中(A)列为单独反演模型结果，(B)列为归一化联合反演模型结果，自上而下分别为电阻率、密度、磁化率和速度反演模型结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，但并不用于限定本发明。

如图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为交叉梯度函数τ的导数；

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

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

其中，Γ为拉格朗日算子，

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

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

其中，

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

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

S5：计算获得迭代模型：

式(9)得到的拉格朗日算子反代入到式(8)中，得到最终的迭代模型改变量式(10)：

最后获得迭代模型如式(11)：

m＝m₀+Δm (11)

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

拟合差如式(12)

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

该基于交叉梯度约束的重磁电震二维联合反演方法不依赖于不同物性参数之间的经验关系，在复杂不确定物性参数关系的地区也可以实用。另外，本发明可以同时对重力、磁法、大地电磁和地震初至波走时四种不同物性方法进行同步联合反演，在多交叉梯度约束下，获得多参数多约束的联合反演结果，相比于单独反演和传统联合反演，可以有效地降低反演多解性，提高反演结果的分辨率，同时还可以有效的避免由于数量级相差较大的物性参数直接耦合反演，所导致的错误反演结果，为深部资源勘查提供了有力的技术支持。

下面对本发明的基于结构约束的归一化重磁电震联合反演方法的准确性和有效性进行验证。

设计一个如图2所示的复杂理论模型，相比于简单的模型而言，类似于真实地下结构的复杂理论模型更具有验证算法的说服力。复杂理论模型中包含了一些典型的地质单元体和断层结构，在均匀半空间D中镶嵌了三个异常体，异常体A出露地表，下部为尖锐的阶梯状断裂结构。异常体B为孤立的地壳岩体；在模型深部区域，高阻异常体A和异常体B下面埋藏着低阻异常体C，异常体C左边出现了阶梯状的断裂构造。

接下来，根据上述设计的模型检验不同反演方法是否具有划分横纵界面和断层、识别孤立的地层岩体和深部基底结构的能力。

上述复杂理论模型正演计算可以获得大地电磁(MT)的视相位和视电阻率数据，重力异常和磁异常数据，地震初至波走时数据。将MT的频点人为设定为10个，频率范围在1～1000Hz之间，人为设定9个观测点等距离分布在测线0～6km之间(其中，频率的个数和观测点个数可以根据需要修改)，采用Wanamaker(1987)提出的三角形剖分有限元算法进行正演计算；重力、磁法地表观测点各有30个，重磁正演计算采用了Singh(2002)提出的算法，该算法可以对任意多边形棱柱进行重磁异常计算；地震震源位于地下50m处，共17个，接收点位于两口井中，水平位置分别为1.5km和4.5km，每口井内有10个接收点，间隔为0.2km，采用Colin Zelt(1998)改进的有限差分方法求解Eikonal方程进行正演地震走时计算；重力、磁法、地震走时地下二维等间隔网格剖分都为140×60，MT地下二维网格剖分为174×74，需要在重磁震等间距网格剖分的基础上向外延展，延展网格以不等间隔剖分，重磁震网格剖分区域即为联合反演区域，反演网格剖分为70×30。

对大地电磁、重力、磁法和地震初至波走时法分别进行单独反演和归一化联合反演。所有的反演方法初始模型均采用均匀半空间模型，电阻率值为100Ω·m，剩余密度为0g/cm3，磁化率为0.0001SI，速度为4000m/s。单独反演经历了6次迭代最终达到收敛，大地电磁法、重力、磁法和地震初至波走时法的拟合差分别为0.425442，0.317181，0.397370，0.421690；归一化联合反演同样经历了6次迭代达到收敛，最终的拟合差分别为0.459097，0.378324，0.50919，0.370272。二种反演拟合差均达到了阈值以下。图3为上述复杂理论模型的单独反演结果(A)和归一化联合反演结果(B)，反演结果由上到下分别为电阻率、密度、磁化率、速度模型，在单独反演中，四种地球物理方法都有各自的优势，大地电磁方法相比于重磁法具有较高的纵向分辨能力，单独反演图3(A)中，电阻率模型基本上恢复了真实模型结构，特别在纵向方面，可以不受高阻层的屏蔽，灵敏的反映出基底深部的低阻层的位置，但是在横向上，MT反演结果没有识别深部异常体C左侧阶梯状的断裂边界；重磁方法横向分辨率要高于MT和地震方法，在重磁单独反演结果中，可以分辨出水平方向的异常界面，但是纵向方面分辨能力很差，容易将深部区域的高密度、高磁化率异常体恢复在浅地表区域得到假的地下反演模型，从而对地质解释产生错误；地震方法相比于非震方法具有更好的分辨能力，但是地震射线不可能覆盖地下全部的勘探区域，深部射线分布相比于浅部区域较少，通过得到的反演结果可知，不能识别出异常体C左侧的阶梯结构边界。以上所有方法都存在各自的优缺点，如果对差异很大的单独反演结果进行地质解释，研究者很难准确的划分地质结构和判断地下地层单元。

图3(B)为归一化联合反演结果，该反演结果和单独反演结果之间具有很大的结构差异，由于结构相似性的约束，联合反演得到了结构一致性更高的反演模型结果。其中，MT反演结果在异常体C与围岩的界面处划分清晰明显，提高了深部区域横向分辨率；同时，在射线分布少的深部区域，地震反演结果也明显改善，准确的识别了深部异常体C与围岩的界面，得到了与真实模型基本吻合的地震模型结果；在重磁反演中，交叉梯度约束后的重磁反演很大程度改善了纵向分辨率不高和深部异常体上移的现象，同时可以准确的识别出孤立岩体的位置和大小.综上所述，相比于单独反演，归一化联合反演得到的结果更接近于真实模型，无论在模型几何形态上，还是在异常体物性参数值上，都得到了有效的改善和提高。

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为交叉梯度函数τ的导数；

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

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

其中，Γ为拉格朗日算子，

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

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

其中，

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

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

S5：计算获得迭代模型：

式(9)得到的拉格朗日算子反代入到式(8)中，得到最终的迭代模型改变量式(10)：

最后获得迭代模型如式(11)：

m＝m₀+Δm (11)

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

拟合差如式(12)

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