CN104537714A - 磁共振与瞬变电磁空间约束联合反演方法 - Google Patents

Abstract

本发明涉及一种高精度的磁共振与瞬变电磁空间约束联合反演方法,基于连分式展开数值积分等价变换法与数字滤波法相结合的正演方法，兼顾磁场计算精度与速度，提出初始模型与最简模型，设计旋转系数矩阵方程，仿真不同坡体倾角与探测剖面角度变化时激发场分布。依据倾角参数变化下的正演结果，结合分层全采集数据反演方法，构建联合反演目标函数与迭代方程和三维大地模型，引入先验信息约束矩阵、粗糙度矩阵及层厚度约束矩阵等空间约束矩阵构建反演迭代方程组，基于预处理共轭梯度法的大型矩阵，实现了电阻率、层厚度、含水量、弛豫时间的光滑连续约束，解决了反演非唯一性问题，解决了海量数据参与空间约束反演大型矩阵优化与计算的难题。

Description

磁共振与瞬变电磁空间约束联合反演方法

技术领域：

本发明涉及一种地球物理探测方法及其数据解释，在地势起伏、地层电性以及大地模型连续变化的条件下，利用磁共振与瞬变电磁空间约束联合反演的方法进行数据反演解释。

背景技术：

水资源短缺严重制约我国西部经济发展，而复杂的地质环境中赋存的裂隙水、岩溶水增加了勘探难度。国际上采用的地震、电磁法和高密度电法等技术能够根据空间构造和电阻率信息间接判断水源地，但不具备定性和定量分析能力。

磁共振测深(MRS)是世界上唯一直接探测地下水的地球物理新方法，具有定性、定量分析水含量的特点，但MRS相位信息数据量少，易受硬件系统及环境干扰影响，导致实际数据反演得到的电阻率分布并不准确。瞬变电磁法(TEM)利用不接地的回线发射一次磁场，在一次磁场的间歇期间接收由地下良导体受激励引起的涡流所产生的二次感应磁场，能够获取地下几百米范围内的电阻率信息。

Behroozmand在GEOPHYSICS[2012,77(4),191-200]“Improvementin MRS parameter estimation by joint and laterally constrainedinversion of MRS and TEM data”中采用MRS与TEM联合反演，引入电阻率信息，适用于电性复杂的地质环境解释，将联合反演与LCI横向约束结合进行地下水探测，判定含水层位置及分布情况。

万玲在Chinese Journal of Geophysics[2013,56(11)]“基于自适应遗传算法的MRS-TEM联合反演方法研究”中通过电阻率分布信息对含水量反演过程的实时修正，反演算法采用自适应遗传算法(AGA)进行，基于繁殖规则，动态调整交叉概率和变异概率，解决了标准遗传算法易未成熟收敛而难以得到全局最优解问题。但以上方法没有考虑相邻测线与测点在空间范围内的相关性，也没有将探测剖面倾斜角度等因素引入反演，致使误差传递，而且将反演扩展至二维、三维也面临着数据矩阵过于庞大的计算难题。

发明内容：

本发明的目的在于针对上述现有技术的不足，提供一种适用于磁共振与瞬变电磁空间约束联合反演数据解释的方法。主要包括在地磁场环境下，设计旋转系数矩阵建立探测倾角模型，将地磁场方向参数引入磁共振反演；采用三维大地模型中磁共振与瞬变电磁空间约束联合反演方案，实现电阻率、层厚度、含水量、弛豫时间的光滑连续约束，解决反演的非唯一性；基于相邻测点空间范围内的关联性，同时反演测区全部测线测点的磁共振与瞬变电磁数据，提出预处理共轭梯度大型矩阵计算技术，实现联合反演空间约束。

本发明的目的是通过以下技术方案方式实现的：

用等价变换法与数字滤波法相结合进行回线源任意接收位置磁共振激发场及瞬变电磁接收场快速正演；通过设计旋转系数矩阵建立探测倾角模型，将参数引入磁共振反演；由于数据噪声大小与反演权系数成反比，即1/噪声，根据实测MRS和TEM数据确定MRS-TEM联合反演权系数，自适应调整正则化参数，构建MRS-TEM全数据联合反演目标函数；结合地下空间网格化，建立三维大地模型；引入先验信息约束矩阵、粗糙度矩阵及层厚度空间约束矩阵构建反演迭代方程组，实现模型三维方向约束矩阵。

磁共振与瞬变电磁空间约束联合反演方法，包括以下步骤：

a、正演计算三维核函数K_3D(q；r)，代表了位置r处MRS信号的灵敏度；

b、利用等价变换法和数字滤波法相结合的方法进行回线源任意接收位置磁共振激发场及瞬变电磁接收场快速正演计算；

c、考虑地磁场方向和天线方向，建立初始模型，在地磁倾角D与地磁偏角I影响下，激发场的表达式为：

${\stackrel{\→}{B}}_T^{\′\′}=R_I{R}_D{\stackrel{\→}{B}}_T$

其中，为任意方向发射的激发场矢量，R_I，R_D分别为角度矩阵；

简化MRS激发场计算方法，定义地磁场方向指向正北，地磁倾角D＝0°，地磁偏角I＝0°时为初始模型，首先以z轴为旋转轴，x轴和y轴沿顺时针旋转角度D，记为坐标系x′y′z′，地磁场方向转为其次以y′轴为旋转轴，x′轴和z′轴沿顺时针旋转角度I，记为坐标系x″y″z″，此时地磁场方向为

d、考虑剖面角度与坡体倾角，建立最简模型，天线的法向倾角α(剖面角度)和偏角β(坡体倾角)，分别表示北向东偏离的角度和水平向地下倾斜的角度，激发场的表达式为：

${\stackrel{\→}{B}}_T^{\′\′\′}=R_{\mathrm{\β}}{R}_{\mathrm{\α}}{\stackrel{\→}{B}}_T$

其中，为任意方向发射的激发场矢量，R_β，R_α分别为角度矩阵；

天线的法向倾角α(剖面角度)和偏角β(坡体倾角)，分别表示北向东偏离的角度和水平向地下倾斜的角度，当α＝0，β＝0时，天线垂向放置，法向方向指向x轴，定义为最简模型，采用步骤c中的旋转方式进行坐标旋转。

e、由于数据噪声大小与反演权系数成反比，即1/噪声，根据实测MRS和TEM数据确定MRS-TEM联合反演权系数，自适应调整正则化参数，构建MRS-TEM全数据联合反演目标函数；

f、根据实际地质情况及岩石电性特征确定模型参数搜索空间，构建三维光滑大地模型；

g、将测区地磁倾角、地磁偏角、坡体倾角和剖面倾角引入旋转系数矩阵进行计算，得到测区的激发场表达式：

h、抽取测区大地模型，进行MRS、TEM快速正演计算，获得模型数据d(m_mod)，判断数据拟合误差是否小于设定误差，若是，则输出反演结果并快速成像；若不是，则进行空间约束矩阵计算；

i、进行空间约束矩阵计算，包括先验信息约束矩阵计算、粗糙度矩阵计算及层厚度约束矩阵计算；

j、构建反演迭代方程组，并运用预处理共轭梯度法解方程，得到新模型参数；

k、计算新模型参数的磁共振、瞬变电磁响应d(m_n+1)，重复执行上述反演过程，直到数据拟合误差小于给定反演拟合误差，然后输出反演结果并快速成像。

有益效果：本发明公开的磁共振与瞬变电磁空间约束联合反演方法适用于地势起伏、地层电性复杂多变的二维、三维地下水探测，设计旋转系数矩阵建立探测倾角模型，将地磁场方向参数引入磁共振反演，为复杂地形的含水层准确解释提供理论依据；采用三维大地模型中磁共振与瞬变电磁空间约束联合反演方案，实现电阻率、层厚度、含水量、弛豫时间的光滑连续约束，解决反演的非唯一性问题；由于全体采集数据集参与反演，提出基于预处理共轭梯度法的大型矩阵计算技术，解决海量数据参与空间约束反演的大型矩阵优化与计算难题，现实意义巨大。

附图说明：

图1为初始模型、最简模型及坐标旋转

  (a)为地磁场方向坐标旋转

  (b)为天线方向坐标旋转

X(North)正北方向

Y(East)正东方向

Z(Down)垂直向下方向

I：地磁偏角

D：地磁倾角

α：天线的法向倾角

β：天线的法向偏角

n：探测线圈法线单位方向向量

b：地磁场单位方向向量

B_T：激发场场强矢量

图2为测线方向大地模型示意图

x代表所分层数，ρ、thk、W和分别代表各层介质电阻率、层厚度、含水量和平均弛豫时间。

图3为磁共振与瞬变电磁空间约束联合反演流程图。

具体实施方式：

磁共振与瞬变电磁空间约束联合反演是一种适用于地势起伏、地层电性复杂多变的二维、三维地下水探测的方法，其将地磁场方向参数引入磁共振反演，准确计算激发(接收)磁场分量以及地磁场方向向量，使结果更加准确、可靠；由于数据噪声大小与反演权系数成反比，即1/噪声，根据实测MRS和TEM数据确定MRS-TEM联合反演权系数，构建MRS-TEM全数据联合反演目标函数，引入先验信息约束矩阵、粗糙度矩阵及层厚度约束矩阵等空间约束矩阵构建反演迭代方程组，并基于预处理共轭梯度法的大型矩阵计算技术，解决海量数据参与空间约束反演的大型矩阵优化与计算难题，对非层状含水层精确反演和准确解释水文地质参数，解决地势起伏、地层电性复杂多变地区水资源的勘察具有重要意义。

下面结合附图和实施例对本发明作进一步详细说明：

磁共振与瞬变电磁空间约束联合反演方法，包括下列顺序和步骤：

a、正演计算三维核函数K_3D(q；r)，代表了位置r处MRS信号的灵敏度。

b、利用等价变换法和数字滤波法相结合的方法进行回线源任意接收位置磁共振激发场及瞬变电磁接收场快速正演计算，兼顾计算速度与精度，是联合反演顺利进行的基础。

c、如图1(a)所示，简化MRS激发场计算方法，定义地磁场方向指向正北，地磁倾角D＝0°，地磁偏角I＝0°时为初始模型，首先以z轴为旋转轴，x轴和y轴沿顺时针旋转角度D，记为坐标系x′y′z′，地磁场方向转为其次以y′轴为旋转轴，x′轴和z′轴沿顺时针旋转角度I，记为坐标系x″y″z″，此时地磁场方向为

考虑地磁场方向和天线方向，建立初始模型，在地磁倾角D与地磁偏角I影响下，激发场的表达式为：

${\stackrel{\→}{B}}_T^{\′\′}=R_I{R}_D{\stackrel{\→}{B}}_T$

其中，为任意方向发射的激发场矢量，R_I，R_D分别为角度矩阵；

d、如图1(b)所示，天线的法向倾角α(剖面角度)和偏角β(坡体倾角)，分别表示北向东偏离的角度和水平向地下倾斜的角度，当α＝0，β＝0时，天线垂向放置，法向方向指向x轴，定义为最简模型，采用步骤c中的旋转方式进行坐标旋转。

考虑剖面角度与坡体倾角，建立最简模型，天线的法向倾角α(剖面角度)和偏角β(坡体倾角)，分别表示北向东偏离的角度和水平向地下倾斜的角度，激发场的表达式为：

${\stackrel{\→}{B}}_T^{\′\′\′}=R_{\mathrm{\β}}{R}_{\mathrm{\α}}{\stackrel{\→}{B}}_T$

其中，为任意方向发射的激发场矢量，R_β，R_α分别为角度矩阵；

e、分别进行MRS和TEM野外数据全波采集，由于数据噪声大小与反演权系数成反比，即1/噪声，根据实测MRS和TEM数据确定MRS-TEM联合反演权系数，自适应调整正则化参数，构建MRS-TEM全数据联合反演目标函数。

f、根据实际地质情况及岩石电性特征确定模型参数搜索空间，构建三维光滑大地模型。

g、将测区地磁倾角、地磁偏角、坡体倾角和剖面倾角引入旋转系数矩阵进行计算，得到测区的激发场表达式。

h、抽取测区大地模型，进行MRS、TEM快速正演计算，获得模型数据d(m_mod)，判断数据拟合误差是否小于设定误差，若是，则输出反演结果并快速成像；若不是，则进行空间约束矩阵计算。

i、进行空间约束矩阵计算，包括先验信息约束矩阵计算、粗糙度矩阵计算及层厚度约束矩阵计算。

j、构建反演迭代方程组，并运用预处理共轭梯度法解方程，得到新模型参数。

k、计算新模型参数的磁共振、瞬变电磁响应d(m_n+1)，重复执行上述反演过程，直到数据拟合误差小于给定设定误差，然后输出反演结果并快速成像。

整个反演流程如图3所示。

步骤b中所述方法：

位于地表半径为a的水平圆回线源在第i层介质中的点(r,z)处磁场垂直分量、径向分量表达式为：

$\begin{array}{c}{B}_{\mathrm{iz}}(r,z)=\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}[a_i{e}^{-u_{i}z}+b_i{e}^{u_{i}z}]u_i{J}_0\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\end{array}$

$\begin{array}{c}{B}_{\mathrm{ir}}(r,z)=\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}[a_i{e}^{-u_{i}z}-b_i{e}^{u_{i}z}]u_i{J}_1\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\end{array}$

其中I为回线源中的电流强度，λ是积分变量，k_i是复波数且ω是线圈中的谐变电流角频率，μ₀为真空磁导率，ρ_i为第i层电阻率，r为收发距，J₀(x)与J₁(x)为第一类的0阶和1阶Bessel函数，a_i和b_i是与整个地电断面有关的待定系数，可采用循环递推算法快速获取。

由于公式中积分项震荡且不收敛，采用“等价变换法”改写为如下形式：

$\begin{array}{c}{B}_{\mathrm{iz}}(r,z)=\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}[a_i{e}^{-u_{i}z}+b_i{e}^{u_{i}z}]u_i{J}_0\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\\ =\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}\left\{\right[a_i{e}^{-u_{i}z}+b_i{e}^{u_{i}z}]\mathrm{\λ}-\mathrm{\λ}[a_i{e}^{-\mathrm{\λz}}+b_i{e}^{\mathrm{\λz}}\left]\right\}{J}_0\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\\ +\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}\left\{\mathrm{\λ}\right[a_i{e}^{-\mathrm{\λz}}+b_i{e}^{\mathrm{\λz}}\left]\right\}{J}_0\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\end{array}$

$\begin{array}{c}{B}_{\mathrm{ir}}(r,z)=\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}[a_i{e}^{-u_{i}z}-b_i{e}^{u_{i}z}]u_i{J}_1\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\\ =\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}\left\{\right[a_i{e}^{-u_{i}z}-b_i{e}^{u_{i}z}]u_i-\mathrm{\λ}[a_i{e}^{-\mathrm{\λz}}-b_i{e}^{\mathrm{\λz}}\left]\right\}{J}_1\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\\ +\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}\left\{\mathrm{\λ}\right[a_i{e}^{-\mathrm{\λz}}-b_i{e}^{\mathrm{\λz}}\left]\right\}{J}_1\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\end{array}$

上两式前一项均采用数字滤波法计算，此处不再赘述。后一项均采用等价关系计算：

${\∫}_0^{\∞}{e}^{-\mathrm{at}}{t}^{\mathrm{\μ}-v}{J}_{\mathrm{\μ}}\left(\mathrm{bt}\right)J_v\left(\mathrm{ct}\right)\mathrm{dt}=\frac{{\left(\frac{1}{2}b\right)}^{\mathrm{\μ}}{\left(\frac{1}{2}c\right)}^v\mathrm{\Γ}(2\mathrm{\μ}+1)}{\mathrm{\Γ}(v+\frac{1}{2})\mathrm{\Γ}\left(\frac{1}{2}\right)}{\∫}_0^{\mathrm{\π}}\frac{{\sin }^{2v}\mathrm{\φd\φ}}{{(a^2+2\mathrm{iac}{\cos }^2\mathrm{\φ}+b^2)}^{\mathrm{\μ}+\frac{1}{2}}}$

${\∫}_0^{\∞}{e}^{-\mathrm{at}}{t}^{\mathrm{\μ}-1}{J}_{\mathrm{\μ}}\left(\mathrm{bt}\right)J_v\left(\mathrm{ct}\right)\mathrm{dt}=\frac{{\left(\mathrm{bc}\right)}^v\mathrm{\Γ}(2\mathrm{\μ}+1)}{\mathrm{\Γ}(v+\frac{1}{2})\mathrm{\Γ}\left(\frac{1}{2}\right)}{\∫}_0^{\mathrm{\π}}\frac{{\sin }^{2v}\mathrm{\φd\φ}}{{(a^2+2\mathrm{iac}{\cos }^2\mathrm{\φ}+b^2)}^{\mathrm{\μ}+\frac{1}{2}}}$

将层状大地上方，地面任意接收位置磁场垂直分量计算公式改写为下式，再利用等价变换法和数字滤波法相结合的方法计算。

$\begin{array}{c}{B}_z=\frac{{\mathrm{\μ}}_0\mathrm{Ia}}{2}{\∫}_0^{\∞}(1+r_{\mathrm{TE}})\mathrm{\λ}{J}_0\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\\ =\frac{\mathrm{Ia}}{2}{\∫}_2^{\∞}[(1+r_{\mathrm{TE}})\mathrm{\λ}-\mathrm{\λ}]J_0\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}+\frac{\mathrm{Ia}}{2}{\∫}_0^{\∞}\mathrm{\λ}{J}_0\left(\mathrm{\λr}\right)J_1\left(\mathrm{\λa}\right)\mathrm{d\λ}\end{array}$

其中r_TE为反射系数。

步骤c、d中所述方法：

激发场B_T在xyz坐标系下分解为：

${\stackrel{\→}{B}}_T=B_x{\stackrel{\→}{e}}_x+B_y{\stackrel{\→}{e}}_y+B_z{\stackrel{\→}{e}}_z$

经过步骤c的两次旋转后激发场变化为:

${\stackrel{\→}{B}}_T^{\′\′}=B_x^{\′\′}{\stackrel{\→}{e}}_x+B_y^{\′\′}{\stackrel{\→}{e}}_y+B_z^{\′\′}{\stackrel{\→}{e}}_z$

仍然采用坐标系xyz，只是各分量的幅度发生变化。B_T″与B_T的关系可用下式表示为：

${\stackrel{\→}{B}}_T^{\′\′}=R_I{R}_D{\stackrel{\→}{B}}_T---\left(1\right)$

其中，

$R_D=\left[\begin{array}{ccc}\cos D& \sin D& 0\\ -\sin D& \cos D& 0\\ 0& 0& 1\end{array}\right]$

$R_I=\left[\begin{array}{ccc}\cos I& 0& \sin I\\ 0& 1& 0\\ -\sin I& 0& \cos I\end{array}\right]$

经过步骤d的两次旋转后激发场变化为:

${\stackrel{\→}{B}}_T^{\′\′\′}=R_{\mathrm{\β}}{R}_{\mathrm{\α}}{\stackrel{\→}{B}}_T$

其中，

$R_{\mathrm{\α}}=\left[\begin{array}{ccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 1\end{array}\right]$

$R_{\mathrm{\β}}=\left[\begin{array}{ccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}\\ 0& 1& 0\\ -\sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right]$

经过逆旋转，得到xyz坐标系下激发场表达式：

${\stackrel{\→}{B}}_T=R_{\mathrm{\α}}^{-1}{R}_{\mathrm{\β}}^{-1}{\stackrel{\→}{B}}_T^{\′\′\′}---\left(2\right)$

综合表达式(1)和(2)，得到任意地磁场方向和天线方向的激发场表达式：

${\stackrel{\→}{B}}_T^{\mathrm{sim}}=R_I{R}_D{R}_{\mathrm{\α}}^T{R}_{\mathrm{\β}}^T{\stackrel{\→}{B}}_T^{\mathrm{ini}}$

步骤e中所述设计目标函数方法：

磁共振-瞬变电磁联合反演数据空间和模型空间分别为：

d＝{[log(V_i,j),log(ρ_ak)]_L}^T,i＝1,N_q；j＝1,N_MRS；k＝1,N_TEM；L＝1,N_soundings

$m={\left\{{[\log \left({\mathrm{\ρ}}_i\right),\log \left({\mathrm{thk}}_j\right),\log \left(W_i\right),\log \left(T_{2i}^{*}\right)]}_L\right\}}^T,i=1,N_{\mathrm{layers}}-1;L=1,N_{\mathrm{soundings}}$

其中V为磁共振响应，是激发脉冲矩的函数，ρ_a为瞬变电磁响应计算的视电阻率，是取样时间的函数，N_q、N_MRS、N_TEM、N_soundings分别代表激发脉冲矩数，磁共振数据个数，瞬变电磁抽道数及测点数，ρ、thk、W和分别代表各层介质电阻率、层厚度、含水量和平均弛豫时间。

根据磁共振、瞬变电磁联合反演实际情况，设计目标函数：

$\mathrm{\Φ}\left(m\right)\≅{\mathrm{\ω}}_{\mathrm{MRS}}{\left|\right|V_{\mathrm{obs}}-V(W,\mathrm{thk},T_2^{*})\left|\right|}_2+{\mathrm{\ω}}_{\mathrm{TEM}}{\left|\right|{\mathrm{\ρ}}_{\mathrm{aobs}}-{\mathrm{\ρ}}_a(\mathrm{\ρ},\mathrm{thk})\left|\right|}_2+{\mathrm{\η}}^2{\left|\right|\mathrm{Cm}\left|\right|}_2$

其中V_obs为MRS实测数据，为理论模型MRS正演计算结果，ρ_aobs为TEM实测数据计算的视电阻率，ρ_a(ρ,thk)为理论模型计算的视电阻率，η²为正则化参数，C为平滑度矩阵。为解决这一优化问题，将目标函数表示成迭代格式：

m_n+1＝m_n+δm

其中，n是当前迭代次数，新的模型增量δm可按下述方法获得：

磁共振、瞬变电磁响应均为模型参数的非线性函数，在初始模型m_mod邻域内进行泰勒展开，忽略高次项，得到：

d_obs≌d(m_mod)+G(m_true-m_mod)+e_obs

简化为：

Gδm_true＝δd_obs+e_obs    (3)

其中d_obs为实测数据，d(m_mod)为模型m_mod的正演计算结果，为雅克比矩阵，矩阵中元素为数据对模型参数的导数，m_true为真实模型，e_obs代表数据观测误差。

步骤i中所述方法：

首先，引入模型先验信息约束矩阵：

Iδm_true＝δm_prior+e_prior    (4)

其中δm_prior＝m_prior-m_ref，Im_true＝m_prior+e_prior，e_prior代表模型先验信息误差，I表示单位矩阵，对角线元素为模型参数。

其次，引入模型粗糙度矩阵：

R_pδm_true＝δr_p+e_rp    (5)

其中p代表初始模型参数，e_rp为约束矩阵误差，δr_p＝-R_pm_ref。

$R_p=\left[\begin{array}{cccccccccc}1& 0& ...& 0& -1& 0& ...& 0& 0& 0\\ 0& 1& 0& ...& 0& -1& 0& ...& 0& 0\\ & & ......& & & & ......& & & \\ 0& 0& 0& ...& 0& 1& 0& ...& 0& 0\end{array}\right]$

最后，引入层厚度约束矩阵：

R_hδm_true＝δr_h+e_rh    (6)

其中h代表各层深度，e_rh为约束矩阵误差，δr_h＝-R_hm_ref。

采用一阶偏导数计算层厚度约束矩阵，如图2所示，以沿测线方向为例，得到R_h矩阵元素为：

步骤j中所述构建联合反演迭代方程组：

综合(3)、(4)、(5)和(6)式得到空间约束磁共振-瞬变电磁联合反演迭代方程组，

$\left[\begin{array}{c}G\\ I\\ R_P\\ R_h\end{array}\right]\·\mathrm{\δ}{m}_{\mathrm{true}}=\left[\begin{array}{c}\mathrm{\δ}{d}_{\mathrm{obs}}\\ \mathrm{\δ}{m}_{\mathrm{prior}}\\ \mathrm{\δ}{r}_p\\ \mathrm{\δ}{r}_h\end{array}\right]+\left[\begin{array}{c}{e}_{\mathrm{obs}}\\ {}^{}{e}_{\mathrm{prior}}\\ e_{\mathrm{rp}}\\ e_{\mathrm{rh}}\end{array}\right]$

简化为：

G'δm_true＝δd+e'    (7)

模型参数更新为：

m_n+1＝m_n+δm_true

其中λ为马奎特阻尼因子。得到新模型参数后，计算新模型参数的磁共振、瞬变电磁响应d(m_n+1)，重复执行上述反演过程，直到数据拟合误差小于给定反演拟合误差。

步骤j中所述预处理共轭梯度法：

对于空间约束全数据磁共振与瞬变电磁联合反演方法，在反演迭代过程中需要求解大型矩阵方程式(式(7))，尤其是雅克比矩阵G的计算。设其中G_ρ，G_w和分别代表数据对电阻率，含水量和弛豫时间的雅克比矩阵，则

$G_w(n,l;m)=\frac{\∂V(n,l)}{\∂w\left(m\right)}{=K}_{3D}(n,m)e^{-t\left(l\right)/T_2^{*}\left(m\right)}$

$G_{T_2^{*}}(n,l;m)=\frac{\∂V(n,l)}{\∂T_2^{*}\left(m\right)}=\frac{t\left(l\right)}{T_2^{*}{\left(m\right)}^2}{K}_{3D}(n,m)e^{-t\left(l\right)/T_2^{*}\left(m\right)}w\left(m\right)$

当大型矩阵条件数很大时，共轭梯度法收敛速度很慢。预处理共轭梯度法(PCG)基于共轭梯度法，对系数矩阵做预处理，降低系数矩阵的条件数(条件数越小收敛速度越快)，具有更高的迭代效率和收敛性。PCG法引入预处理因子矩阵M来降低系数矩阵的条件数，加快收敛速度。将公式(7)写成方程Ax＝b形式，设预处理因子矩阵M＝W^TW，则方程变为：

(W^-TAW^-1)Wx＝W^-Tb

或A'x'＝b'    (8)当矩阵M与A结构近似时，W^-TAW^-1或A'的条件数较小，收敛速度加快。求解方程组(8)的预处理共轭梯度法可描述如下：

1)取x的初始估计x₀，计算r₀＝b'-A'x₀，p₀＝M^-1r₀；

2)对迭代步数k＝1，2，…，进行如下计算

① ${\mathrm{\&alpha;}}_k=\frac{<r_k,M^{-1}{r}_k>}{<p_k,A^{\&prime;}{p}_k>}$

②x_k+1＝x_k+α_kp_k

③r_k+1＝r_k-α_kA'p_k

④ ${\mathrm{\&beta;}}_k=\frac{<r_{k+1},M^{-1}{r}_{k+1}>}{<r_k,M^{-1}{r}_k>}$

⑤p_k+1＝M^-1r_k+1+β_kp_k

3)当满足||r_k+1||/||b||＜ε时，迭代停止。

其中，x_k为每次迭代得到的方程的解；r_k为每次迭代解的残差；ε为迭代误差的判断值；＜h_k,r_k＞表示h_k,r_k的内积，其余变量为每次迭代的过程变量。

Claims (7)

1.一种磁共振与瞬变电磁空间约束联合反演方法，其特征在于，用等价变换法与数字滤波法相结合进行回线源任意接收位置磁共振激发场及瞬变电磁接收场快速正演；通过设计旋转系数矩阵建立探测倾角模型，将参数引入磁共振反演；由于数据噪声大小与反演权系数成反比，即1/噪声，根据实测MRS和TEM数据确定MRS-TEM联合反演权系数，自适应调整正则化参数，构建MRS-TEM全数据联合反演目标函数；结合地下空间网格化，建立三维大地模型；引入先验信息约束矩阵、粗糙度矩阵及层厚度空间约束矩阵构建反演迭代方程组，实现模型三维方向约束矩阵。

2.按照权利要求1所述的磁共振与瞬变电磁空间约束联合反演方法，其特征在于，将层状大地上方，地面任意接收位置磁场垂直分量计算公式改写为下式，

其中r_TE为反射系数；

再利用等价变换法和数字滤波法相结合的方法进行快速正演计算。

3.按照权利要求1所述的磁共振与瞬变电磁空间约束联合反演方法，其特征在于，将地磁场方向参数引入磁共振反演，设计旋转系数矩阵方程，得到任意地磁场方向和天线方向的激发场表达式：

实现磁共振激发场精确正演计算。

4.按照权利要求1所述的磁共振与瞬变电磁空间约束联合反演方法，其特征在于，将磁共振探测全部数据引入准三维反演，由于数据噪声大小与反演权系数成反比，即1/噪声，根据实测MRS和TEM数据确定MRS-TEM联合反演权系数，构建MRS-TEM全数据联合反演目标函数：

。

5.按照权利要求1所述的磁共振与瞬变电磁空间约束联合反演方法，其特征在于，引入先验信息约束矩阵、粗糙度矩阵及层厚度约束矩阵等空间约束矩阵构建反演迭代方程组。

6.按照权利要求1所述的磁共振与瞬变电磁空间约束联合反演方法，其特征在于，应用预处理共轭梯度法PCG计算反演迭代方程组中的大型矩阵方程式， 使算法具有更高的迭代效率和收敛性。

7.按照权利要求1所述的磁共振与瞬变电磁空间约束联合反演方法，其特征在于，包括以下步骤：

a、正演计算三维核函数K_3D(q；r)，代表了位置r处MRS信号的灵敏度；

b、利用等价变换法和数字滤波法相结合的方法进行回线源任意接收位置磁共振激发场及瞬变电磁接收场快速正演计算，兼顾计算速度与精度，是联合反演顺利进行的基础；

c、考虑地磁场方向和天线方向，建立初始模型，在地磁倾角D与地磁偏角I影响下，激发场的表达式为：

其中，为任意方向发射的激发场矢量，R_I，R_D分别为角度矩阵；

d、考虑剖面角度与坡体倾角，建立最简模型，天线的法向倾角α(剖面角度)和偏角β(坡体倾角)，分别表示北向东偏离的角度和水平向地下倾斜的角度，激发场的表达式为：

其中，为任意方向发射的激发场矢量，R_β，R_α分别为角度矩阵；

e、分别进行MRS和TEM野外数据全波采集，由于数据噪声大小与反演权系数成反比，即1/噪声，根据实测MRS和TEM数据确定MRS-TEM联合反演权系数，自适应调整正则化参数，构建MRS-TEM全数据联合反演目标函数；

f、根据实际地质情况及岩石电性特征确定模型参数搜索空间，构建三维光滑大地模型；

g、将测区地磁倾角、地磁偏角、坡体倾角和剖面倾角引入旋转系数矩阵进行计算，得到测区的激发场表达式：

h、抽取测区大地模型，进行MRS、TEM快速正演计算，获得模型数据d(m_mod)，判断数据拟合误差是否小于设定误差，若是，则输出反演结果并快速成像；若不是，则进行空间约束矩阵计算；

i、进行空间约束矩阵计算，包括先验信息约束矩阵计算、粗糙度矩阵计算及层厚度约束矩阵计算；

j、构建反演迭代方程组，并运用预处理共轭梯度法解方程，得到新模型参数；

k、计算新模型参数的磁共振、瞬变电磁响应d(m_n+1)，重复执行上述反演过程，直到数据拟合误差小于设定误差，然后输出反演结果并快速成像。