中南大学范文.doc

1、2010 高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料） ，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从 A/B/C/D 中选择一项填写）： C 我们的参赛报名号为（如果

2、赛区设置报名号的话）： 80 所属学校（请填写完整的全名）： 中南大学 参赛队员 (打印并签名) ：1. 欧阳彬 2. 韩帅 3. 王佳 指导教师或指导教师组负责人 (打印并签名)： 日期： 2010 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2010 高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：1核磁共振驰豫信号反演问题摘 要针对岩石核磁共振分析中横向弛豫时间 谱的反演问题，本文分驰豫时间2T分布

3、预先给定与否两种情况进行探究，分别得到了相应的结果，最后考虑误差影响提出了修正模型算法。预先给定驰豫时间分布时对 谱的求解。在 分布给定下，欲求 谱，只22j 2T需确定相应的 。考虑到实际中解 应该是非负值，对磁化强度信号 和jf jf ()yt进行非负最小二乘法拟合，建立非负多元线性拟合模型（） 。通过tMATLAB7.10 对不同的布点个数进行试算求解，得到三种布点方式 (对数均匀布点、幂指数布点、线性均匀布点)的 谱，并对结果分析得到：对数均匀布点2T在该实例中布点最为合理，能区分 10 种类型弛豫时间对应的孔隙。驰豫时间分布未知时对 谱的求解。基于模型，可以得到较为准确的驰2豫时间

4、的种类数 ，以及 和 的大致范围。考虑到目标函数为非线性函2jTkjfj数，并且变量具有边界值，为典型的二次约束规划问题，由于遗传算法求解这种问题不会陷入局部最优解且不会受到目标表达式形式影响的优点，因此采用遗传算法（）对其进行求解，通过 MATLAB7.10 求得 谱（表 1） 。为了验2T证解的精度，用蒙特卡洛算法在最优解邻域内进行了 1000 次随机检验，结果表明该最优解具有可靠性。考察随机误差对反演结果产生的影响，首先求得信噪比 SNR 约为 92，然后将以原始数据求解得到的 谱与理论结果对比，发现：误差对模型有较大2T的影响，使得模型对空隙种类的区分能力下降，驰豫时间种类 k 由 1

5、0 种变为 3种；误差对模有较小的影响， 谱 值 10 个中有 4 个产生了波动，其中 3jf的波动都不超过 20%，表明遗传算法有一定的抗噪能力。为克服误差影响，提出改进规划模型（） ，在目标模型中加入惩罚函数，考虑到曲率平滑修正方法比较稳定，即使原始回波串质量较差时，也能得出较好的结果，采用曲率平滑的方式减小模型算法对噪声的敏感性。以原始数据求解验证，将所得的 谱与理论结果进行误差分析（表 4） ，发现二者较为接近2T（图 10） ，说明该模型对信噪比较低的情况也可适用。最后，对三个模型进行了评价。关键词：非负最小二乘 遗传算法 蒙特卡洛 曲率平滑21 问题重述在核磁共振测井中, 一般采用

6、 CMPG 方法测量自旋回波串，纵向驰豫时间和横向驰豫时间 都是用来描述自旋回波信号的驰豫特征。由于 谱能够提1T2T 2T供许多岩石物性和流体特性的信息，越来越受到人们的关注。根据核磁共振理论，单个空隙中的磁化强度信号的衰减满足单指数衰减规律，但由于岩石内部是一系列大小不等的孔隙群体组成，所以在岩石核磁共振中测得的中磁化强度信号 为一系列单个孔隙磁化强度信号的叠加，因此 可描述为()yt ()yt21(),jtmTjytfen（其中 为第 类孔隙在总孔隙中所占的份额， 为第 类孔隙的 驰豫jf 2j 2T时间，通常范围为 ， 为回波间隔时间， 为采集时间）20.10jsTst岩石核磁共振分析

7、的关键是如何从上式中解出各类孔隙的 驰豫时间 以2j及各类孔隙在总孔隙中所占的份额 ，称之为 谱。jf2T分析题意可知，欲得到横向驰豫时间的 T2 谱，有以下 3 个问题需要解决：（1） 预先给定 驰豫时间分布 ，典型取法为在 内对数均匀选取22j min2ax(,)个点，称之为驰豫时间布点，也可考虑用 2 的幂指数布点、线性均匀布点等方m式。基于上述几种布点方式建立求解 谱的数学模型。（2） 假设不预先给定 驰豫时间分布 ，建立求解 谱的数学模型。2T2jT2T（3） 考虑到测量数据的误差，通过算例分析误差给前面两个问题带来的影响；并考虑如何克服误差带来的不良影响。2 问题分析2.1 问题背

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

9、测量等方面都存在随机误差，表现为反演求解过程中的信噪比（SNR，其值越小表示误差影响越大） ，其对反演求解的稳定性有着很大的影响，在算法实现的过程中需要对其进行降低和克服处理。32.2 问题具体分析根据分析可知，问题的解决在于如下几个方面：（1） 怎样人为地指定 系列？2jT通过查阅相关文献，发现文献中均提到以下三种最为常见的布点方法 1,2,3：对数均匀布点;2 的幂指数方式布点;线性均匀布点。因此，本文主要考虑以上三种布点方法。（2） 当人为的指定一系列 的值后，怎样求解各类孔隙在总孔隙中所占份2j额 ，从而确定 谱？jf2T对于给定的横向驰豫时间 分布 (对数均匀布点、幂指数布点、线性均

10、jT匀布点)，欲确定 谱，则需对磁化强度信号 和 进行拟合，但考虑到在此)yt实际问题中，解 应该是非负的，可以考虑非负最小二乘拟合。考虑到 谱中jf 2T每一个峰对应一个种类孔隙的特征（包含弛豫时间和该类孔隙所占份额） ，因此谱中峰的个数就等于驰豫时间 的种类数 。对于 的选取，可以根据峰值2T2jTk越多和 谱连续的原则来选取。（3） 若不考虑 分布，如何直接求解各类孔隙的 驰豫时间 以及各类2 2T2j孔隙在总孔隙中所占的份额 ？jf在上述采用非负最小二乘法可以得到驰豫时间 的种类数 ，同时还可以jk得到第 类孔隙在总孔隙中所占的份额 和第 类孔隙的 驰豫时间 的大致j jf22jT范围

11、，考虑到目标函数为非线性函数，基于遗传算法具有对可行解的广泛性，采用遗传算法对其进行求解。但由于遗传算法往往只能找到次优解，因此为了验证结果的可靠性和精确性，可以采用蒙特卡洛算法对其进行验证。（4） 如何通过算例分析误差给前面两个问题带来的影响？可以将题附表中原始数据代入前面两个问题所建立的模型进行求解。将两者的得到的结果进行比较，分析误差给前面两个问题带来的影响，由于测量数据会有一定的误差，因此，误差的影响是可以体现的。（5） 如何克服误差带来的不良影响？考虑在目标模型中加入误差影响因素，采用惩罚函数的方法，可以采用模平滑和曲率平滑的方式对模型算法进行修正。3 问题假设（1） 岩石核磁共振测

12、得磁化强度信号为一系列单个孔隙磁化强度信号的叠加；（2） 岩石中孔隙的 驰豫时间，通常范围为 ；2T20.10jmsTs（3） 谱仅与各类孔隙的驰豫时间 和各类孔隙在总孔隙中所占的份额2T2jTjf有关。44 名词解释4.1 核磁共振现象 2对于被磁化后的核自旋系统，如果在垂直于静磁场的方向加一个射频场，根据量子力学原理，核自旋系统将发生共振吸收现象，即处于低能态的核自旋将通过吸收射频场提供的能量，跃迁到高能态。4.2 核磁共振 2一种有效的岩心分析方法。利用核磁共振直接观测到对应孔隙流体含量的信号幅度和反映流体自身性质及其所处孔隙环境的弛豫特征与扩散行为。通过对这些观测信息的分析，可以得到样

13、品总孔隙度、有效孔隙度、孔径分布、渗透率，直至流体饱和度。4.3 弛豫 3在射频场施加以前，系统处于平衡状态，射频场作用期间，磁化矢量偏离静磁场方向;射频场作用结束后，核自旋从高能级的非平衡状态恢复到低能级的平衡状态，宏观磁化矢量恢复到平衡状态的过程。4.4 弛豫时间谱 3在实际测量过程中，获取的衰减曲线由许多不同孔隙中流体衰减信号的叠加而成的。采用数学反演技术，可以计算出不同大小孔隙中的流体所占的份额。4.5 信噪比 SNR2从测量数据中估计出的信号强度与噪音之比，定义为第一个回波的幅度值除以误差矢量的标准差。5 符号说明： 采集时间；t： 回波间隔时间；： 横向驰豫时间；2T： 磁化强度信

14、号；()yt： 测量的到的回波个数；n： 驰豫时间 （孔隙）的种类数；k2jT： 第 类孔隙在总孔隙中所占的份额；jfj： 在 内，横向弛豫时间 假定布点数；mminax(,)2T： 第 类孔隙的 驰豫时间，通常范围为 。2jTj2 20.10jmsTs6 问题一：非负多元线性拟合模型（）6.1 模型分析岩石核磁共振分析的关键是如何从磁化强度信号 表达式中解出各类孔()yt隙的 驰豫时间 以及各类孔隙在总孔隙中所占的份额 ，称之为 谱。对2T2j jf2T5于给定的横向驰豫时间 分布 (对数均匀布点、幂指数布点、线性均匀布点)，2Tj欲确定 谱，则需对磁化强度信号 和 进行拟合，但考虑到在此实

15、际问题中，2Tyt解 应该是非负的，于是采用非负最小二乘法进行拟合。对于 谱，其中每一jf 2T个峰对应一个种类孔隙的特征（包含弛豫时间和该类孔隙所占份额） ，因此 谱2中峰的个数就等于驰豫时间 的种类数 。对于 的选取，可以根据峰值越多2jTk和 谱连续的原则来选取。2T6.2 模型假设（1） 不考虑测量数据的误差；（2） 假设预先给定 驰豫时间分布 ；2T2jT6.3 模型建立根据核磁共振理论，单个空隙中的磁化强度信号的衰减满足单指数衰减规律，但由于岩石内部是一系列大小不等的孔隙群体组成，所以在岩石核磁共振中测得的中磁化强度信号 为一系列单个孔隙磁化强度信号的叠加，描述为：()yt* ME

16、RGEFORMAT (1)21(),jtmTjtfen其中 ， 分别为测量到的回波个数和横向弛豫时间 分布个数， 为第nm2Tjf类孔隙在总孔隙中所占的份额， 为第 类孔隙的 驰豫时间（通常范围为j 2j） ， 为回波间隔时间， 为采集时间。20.10jsTst如果预先给定 驰豫时间分布 ，从（1）式可知，只需求出 ，便可求2jTjf得 谱。考虑将 的全部采样值构成的向量 可以写成如下矩阵形式：()yt Y* MERGEFORMAT (2)AF其中，磁化强度信号矩阵 ，第 类孔隙在总孔隙中所占1(),()TnYytt j的份额矩阵 ，矩阵 表示如下：1(,)mFf 1122221mnnmttT

17、TttTTeA 欲求 谱，只需求得 。从物理意义上讲,方程（ 1）的解 应该是非负的，2TFjf可以建立非负约束：* MERGEFORMAT (3)0jf基于上述分析，需要反演出 和 ，得到横向弛豫时间 和 的关系图, 2Tj 2Tjf即 谱图。于是将方程（ 1）的求解原问题转化为如下带非负约束的优化问题。2T6* MERGEFORMAT (4)21min().0niiZyAFst其中， 驰豫时间分布 由布点方法给定。2T2jT6.4 模型求解6.4.1. 驰豫时间分布22j（1） 对 驰豫时间采用对数均匀布点考虑在 上对数均匀取 个点，以 为底，有如下约minax(,)Tm(10)a一 般

18、取束：* MERGEFORMAT (5)2max2in(logl)/2max2 in()jTjj T（2） 对 驰豫时间采用幂指数布点2T考虑在 上采用 幂指数布点，有如下约束：min2ax(,)()b一 般 取* MERGEFORMAT (6)21min2ax1,jjTjT（3） 对 驰豫时间采用线性均匀布点2T考虑在 上采用线性均匀布点，有如下约束：min2ax(,)* MERGEFORMAT (7)maxin1,jTjm6.4.2. 不同驰豫时间分布下 的模型求解jf首先，横向弛豫时间分布点数 的确定。当布点个数较少时，弛豫时间谱的分辨率较低，不能精细的反映储层孔隙结构；当布点个数较多时

19、，参与反演的弛豫分量的个数较多，可以提高弛豫谱对孔隙结构的分辨率。但是，布点数越多需要反演的矩阵的列数就越大，算法的计算速度就越慢。因此为了保证弛豫时间谱的真实性,同时为了反演的速度,在反演过程中应该保证一定的分布点数和覆盖范围。通过查阅相关文献得到：分布点数 为 3264 时较为合适 ,布点区间应该跨越 3 个数量级 5。m其次，驰豫时间 的种类数 的确定。2jTk如果 选取的布点数 取值不同，对应的谱线将会不一样，从而图中得到2m的峰的个数将会不一样。由于 谱中每一个峰对应一个种类孔隙的特征（包含2弛豫时间和该类孔隙所占份额） ，因此 谱中峰的个数就等于驰豫时间 的种2T2jT类数 。对于

20、 的选取，可以根据峰值越多和 谱连续的原则来选取。k 2然后，对三种不同弛豫时间分布，采用 MATLAB7.10 优化工具箱分别进行求解，求解如下：（1） 对数（以 10 为底）均匀布点的求解7通过 MATLAB7.10 编程（源程序见附录 1） ，选取不同的横向弛豫时间分布个数 值，其中 ，得到结果如下图所示：2Tm8,204,6图 1 对数均匀布点反演得到的 T2 谱从上图分析可知： 取值不同，得到的峰值个数不一定相同，每个峰值对m应的 也不同。根据以上分析，从图中看出，最多的峰值个数为 10 个，于是得f到驰豫时间 的种类数 的最优值： 。为了更好的体现第 类的孔隙对2jTk10kj应的

21、弛豫时间 以及其在总孔隙中所占的份额，由以上分析可知分布点数为3264 时较为合适，不失一般性，选取 对应的 个峰值数据画图如36T10k下：图 2 对数均匀布点反演得到孔隙的弛豫时间与其所占份额从上图可以看出：第 9 类孔隙在总孔隙中所占的份额最大，第 1 类孔隙在总孔隙中所占的份额最小。（2） 幂指数（以 2 为例）布点的求解8通过 MATLAB7.10 编程（源程序见附录 1） ，考虑到 ，单位2(0.1,)T为毫秒，于是取分布点数为 14，得到 谱如下，可以看出共四个峰值，因此2T。4k图 3 幂指数布点所得 T2 谱（3） 线性均匀布点的求解通过 MATLAB7.10 编程（源程序见

22、附录 1） ，选取不同的 值，其中m，得到结果如下图所示。16,24m图 4 幂指数布点反演得到的 T2 谱同理，从上图分析可知：图中只有一个峰值，因此 的最优值： 。这k1k是由于线性均匀布点在 时，取到的第一个点便已超过前 9 个峰值所16,32m对应的弛豫时间。9最后，对三种布点方法进行比较。根据以上三种不同布点方式，通过相应的结果分析得到：对数均匀布点在该实例中布点最为合理，且可以得到横向弛豫时间 分布个数 、第 类孔隙2Tmj在总孔隙中所占的份额 和第 类孔隙的 驰豫时间 的大致范围。后两者得jf2j到的结果不是很合理。其中，对于线性均匀布点，若想得到完整的峰的个数，则需要的布点数非

23、常多 , 但是可能会造成计算量太大 , 很难实现。6.5 模型评价6.5.1 优点（1） NNLS 方法数字稳定性好、 容易实现；（2） 对于解出现负数的情况，NNLS 解比方程的精确解更有意义；（3） NNLS 方法具有有限步收敛、计算速度快的优点，虽然得不到准确的,但从反演结果，可以确定弛豫数据的组分数，而且还可以得到 的f f大致范围。6.5.2 缺点（1） 对于线性均匀布点，要求布点数要非常多,计算量太大 ,很难实现（2） 对于对数均匀布点和幂指数布点，在 分布不连续 , 尤其是 组分2T2T离散且分布较宽时 , 反演得到的 谱分辨率较低。7 问题二：遗传算法 （）7.1 模型分析在不

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

25、 类孔隙的 驰豫时间下限；2injTj2T： 第 类孔隙的 驰豫时间上限；max10： 第 类孔隙在总孔隙中所占的份额下限；minjfj： 第 类孔隙在总孔隙中所占的份额上限；ax7.3 模型假设（1） 不考虑测量数据的误差；（2） 假设不预先给定 驰豫时间分布 ；2T2jT7.4 模型算法的建立7.4.1 目标函数的建立基于模型结果分析，可以得到驰豫时间 的种类数 的值、第 类孔隙2jkj在总孔隙中所占的份额 和第 类孔隙的 驰豫时间 的大致范围。因此 ,欲jfTj同时反演出 和 ，得到横向弛豫时间 和 的关系图（ 谱） ，可以将原2Tj jf2T问题转化为带非负约束的优化问题。从而需要在模

26、型的基础上加上以下约束：（1） 第 类孔隙在总孔隙中所占的份额 取值约束：j jfminaxjjjff（2） 第 类孔隙的 驰豫时间 的取值约束：2T2j 22mT由此建立的改进后的二次规划模型（共有 个未知数）如下：* 211minax22m()(). 1ijtnkTijijjjZyfefst jkT A份 额 取 值 约 束（ 弛 豫 时 间 约 束 ）MERGEFORMAT (8)7.4.2 遗传算法的基本思想遗传算法（Genetic Algorithms，简称 GA）是一种基于自然选择原理和自然遗传机制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目标的优化。

27、遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终得到最优解或准最优解。它必须做以下操作：初始群体的产生、求每一个体的适应度、根据适者生存的原则选择优良个体、被选出的优良个体两两配对，通过随机交叉其染色体的基因并随机变异某些染色体的基因后生成下一代群体，按此方法使群体逐代进化，直到满足进化终止条件。7.4.3 遗传算法实现方法-step1 随机选取一组符合变量边界约束的 和 的值，对初始值进行编码，得jf2jT到初始染色体。将初始种群对应的目标函数值作为种群初始适应度；step2 将初始的种群进行交叉操作和变异操作，得到子代的染色体；step3 通过比较子代染色体的适应度，即所

28、对应的目标函数值大小，来筛选出较 优的子代染色体。分别作为目前的较优染色体和最大适应度；step4 如果没有达到最大遗传代数，以目前的较优染色体和适应度重复 step211和 step3 过程；step5 如果达到了最大的遗传代数，或者个父代与子代染色体的变化值已经达到 了最大精度，退出循环；step6 对最优染色体进行反编码，输出最优的 和 值，和对应的目标函数值jf2jT即最大适应度。7.4.4 遗传算法的流程图选取初始种群 ， 编码交叉 、 变异操作选择是否达到精度要求是否达到最大代数得到符合要求的种群是是输出结果 ，退出否否图 5 遗传算法流程图7.5 模型求解采用 MATLAB7.1

29、0 优化工具箱中遗传算法求解器进行求解（源程序见 2） 。为了得到更精确的结果，希望得最优解的搜索能够充分的进行，于是将遗传的代数设为 200 代，总群数量为 30。利用 MATLAB7.10 将数据及目标函数进行处理之后，编程，调用优化工具箱求得在分类数 的时候，得到第 类孔隙在总孔隙中所占的份额 和第 类10kj jf孔隙的 驰豫时间 如下表：2T2j表 1 去噪信号各类孔隙弛豫时间与其所占份额1 2 3 4 5 6 7 8 9 102jT0.62 0.63 2.25 2.26 4.16 7.86 18.36 36.59 100.275505.33（10 4jf）0.86 0.5 0.8

30、1.2 1.5 1.7 3 2.5 7.8 1.312为了更直观的体现第 类孔隙在总孔隙中所占的份额 和第 类孔隙的j jf2T持豫时间 之间的关系，将 谱画图如下：2jT2T7.6 模型结果验证由于遗传算法可能会收敛过快，而且对其解得精度和可信度目前还没有较为有效的定量分析方法，为了验证结果的可靠性和精确性，采用蒙特卡洛算法随机搜索遗传算法求解出来的解的领域对应的目标函数值，从而实现对遗传算法的验证。7.6.1 蒙特卡洛算法设计蒙特卡洛算法具体过程如下：-step1 将遗传算法的求解结果作为最优解，其求解到的目标函数的值为最优目标值；设定随机检测的次数；step2 产生一组符合要求的随机数，

31、加在最优解上，同时保证得到的解向量在约束范围内；step3 利用新产生解向量来求解目标函数值；step4 若 step2 中得到的目标函数值小于最优目标值, 则取代原来值；否则进行 step5；step5 重复 step2，step3 直到随机检测完成；step6 输出搜索得到的目标函数值。图 6 遗传算法得到的 T2 谱13-7.6.2 蒙特卡洛算法实现和求解通过 MATLAB7.10 编程（源程序见附录 3）实现蒙特卡洛算法，对遗传算法得出来的最优解进行检验，随机检验的次数设为 1000，结果并没有发现该解的周围还有更优的解向量，于是说明此次求解结果具有一定的合理性和精确性。7.7 模型评

32、价7.7.1 优点遗传算法对该类复杂非线性规划问题的求解具有一定的适用性，在求解过程中不易于陷入局部最优解，也不会受到目标函数表达形式的限制，对该类问题的求解具有通用性。7.7.2 缺点遗传算法容易出现过早收敛的现象，同时对遗传算法的解往往没有一个有效的精度和可行度的定量分析方法。8 问题三：改进规划模型（）8.1 模型分析在实际测量过程中，测量仪器受到电子部件和环境的影响，不可避免的在测量数据中会存在误差，产生随机噪声，噪声的大小对反演结果会带来很大的影响。本文中，随机误差的影响可以表现为信噪比（SNR） 。其中，信噪比的值越大表明误差影响越小，相反，信噪比的值越小表明误差影响越大；无噪声(

33、SNR= )为理想回波。NNLS 等一般算法在没有误差，或者在信噪比较高的情况下计算较为准确，但在有测量误差比较大的情况下，反演得到的结果振荡比较厉害；对于低信噪比的情况，会出现基线偏离的情况，导致孔隙度估计不准。考虑测量数据的误差，通过算例分析误差给问题一和问题二的影响，未经去噪处理的原始数据代入以上模型算法中求解，分别得到结果画图如下。8.1.1 对于问题一14图 7 NNLS 算法对原始数据处理结果上图为 NNLS 算法对原始数据的处理结果，通过与图 1 的对比，发现在不同布点情况下， 的峰值个数只有 3 个且被相互分割开来，同时 的连续性较jf jf差，在分布的中间部分会连续地取到 0

34、 值。可见岩石空隙的种类数没有被很好的区分；在预先假定 分布、 和 未知的情况下，增加了 被忽略的可能性，2jTmkk误差作用对 NNLS 影响较大。8.1.2 对于问题二通过 MATLAB 得到原始信号各类孔隙弛豫时间与其所占份额数据如下表，谱如下图。2T表 2 原始信号各类孔隙弛豫时间与其所占份额1 2 3 4 5 6 7 8 9 102j0.54 0.63 1.23 2.26 4.16 7.86 18.36 36.59 100.27 5505.33（10 4jf） 0.43 0.43 0.75 1.17 1.43 1.69 3.00 2.43 7.76 1.3015图 8 遗传算法对原始

35、数据处理结果上图为布点数采用理论值 等于 10 时，通过遗传算法求解得到的各个k和对应的 的谱相图，与图 6 对比发现只有部分 和 的取值有轻微变动，2jTjf 2jTjf说明误差对遗传算法计算的影响较小，而且还可能是因为在 取定、 分布未k2jT定的情况下，本文的信噪比 SNR 约为 92，为较大值，误差作用本身不明显。8.2 模型建立通过上述对比分析，发现误差的影响（特别是当信噪比 SNR 较小时）是不可忽略的，为此需要对模型算法进行改进，在目标模型中加入误差影响因素，采用惩罚函数的方法，从而使使目标函数变为：* MERGEFORMAT (9)2211min()()ijtmTij jiZy

36、feBfA其中， 为惩罚项（也叫平滑项） ， 为平滑参数，平滑项的构造有多()jBf种，常用的有 1,4：模平滑为:* MERGEFORMAT (10)22111min()ijtmmTij jiZyfefA曲率平滑为:* 2 21111in()ijtmTij jjji jyfeffAMERGEFORMAT (11)对于此模型的求解，可以利用各种算法求出其控制方程，最后通过反复迭代可求解。168.3 模型求解8.3.1 平滑项的构造方法的选取考虑到模平滑对原始回波中的噪声比较敏感，而曲率平滑的算法很稳定，原始回波串质量较差时，也能得出较好的结果。于是在求解过程中选用曲率平滑。8.3.2 平滑参数

37、 的选取首先，一般 值越大，解的结果越稳定。大的 取值抑制了原始数据对 T2弛豫谱的贡献，所以从这个角度来看， 值要取得尽可能小；其次 值的选取应该和噪声水平相匹配。比较合理的做法是，使 值尽可能的小，并同时保证同一样品的两次弛豫实验数据的反演有相同结果。经过计算可得本文中的数据的信噪比为约为 92，这里根据参考文献，选取等于 10-81。8.4 模型结果分析将原始数据再代入修正模型，通过 MATLAB 求解，数据见下表，得到 谱2T如下图：表 3 改进算法后各类孔隙弛豫时间与其所占份额1 2 3 4 5 6 7 8 9 102jT0.60 0.85 1.23 2.26 4.16 7.86 1

38、8.36 36.59 100.27 5505.33（10 4jf） 0.33 0.50 0.75 1.13 1.43 1.62 2.97 2.39 7.73 1.24图 9 考虑误差影响后得到结果通过图 9 与图 6 对比，画图如下：17图 10 去噪修正与理论结果对比可以发现，经过修正后的模型算法的求解，误差的影响大大减小，实际所得值与理论值相差很小，误差见下表，可见一定程度上较好的克服了误差影响。表 4 误差表j1 2 3 4 5 6 7 8 9 10误差jf0.11 0.79 0.30 0.03 0.06 0.02 0.01 0.04 0.01 0.03从上表可以看出，除了第二个数据出现

39、异常外，其他的误差均比较小。8.4.1 优点优化算法引入了正则化项，抗噪能力比较强，可用于低信噪比驰豫谱反演.8.4.2 缺点曲率平滑对正则化参数 的选取不是很敏感。19 模型的评价9.1 优点（1） 对于解出现负数的情况，NNLS 解比方程的精确解更有意义，数字稳定性好、 容易实现；（2） NNLS 方法具有有限步收敛、计算速度快的优点，虽然得不到准确的 ,f但从反演结果，可以确定弛豫数据的组分数，而且还可以得到 的大致范围；（3） 遗传算法对该类复杂非线性规划问题的求解具有一定的适用性，在求解过程中不易于陷入局部最优解，也不会受到目标函数表达形式的限制，对该类问题的求解具有通用性；（4）

40、优化算法引入了正则化项，抗噪能力比较强，可用于低信噪比驰豫谱反演。9.2 缺点（1） 对于线性均匀布点，要求布点数要非常多,计算量太大 ,很难实现；（2） 对于对数均匀布点和幂指数布点，在 分布不连续 , 尤其是 组分 离2T2T散且分布较宽时 , 反演得到的 谱分辨率较低；（3） 遗传算法容易出现过早收敛的现象，同时对遗传算法的解往往没有一个有效的精度和可行度的定量分析方法；（4） 曲率平滑对正则化参数的选取不是很敏感。1参考文献1 潘克家，石油地球物理勘探中若干反问题研究，复旦大学，博士论文，2009。2 李艳，复杂储层岩石核磁共振特性实验分析与应用研究，中国石油大学，硕士论文，2007。

41、3 陈权，岩石核磁共振及其在渗流力学和油田开发中的应用研究，中国科学院，2001。4 潘克家 陈华 谭永基，基于差分进化算法的核磁共振 T2 谱多指数反演，物理学报，第 57 卷 第 9 期 2008 年 9 月。5 童茂松 李莉 王伟男 范清华 姜亦忠 张加举 丁柱 崔杰,岩石激发极化弛豫时间谱与孔隙结构、渗透率的关系, 地球物理学报, 第 48 卷第 3 期,2005。6 杨海滨 韩 朋,，遗传算法浅析，文化教育， 2006。2附录1.非负最小二乘拟合源程序% 数据处理load 信号强度采样数据 1.matm=14;%横向弛豫时间的分布个数% color=rand(1,3);% lin_d

42、ot=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)-m

43、in(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:mA(i,j)=exp(-t1(i)./log2_dot(j);endend%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_do

44、t3(j);% end% end% %b2,bint2,r2,rint2,stats2=regress(y,A);% x3=lsqnonneg(A,y);semilogx(log2_dot,x,+-);3hold on;axis(0.1 10000 -4000 90000);% clear;clc;2.遗传算法 MATLAB 源程序function f=fitness2(x)load 信号强度采样数据 1.matnum1=size(t1,1);m=10;%定义空隙的种类for i=1:size(y_an,1)for j=1:mA(i,j)=exp(-t1(i)./x(j);endendf=0;

45、for i=1:num1for j=1:mf=f+A(i,j)*x(m+j);endend% load 信号强度采样数据 1.matlb=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=gaoptims

46、et;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 为随机模拟次数4x_best=x;for

47、 i=1:loopx_rand=x_best+rand_num(x_best);if (fitness2(x_rand)fitness2(x_best)x_best=x_rand;endenddisp(检验优化结果如下：)if(isequal(x_best,x)disp(x_best=x);disp(显然，原来的结果还是比较优的);disp(fitness2(x_best);elsedisp(优化后的结果为：);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;endend% 产生符合要求的随机数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