CN104007477A - 一种地面核磁共振三维反演方法 - Google Patents

Abstract

本发明公开了一种地面核磁共振三维反演方法，首先通过拉直变换将三维正演模型进行降维处理，将三维反演问题抽象为矩阵方程求解模型，然而针对SNMR信号三维反演问题，提出了约束总体GMERR算法，该算法首先通过总体GMERR算法求解含水量的中间解，之后通过约束条件将中间解的取值范围强制转化到0～100％之间，经多次迭代运算求解最优解。本发明在不同信噪比下反演精度均较高，且稳定性好、抗噪能力强。

Description

一种地面核磁共振三维反演方法

技术领域

本发明涉及地面核磁共振技术领域，具体涉及一种地面核磁共振三维反演方法。 

背景技术

地面核磁共振SNMR技术是目前世界上唯一的一种直接找水的物探方法，该项技术已在探测地下水、考古、地下水污染检测等领域得到了一定的应用。近年来，随着该领域专家和学者们的逐渐深入研究，SNMR技术得到了进一步的完善。反演计算含水率是该技术研究过程中的关键环节，而反演准确度和分辨率是衡量反演算法性能的关键指标。其中，一维正反演理论较为成熟，已经相继刊登出多种有效算法，如：改进的模拟退火算法反演，提高了现有反演算法的稳定度和收敛速度；QT反演算法，利用各个激发脉冲矩对应的全部采样点数据进行反演，充分挖掘了接收信号信息，在一定程度上提高了反演精度。但是，由于接收信号呈现近似指数衰减，晚期信号信噪比很低，该方法只适用于高信噪比环境；有的采用了积分门技术接收信号，提高了各个采样点数据的精度，并进行全衰减反演，是对QT反演的一种改进。在二维反演方面，Boucher、Girard和Legchenko等研究了在二维剖面方向上E0-q曲线随地下含水构造的变化趋势，但他们只对二维反演做了定性研究，没有给出具体的二维反演公式。Legchenko等对三维反演做了一定的研究，虽然能在三维空间反演出模型的含水构造，但是由于在三维空间设定的网格尺寸较大，只能粗略的估计出地下含水构造，其反演分辨率有待提高。由于二维、三维反演算法存在运算量大、待求解变量数多、非线性等问题，目前世界上唯一商业版反演软件NUMISPLUS仍采用一维反演，而二维、三维正反演研究仍处于起步阶段。 

发明内容

本发明所要解决的技术问题是现有三维空间反演方法只能粗略的估计出地下含水构造，反演分辨率不高的问题，提供一种地面核磁共振三维反演方法。 

为解决上述问题，本发明是通过以下技术方案实现的： 

一种地面核磁共振三维反演方法，包括如下步骤： 

步骤1，用拉直变换方法将三维正演模型进行降维处理，将其抽象为矩阵方程求解模型，即 

E^*＝K^*·n   ① 

式中，E^*为地面核磁共振信号初始振幅值矩阵，K^*为核函数矩阵，n为待求解的含水量空间分布值； 

步骤2，通过补0使得矩阵方程求解模型中的核函数矩阵K^*变为核函数方阵K，相应的地面核磁共振信号初始振幅值矩阵E^*变为初始振幅值矩阵E； 

步骤3，对核函数方阵K做奇异值分解，以获得核函数方阵K的奇异值σ_i； 

步骤4，根据所得奇异值σ_i计算方阵K的有效秩r^*，并令m＝r^*； 

步骤5，根据下式计算初始标准正交基V₁， 

$\left\{\begin{array}{c}{R}_0=E-K\·n_0\\ V_1=R_0/{\left|\right|R_0\left|\right|}_F\end{array}\right.$    ② 

式中，V₁为初始标准正交基，R₀是待求变量矩阵中第一个元素值，K为核函数方阵，E为初始振幅值矩阵，n₀为设定的初始化含水量值，||·||_F为矩阵的弗罗伯尼范数算子； 

步骤6，根据所得初始标准正交基V₁和核函数方阵K的转置广义克雷洛夫子空间F_m(K^T,V₁)，并对所构建的F_m(K^T,V₁)利用总体阿诺尔迪方法，在弗罗伯尼范数下，生成一组标准正交基V₁,V₂,…,V_m，并令标准正交基集合U_m＝[V₁,V₂,…,V_m]； 

步骤7，求解向量y，其中向量y满足下式 

   ③ 

式中，y为待求解的向量，V_i是标准正交基集合U_m里面的元素，V_i ^T为V_i的转置，i＝1,2,3……,m，K为核函数方阵，K^T为K的转置，R₀是待求变量矩阵中第一个元素值； 

步骤8，求解含水量变量矩阵n_i，其中含水量变量矩阵n_i满足下式 

$n_i\mathrm{abs}(n_0+K^T{U}_m\⊗y)$    ④ 

式中，n_i为待求解的含水量变量矩阵，i＝1,2,3……,m，abs(·)为复数取模算子，n₀为设定的初始化含水量值，K^T为核函数方阵K的转置，U_m为标准正交基集合，y为步骤7所得的向量；y_i为向量y中的元素，V_i是标准正交基集合U_m里面的元素； 

步骤9，在矩阵n_i中，将值大于1的元素置为1，值小于0的元素置为0； 

步骤10，求解矩阵R_i的弗罗伯尼范数算子||R_i||_F，其中矩阵R_i的弗罗伯尼范数算子||R_i||_F满足下式 

||R_i||_F＝||E-K·n_i||   ⑤ 

式中，||R_i||_F为待求解的弗罗伯尼范数算子，K为核函数方阵，E为初始振幅值矩阵，n_i为步骤8所得的矩阵，i＝1,2,3……,m； 

步骤11：如果实际迭代次数N_s达到设定的最大迭代次数N_max时即N_s≥N_max时，或者所得弗罗伯尼范数算子||R_m||_F小于设定的误差阈值θ时即||R_m||_F<θ时，则停止迭代；否则令将实际迭代次数N_s加1即N_s＝N_s+1，并将含水量值n₀设为本次迭代所得含水量n_s即n₀＝n_s，返回步骤5；其中s＝1,2,3……,r，r为方阵K奇异值的个数。 

上述步骤3中，采用下式对核函数方阵K做奇异值分解，即 

$K=U_{l\&times;l}\left[\begin{array}{cc}{\mathrm{\&Lambda;}}_r& O_3\\ O_4& O_5\end{array}\right]V_{l\&times;l}^H$

式中，σ_i为方阵K的奇异值σ₁≥σ₂≥…≥σ_r，r为方阵K的奇异值的个数；O₃、O₄和O₅是0矩阵；U_l×l和是方阵K的l阶酉矩阵且是正交矩阵。 

上述步骤4中，采用下式计算方阵K的有效秩r^*即令 

$\mathrm{\&psi;}\left(s\right)=\frac{\sqrt{{\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2+...+{\mathrm{\&sigma;}}_s}}{\sqrt{{\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2+...+{\mathrm{\&sigma;}}_r}},s=\mathrm{1,2},...,r$     ⑥ 

则r^*为第一个满足ψ(s)>η的s值； 

式中，ψ(s)为求取有效秩r^*的表达式；σ_i为方阵K的奇异值；r为方阵K奇异值的个数；s为方阵K的第s个奇异值；η为设定值。 

上述步骤11中，最大迭代次数N_max设定为方阵K奇异值的个数r，即N_max＝r。 

与现有技术相比，本发明具有如下特点： 

1、通过拉直变换将三维正演模型进行降维处理，将三维反演问题抽象为矩阵方程求解模型，可参考现有的矩阵方程求解算法对其求解，降低了三维反演求解问题的复杂度； 

2、针对SNMR信号三维反演问题，提出了约束总体GMERR算法，该算法首先通过总体GMERR算法求解含水量的中间解，之后通过约束条件将中间解的取值范围强制转化到0～100％之间，经多次迭代运算求解最优解； 

3、在不同信噪比下反演精度均较高，且稳定性好、抗噪能力强，随着信噪比的降低反演结果的方均根值略有增加，在信噪比为-5dB时，其反演精度仍然很高，其反演得到的含水量值的方均根为11.13％，当信噪比降到-10dB 时，反演结果中的含水体与背景区域的边界变得比较模糊，已经不能较好的定位出含水体的空间分布，而世界上唯一的商业版反演软件NUMISPLUS软件只能对信噪比大于2dB的数据进行有效的一维反演。 

附图说明

图1为仿真模型的三维空间模型示意图。 

图2为f_SNR＝5dB时，含水量反演切片图。 

图3为f_SNR＝0dB时，含水量反演切片图。 

图4为f_SNR＝-5dB时，含水量反演切片图。 

图5为f_SNR＝-10dB时，含水量反演切片图。 

具体实施方式

本发明地面核磁共振三维反演方法，包括如下步骤： 

步骤I：在收/发天线共圈模式SNMR三维正反演研究中，SNMR三维反演算法首先需要建立三维水文地质模型：含水量在X-Y-Z三维空间内的分布是非均匀的，p测点处接收信号的初始振幅E₀(p,q)公式为： 

$E_0(p,q)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{K}_{3D}(p,q;x,y,z)\&CenterDot;n(p;x,y,z)\mathrm{dxdydz}---\left(1\right)$

其中，q为激发脉冲矩；K_3D(p,q；x,y,z)为p点处三维核函数；n(p；x,y,z)为p点处三维含水量分布函数。将探测空间区域含水量n(x,y,z)沿着X-Y-Z三维空间方向做离散化处理为： 

$\left\{\begin{array}{c}{x}_{k1-1}\&le;x\&le;x_{k1},k1=\mathrm{1,2},...,N_x\\ y_{k2-1}\&le;y\&le;y_{k2},k2=\mathrm{1,2},...,N_y\\ z_{k3-1}\&le;z\&le;z_{k3},k3=\mathrm{1,2},...,N_z\end{array}\right.---\left(2\right)$

得到L＝N_z×N_x×N_y个含水微元。在探测空间区域内有N_p个测点，每个测点处发射N_q个激发脉冲矩，每个测点覆盖空间范围是l＝n_z×n_x×n_y。测点p_i处的第q_h个激发脉冲矩对应的初始振幅为： 

$\begin{array}{c}{E}_0=(p_i,q_h)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{K}_{3D}(p_i,q_h;x,y,z)\&CenterDot;n(p_i,q_h;x,y,z)\mathrm{dxdydz}\\ =\underset{k1=1}{\overset{n_x}{\mathrm{\&Sigma;}}}\underset{k2=1}{\overset{n_y}{\mathrm{\&Sigma;}}}\underset{k3=1}{\overset{n_z}{\mathrm{\&Sigma;}}}{K}_{3D}(p_i,q_h;x_{k1},y_{k2},z_{k3})\&CenterDot;n(p_i,q_h;x_{k1},y_{k2},z_{k3}){\mathrm{\&Delta;z}}_{k3}{\mathrm{\&Delta;y}}_{k2}{\mathrm{\&Delta;x}}_{k1}\end{array}---\left(3\right)$

对三维矩阵和进行拉直变换，将其转化为一维向量形式，即： 

$\left\{\begin{array}{c}\mathrm{vec}\left(K_{3D}(p_i,q_h)\right)=K_h=[K_{h1},K_{h2},...K_{\mathrm{hl}}]\\ \mathrm{vec}\left(n(p_i,q_h)\right)=n_i=[n_1,n_2,...,n_l]\end{array}\right.---\left(4\right)$

其中，vec(·)为拉直变换算子，则(3)式转化为： 

$E_{\mathrm{hi}}=E_0(p_i,q_h)=K_h\&CenterDot;n_i^T.---\left(5\right)$

因此，SNMR三维反演问题可以抽象为求解矩阵方程E^*＝K^*·n，即 

式中，E为测点处初始振幅值，K为三维核函数经拉直变换后的系数，n为待求解的含水量空间分布值，N_p为探测空间总测点数目，N_q为每个测点处发射的激发脉冲矩数目。 

步骤II：设定起始步数为m，最大迭代次数为N_max，初始迭代次数为N_s＝0，初始化含水量值n₀，通过补0使得矩阵方程中的核函数矩阵K^*变为方阵K，相应的SNMR信号初始振幅值矩阵E^*变为矩阵E，计算R₀＝E-K·n₀，V₁＝R₀/||R₀||_F。其中，n₀为0～1之间的随机数；K^*是N_q行l列矩阵，且l>N_q，所以满足 $K=\left[\begin{array}{c}{K}^{*}\\ O_1\end{array}\right],$ O₁为l-N_q行l列的0矩阵；K^*是N_q行N_P列矩阵，所以满足 $E=\left[\begin{array}{c}{E}^{*}\\ O_2\end{array}\right],$ O₂为l-N_q行N_P列的0矩阵；||·||_F为矩阵的Frobenius范数算子。 

对核函数矩阵K做奇异值分解 

$K=U_{l\&times;l}\left[\begin{array}{cc}{\mathrm{\&Lambda;}}_r& O_3\\ O_4& O_5\end{array}\right]V_{l\&times;l}^H---\left(7\right)$

其中，U_l×l和是正交矩阵，O₃、O₄和O₅是0矩阵， 

,σ₁≥σ₂≥…≥σ_r。 

计算K的有效秩r^*，令 

$\mathrm{\&psi;}\left(s\right)=\frac{\sqrt{{\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2+...+{\mathrm{\&sigma;}}_s}}{\sqrt{{\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_2+...+{\mathrm{\&sigma;}}_r}},s=\mathrm{1,2},...,r---\left(8\right)$

则r^*为第一个满足ψ(s)>η的s值，通常η取接近于1的值，在本发明优选实施例中η＝0.85。令m＝r^*，N_max＝r。 

步骤III：利用总体Arnoldi算法生成广义Krylov子空间F_m(K^T,V₁)在 Frobenius范数下的一组标准正交基V₁,V₂,…,V_m，令U_m＝[V₁,V₂,…,V_m]。 

步骤IV：求解向量y∈C^m×1。y满足方程： 

步骤V：求解矩阵n_m。n_m满足关系式： 

$n_m=\mathrm{abs}(n_0+K^T{U}_m\&CircleTimes;y)---\left(10\right)$

其中，abs(·)为复数取模算子，在矩阵n_m中，将值大于1的元素置1，小于0的元素置0。 

步骤VI：求解||R_m||_F＝||E-K·n_m||，如果||R_m||_F<θ或者N_max≤N_s，则停止迭代；否则令n₀＝n_m，N_s＝N_s+1，返回执行步骤1。其中，θ为误差阈值，文中取值为10^-8。 

本发明首先建立三维水文地质模型：含水量在X-Y-Z三维空间内的分布是非均匀的，p测点处接收信号的初始振幅E₀(p,q)公式为： 

$E_0(p,q)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{K}_{3D}(p,q;x,y,z)\&CenterDot;n(p;x,y,z)\mathrm{dxdydz}---\left(11\right)$

其中q为激发脉冲矩，K_3D(p,q；x,y,z)为p点处三维核函数，n(p；x,y,z)为p点处三维含水量分布函数。将探测空间区域含水量n(x,y,z)沿着X-Y-Z三维空间方向做离散化处理，得到L＝N_z×N_x×N_y个含水微元。在探测空间区域内有N_p个测点，每个测点处发射N_q个激发脉冲矩，每个测点覆盖空间范围是l＝n_z×n_x×n_y。对三维矩阵和进行拉直变换，将其转化为一维向量形式，SNMR三维反演问题抽象为求解矩阵方程。为了高精度快速求解矩阵方程，根据SNMR信号的特性设计了带约束的总体GMERR算法，该算法将初始残差矩阵投影到广义Krylov子空间上，并通过迭代计算以求解误差范数最小解，迭代过程中的含水量中间解取值范围是0～100％。 

下面结合图表对本发明进一步说明。首先建立三维水文地质模型；然后，对其进行正演计算，求出各个测点的E₀-q曲线；之后，对E₀加入一定信噪比的噪声，得到含噪信号E^*；最后，用本发明对上述含噪信号进行三维反演，并用反演得到的含水量值的方均根f_RMS评价反演结果性能，f_RMS公式为： 

$f_{\mathrm{RMS}}=\sqrt{\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{(n_i-n_i^{\#})}^2}\&times;100\%---\left(12\right)$

其中，n_i和分别为理论模型和反演结果中各地质微元的含水量值。信噪比f_SNR公式为： 

$f_{\mathrm{SNR}}=10\lg \left(\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{E}_i^2}{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{(E_i-E_{0i})}^2}\right)---\left(13\right)$

其中，E_i为信号序列，E_0i为噪声序列。 

为了保证实验结果的实用性，仿真模型尽可能的模拟实测环境，模型在X-Y-Z三维空间内的尺寸为400m×400m×90m，剖分的微元大小为20m×20m×5m，三维空间内共有20×20×18个含水微元，该空间内有两个含水体，一个是大小为100m×100m×10m，含水量为40％，位于X：50m-150m、Y：50m-150m、Z：15m-25m；另一个是大小为100m×100m×20m，含水量为70％，位于X：250m-350m、Y：250m-350m、Z：30m-50m，其三维空间模型如图1所示。共铺设121个测点，测点坐标分布矩阵为 

由于地面核磁共振设备接收灵敏度为10nV，每个测点p_i(x,y)处只能接收X-Y-Z空间x-100m至x+100m、y-100m至y+100m、垂深90m以内区域的含水量信息，区域以外的核函数微元为0。每个测点发送20个激发脉冲矩，其最大值为10As。地磁场强度为47000nT，磁偏角为40°，收/发线圈边长为100m，背景电阻率为500Ω·m。在不同信噪比下，反演结果对比图如图2-图5所示，表1中列出了其对应的方均根值。 

表1不同信噪比下，反演结果对比表 

fSNR/dB	5	0	-5	-10
					fRMS/％	11.09	11.12	11.13	11.35

从图2-图5和表1的反演结果中可以看出，本发明的反演算法在不同信噪比下反演精度均较高，算法稳定性好，抗噪能力强，随着信噪比的降低，反演结果的方均根值略有增加，在f_SNR＝-5dB时，其反演精度仍然很高，当信噪f_SNR比降到-10dB时，反演结果中的含水体与背景区域的边界变得比较模糊，已经不能较好的定位出含水体的空间分布。因此，本发明适用于信噪比大于-5dB的SNMR信号三维反演。 

Claims (4)

1.一种地面核磁共振三维反演方法，其特征是包括如下步骤： 

步骤1，用拉直变换方法将三维正演模型进行降维处理，将其抽象为矩阵方程求解模型，即 

E^*＝K^*·n   ① 

式中，E^*为地面核磁共振信号初始振幅值矩阵，K^*为核函数矩阵，n为待求解的含水量空间分布值； 

步骤2，通过补0使得矩阵方程求解模型中的核函数矩阵K^*变为核函数方阵K，相应的地面核磁共振信号初始振幅值矩阵E^*变为初始振幅值矩阵E； 

步骤3，对核函数方阵K做奇异值分解，以获得核函数方阵K的奇异值σ_i； 

步骤4，根据所得奇异值σ_i计算方阵K的有效秩r^*，并令m＝r^*； 

步骤5，根据下式计算初始标准正交基V₁， 

   ② 

式中，V₁为初始标准正交基，R₀是待求变量矩阵中第一个元素值，K为核函数方阵，E为初始振幅值矩阵，n₀为设定的初始化含水量值，||·||_F为矩阵的弗罗伯尼范数算子； 

步骤6，根据所得初始标准正交基V₁和核函数方阵K的转置广义克雷洛夫子空间F_m(K^T,V₁)，并对所构建的F_m(K^T,V₁)利用总体阿诺尔迪方法，在弗罗伯尼范数下，生成一组标准正交基V₁,V₂,…,V_m，并令标准正交基集合U_m＝[V₁,V₂,…,V_m]； 

步骤7，求解向量y，其中向量y满足下式 

   ③ 

式中，y为待求解的向量，V_i是标准正交基集合U_m里面的元素，V_i ^T为V_i的转置，i＝1,2,3……,m，K为核函数方阵，K^T为K的转置，R₀是待求变量矩阵中第一个元素值； 

步骤8，求解含水量变量矩阵n_i，其中含水量变量矩阵n_i满足下式 

   ④ 

式中，n_i为待求解的含水量变量矩阵，i＝1,2,3……,m，abs(·)为复数取模算子，n₀为设定的初始化含水量值，K^T为核函数方阵K的转置， U_m为标准正交基集合，y为步骤7所得的向量；y_i为向量y中的元素，V_i是标准正交基集合U_m里面的元素； 

步骤9，在矩阵n_i中，将值大于1的元素置为1，值小于0的元素置为0； 

步骤10，求解矩阵R_i的弗罗伯尼范数算子||R_i||_F，其中矩阵R_i的弗罗伯尼范数算子||R_i||_F满足下式 

||R_i||_F＝||E-K·n_i||   ⑤ 

式中，||R_i||_F为待求解的弗罗伯尼范数算子，K为核函数方阵，E为初始振幅值矩阵，n_i为步骤8所得的矩阵，i＝1,2,3……,m； 

步骤11：如果实际迭代次数N_s达到设定的最大迭代次数N_max时即N_s≥N_max时，或者所得弗罗伯尼范数算子||R_m||_F小于设定的误差阈值θ时即||R_m||_F<θ时，则停止迭代；否则令将实际迭代次数N_s加1即N_s＝N_s+1，并将含水量值n₀设为本次迭代所得含水量n_s即n₀＝n_s，返回步骤5；其中s＝1,2,3……,r，r为方阵K奇异值的个数。 

2.根据权利要求1所述的一种地面核磁共振三维反演方法，其特征是步骤3中，采用下式对核函数方阵K做奇异值分解，即 

式中σ_i为方阵K的奇异值σ₁≥σ₂≥…≥σ_r，r为方阵K的奇异值的个数；O₃、O₄和O₅是0矩阵；U_l×l和是方阵K的l阶酉矩阵且是正交矩阵。 

3.根据权利要求1所述的一种地面核磁共振三维反演方法，其特征是步骤4中，采用下式计算方阵K的有效秩r^*即令 

    ⑥ 

则r^*为第一个满足ψ(s)>η的s值； 

式中，ψ(s)为求取有效秩r^*的表达式；σ_i为方阵K的奇异值；r为方 阵K奇异值的个数；s为方阵K的第s个奇异值；η为设定值。 

4.据权利要求1所述的一种地面核磁共振三维反演方法，其特征是在步骤11中，最大迭代次数N_max设定为方阵K奇异值的个数r，即N_max＝r。