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

Abstract

本发明公开了一种地面核磁共振二维反演方法，其用拉直变换方法将二维正演模型进行降维处理，将其抽象为矩阵方程求解模型，并用最小二乘奇异值分解(LS-SVD)与改进的随机梯度下降法(ISGD)相结合的方法进行反演求解，采用LS-SVD求取矩阵方程的粗略解，在该粗略解的基础上，用ISGD求取其精细解。在不同信噪比的条件下，本发明的反演结果均与模型中含水构造分布相吻合，即使在信噪比为0dB时，其反演结果仍能分辨出地下水文地质构造，其反演得到的含水量值的方均根为8.26％，而此时LS-SVD和ISGD两种方法的反演结果均无效，其方均根值分别为30.14％和15.35％。

Description

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

技术领域

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

背景技术

地面核磁共振(Surface Nuclear Magnetic Resonance，简称SNMR)技术是目前世界上唯一的一种直接找水的物探方法，该项技术已在探测地下水、考古、地下水污染检测等领域得到了一定的应用。近年来，随着专家和学者们的逐渐深入研究，SNMR技术得到了进一步的完善。反演计算含水率是该技术研究过程中的关键环节，而反演准确度和分辨率是衡量反演算法性能的关键指标。其中，一维正反演理论较为成熟，已经相继刊登出多种有效算法，如：文献1[DAI Miao,HU Xiangyun,WU Haibo,et al.“Inversion of surfacenuclear magnetic resonance for groundwater exploration,”ChineseJournal of Geophysics,2009,52(5):1166-1173.]提出了改进的模拟退火算法反演，提高了现有反演算法的稳定度和收敛速度；文献2[Mueller-PetkeM.,Yaramanci U..QT inversion-comprehensive use of the completesurface NMR data set[J].Geophysics,2010,75:199–209.]提出了QT反演算法，利用各个激发脉冲矩对应的全部采样点数据进行反演，充分挖掘了接收信号信息，在一定程度上提高了反演精度；文献3[Ahmad A.Behroozmand,Esben Auken,Gianluca Fiandaca,et al.Efficient fulldecay inversion of MRS data with a stretched-exponential approximationof the T2*distribution[J].Geophysical Journal International,2012,190:900–912.]采用了积分门技术接收信号，提高各个采样点数据的精度，并进行全衰减反演，是对QT反演的一种改进。在二维反演方面，Boucher、Girard和Legchenko等研究了在二维剖面方向上E0-q曲线随地下含水构造的变化趋势，但他们只对二维反演做了定性研究，没有给出具体的二维反演公式。Legchenko等对三维反演做了一定的研究，虽然能在三维空间反演出模型的含水构造，但是由于在三维空间设定的网格尺寸较大，只能粗略的估计出地下含水构造，其反演分辨率有待提高。由于二维、三维反演算法存在运算量大、待求解变量数多、非线性等问题，目前世界上唯一商业版反演软件NUMISPLUS仍采用一维反演，而二维、三维正反演研究仍处于起步阶段。

发明内容

本发明所要解决的技术问题是现有地面核磁共振技术的实用性不强，一维反演算法横向分辨率低的不足，提出一种地面核磁共振二维反演方法。

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

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

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

步骤2、采用最小二乘奇异值分解法求取步骤1所抽象出的矩阵方程求解模型的粗略解n_LS-SVD；

步骤2.1、对步骤1所抽象出的矩阵方程求解模型中的核函数矩阵K做奇异值分解，以获得核函数矩阵K的奇异值σ；

步骤2.2、根据步骤2.1所获得的核函数矩阵K的奇异值σ计算核函数矩阵K的有效秩r^*；

步骤2.3、根据步骤2.1所获得的核函数矩阵K的奇异值σ和步骤2.3所计算出的核函数矩阵K的有效秩r^*，求取步骤1所抽象出的矩阵方程求解模型的最小二乘解；

步骤3、采用随机梯度下降法求取步骤1所抽象出的矩阵方程求解模型的精细解n^#；

步骤3.1、采用吉洪诺正则化方法构建模型的适应度函数，并将该适应度函数的当前一次迭代的解初始化为n_LS-SVD，初始迭代次数初始化为0；

步骤3.2、根据搜索路径更新公式计算模型的适应度函数的下一次迭代的搜索路径，并将更新的搜索路径视为当前一次迭代的搜索路径；

步骤3.3、将当前一次迭代的搜索路径中的每一个个体n_h+1作为当前值n_cur带入步骤3.1所构建的模型的适应度函数中，计算适应度函数的当前一次迭代的最优解n^#；

步骤3.4、如果迭代次数达到设定的最大迭代次数N_max或者当前最优适应度函数值小于设定的反演精度阈值φ，则停止迭代，反演结果为当前最优解n^#；否则，返回步骤3.2。

上述步骤1所抽象出的矩阵方程求解模型为：

  ①

式中，E为初始振幅，K为核函数矩阵，n为所求向量，M为剖面方向总共发射的激发脉冲矩的个数，L为核函数矩阵K的列向量的个数。

上述步骤2.1，采用式②对矩阵方程求解模型中的核函数矩阵K做奇异值分解

$K_{M\×L}=U_{M\×M}{\mathrm{\ΛV}}_{L\×L}^H$   ②

式中，U_M×M和是正交矩阵，为V_L×L的复共轭转置，M为剖面方向总共发射的激发脉冲矩的个数，L为核函数矩阵K的列向量的个数，Λ为对角阵，

  ③

式中，Λ为对角阵，M为剖面方向总共发射的激发脉冲矩的个数，L为矩阵K的列向量的个数，r为核函数矩阵K奇异值的个数，σ为核函数矩阵K的奇异值。

上述步骤2.2具体为：

令

$\mathrm{\ψ}\left(s\right)=\frac{\sqrt{{\mathrm{\σ}}_1+{\mathrm{\σ}}_2+...+{\mathrm{\σ}}_s}}{\sqrt{{\mathrm{\σ}}_1+{\mathrm{\σ}}_2+...+{\mathrm{\σ}}_r}},s=\mathrm{1,2},...,r$   ④

则有效秩r^*为第一个满足ψ(s)>θ的s值；

式中，ψ(s)为求取有效秩r^*的表达式，σ为核函数矩阵K的奇异值，r为核函数矩阵K奇异值的个数，s为核函数矩阵K的第s个奇异值，θ为设定值。

上述步骤2.3具体为：

$n_{\mathrm{LS}-\mathrm{SVD}}=\underset{f=1}{\overset{r^{*}}{\mathrm{\Σ}}}\frac{u_f^{T}E}{{\mathrm{\σ}}_f}{v}_f$   ⑤

式中，n_LS-SVD为用最小二乘奇异值分解求取的解，E为初始振幅，σ为核函数矩阵K的奇异值，u_f和v_f分别为U_M×M和V_L×L的第f个列向量(f＝1,2,…r^*)，u_f ^T为u_f的转置，M为剖面方向总共发射的激发脉冲矩的个数，r^*为核函数矩阵K的有效秩。

上述步骤3.1构建出的模型的适应度函数为：

  ⑥

式中，Η(n_cur)为模型的适应度函数，K为核函数矩阵，n_cur为当前值，E为初始振幅，R为正则化因子，函数为用于求解n_cur的1次偏导数。

上述步骤3.2中，搜索路径更新公式为：

$\left\{\begin{array}{c}{n}_{h+1}^{*}=n_h+w{\mathrm{\&Delta;v}}_h\\ n_{h+1}=\left\{\begin{array}{c}0,n_{h+1}^{*}<0\\ n_{h+1}^{*},0\&le;n_{h+1}^{*}\&le;1\\ 1,n_{h+1}^{*}>1\end{array}\right.\end{array}\right.$   ⑦

式中，n_h+1为当前一次迭代的适应度函数中n_cur的值，n_h为上一次迭代的适应度函数中n_cur的值，w为移动权重，△v_h为移动方向向量，△v_h为移动方向矩阵△_v中的第h列即△v_h＝(△v_1h,△v_2h,...,△v_Lh)^T，其中

  ⑧

式中，△_v为移动方向矩阵，△_v的子个体△v_ij∈{x|-1,0,1}，且每次迭代△v均随机生成，M为每一次迭代的子个体数，亦为剖面方向总共发射的激发脉冲矩的个数，L为核函数矩阵K的列向量的个数。

与现有技术相比，本发明的有益效果是：

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

2)提出了将成熟的一维反演算法LS-SVD与ISGD相结合的二维反演算法。首先，采用LS-SVD求取矩阵方程的粗略解n_LS-SVD；然后，以n_LS-SVD作为模型的初始化值，以Tikhonov正则化方法构建模型的适应度函数，用ISGD求取矩阵方程的精细解。

3)对现有的随机梯度下降法进行改进，采用了变步长搜索，这既能加快算法的收敛速度，又能保证计算精度。

4)在不同信噪比下性能均优于LS-SVD和ISGD，其反演结果均与模型中含水构造分布比较吻合，即使在信噪比为0dB时，其反演结果仍能分辨出地下水文地质构造，其反演得到的含水量值的方均根为8.26％，而此时LS-SVD和ISGD反演结果无效，其方均根值分别为30.14％和15.35％。

附图说明

图1为二维剖面含水量模型；

图2为f_SNR＝15dB时，S-SVD反演(图2a)、ISGD反演(图2b)、本发明反演(图2c)三种算法反演结果对比图。

图3为f_SNR＝10dB时，LS-SVD反演(图3a)、ISGD反演(图3b)、本发明反演(图3c)三种算法反演结果对比图。

图4为f_SNR＝5dB时，LS-SVD反演(图4a)、ISGD反演(图4b)、本发明反演(图4c)三种算法反演结果对比图。

图5为f_SNR＝0dB时，LS-SVD反演(图5a)、ISGD反演(图5b)、本发明反演(图5c)三种算法反演结果对比图。

具体实施方式

本发明所设计的地面核磁共振二维反演方法，采用最小二乘奇异值分解与改进的随机梯度下降法相结合的方法进行反演求解。

在收/发天线共圈模式地面核磁共振(简称SNMR)二维正反演研究中，含水量沿着深度和剖面方向变化，测点p处接收信号的初始振幅E₀(p,q)公式为：

$\left\{\begin{array}{c}{E}_0(p,q)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{K}_{2D}(p,q;x,z)\&CenterDot;n(p;x,z)\mathrm{dxdz}\\ K_{2D}(p,q;x,z)={\&Integral;}_{-\&infin;}^{\&infin;}{K}_{3D}(p,q;x,y,z)\mathrm{dy}\\ \&PartialD;n(p;x,y,z)/\&PartialD;y=0\end{array}\right.---\left(1\right)$

其中，q为激发脉冲矩；K_2D(p,q；x,z)和K_3D(p,q；x,y,z)分别为p点处二维和三维核函数；n(p；x,z)和n(p；x,y,z)分别为p点处二维和三维含水量分布函数。

将探测剖面区域二维含水量n(x,z)沿着深度和剖面方向做离散化处理为：

$\left\{\begin{array}{c}{z}_{k-1}\&le;z\&le;z_k,k=\mathrm{1,2},...,N_z\\ x_{j-1}\&le;x\&le;x_j,j=\mathrm{1,2},...,N_x\end{array}\right.$

得到L＝N_z×N_x个含水微元。在剖面方向上有N_p个测点，每个测点处发射N_q个激发脉冲矩，则剖面方向上总共发射M＝N_p×N_q个激发脉冲矩。测点p_i处的第q_h个激发脉冲矩对应的初始振幅为：

$\begin{array}{c}{E}_0(p_i,q_h)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{K}_{2D}(p_i,q_h;x,z)\&CenterDot;n(p_i;x,z)\mathrm{dxdz}\\ =\underset{j=1}{\overset{N_x}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_z}{\mathrm{\&Sigma;}}}{K}_{2D}(p_i,q_h;x_j,z_k)\&CenterDot;n_{\mathrm{kj}}{\mathrm{\&Delta;z}}_k{\mathrm{\&Delta;x}}_j\end{array}---\left(2\right)$

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

$\left\{\begin{array}{c}{K}_{2D}\&RightArrow;K_m=[K_{m1},K_{m2},...K_{\mathrm{mL}}\\ n\&RightArrow;n^{*}={[n_1,n_2,...,n_L]}^T\end{array}\right.$

则(2)式转化为：

E_m＝E₀(p_i,q_h)＝K_m·n^*  (3)

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

SNMR二维反演是一个带约束的高维欠定非线性方程求解问题，且核函数具有病态性，为了高精度快速求解方程(4)，本发明提出了基于LS-SVD的随机梯度下降法，首先用LS-SVD求取粗略解，之后用该解作为ISGD的初始解，以在该粗略解周围搜索精细解。

步骤1：采用LS-SVD求取方程(4)的粗略解n_LS-SVD。

对方程(4)中的核函数矩阵K做奇异值分解

$K_{M\&times;L}=U_{M\&times;M}{\mathrm{\&Lambda;V}}_{L\&times;L}^H---\left(5\right)$

其中，U_M×M和是正交矩阵，Λ为对角阵，即：

且σ₁≥σ₂≥…≥σ_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(6\right)$

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

方程(4)的最小二乘解为：

$n_{\mathrm{LS}-\mathrm{SVD}}=\underset{f=1}{\overset{r^{*}}{\mathrm{\&Sigma;}}}\frac{u_f^{T}E}{{\mathrm{\&sigma;}}_f}{v}_f---\left(7\right)$

其中，u_f和v_f分别为U_M×M和V_L×L的第f个列向量(f＝1,2,…r^*)。

步骤2：采用ISGD求取方程(4)的精细解n^#。

2.1)算法初始化。解向量初始化为n₀＝n_LS-SVD；最大迭代次数为N_max；初始迭代次数为N_num＝0；每一代子个体数为M；反演精度阈值为φ，在本发明优选实施例中，φ的取值为10^-8；采用Tikhonov正则化方法构建模型的适应度函数为：

其中，R为正则化因子，在本发明优选实施例中，R的取值为10^-8；函数用于求解n_cur的1次偏导数。

2.2)计算下一代搜索路径。下一代搜索路径更新公式为：

$\left\{\begin{array}{c}{n}_{h+1}^{*}=n_h+w{\mathrm{\&Delta;v}}_h\\ n_{h+1}=\left\{\begin{array}{c}0,n_{h+1}^{*}<0\\ n_{h+1}^{*},0\&le;n_{h+1}^{*}\&le;1\\ 1,n_{h+1}^{*}>1\end{array}\right.\end{array}\right.---\left(9\right)$

其中，w为移动权重，△v_h为移动方向向量。

移动方向矩阵为其中，△v_ij∈{x|-1,0,1}，且每次迭代△v均随机生成。将△v向量化为△v＝(△v₁,△v₂,…,△v_M)，其中，△v_h＝(△v_1h,△v_2h,...,△v_Lh)^T。

为了加快收敛速度，同时保证计算精度，本发明对现有的随机梯度下降法进行改进，采用变步长搜索，即：w＝(1-N_num/N_max)△w，其中，△w为固定步长，在本发明优选实施例中，△w的取值为0.05。

w＝(1-N_num/N_max)△w,

2.3)计算适应度函数，更新最优解。将步骤2.2)中所得当前一代中每一个体n_h+1分别带入公式(8)，计算当前一代最优解。

2.4)判断迭代是否停止。如果迭代次数达到N_max，或者当前最优适应度函数值小于φ，则停止迭代，反演结果为当前最优解n^#；否则，返回步骤2.2)，重复执行步骤2.2)-2.4)。

本发明提出了采用拉直变换方法将二维正演模型进行降维处理，将二维反演模型转化为矩阵方程求解的问题；由于转化后的矩阵方程具有待求解变量数多、非线性、病态等特性，为了提高对其的求解速度和精度，并提出了用最小二乘奇异值分解(the least-squares singular value decomposition,简称LS-SVD)与改进的随机梯度下降法(the improved stochastic gradientdescent，简称ISGD)相结合的方法进行反演求解，采用LS-SVD求取矩阵方程的粗略解，在该粗略解的基础上，用ISGD求取其精细解。最后，通过仿真实验在不同信噪比下对本文研究内容进行了验证。

为了验证本发明所提出的SNMR二维反演算法的有效性，首先建立二维水文地质模型；然后，对其进行正演计算，求出各个测点的E0-q曲线；之后，对E0加入一定信噪比的噪声，得到含噪信号E；最后，分别用LS-SVD、ISGD和本文算法对上述含噪信号进行二维反演，并用反演得到的含水量值的方均根fRMS评价反演结果性能，fRMS公式为：

$f_{\mathrm{RMS}}=\sqrt{\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{(n_i-n_i^{\#})}^2}\&times;100\%---\left(10\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(11\right)$

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

为了保证实验结果的实用性，仿真模型尽可能的模拟实测环境，模型在X-Z剖面的尺寸为600m×90m，剖分的微元大小为20m×5m，剖面内共有30×18个含水微元，剖面270m至380m之间垂深25m至40m的区域含水量为50％，其他区域含水量为5％，其二维模型如图1所示。共铺设21个测点，测点坐标分别为p₁(100m,0)、p₂(120m,0)、……、p₂₁(500m,0)。由于地面核磁共振设备接收灵敏度为10nV，每个测点p_i(x,0)处只能接收X-Z剖面上x-100m至x+100m之间垂深90m以内区域的含水量信息，区域以外的核函数微元K_ml为0。每个测点发送20个激发脉冲矩，其最大值为10As。地磁场强度为44630nT，磁偏角为24°，收/发线圈边长为100m，背景电阻率为100Ω·m。在不同信噪比下，三种算法反演结果对比图如图2-图5所示，表1中列出了其对应的方均根值，其中，ISGD和本发明的迭代次数都是1000次。

表1不同信噪比下，三种算法反演结果对比表

fSNR/dB	fRMS(LS-SVD)/％	fRMS(ISGD)/％	fRMS(the new algorithm)/％
				15	6.13	18.58	5.19
10	9.24	18.71	5.85
				5	16.17	17.28	7.06
0	30.14	15.35	8.26

从图2-图5和表1的反演结果对比中，可以看出LS-SVD反演在高信噪比下，其反演精度较高，但是随着信噪比的降低其精度快速降低，在0dB时，已经不能从反演结果中分辨出目标含水构造；2D SNMR反演是一个高维欠定非线性方程求解问题，由于求解过程中待求解变量数较多，ISGD在不同信噪比下反演精度均较差，虽然其反演结果能分辨出目标含水构造，但是剖面图的边界处和垂深70m以下区域存在虚假含水信息；本发明的反演算法在不同信噪比下反演精度均优于LS-SVD和ISGD，算法稳定性好，抗噪能力强，随着信噪比的降低反演结果的方均根值略有增加，在fSNR＝0dB时，其反演精度仍然很高。

Claims (7)

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

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

步骤2、采用最小二乘奇异值分解法求取步骤1所抽象出的矩阵方程求解模型的粗略解n_LS-SVD；

步骤2.1、对步骤1所抽象出的矩阵方程求解模型中的核函数矩阵K做奇异值分解，以获得核函数矩阵K的奇异值σ；

步骤2.2、根据步骤2.1所获得的核函数矩阵K的奇异值σ计算核函数矩阵K的有效秩r^*；

步骤2.3、根据步骤2.1所获得的核函数矩阵K的奇异值σ和步骤2.3所计算出的核函数矩阵K的有效秩r^*，求取步骤1所抽象出的矩阵方程求解模型的最小二乘解；

步骤3、采用随机梯度下降法求取步骤1所抽象出的矩阵方程求解模型的精细解n^#；

步骤3.1、采用吉洪诺正则化方法构建模型的适应度函数，并将该适应度函数的当前一次迭代的解初始化为n_LS-SVD，初始迭代次数初始化为0；

步骤3.2、根据搜索路径更新公式计算模型的适应度函数的下一次迭代的搜索路径，并将更新的搜索路径视为当前一次迭代的搜索路径；

步骤3.3、将当前一次迭代的搜索路径中的每一个个体n_h+1作为当前值n_cur带入步骤3.1所构建的模型的适应度函数中，计算适应度函数的当前一次迭代的最优解n^#；

步骤3.4、如果迭代次数达到设定的最大迭代次数N_max或者当前最优适应度函数值小于设定的反演精度阈值φ，则停止迭代，反演结果为当前最优解n^#；否则，返回步骤3.2。

2.根据权利要求1所述的一种地面核磁共振二维反演方法，其特征是步骤1所抽象出的矩阵方程求解模型为：

  ①

式中，E为初始振幅，K为核函数矩阵，n为所求向量，M为剖面方向总共发射的激发脉冲矩的个数，L为核函数矩阵K的列向量的个数。

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

$K_{M\&times;L}=U_{M\&times;M}{\mathrm{\&Lambda;V}}_{L\&times;L}^H$   ②

式中，U_M×M和是正交矩阵，为V_L×L的复共轭转置，M为剖面方向总共发射的激发脉冲矩的个数，L为核函数矩阵K的列向量的个数，Λ为对角阵，

  ③

式中，Λ为对角阵，M为剖面方向总共发射的激发脉冲矩的个数，L为矩阵K的列向量的个数，r为核函数矩阵K奇异值的个数，σ为核函数矩阵K的奇异值。

4.根据权利要求1所述的一种地面核磁共振二维反演方法，其特征是步骤2.2具体为：

令

$\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^*的表达式，σ为核函数矩阵K的奇异值，r为核函数矩阵K奇异值的个数，s为核函数矩阵K的第s个奇异值，θ为设定值。

5.根据权利要求1所述的一种地面核磁共振二维反演方法，其特征是步骤2.3具体为：

$n_{\mathrm{LS}-\mathrm{SVD}}=\underset{f=1}{\overset{r^{*}}{\mathrm{\&Sigma;}}}\frac{u_f^{T}E}{{\mathrm{\&sigma;}}_f}{v}_f$   ⑤

式中，n_LS-SVD为用最小二乘奇异值分解求取的解，E为初始振幅，σ为核函数矩阵K的奇异值，u_f和v_f分别为U_M×M和V_L×L的第f个列向量(f＝1,2,…r^*)，u_f ^T为u_f的转置，M为剖面方向总共发射的激发脉冲矩的个数，r^*为核函数矩阵K的有效秩。

6.根据权利要求1所述的一种地面核磁共振二维反演方法，其特征是步骤3.1构建出的模型的适应度函数为：

  ⑥

式中，Η(n_cur)为模型的适应度函数，K为核函数矩阵，n_cur为当前值，E为初始振幅，R为正则化因子，函数为用于求解n_cur的1次偏导数。

7.根据权利要求1所述的一种地面核磁共振二维反演方法，其特征是步骤3.2中，搜索路径更新公式为：

$\left\{\begin{array}{c}{n}_{h+1}^{*}=n_h+w{\mathrm{\&Delta;v}}_h\\ n_{h+1}=\left\{\begin{array}{c}0,n_{h+1}^{*}<0\\ n_{h+1}^{*},0\&le;n_{h+1}^{*}\&le;1\\ 1,n_{h+1}^{*}>1\end{array}\right.\end{array}\right.$   ⑦

式中，n_h+1为当前一次迭代的适应度函数中n_cur的值，n_h为上一次迭代的适应度函数中n_cur的值，w为移动权重，△v_h为移动方向向量，△v_h为移动方向矩阵△v中的第h列即△v_h＝(△v_1h,△v_2h,...,△v_Lh)^T，其中

  ⑧

式中，△_v为移动方向矩阵，△_v的子个体△v_ij∈{x|-1,0,1}，且每次迭代△v均随机生成，M为每一次迭代的子个体数，亦为剖面方向总共发射的激发脉冲矩的个数，L为核函数矩阵K的列向量的个数。