中南大学范文

2010高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C

我们的参赛报名号为（如果赛区设置报名号的话）： 80

所属学校（请填写完整的全名）： 中南大学

参赛队员 (打印并签名) ：1. 欧阳彬 2. 韩帅 3. 王佳 指导教师或指导教师组负责人 (打印并签名)：

日期： 2010 年 月 日

赛区评阅编号（由赛区组委会评阅前进行编号）：

2010高教社杯全国大学生数学建模竞赛

编 号 专 用 页

评阅人 评分 备注

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

核磁共振驰豫信号反演问题

摘 要

针对岩石核磁共振分析中横向弛豫时间T2谱的反演问题，本文分驰豫时间分布预先给定与否两种情况进行探究，分别得到了相应的结果，最后考虑误差影响提出了修正模型算法。

预先给定驰豫时间分布时对T2谱的求解。在T2j分布给定下，欲求T2谱，只需确定相应的fj。考虑到实际中解fj应该是非负值，对磁化强度信号y(t)和t进行非负最小二乘法拟合，建立非负多元线性拟合模型（Ⅰ）。通过MATLAB7.10对不同的布点个数进行试算求解，得到三种布点方式(对数均匀布点、幂指数布点、线性均匀布点)的T2谱，并对结果分析得到：对数均匀布点在该实例中布点最为合理，能区分10种类型弛豫时间对应的孔隙。

驰豫时间分布未知时对T2谱的求解。基于模型Ⅰ，可以得到较为准确的驰豫时间T2j的种类数k，以及fj和T2j的大致范围。考虑到目标函数为非线性函数，并且变量具有边界值，为典型的二次约束规划问题，由于遗传算法求解这种问题不会陷入局部最优解且不会受到目标表达式形式影响的优点，因此采用遗传算法（Ⅱ）对其进行求解，通过MATLAB7.10求得T2谱（表1）。为了验证解的精度，用蒙特卡洛算法在最优解邻域内进行了1000次随机检验，结果表明该最优解具有可靠性。

考察随机误差对反演结果产生的影响，首先求得信噪比SNR约为92，然后将以原始数据求解得到的T2谱与理论结果对比，发现：误差对模型Ⅰ有较大的影响，使得模型对空隙种类的区分能力下降，驰豫时间种类k由10种变为3种；误差对模Ⅱ有较小的影响，T2谱fj值10个中有4个产生了波动，其中3的波动都不超过20%，表明遗传算法有一定的抗噪能力。

为克服误差影响，提出改进规划模型（Ⅲ），在目标模型中加入惩罚函数，考虑到曲率平滑修正方法比较稳定，即使原始回波串质量较差时，也能得出较好的结果，采用曲率平滑的方式减小模型算法对噪声的敏感性。以原始数据求解验证，将所得的T2谱与理论结果进行误差分析（表4），发现二者较为接近（图10），说明该模型对信噪比较低的情况也可适用。

最后，对三个模型进行了评价。

关键词：非负最小二乘 遗传算法 蒙特卡洛 曲率平滑

1

1 问题重述

在核磁共振测井中, 一般采用CMPG方法测量自旋回波串，纵向驰豫时间T1和横向驰豫时间T2都是用来描述自旋回波信号的驰豫特征。由于T2谱能够提供许多岩石物性和流体特性的信息，越来越受到人们的关注。根据核磁共振理论，单个空隙中的磁化强度信号的衰减满足单指数衰减规律，但由于岩石内部是一系列大小不等的孔隙群体组成，所以在岩石核磁共振中测得的中磁化强度信号y(t)为一系列单个孔隙磁化强度信号的叠加，因此y(t)可描述为

mtT2j

y(t)j1fje,tn

（其中fj为第j类孔隙在总孔隙中所占的份额，T2j为第j类孔隙的T2驰豫时间，通常范围为0.1msT2j10000ms，为回波间隔时间，t为采集时间）

岩石核磁共振分析的关键是如何从上式中解出各类孔隙的T2驰豫时间T2j以及各类孔隙在总孔隙中所占的份额fj，称之为T2谱。

分析题意可知，欲得到横向驰豫时间的T2谱，有以下3个问题需要解决： （1） 预先给定T2驰豫时间分布T2j，典型取法为在(T2min,T2max)内对数均匀选取m个点，称之为驰豫时间布点，也可考虑用2的幂指数布点、线性均匀布点等方式。基于上述几种布点方式建立求解T2谱的数学模型。

（2） 假设不预先给定T2驰豫时间分布T2j，建立求解T2谱的数学模型。

（3） 考虑到测量数据的误差，通过算例分析误差给前面两个问题带来的影响；并考虑如何克服误差带来的不良影响。

2 问题分析

2.1 问题背景分析

对岩石流体系统的弛豫特性进行研究，一般需从弛豫数据反解出其弛豫谱，对这样一个问题的求解实际上是解第一类Ferdhoml积分方程，这是一个非适定问题[1]，也就是说，存在许多不同的T2横向弛豫时间分布都能相当好的拟合到原始回波衰减曲线，这种相似性增大了反演过程的难度。处理此问题的做法一般可以分为两步：（1）预先给定一系列T2j的分布，使其尽可能多的均匀覆盖测试的驰豫时间范围(T2min,T2max)，进而经过试算，来确定布驰豫时间T2j的种类数k；（2）确定k后，不给定T2j的分布，通过算法直接求解各类孔隙的T2驰豫时间T2j以及各类孔隙在总孔隙中所占的份额fj。

在实际问题中，数据的测量等方面都存在随机误差，表现为反演求解过程中的信噪比（SNR，其值越小表示误差影响越大），其对反演求解的稳定性有着很大的影响，在算法实现的过程中需要对其进行降低和克服处理。

2

2.2 问题具体分析

根据分析可知，问题的解决在于如下几个方面：

（1） 怎样人为地指定T2j系列？

通过查阅相关文献，发现文献中均提到以下三种最为常见的布点方法[1,2,3]：对数均匀布点;2的幂指数方式布点;线性均匀布点。因此，本文主要考虑以上三种布点方法。

（2） 当人为的指定一系列T2j的值后，怎样求解各类孔隙在总孔隙中所占份额fj，从而确定T2谱？

对于给定的横向驰豫时间T2分布T2j(对数均匀布点、幂指数布点、线性均匀布点)，欲确定T2谱，则需对磁化强度信号y(t)和t进行拟合，但考虑到在此实际问题中，解fj应该是非负的，可以考虑非负最小二乘拟合。考虑到T2谱中每一个峰对应一个种类孔隙的特征（包含弛豫时间和该类孔隙所占份额），因此T2谱中峰的个数就等于驰豫时间T2j的种类数k。对于k的选取，可以根据峰值越多和T2谱连续的原则来选取。

（3） 若不考虑T2分布，如何直接求解各类孔隙的T2驰豫时间T2j以及各类孔隙在总孔隙中所占的份额fj？

在上述采用非负最小二乘法可以得到驰豫时间T2j的种类数k，同时还可以得到第j类孔隙在总孔隙中所占的份额fj和第j类孔隙的T2驰豫时间T2j的大致范围，考虑到目标函数为非线性函数，基于遗传算法具有对可行解的广泛性，采用遗传算法对其进行求解。但由于遗传算法往往只能找到次优解，因此为了验证结果的可靠性和精确性，可以采用蒙特卡洛算法对其进行验证。 （4） 如何通过算例分析误差给前面两个问题带来的影响？

可以将题附表中原始数据代入前面两个问题所建立的模型进行求解。将两者的得到的结果进行比较，分析误差给前面两个问题带来的影响，由于测量数据会有一定的误差，因此，误差的影响是可以体现的。 （5） 如何克服误差带来的不良影响？

考虑在目标模型中加入误差影响因素，采用惩罚函数的方法，可以采用模 平滑和曲率平滑的方式对模型算法进行修正。

3 问题假设

（1） 岩石核磁共振测得磁化强度信号为一系列单个孔隙磁化强度信号的叠加； （2） 岩石中孔隙的T2驰豫时间，通常范围为0.1msT2j10000ms；

（3） T2谱仅与各类孔隙的驰豫时间T2j和各类孔隙在总孔隙中所占的份额fj

有关。

3

4 名词解释

4.1 核磁共振现象[2] 对于被磁化后的核自旋系统，如果在垂直于静磁场的方向加一个射频场，根据量子力学原理，核自旋系统将发生共振吸收现象，即处于低能态的核自旋将通过吸收射频场提供的能量，跃迁到高能态。 4.2 核磁共振[2] 一种有效的岩心分析方法。利用核磁共振直接观测到对应孔隙流体含量的信号幅度和反映流体自身性质及其所处孔隙环境的弛豫特征与扩散行为。通过对这些观测信息的分析，可以得到样品总孔隙度、有效孔隙度、孔径分布、渗透率，直至流体饱和度。 4.3 弛豫[3] 在射频场施加以前，系统处于平衡状态，射频场作用期间，磁化矢量偏离静磁场方向;射频场作用结束后，核自旋从高能级的非平衡状态恢复到低能级的平衡状态，宏观磁化矢量恢复到平衡状态的过程。 4.4 弛豫时间谱[3] 在实际测量过程中，获取的衰减曲线由许多不同孔隙中流体衰减信号的叠加而成的。采用数学反演技术，可以计算出不同大小孔隙中的流体所占的份额。 4.5 信噪比SNR[2]

从测量数据中估计出的信号强度与噪音之比，定义为第一个回波的幅度值除以误差矢量的标准差。

5 符号说明

： T2： y(t)： n：

k：

fj： m：

T2jt：

采集时间；

回波间隔时间； 横向驰豫时间； 磁化强度信号；

测量的到的回波个数；

驰豫时间T2j（孔隙）的种类数； 第j类孔隙在总孔隙中所占的份额；

在(T2min,T2max)内，横向弛豫时间T2假定布点数； 第j类孔隙的T2驰豫时间，通常范围为0.1msT2j10000ms ： 。

6 问题一：非负多元线性拟合模型（Ⅰ）

6.1 模型Ⅰ分析

岩石核磁共振分析的关键是如何从磁化强度信号y(t)表达式中解出各类孔隙的T2驰豫时间T2j以及各类孔隙在总孔隙中所占的份额fj，称之为T2谱。对于给定的横向驰豫时间T2分布T2j(对数均匀布点、幂指数布点、线性均匀布点)，

4

欲确定T2谱，则需对磁化强度信号y(t)和t进行拟合，但考虑到在此实际问题中，解fj应该是非负的，于是采用非负最小二乘法进行拟合。对于T2谱，其中每一个峰对应一个种类孔隙的特征（包含弛豫时间和该类孔隙所占份额），因此T2谱中峰的个数就等于驰豫时间T2j的种类数k。对于k的选取，可以根据峰值越多和T2谱连续的原则来选取。

6.2 模型Ⅰ假设

（1） 不考虑测量数据的误差；

（2） 假设预先给定T2驰豫时间分布T2j；

6.3 模型Ⅰ建立

根据核磁共振理论，单个空隙中的磁化强度信号的衰减满足单指数衰减规律，但由于岩石内部是一系列大小不等的孔隙群体组成，所以在岩石核磁共振中测得的中磁化强度信号y(t)为一系列单个孔隙磁化强度信号的叠加，描述为：

mtT2j

y(t)j1fje,tn (1)

其中n，m分别为测量到的回波个数和横向弛豫时间T2分布个数，fj为第j类孔隙在总孔隙中所占的份额，T2j为第j类孔隙的T2驰豫时间（通常范围为0.1msT2j10000ms），为回波间隔时间，t为采集时间。

如果预先给定T2驰豫时间分布T2j，从（1）式可知，只需求出fj，便可求得T2谱。考虑将y(t)的全部采样值构成的向量

YAF

Y可以写成如下矩阵形式：

(2)

其中，磁化强度信号矩阵Y(y(t1),,y(tn))T，第j类孔隙在总孔隙中所占的份额矩阵F(f1,,fm,)，矩阵A表示如下：

1T21eAtnT21ett1T2me

tnT2me

欲求T2谱，只需求得F。从物理意义上讲,方程（1）的解

可以建立非负约束：

fjfj应该是非负的，

(3)

0

基于上述分析，需要反演出T2和fj，得到横向弛豫时间T2和fj的关系图, 即T2谱图。于是将方程（1）的求解原问题转化为如下带非负约束的优化问题。

n

mins.t.Zi1(yiAF)2 (4)

F05

其中，T2驰豫时间分布T2j由布点方法给定。

6.4 模型Ⅰ求解

6.4.1. T2驰豫时间分布T2j

（1） 对T2驰豫时间采用对数均匀布点

考虑在(T2min,T2max)上对数均匀取m个点，以a(一般取10)为底，有如下约束：

T2ja(logaT2maxlogaT2min)j/m(T2maxT2minj)m (5)

（2） 对T2驰豫时间采用幂指数布点

考虑在(T2min,T2max)上采用b(一般取2)幂指数布点，有如下约束：

T2jb1jj1,2,,mm (6)

bT2min,bT2max（3） 对T2驰豫时间采用线性均匀布点

考虑在(T2min,T2max)上采用线性均匀布点，有如下约束：

6.4.2. 不同驰豫时间分布下fj的模型求解

首先，横向弛豫时间分布点数m的确定。 当布点个数较少时，弛豫时间谱的分辨率较低，不能精细的反映储层孔隙结构；当布点个数较多时，参与反演的弛豫分量的个数较多，可以提高弛豫谱对孔隙结构的分辨率。但是，布点数越多需要反演的矩阵的列数就越大，算法的计算速度就越慢。因此为了保证弛豫时间谱的真实性,同时为了反演的速度,在反演过程中应该保证一定的分布点数和覆盖范围。通过查阅相关文献得到：分布点数m为32～ 时较为合适,布点区间应该跨越3个数量级[5]。

其次，驰豫时间T2j的种类数k的确定。

如果T2选取的布点数m取值不同，对应的谱线将会不一样，从而图中得到的峰的个数将会不一样。由于T2谱中每一个峰对应一个种类孔隙的特征（包含弛豫时间和该类孔隙所占份额），因此T2谱中峰的个数就等于驰豫时间T2j的种类数k。对于k的选取，可以根据峰值越多和T2谱连续的原则来选取。

然后，对三种不同弛豫时间分布，采用MATLAB7.10优化工具箱分别进行求解，求解如下：

（1） 对数（以10为底）均匀布点的求解

通过MATLAB7.10编程（源程序见附录1），选取不同的横向弛豫时间T2分布个数m值，其中m8,20,40,，得到结果如下图所示：

T2jTmaxTminmjj1,2,m (7)

6

图 1对数均匀布点反演得到的T2谱

从上图分析可知：m取值不同，得到的峰值个数不一定相同，每个峰值对应的f也不同。根据以上分析，从图中看出，最多的峰值个数为10个，于是得到驰豫时间T2j的种类数k的最优值：k10。为了更好的体现第j类的孔隙对应的弛豫时间T2j以及其在总孔隙中所占的份额，由以上分析可知分布点数为32～ 时较为合适，不失一般性，选取T36对应的k10个峰值数据画图如下：

图 2对数均匀布点反演得到孔隙的弛豫时间与其所占份额

从上图可以看出：第9类孔隙在总孔隙中所占的份额最大，第1类孔隙在总孔隙中所占的份额最小。

（2） 幂指数（以2为例）布点的求解

通过MATLAB7.10编程（源程序见附录1），考虑到T2(0.1,10000)，单位为毫秒，于是取分布点数为14，得到T2谱如下，可以看出共四个峰值，因此k4。

7

图 3幂指数布点所得T2谱

（3） 线性均匀布点的求解

m通过MATLAB7.10编程（源程序见附录1），选取不同的m值，其中

，得到结果如下图所示。 16,32,

图 4幂指数布点反演得到的T2谱

同理，从上图分析可知：图中只有一个峰值，因此k的最优值：k1。这是

由于线性均匀布点在m16,32,时，取到的第一个点便已超过前9个峰值所对应的弛豫时间。

最后，对三种布点方法进行比较。

根据以上三种不同布点方式，通过相应的结果分析得到：对数均匀布点在该实例中布点最为合理，且可以得到横向弛豫时间T2分布个数m、第j类孔隙在总

8

孔隙中所占的份额

fj和第j类孔隙的T2驰豫时间T2j的大致范围。后两者得到的

结果不是很合理。其中，对于线性均匀布点，若想得到完整的峰的个数，则需要的布点数非常多 , 但是可能会造成计算量太大 , 很难实现。

6.5 模型Ⅰ评价

6.5.1 优点

（1） NNLS方法数字稳定性好、 容易实现；

（2） 对于解出现负数的情况，NNLS解比方程的精确解更有意义； （3） NNLS方法具有有限步收敛、计算速度快的优点，虽然得不到准确的f,

但从反演结果，可以确定弛豫数据的组分数，而且还可以得到f的大致范围。 6.5.2 缺点

（1） 对于线性均匀布点，要求布点数要非常多,计算量太大 ,很难实现

（2） 对于对数均匀布点和幂指数布点，在T2分布不连续 , 尤其是T2组分

离散且分布较宽时 , 反演得到的T2谱分辨率较低。

7 问题二：遗传算法 （Ⅱ）

7.1 模型Ⅱ分析

在不预先给定T2驰豫时间分布T2j的情况下，考虑到上述目标函数为非线性函数，相应的变量具有边界值，于是就转化为一个约束非线性规划问题的求解，由于问题具有一定的复杂性，所以传统的算法可能会陷入局部最优或者无法求解。于是考虑到使用遗传算法来进行求解，该算法采用群体启发式随机搜索的方式，且对目标函数的可微性等不具备特殊要求，比较适合此个问题的求解[6]。 然而遗传算法容易出现过早收敛的现象，同时对遗传算法的解往往没有一个有效的精度和可行度的定量分析方法。因此为了验证结果的可靠性和精确性，可以采用蒙特卡洛算法对其进行验证。

7.2 模型Ⅱ符号说明

Pm： 变异概率；

交叉概率；

进化参数群体规模；

第j类孔隙的T2驰豫时间下限； 第j类孔隙的T2驰豫时间上限；

第j类孔隙在总孔隙中所占的份额下限； 第j类孔隙在总孔隙中所占的份额上限；

： M: T2jmin： T2jmax： fjmin： fjmax：

Pc9

7.3 模型Ⅱ假设

（1） 不考虑测量数据的误差；

（2） 假设不预先给定T2驰豫时间分布T2j；

7.4 模型Ⅱ算法的建立

7.4.1 目标函数的建立

基于模型Ⅰ结果分析，可以得到驰豫时间T2j的种类数k的值、第j类孔隙在总孔隙中所占的份额fj和第j类孔隙的T2驰豫时间T2j的大致范围。因此 ,欲同时反演出 T2和fj，得到横向弛豫时间 T2和fj的关系图（T2谱），可以将原问题转化为带非负约束的优化问题。从而需要在模型Ⅰ的基础上加上以下约束： （1） 第j类孔隙在总孔隙中所占的份额fj取值约束：fjminfjfjmax

T2jminT2jT2jmax （2） 第j类孔隙的T2驰豫时间T2j的取值约束：

由此建立的改进后的二次规划模型（共有2m个未知数）如下：

nktiT2jminZi1(yij1fje)2

(份额取值约束)1jk（弛豫时间约束）(8)

fjminfjfjmaxs.t.TT2jT2jmax2jmin7.4.2 遗传算法的基本思想

遗传算法（Genetic Algorithms，简称 GA）是一种基于自然选择原理和自然遗传机制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目标的优化。遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终得到最优解或准最优解。它必须做以下操作：初始群体的产生、求每一个体的适应度、根据适者生存的原则选择优良个体、被选出的优良个体两两配对，通过随机交叉其染色体的基因并随机变异某些染色体的基因后生成下一代群体，按此方法使群体逐代进化，直到满足进化终止条件。 7.4.3 遗传算法实现方法

---------------------------------------------------------------------------------------------------------------- step1 随机选取一组符合变量边界约束的fj和T2j的值，对初始值进行编码，得

到初始染色体。将初始种群对应的目标函数值作为种群初始适应度； step2 将初始的种群进行交叉操作和变异操作，得到子代的染色体； step3

step4

step5

step6

通过比较子代染色体的适应度，即所对应的目标函数值大小，来筛选出较优的子代染色体。分别作为目前的较优染色体和最大适应度；

如果没有达到最大遗传代数，以目前的较优染色体和适应度重复step2和step3过程；

如果达到了最大的遗传代数，或者个父代与子代染色体的变化值已经达到了最大精度，退出循环；

对最优染色体进行反编码，输出最优的fj和T2j值，和对应的目标函数值即最大适应度。

10

7.4.4 遗传算法的流程图

选取初始种群，编码交叉、变异操作选择得到符合要求的种群是是否达到精度要求否是否达到最大代数是输出结果，退出否

图 5遗传算法流程图

7.5 模型Ⅱ求解

采用MATLAB7.10优化工具箱中遗传算法求解器进行求解（源程序见2）。为了得到更精确的结果，希望得最优解的搜索能够充分的进行，于是将遗传的代数设为200代，总群数量为30。

利用MATLAB7.10将数据及目标函数进行处理之后，编程，调用优化工具箱求得在分类数k10的时候，得到第j类孔隙在总孔隙中所占的份额fj和第j类孔隙的T2驰豫时间T2j如下表：

表 1 去噪信号各类孔隙弛豫时间与其所占份额

T2j1 2 3 4 5 6 7 8 9 10

0.62 0.63 2.25 2.26 4.16 4

fj（10） 0.86 0.5 0.8 1.2 1.5

7.86 18.36 36.59 100.27 5505.33

1.7 3 2.5 7.8 1.3

f

为了更直观的体现第j类孔隙在总孔隙中所占的份额持豫时间T2j之间的关系，将T2谱画图如下：

11

j和第j类孔隙的T2

图 6 遗传算法得到的T2谱

7.6 模型Ⅱ结果验证

由于遗传算法可能会收敛过快，而且对其解得精度和可信度目前还没有较为有效的定量分析方法，为了验证结果的可靠性和精确性，采用蒙特卡洛算法随机搜索遗传算法求解出来的解的领域对应的目标函数值，从而实现对遗传算法的验证。

7.6.1 蒙特卡洛算法设计 蒙特卡洛算法具体过程如下：

-------------------------------------------------------------------------------------------------------------- step1 将遗传算法的求解结果作为最优解，其求解到的目标函数的值为

最优目标值；设定随机检测的次数； step2 产生一组符合要求的随机数，加在最优解上，同时保证得到的解

向量在约束范围内； step3 利用新产生解向量来求解目标函数值； step4 若step2中得到的目标函数值小于最优目标值,则取代原来值；否

则进行step5； step5 重复step2，step3直到随机检测完成； step6 输出搜索得到的目标函数值。

--------------------------------------------------------------------------------------------------------------

7.6.2 蒙特卡洛算法实现和求解

通过MATLAB7.10编程（源程序见附录3）实现蒙特卡洛算法，对遗传算法得出来的最优解进行检验，随机检验的次数设为1000，结果并没有发现该解的周围还有更优的解向量，于是说明此次求解结果具有一定的合理性和精确性。

12

7.7 模型Ⅱ评价

7.7.1 优点

遗传算法对该类复杂非线性规划问题的求解具有一定的适用性，在求解过程中不易于陷入局部最优解，也不会受到目标函数表达形式的，对该类问题的求解具有通用性。 7.7.2 缺点

遗传算法容易出现过早收敛的现象，同时对遗传算法的解往往没有一个有效的精度和可行度的定量分析方法。

8 问题三：改进规划模型（Ⅲ）

8.1 模型Ⅲ分析

在实际测量过程中，测量仪器受到电子部件和环境的影响，不可避免的在测量数据中会存在误差，产生随机噪声，噪声的大小对反演结果会带来很大的影响。本文中，随机误差的影响可以表现为信噪比（SNR）。其中，信噪比的值越大表明误差影响越小，相反，信噪比的值越小表明误差影响越大；无噪声(SNR=∞)为理想回波。

NNLS等一般算法在没有误差，或者在信噪比较高的情况下计算较为准确，但在有测量误差比较大的情况下，反演得到的结果振荡比较厉害；对于低信噪比的情况，会出现基线偏离的情况，导致孔隙度估计不准。

考虑测量数据的误差，通过算例分析误差给问题一和问题二的影响，未经去噪处理的原始数据代入以上模型算法中求解，分别得到结果画图如下。 8.1.1 对于问题一

图 7 NNLS算法对原始数据处理结果

上图为NNLS算法对原始数据的处理结果，通过与图1的对比，发现在不

13

同布点情况下，

fj的峰值个数只有3个且被相互分割开来，同时

fj的连续性较

差，在分布的中间部分会连续地取到0值。可见岩石空隙的种类数没有被很好的区分；在预先假定T2j分布、m和k未知的情况下，增加了k被忽略的可能性，误差作用对NNLS影响较大。 8.1.2 对于问题二

通过MATLAB得到原始信号各类孔隙弛豫时间与其所占份额数据如下表，T2谱如下图。

表 2 原始信号各类孔隙弛豫时间与其所占份额

T2fj

1 2 3 4 5 6 7 8 9 10 0. 0.63 1.23 2.26 4.16 7.86 18.36 36.59 100.27 5505.33

3.00

2.43

7.76

1.30

4（10） 0.43 0.43 0.75 1.17 1.43 1.69 j

图 8 遗传算法对原始数据处理结果

上图为布点数采用理论值k等于10时，通过遗传算法求解得到的各个T2j和对应的fj的谱相图，与图6对比发现只有部分T2j和fj的取值有轻微变动，说明误差对遗传算法计算的影响较小，而且还可能是因为在k取定、T2j分布未定的情况下，本文的信噪比SNR约为92，为较大值，误差作用本身不明显。

8.2 模型Ⅲ建立

通过上述对比分析，发现误差的影响（特别是当信噪比SNR较小时）是不可忽略的，为此需要对模型算法进行改进，在目标模型中加入误差影响因素，采用惩罚函数的方法，从而使使目标函数变为：

nmtiT2j

minZi1(yij1fje)B(fj)22 (9)

14

其中，B(fj)为惩罚项（也叫平滑项），为平滑参数，平滑项的构造有多种，

常用的有[1,4]：

模平滑为:

nmtiT2jm

曲率平滑为:

minZi1(yij1fje)2j1fj2 (10)

nmtiT2jm

minZi1(yij1fje)2j12fjfj1fj12 (11)

对于此模型的求解，可以利用各种算法求出其控制方程，最后通过反复迭代可求解。

8.3 模型Ⅲ求解

8.3.1 平滑项的构造方法的选取

考虑到模平滑对原始回波中的噪声比较敏感，而曲率平滑的算法很稳定，原始回波串质量较差时，也能得出较好的结果。于是在求解过程中选用曲率平滑。

8.3.2 平滑参数的选取

首先，一般值越大，解的结果越稳定。大的取值抑制了原始数据对T2弛豫谱的贡献，所以从这个角度来看，值要取得尽可能小；其次值的选取应该和噪声水平相匹配。比较合理的做法是，使值尽可能的小，并同时保证同一样品的两次弛豫实验数据的反演有相同结果。

经过计算可得本文中的数据的信噪比为约为92，这里根据参考文献，选取等于10-8[1]。

8.4 模型Ⅲ结果分析

将原始数据再代入修正模型，通过MATLAB求解，数据见下表，得到T2谱如下图：

表 3改进算法后各类孔隙弛豫时间与其所占份额

T2fjj

1 2 3 4 5 6 7 8 9 10 0.60 0.85 1.23 2.26 4.16 7.86 18.36 36.59 100.27 5505.33

2.97

2.39

7.73

1.24

（104） 0.33 0.50 0.75 1.13 1.43 1.62

15

图 9 考虑误差影响后得到结果

通过图9与图6对比，画图如下：

图 10 去噪修正与理论结果对比

可以发现，经过修正后的模型算法的求解，误差的影响大大减小，实际所得值与理论值相差很小，误差见下表，可见一定程度上较好的克服了误差影响。

表 4 误差表

jf

1 2 3 4 5 0.06

6 0.02

7 0.01

8 0.04

9 0.01

10 0.03

j误差

0.11 0.79 0.30 0.03

从上表可以看出，除了第二个数据出现异常外，其他的误差均比较小。 8.4.1 优点

优化算法引入了正则化项，抗噪能力比较强，可用于低信噪比驰豫谱反演. 8.4.2 缺点 曲率平滑对正则化参数的选取不是很敏感。

16

9 模型的评价

9.1

优点

（1） 对于解出现负数的情况，NNLS解比方程的精确解更有意义，数字稳定性

好、 容易实现；

（2） NNLS方法具有有限步收敛、计算速度快的优点，虽然得不到准确的f,但

从反演结果，可以确定弛豫数据的组分数，而且还可以得到f的大致范围； （3） 遗传算法对该类复杂非线性规划问题的求解具有一定的适用性，在求解过

程中不易于陷入局部最优解，也不会受到目标函数表达形式的，对该类问题的求解具有通用性； （4） 优化算法引入了正则化项，抗噪能力比较强，可用于低信噪比驰豫谱反演。

9.2 缺点

（1） 对于线性均匀布点，要求布点数要非常多,计算量太大 ,很难实现； （2） 对于对数均匀布点和幂指数布点，在T2分布不连续 , 尤其是T2组分 离

散且分布较宽时 , 反演得到的T2谱分辨率较低； （3） 遗传算法容易出现过早收敛的现象，同时对遗传算法的解往往没有一个有

效的精度和可行度的定量分析方法；

（4） 曲率平滑对正则化参数的选取不是很敏感。

1

参考文献

［1］ 潘克家，石油地球物理勘探中若干反问题研究，复旦大学，博士论文，2009。 ［2］ 李艳，复杂储层岩石核磁共振特性实验分析与应用研究，中国石油大学，硕

士论文，2007。 ［3］ 陈权，岩石核磁共振及其在渗流力学和油田开发中的应用研究，中国科学院，

2001。

［4］ 潘克家 陈华 谭永基，基于差分进化算法的核磁共振T2谱多指数反演，物

理学报，第57卷 第9期 2008年9月。

［5］ 童茂松 李莉 王伟男 范清华 姜亦忠 张加举 丁柱 崔杰,岩石激发极化弛

豫时间谱与孔隙结构、渗透率的关系, 地球物理学报, 第48卷第3期,2005。 ［6］ 杨海滨 韩 朋,，遗传算法浅析，文化教育，2006。

1

附录

1. 非负最小二乘拟合源程序 %% 数据处理

load 信号强度采样数据1.mat m=14;%横向弛豫时间的分布个数 % color=rand(1,3);

% lin_dot=linspace(0.1,10000,m+2); % lin_dot=lin_dot(2:m+1);

log2_dot=log2space(0,14,m+2); log2_dot=log2_dot(2:m+1);

% log_dot3=logspace(-1,4,m+2); % log_dot3=log_dot3(2:m+1); %y_test=y(7:26);

% y_test=y_test/max(y_test);

%y_test=(max(y_test)-y_test)/(max(y_test)-min(y_test)); %t_test=t(7:26);

% y=(y-min(y))/(max(y)-min(y)); %线性均匀选取弛豫时间分布 % for i=1:size(t,1) % for j=1:30

% A(i,j)=exp(-t(i)./test_dot(j)); % end % end %

% lsqnonneg(A,y) for i=1:size(t1,1) for j=1:m

A(i,j)=exp(-t1(i)./log2_dot(j)); end end

%[b1,bint1,r1,rint1,stats1]=regress(y',A); x=lsqnonneg(A,y_an);

%对数均匀选取弛豫时间分布 % for i=1:size(y,1) % for j=1:m

% A(i,j)=exp(-t(i)./log_dot3(j)); % end % end

% %[b2,bint2,r2,rint2,stats2]=regress(y',A); % x3=lsqnonneg(A,y);

semilogx(log2_dot,x','+-');

2

hold on;

axis([0.1 10000 -4000 90000]); % clear;clc;

2. 遗传算法MATLAB源程序 function f=fitness2(x)

load 信号强度采样数据1.mat num1=size(t1,1); m=10;%定义空隙的种类 for i=1:size(y_an,1) for j=1:m

A(i,j)=exp(-t1(i)./x(j)); end end f=0;

for i=1:num1 for j=1:m

f=f+A(i,j)*x(m+j); end end

% load 信号强度采样数据1.mat

lb=[0.34,0.63,1.23,2.26,4.16,7.86,18.36,36.59,100.27,5505.33]; ff=[.3 .5 .8 1.2 1.5 1.7 3 2.5 7.8 1.3]; lb=[lb ff];

ub=[0.63,1.23,2.26,4.16,7.86,18.36,36.59,67.35,165.98,10000]; fb=[inf inf inf inf inf inf inf inf inf inf]; ub=[ub fb];

options=gaoptimset;

options.Generations=200; options.PopulationSize=30;

[x,fval,exitflag]=ga(@fitness2,20,[],[],[],[],lb,ub,[],options); % m=5;

% lin_dot1=linspace(-1,4,m+2); % lin_dot1=lin_dot1(2:m+1); % semilogx(log_dot1,x1,'-*');

3. 蒙特卡洛算法对遗传算法检验的源程序 %% 利用蒙特卡洛对结果进行改进和优化 function x_best=mentekaro(x,loop) %x为待优化的变量,为列向量 %loop为随机模拟次数

3

x_best=x; for i=1:loop

x_rand=x_best+rand_num(x_best);

if (fitness2(x_rand)>fitness2(x_best)) x_best=x_rand; end end

disp('检验优化结果如下：') if(isequal(x_best,x)) disp('x_best==x');

disp('显然，原来的结果还是比较优的'); disp(fitness2(x_best)); else

disp('优化后的结果为：'); disp(x_best);

disp(fitness2(x_best)); end

%% 判断函数

function flag=isequal(x1,x2) flag=1;

for i=1:size(x1,1) if(x1(i)~=x2(i)) flag=0; end end

%% 产生符合要求的随机数 function delta=rand_num(x)

lb=[0.34,0.63,1.23,2.26,4.16,7.86,18.36,36.59,100.27,5505.33,.3,.5,.8,1.2,1.5,1.7,3,2.5,7.8,1.3]; delta=abs(x-lb); for i=1:size(x,1)

delta(i)=rand*delta(i); end

%% 赋值函数

% function f=fval(x) % %初始化数据

% % load 信号强度采样数据.mat % % num1=size(t,1); % % num2=size(x,1);

% % log_dot=logspace(-1,4,num2+2); % % log_dot=log_dot(2:num2+1); % % for i=1:size(y,1)

4

% % for j=1:num2

% % A(i,j)=exp(-t(i)./log_dot(j)); % % end % % end

% % delta=sum((A*x-y).^2); % load 信号强度采样数据1.mat % num1=size(t1,1); % m=10;%定义空隙的种类 % for i=1:size(y_an,1) % for j=1:m

% A(i,j)=exp(-t1(i)./x(j)); % end % end % f=0;

% for i=1:num1 % for j=1:m

% f=f+A(i,j)*x(m+j); % end % end

5