CN111190233A - 一种基于展宽指数c的预极化场磁共振正反演方法 - Google Patents

Abstract

本发明涉及一种地球物理探测方法及其数据的反演解释，能够提高自由衰减信号单指数拟合的精确度，并且避免了多指数拟合数据量过大的问题，通过令地磁场方向为x轴，水平面上与地磁场方向垂直为y轴，垂直地面向下为z轴，由初始振幅数据计算地下空间位置的灵敏度核函数K，利用核磁共振全波信号进行正反演，通过吉洪诺夫法搜索正则化参数，对目标函数进行高斯牛顿迭代法求解，利用共轭梯度方法来求取每次迭代的模型增量，利用线性搜索来获取最优搜索步长。

Description

一种基于展宽指数C的预极化场磁共振正反演方法

技术领域：

本发明涉及一种地球物理探测方法及其数据的反演解释，能够提高自由衰减信号单指数拟合的精确度，并且避免了多指数拟合数据量过大的问题，利用磁共振技术(Magnetic Resonance Sounding，MRS)可以高效精确的进行层状地下水结构的探测。

背景技术：

地面核磁共振(Surface Nuclear Magnetic Resonance，简称SNMR)是国际上发展起来的一种新的地球物理直接探测地下水的方法，这种在地面直接探测地下介质中氢核丰度的技术，不仅可以用于缺水地区的地下水资源勘查与评价，还可以在地下水引起的堤坝渗漏、矿井突水、隧道涌水、滑坡等地质灾害水源的探测预警中发挥独特的作用。

Mike Mueller-Petke和Ugur Yaramanci在论文QT inversion-Comprehensiveuse of the complete surface NMR data set(GEOPHYSICS卷：75期：4页WA199-WA209)中提出了一种新的反演方式，即QT反演，是利用全部核磁共振包络信号一次性带入反演算法当中，提取含水构造的局部含水量分布信息和平均弛豫时间T₂*分布信息，再进行综合，得到地下不同位置、不同弛豫时间下的总含水量赋存情况，实现含水量w和弛豫时间T₂*两个参数的高精度二维成像，提高了地下弛豫时间分布的空间分辨率和反演问题的稳定性。O.Mohnke和U.Yaramanci在论文Pore size distributions and hydraulicconductivities of rocks derived from Magnetic Resonance Sounding relaxationdata using multi-exponential decay time inversion(JOURNAL OF APPLIEDGEOPHYSICS卷：66期：3-4页：73-81特刊：SI)提出了用多指数衰减时间反演磁共振测深弛豫资料反演岩石的孔径分布和磁共振孔隙度以及衰减时间的方法，比较和讨论了目前水文地质学和核磁共振实验中用于估算磁共振孔隙度以及衰减时间的关系式及其在磁共振波谱中的适用性，并考虑了各种实验室核磁共振和磁共振波谱实验的结果，即岩石比定标系数。

上述论文QT inversion-Comprehensive use of the complete surface NMRdata set中QT(GEOPHYSICS卷：75期：4页WA199-WA209)反演解决了一维地面核磁共振反演问题，并且是目前稳定性最好的一维核磁共振反演方法，但是其为全波反演，在单指数反演时其数据量较大，若采用多指数反演，计算难度会大大提高，在三维二维反演中很难实现。上述论文Pore size distributions and hydraulic conductivities of rocks derivedfrom Magnetic Resonance Sounding relaxation data using multi-exponentialdecay time inversion(JOURNAL OF APPLIED GEOPHYSICS卷：66期：3-4页：73-81特刊：SI)采用多指数拟合的方式，反演数据量过大，计算速度慢，在实际应用中并不实用。

发明内容

本发明的目的在于针对上述现有技术的不足，提供了一种适用于层状水体的高效率高分辨率探测，既能节省反演时间，又能得到高精度平滑反演结果，针对复杂孔隙情况产生多个T₂*的FID信号时能够得到更加精确地孔隙度解释结果，是对层状水的一种基于展宽指数C的预极化场磁共振正反演方法。

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

一种基于展宽指数C的预极化场磁共振正反演方法，该方法包括：

1)根据预探测地实际情况设置材料电阻率，建立三维可视化模型；

2)先向极化线圈通入极化电流，产生预极化场，计算预极化场的氢质子净磁化强度；

3)根据步骤2)中的氢质子净磁化强度计算灵敏度核函数K(q,r)，并根据灵敏度核函数K(q,r)计算z方向的一维核函数；

4)利用步骤3)中z方向的一维核函数计算核磁共振包络曲线，将核磁共振包络曲线表达式离散为矩阵形式，并采用共轭梯度求解器缩小数据矩阵；

5)利用步骤4)的离散为矩阵形式的核磁共振包络曲线构建基于展宽指数C的总体目标函数；

6)利用吉洪诺夫法搜索最优正则化参数；

7)将步骤5)中的总体目标函数重新表示成迭代格式，采用步骤6)中的最优正则化参数，对总体目标函数进行高斯牛顿迭代求解，利用共轭梯度法来求取每次迭代的模型增量，利用线性搜索获得搜索步长；

8)对搜索步长进行误差判断，若大于设定误差则返回步骤7)；若是，则进行下一步；

9)输出最佳搜索步长，根据最佳搜索步长得到反演结果，快速成像。

进一步地，步骤2)中，氢质子净磁化强度为：

其中B_p0为氢质子的静态磁场的强度大小，N为单位体积氢原子数量，N＝6.692×10²⁸；B_p0为预极化场的磁场强度；γ为氢原子旋磁比,大小为0.267518×10⁹；

为约化普朗克常数；K_B为玻尔兹曼常数，其大小为1.3805×10^-23；T为开氏温度，通常为293K。

进一步地，步骤2)求取氢质子的静态磁场的强度大小的方法为：

21)利用毕奥萨伐定理计算磁场，其中地磁场为记为B₀，由预极化线圈激发的预极化磁场记B_p，预极化磁场按三分量形式表示为：

式中，B_p为预极化磁场；B_px为沿x方向的矢量磁场；B_py为沿y方向的矢量磁场；B_pz为沿在z方向的矢量磁场；

和

为直角坐标系的方向向量；

22)地磁场同时包含地磁倾角I和地磁偏角D，为了简化磁场计算，对坐标系进行旋转，引入旋转坐标系，旋转前的原始坐标系记为坐标系x、y、z，首先，坐标系以z轴为轴，沿顺时针方向水平旋转角度D，水平旋转后的坐标系记为坐标系x′、y′、z′；其次，以y′轴为轴，沿顺时针方向垂直旋转角度I，垂直旋转后的坐标系为坐标系x"、y"、z"，预极化磁场B_p在新坐标系x"、y"、z"之下为：

B_p"与B_p的关系为

B"_p＝R_IR_DB_p

其中，

I为地磁倾角与，D为地磁偏角；

当预极化线圈为矩形时，预极化磁场B_p随预极化线圈法向偏度α_n及法向倾度β_n的变化而变化，原始坐标系下，沿顺时针方向水平旋转角度α_n，再沿顺时针方向垂直旋转角度β_n，坐标旋转后磁场记为B_p"'，坐标旋转后磁场B_p"'与原坐标系下B_p的关系为：

B_p"'＝R_αR_βB_p

其中，

对坐标系进行法向倾角和法向偏角的逆旋转，有：

则由预极化线圈所在坐标系下的磁场B_p"'得到地磁场所在坐标系下的等效磁场：

23)旋转坐标系之后的预极化磁场B_p"的x轴方向与地磁场方向重合，作用于氢质子的静态磁场的强度大小为：

进一步地，步骤5)具体包括：

三维灵敏度核函数记为K(q,r)，计算过程为：

式中，q为激发脉冲矩，r为某一位置，w(r)为剖分单元的含水量，t代表接收时间，

为弛豫时间，C为展宽指数，ω_L为激发频率，M_bp为氢原子核的宏观磁化强度，λ_P为旋磁比，I₀为单位电流强度，

为发射场顺时针分量，

为接收场逆时针分量，ξ_T和ξ_R为发射和接收相位；

将三维核函数K(q,r)写成K(q；x，y，z)的形式计算一维灵敏度核函数，将一维灵敏度核函数记为K_1D(q；z)，其表达式为：

式中，K_1D(q；z)为关于z方向的一维核函数，q为激发脉冲矩。

进一步地，步骤4)中核磁共振包络曲线表达式如下：

式中，t、

分别代表接收时间、弛豫时间，q为激发脉冲矩，z为深度，w(z)为该剖分身深度的含水量，K_1D(q，z)表示核函数，C可以有效的控制弛豫时间的分布。

进一步地，将核磁共振包络曲线表达式离散为矩阵形式：

其中，Δ代表差值；K_w为核函数对含水量的雅克比矩阵；

为核函数对弛豫时间的雅克比矩阵；K_C为核函数展宽指数C的雅克比矩阵，表示为：

式中，I是与C维度相同的单位矩阵。

进一步地，步骤5)构建基于展宽指数C的反演目标函数包括：由数据目标函数Φ_d和模型目标函数Φ_m建立关于实测信号与正演信号差值的最优化问题的反演总体目标函数Φ，则吉洪诺夫正则化标准方程：

Φ＝Φ_d+λΦ_m

式中，λ为正则化参数，数据目标函数为：

式中，D实测信号振幅V的权值，由信号或噪声的不确定度计算得到，模型目标函数为：

式中，C_s是平滑度矩阵，

为平滑限制条件，确保在所有剖分网格内部都能得到稳定的平滑约束解，其中矩阵m为包含含水量和弛豫时间的矩阵，通过下式计算：

进一步地，

总体目标函数重新表示成迭代格式：w_k+1＝w_k+η_kΔw_k

其中，k是当前迭代次数，w_k为当前迭代次数下的含水量值，η_k是搜索步长，Δw_k为含水量的模型增量；

步骤7)中，利用高斯牛顿迭代法建立模型建立包含模型增量Δw_k的方程

式中，T表示矩阵的转置；J为Gker的雅克比矩阵，用共轭梯度最小二乘法法来求取每次迭代的模型增量Δw_k，利用线性搜索来获取最佳搜索步长，进而求得总体目标函数Φ的最小值对应的含水量。

有益效果：本发明公开了一种基于展宽指数C的预极化场磁共振正反演方法，在实际环境中复杂含水构造可能存在同一位置含水体具有多种岩性，或者相邻的两个含水体岩性不同的情况，此时，复杂的含水构造存在多个弛豫时间T₂*。对于复杂含水构造，核磁共振信号以多指数衰减形式存在，单指数拟合并不能准确获得核磁共振信号，但是多指数拟合数据集太大，计算太慢，因此本发明利用展宽指数C的衰减去近似代替多指数拟合，在保持计算速度的同时又提高了反演精度，针对复杂孔隙情况产生多个T₂*的自由衰减信号时能够得到更加精确地孔隙度解释结果，对二维、三维的反演都具有巨大意义。利用吉洪诺夫正则化的方法自动选取最优正则化参数，对目标函数进行高斯牛顿迭代法求解，利用共轭梯度方法来求取每次迭代的模型增量，利用线性搜索来获取最优搜索步长，通过合理的正则化参数和约束条件选择，得到高精度平滑的含水量反演结果，在获得高精度平滑反演结果的同时实现了快速计算。

附图说明：

图1为一维地面核磁共振探测线圈布设位置关系图；

图2为坐标系旋转示意图，(a)原始坐标系；(b)水平旋转后坐标系；(c)垂直旋转后坐标系；

图3为线圈坐标系旋转示意图，(a)原始坐标系；(b)水平旋转后坐标系；(c)垂直旋转后坐标系；

图4基于展宽指数C的核磁共振预极化正反演方法流程图；

具体实施方式：

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

本发明基于展宽指数C的一维地面核磁共振正反演方法是一种适用于层状水体的高效率高分辨率探测，在保持计算时间的同时又能可靠获得高精度平滑的反演结果的反演方法，该方法通过发射线圈向地下发射频率为当地拉莫尔频率的交变电流，产生感应磁场。在感应磁场的作用下，地下水中氢质子产生能级跃迁，大量氢质子跃迁到高能级上。当撤去供电电流，这些高能级氢质子便逐渐回到低能级状态，释放出大量的具有拉莫尔频率的能量子，在接收线圈中感应出核磁共振信号。令地磁场方向为x轴，水平面上与地磁场方向垂直为y轴，垂直地面向下为z轴，由初始振幅数据计算地下空间位置的灵敏度核函数K，利用核磁共振全波信号进行正反演，通过吉洪诺夫法搜索正则化参数，对目标函数进行高斯牛顿迭代法求解，利用共轭梯度方法来求取每次迭代的模型增量，利用线性搜索来获取最优搜索步长，此外，对于地下介质含水量w，其取值应该具有实际意义，所以采用一个正切变换使其取值限制在[0，1]范围内。通过合理的正则化参数和约束条件选择，得到高精度平滑的一维含水量反演结果。

如图1所示，该方法先给预极化线圈通入直流电，产生预极化场，然后通过发射线圈向地下发射频率为当地拉莫尔频率的交变电流，产生感应磁场。在感应磁场的作用下，地下水中氢质子产生能级跃迁，大量氢质子跃迁到高能级上。当撤去供电电流，这些高能级氢质子便逐渐回到低能级状态，释放出大量的具有拉莫尔频率的能量子，在接收线圈中感应出核磁共振信号。其中，x(east)正东方向，y(south)正南方向，z(down)垂直向下方向。

参见图3为一种基于展宽指数C的预极化场磁共振正反演方法，包括

1)根据预探测地实际情况设置材料电阻率，建立三维可视化模型；

2)先向极化线圈通入极化电流，产生预极化场，计算预极化场的氢质子净磁化强度；

3)根据步骤2)中的氢质子净磁化强度计算灵敏度核函数K(q,r)，并根据灵敏度核函数K(q,r)计算z方向的一维核函数；

4)利用步骤3)中z方向的一维核函数计算核磁共振包络曲线，将核磁共振包络曲线表达式离散为矩阵形式，并采用共轭梯度求解器缩小数据矩阵；

5)构建基于展宽指数C的总体目标函数；

6)利用吉洪诺夫法搜索最优正则化参数；

7)将步骤5)中的总体目标函数重新表示成迭代格式，采用步骤6)中的最优正则化参数，对总体目标函数进行高斯牛顿迭代求解，获得搜索步长；

8)对搜索步长进行误差判断，若大于设定误差则返回步骤7)；若是，则进行下一步；

9)输出最佳搜索步长，根据最佳搜索步长得到反演结果，快速成像。

步骤1和2具体包括：如图2和图3所示，当磁场方向同时考虑地磁倾角I与地磁偏角D时，预极化磁场B_p在与B₀方向上的关系非常复杂。针对这种情况，引入旋转坐标系的方法进行计算，从而将磁场合成。旋转前的原始坐标系记为坐标系x、y、z。首先，原始坐标系以z轴为轴，沿顺时针方向水平旋转角度D，水平旋转后的坐标系记为坐标系x′、y′、z′(图2b)；其次，以y′轴为轴，沿顺时针方向垂直旋转角度I，垂直旋转后的坐标系见图2(c)记为坐标系x"、y"、z"。预极化磁场B_p在新坐标系x"、y"、z"之下可表示为：

此时，B_p"与B_p的关系可表示为

B'_p'＝R_IR_DB_p

其中，

参见图3(a)、(b)和(c)所示，当预极化线圈为矩形时，预极化磁场B_p随预极化线圈法向偏度α_n及法向倾度β_n的变化而变化，同理，可通过旋转坐标系的方法对其进行计算。在原始坐标系下，沿顺时针方向水平旋转角度α_n，再沿顺时针方向垂直旋转角度β_n，坐标旋转后磁场记为B_p"'，原坐标系下B_p的关系可表示为：B_p"'＝R_αR_βB_p

其中，

由于实际计算过程中通常在线圈所在的坐标系之下进行计算，因此要对进行了法向偏度α_n及法向倾度β_n旋转的坐标系进行逆旋转，计算预极化磁场B_p，有：

式中R_α和R_β都是正定的，因此取逆计算可用转置T代替。则由线圈所在坐标系下的磁场B_p"'可得到地磁场所在坐标系下的等效磁场：

旋转坐标系之后的预极化磁场B_p"的x轴方向与地磁场方向重合，此时作用于氢质子的静态磁场的强度大小可表示为：

此时氢质子净磁化强度为

步骤3)三维灵敏度核函数记为K(q,r)，其计算过程为：

式中，q为激发脉冲矩，r为某一位置，w(r)为该剖分单元的含水量，t、

分别代表接收时间、弛豫时间，C为展宽指数，ω_L为激发频率，M_bp为氢原子核的宏观磁化强度，λ_P为旋磁比，I₀为单位电流强度，

为发射场顺时针分量，

为接收场逆时针分量，ξ_T和ξ_R为发射和接收相位。

将三维核函数K(q，r)写成K(q；x，y，z)的形式计算一维灵敏度核函数，将一维灵敏度核函数记为K_1D(q；z)，其表达式为：

式中，K_1D为关于z方向的一维核函数，q为激发脉冲矩；

步骤4)在基于核磁共振信号全部包络数据的QT反演(Q-Time Inversion，QTI)方法的基础上，用展宽指数C近似代替核磁自由衰减信号的多指数现象。传统所采用的的单指数拟合方式对于多孔隙情况误差大，所以引入了展宽指数C，它有助于用单个弛豫时间参数来描述孔隙，并且能够更加精确的确定初始振幅和弛豫时间。此外，与核磁共振数据的多指数拟合相比，核磁共振自由衰减信号的数据集可以大量减少，节省了计算时间，提高了反演精度。展宽指数C虽然没有严格的物理基础，但它表征了用单指数公式描述信号衰减的偏差，它是弛豫速率均匀性的一种表示；在核磁共振自由衰减信号的公式中加入参数C，通过参数C的衰减可以更精确的拟合振幅的多指数衰减现象，核磁共振包络曲线表达式如下：

式中，t、

分别代表接收时间、弛豫时间，q为激发脉冲矩，z为深度，w(z)为该剖分身深度的含水量，K_1D(q，z)表示核函数，C可以有效的控制弛豫时间的分布；

将上面的核磁共振包络曲线表达式离散为矩阵形式为：

其中，Δ代表差值；K_w为核函数对含水量的雅克比矩阵；

为核函数对弛豫时间的雅克比矩阵；K_C为核函数展宽指数C的雅克比矩阵，表示为：

式中，I是与C维度相同的单位矩阵。

本发明所用到的数据是整个核磁包络曲线的振幅，数据集会产生大量的数据，不利于后续的指数计算，所以在反演之前利用共轭梯度求解器缩小数据矩阵减少计算数据，加快计算速度；

采用基于展宽指数C的预极化场磁共振反演方法进行反演，即在反演过程中除了对含水量和弛豫时间T₂ ^*的反演外还增加了对展宽指数C的反演，提高了拟合自由衰减信号的精确度；

步骤5)由数据目标函数Φ_d和模型目标函数Φ_m建立关于实测信号与正演信号差值的最优化问题的反演总体目标函数Φ，写作吉洪诺夫正则化标准方程：

Φ＝Φ_d+λΦ_m

式中，λ为正则化参数，数据目标函数表示为：

式中，D为实测信号振幅V的权值，由信号或噪声的不确定度计算得到，模型目标函数表示为：

式中，C_s是平滑度矩阵，

为平滑限制条件，确保在所有剖分网格内部都能得到稳定的平滑约束解，其中矩阵m为包含含水量和弛豫时间的矩阵，可表示为下式：

w是含水量。

步骤6)利用吉洪诺夫正则化法选取最优正则化参数λ；

步骤7)将总体目标函数表达式可以重新表示成迭代格式：

w_k+1＝w_k+η_kΔw_k

其中，k是当前迭代次数，w_k为当前迭代次数下的含水量值，ηk是搜索步长，Δw_k为含水量的模型增量。新的模型增量可用高斯牛顿方法求解。

利用高斯牛顿迭代法建立模型方程

式中，T表示矩阵的转置；J为Gker的雅克比矩阵，用共轭梯度最小二乘法法来求取每次迭代的模型增量Δw_k，利用线性搜索来获取最佳搜索步长，进而求得总体目标函数Φ的最小值对应的含水量；

n、重复上述m过程，直到搜索步长小于设定误差，获得高精度的反演结果，输出并快速成像。

Claims (8)

1.一种基于展宽指数C的预极化场磁共振正反演方法，其特征在于，该方法包括：

1)根据预探测地实际情况设置材料电阻率，建立三维可视化模型；

2)先向极化线圈通入极化电流，产生预极化场，计算预极化场的氢质子净磁化强度；

3)根据步骤2)中的氢质子净磁化强度计算灵敏度核函数K(q,r)，并根据灵敏度核函数K(q,r)计算z方向的一维核函数；

4)利用步骤3)中z方向的一维核函数计算核磁共振包络曲线，将核磁共振包络曲线表达式离散为矩阵形式，并采用共轭梯度求解器缩小数据矩阵；

5)利用步骤4)的离散为矩阵形式的核磁共振包络曲线构建基于展宽指数C的总体目标函数；

6)利用吉洪诺夫法搜索最优正则化参数；

7)将步骤5)中的总体目标函数重新表示成迭代格式，采用步骤6)中的最优正则化参数，对总体目标函数进行高斯牛顿迭代求解，利用共轭梯度法来求取每次迭代的模型增量，利用线性搜索获得搜索步长；

8)对搜索步长进行误差判断，若大于设定误差则返回步骤7)；若是，则进行下一步；

9)输出最佳搜索步长，根据最佳搜索步长得到反演结果，快速成像。

2.按照权利要求1所述的方法，其特征在于，步骤2)中，氢质子净磁化强度为：

其中B_p0为氢质子的静态磁场的强度大小，N为单位体积氢原子数量，N＝6.692×10²⁸；B_p0为预极化场的磁场强度；γ为氢原子旋磁比,大小为0.267518×10⁹；

为约化普朗克常数；K_B为玻尔兹曼常数，其大小为1.3805×10^-23；T为开氏温度，通常为293K。

3.按照权利要求2所述的方法，其特征在于，步骤2)求取氢质子的静态磁场的强度大小的方法为：

21)利用毕奥萨伐定理计算磁场，其中地磁场为记为B₀，由预极化线圈激发的预极化磁场记B_p，预极化磁场按三分量形式表示为：

式中，B_p为预极化磁场；B_px为沿x方向的矢量磁场；B_py为沿y方向的矢量磁场；B_pz为沿在z方向的矢量磁场；

和

为直角坐标系的方向向量；

22)地磁场同时包含地磁倾角I和地磁偏角D，为了简化磁场计算，对坐标系进行旋转，引入旋转坐标系，旋转前的原始坐标系记为坐标系x、y、z，首先，坐标系以z轴为轴，沿顺时针方向水平旋转角度D，水平旋转后的坐标系记为坐标系x′、y′、z′；其次，以y′轴为轴，沿顺时针方向垂直旋转角度I，垂直旋转后的坐标系为坐标系x"、y"、z"，预极化磁场B_p在新坐标系x"、y"、z"之下为：

B_p"与B_p的关系为

B"_p＝R_IR_DB_p

其中，

I为地磁倾角与，D为地磁偏角；

当预极化线圈为矩形时，预极化磁场B_p随预极化线圈法向偏度α_n及法向倾度β_n的变化而变化，原始坐标系下，沿顺时针方向水平旋转角度α_n，再沿顺时针方向垂直旋转角度β_n，坐标旋转后磁场记为B_p"'，坐标旋转后磁场B_p"'与原坐标系下B_p的关系为：

B_p"'＝R_αR_βB_p

其中，

对坐标系进行法向倾角和法向偏角的逆旋转，有：

则由预极化线圈所在坐标系下的磁场B_p"'得到地磁场所在坐标系下的等效磁场：

23)旋转坐标系之后的预极化磁场B_p"的x轴方向与地磁场方向重合，作用于氢质子的静态磁场的强度大小为：

4.按照权利要求1所述的方法，其特征在于，步骤5)具体包括：

三维灵敏度核函数记为K(q,r)，计算过程为：

式中，q为激发脉冲矩，r为某一位置，w(r)为剖分单元的含水量，t代表接收时间，

为弛豫时间，C为展宽指数，ω_L为激发频率，M_bp为氢原子核的宏观磁化强度，λ_P为旋磁比，I₀为单位电流强度，

为发射场顺时针分量，

为接收场逆时针分量，ξ_T和ξ_R为发射和接收相位；

将三维核函数K(q,r)写成K(q；x，y，z)的形式计算一维灵敏度核函数，将一维灵敏度核函数记为K_1D(q；z)，其表达式为：

式中，K_1D(q；z)为关于z方向的一维核函数，q为激发脉冲矩。

5.按照权利要求1所述的方法，其特征在于，步骤4)中核磁共振包络曲线表达式如下：

式中，t、

分别代表接收时间、弛豫时间，q为激发脉冲矩，z为深度，w(z)为该剖分身深度的含水量，K_1D(q，z)表示核函数，C可以有效的控制弛豫时间的分布。

6.按照权利要求5所述的方法，其特征在于，将核磁共振包络曲线表达式离散为矩阵形式：

其中，Δ代表差值；K_w为核函数对含水量的雅克比矩阵；

为核函数对弛豫时间的雅克比矩阵；K_C为核函数展宽指数C的雅克比矩阵，表示为：

式中，I是与C维度相同的单位矩阵。

7.按照权利要求1所述的方法，其特征在于，步骤5)构建基于展宽指数C的反演目标函数包括：由数据目标函数Φ_d和模型目标函数Φ_m建立关于实测信号与正演信号差值的最优化问题的反演总体目标函数Φ，则吉洪诺夫正则化标准方程：

Φ＝Φ_d+λΦ_m

式中，λ为正则化参数，数据目标函数为：

式中，D实测信号振幅V的权值，由信号或噪声的不确定度计算得到，模型目标函数为：

式中，C_s是平滑度矩阵，

为平滑限制条件，确保在所有剖分网格内部都能得到稳定的平滑约束解，其中矩阵m为包含含水量和弛豫时间的矩阵，通过下式计算：

8.按照权利要求1所述的方法，其特征在于，

总体目标函数重新表示成迭代格式：w_k+1＝w_k+η_kΔw_k

其中，k是当前迭代次数，w_k为当前迭代次数下的含水量值，η_k是搜索步长，Δw_k为含水量的模型增量；

步骤7)中，利用高斯牛顿迭代法建立模型建立包含模型增量Δw_k的方程

式中，T表示矩阵的转置；J为Gker的雅克比矩阵，用共轭梯度最小二乘法法来求取每次迭代的模型增量Δw_k，利用线性搜索来获取最佳搜索步长，进而求得总体目标函数Φ的最小值对应的含水量。

Images