第三次模拟赛

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

2、校（请填写完整的全名）： 四川理工学院参赛队员（打印并签名）:1.郭姝霖徐高飞胡婷婷日期：2012 年 08月 29 日评阅编号（由评委团评阅前进行编号）:2012年四川理工学院暑期培训第三次模拟赛编号专用页评阅编号（由评委团评阅前进行编号）:核磁共振弛豫信号反演问题摘要本文针对核磁共振弛豫信号反演问题，在S0反演算法的基础之上，创新性地提出 了变换反演算法，解决了当数据信噪比较低时，反演结果的分辨率比较低，有可能造成 T谱的畸形的情形，建立了非负最小二乘法和DE的优化算法相结合的模型，很好的解 决了问题。针对问题一，本文首先选用对数、2的幕指数、线性均匀布点三种方式进行布点， 对T谱建立数学

3、模型，再将三种模型相互比较，在此基础之上，对模型进行修改，得出 ，2 ，一 了很好的结果，针对问题二，要求在弛豫时间未知的情况下对T谱进行求解，2针对问题三，误差的主要来源是测量误差，而会造成测量数据误差的有等待时间、 回波间隔、轻烃或稠油、顺磁物质以及采集方法等因素。而这些都是在测井的时候就已 经存在，无法避免，只能通过改进现有的计算模型以及相关的算法来使得结果更加精确。 当原始数据的信噪比SNR 80时，可选用SVD反演算法，然而为保证T反演结果的真 实性，本文选用30个T的布点数来建立改进的SVD反演算法，最后将非负最小二乘法 与差分进化的优化算法相结合，很好地排除了不同条件下造成的误差

4、。本文建立了反演 优化模型，并通过差分进化算法求解，最终得出在建模和算法上合理的更加精确的解决 方法。最后，本文对结果进行了分析，对模型的优缺点进行了评价，并对模型进行了推广和改 进。关键词：核磁共振SVD反演算法非负最小二乘法一、问题重述1.1背景核磁共振（NMR）在测井技术和岩心分析中已经得到广泛应用，正在为油气资源的勘 探和开发发挥重要作用。油气井中NMR技术提供的油气藏流体特性和储集层参数，如地 层孔隙度、孔径分布、束缚水与可动流体孔隙体积、渗透率、以及流体的扩散系数和黏 度等信息，都需要经过一个基本的反演处理，即NMR弛豫信号的多指数拟合，得到弛豫 时间的分布。在核磁共振测井中，一般

5、采用CMPG方法测量自旋回波串，纵向弛豫时间广和横向 弛豫时间T都是用来描述自旋回波信号的弛豫特征。由于T谱能够提供许多岩石物性和 流体特性的信息，越来越受到人们的关注。根据核磁共振理论，单个空隙中的磁化强度 信号的衰减满足单指数衰减规律，但由于岩石内部是一系列大小不等的孔隙群体组成， 所以在岩石核磁共振中测得的中磁化强度信号yt）为一系列单个孔隙磁化强度信号的 叠加，因此y（t）可描述为：yG）=Z f - e-T；j, t = n -tj=1 TOC o 1-5 h z 其中f为第j类孔隙在总孔隙中所占的份额，T为第j类孔隙的T弛豫时间，通常范 j一2 j 一2围为0.1ms T 为测量的

6、磁化强度衰减信号，A =A = e*，1 2 nij nxm |_-nxmf = (f, f，f)为弛豫时间T对应各点的幅度值。T (j = 1,2,.,m)为预先指定的T时1 2 m2 j2 j2间分布系列。本文采用在T . , T )区间内对数均匀地选取m个点，即弛豫时间布点，同 时也采用2的幕指数，线性均匀布点的三种方式。为使算法简便化，本文采用变换的反演算法，假设每个点的测量误差是随机独立且 标准差都一样，则可给定如下的目标函数：min x 2 = i=1使得s.t f 0,1 j 为 将(2)式代入(3)式有：则问题的解即为方程的解：3.2.2模型的求解f = AtCn x 1阶矩阵

7、。At C - At )+X I - C = Mt j(a - At +人 I ) C = J(2)(3)运用Matlab软件求解得C的值，再经过线性变换式回代而获得，选择合适的人， 以保证矩阵J - At +人L)的可逆，则求解得方程(1)的最小二乘解为：f = At - Ca - At +人I )】 j则分别可得如下如图4，图5，图6结果，其具体程序见附件1：由图重结果可以看出，f的值有正有负，这在物理意义上来说不符合逻辑，故本文 i对模型改进如下：T谱非负限制的实现，为了克服传统的缺点，不破坏T分布的连续性，防止T谱的.一.222畸形，本文米用迭代法逐渐消去负数的组分，其思想为：通过T域

8、到时域的变换后，将 第k项相邻点的T分布信息叠加到负数负数的T组分中，重新进行反演计算，直到满22 k足所有的f 0（j = 1,2质）为止。具体结果见下图：i由模型一的反演结果分析可知反演的分辨率低的问题，是因为T在运算前就已经4.1问题分析2 j选定，所以本模型中将T也看成未知数。并采用将T和s成为了要求解的未知数，将2 j2 j jT2和s.间的关系由线性转为非线性关系求解。为使求解过程简单，本文采用非线性拟 合求解；求解过程采用L M （高斯一牛顿）法。非线性拟合过程为：实际问题中，匕是非负的，于是采用非负最小二乘法进行拟合，首先设置C的初值，最后对初值进行修正，直到X 2相对于前一次

9、不再增长后者不再明 显增长时为止。针对T谱中每一个峰对应一个种类孔隙的特征，T谱中峰的个数就等于 弛豫时间T.的种类数。对于k的选取，本文在多次拟合的情况下，选用30个2的布 点数来建立改进的SVD反演算法，将非负最小二乘法与差分进化的优化算法相结合，很 好的排除了不同条件下造成的误差。4.2模型的建立及求解4.2.1模型的建立进行如下改进得:y(t)=* fe-tj=iy(t)= C e-ti/c2 + C e-ti/c4 + . . + C e-1./c2k 132 k -1C 即为模型-p的f，C即为T。通过上述变换，T与f间2 k-1i 2 k2j2 j i由问题一可以看出T.在运算前

10、由Matlab软件随机选定，而反演算法分辨率很低， 为提高本文的分辨率，本文对模型：其中，k = m / 2，的关系变为y与C之间的非线性关系，于是本文采用线性拟合求解。i i 的函数值在点G = 1,2, ., m设有实验数据(t, y )(i = 1,2,m)，寻找函数f (t, y)，使得函数在点x (i = 1,2,m)夕卜 )处的函数值与观测数据偏差的平方和达到最小，即：min (f (t, y )- y= 代iii其中， r Q=f (t, y)- y.。 可得系数的计算公式为：i=1i=1 (y - yiR 2 = 1 -i=1 (y- ii=1其中， = 1 lLyi，而R 2

11、越趋于1表明拟合效果越好。.i=1本文采用基于改进的L - M (高斯一牛顿法)算法对其进行求解，将其转化为求解误 差平方和最小。即有：min f (x )= 1 代=12 i 2记：r(x)= r (x), r()., r 1)1, A(x)= (Vr ., Vr (x)。12m1m本文采用改进的L-M算法，其形式为：s = I/t Co)j (w)+ pTL J JU)其中，比例系数p 0为常数，I是单位矩阵。编程求解的结果如下:五、问题三5.1问题分析针对问题三，由于测量方式和测量工具的差异，测量数据必然会存在误差。前两问 中考虑到测量数据的误差，均通过分析去噪信号，求解测量信号的尸频谱

12、。在实际检测 计算过程中，需要考虑测量误差对结果的影响。因此在数据处理时应尽可能的减小测量 误差的程度。本文将测量误差视为测量过程中各种干扰导致的结果。5.2模型的建立及求解信噪比SNR :从测量数据中估计出的信号强度与噪音之比，定义为第一个回波的幅 度值除以误差矢量的标准差。为了衡量测量误差的程度，此处引入信噪比（SNR），信噪比（SNR）按如下方式定义： SNR = 20% 产|2其中y表示有效信号，即去噪信号，表示原信号，即未经去噪处理的信号。SNR 数值越大，表明信号受干扰越少，即测量误差越小。本文采用计算机模拟来分析测量误差对于前面两个问题带来的影响。首先给定一个 具有双峰特征的T谱

13、分布y（t）= 200/16 + 800/1024 （单位：ms ），由频谱可计算出不同 时刻匕，计算出不同信号强度y,再对y加上一定大小的噪音得到，由数据可计算出 受干扰信号的信噪比，模拟信号的噪声由Matlab软件随机生成。对于模拟所得的信号， 在预先给定弛豫时间分布T.和未预先给定弛豫时间分布两种情况下，分别采用非负最 小二乘算法和差分进化算法反演频谱。最终通过对比得出测量误差对于这两种情况的频 谱计算结果的影响。由问题一的结果分析可得，采用对数均匀布点方式曲线拟合相似度较高，因此本问 采用对数均匀布点方式，取L = 2 j（j = 1,2, ,15）（单位：ms ）。分别在信号 y（t

14、）= 200/16 + 800/1024上加以不同程度的噪声，计算模拟信号的信噪比SNR，在利 用问题一中的NNLS算法对模拟信号进行反演，得出反演曲线的残差和对应的f与原信 号f.的方差作为评判反演结果好坏的指标。j（1）加为强度为0-1的随机噪音利用Matlab计算结果如下表1：表1强度为0-1随机噪音信号反演结果信噪比112.2516残差34.5041TT 22T 23T 24T 25T 26T 210T211） 时间（ms ）24163225610242048f.f2ffff6f10f数值12.78644.6748185.03178.47050.6549794.40954.5702J

15、.方差76.2473（2）加为强度为0-5的随机噪音利用Matlab计算结果如下表2：表2强度为0-5随机噪音信号反演结果信噪比86.3802残差561.2619TTTTTT孕 23242529210时间ms )416325121024ffffffj345910数值5709524127.902226.7878224.7314434.0134f j方差38744.9(1)加为强度为0-20的随机噪音 利用Matlab计算结果如下表3：表3强度为0-20随机噪音信号反演结果信噪比64.6336残差8520.4T2jT 24T 29T 210时间ms )165121024fff9f10数值154.8

16、43017.7149786.9547f,方差841.0507结论：通过对比不同程度噪音干扰下的反演结果可知，NNLS算法在信噪比较高的情况 下计算比较准确，但在测量误差比较大的情况下，反演得到的结果与实际信号差距很大， 已经不能正确计算出较好的结果。总结：根据上述分析可知，要减小测量误差带来的不良影响，可从两方面下手。一方面 在进行驰豫信号反演前，对测量数据进行去噪处理。另一方面，在去噪处理后，估算去 噪信号的SNR，如果信号SNR较高，进行驰豫信号反演时采用NNLS算法或者差分进化 算法均可，如果信号SNR较低，则应采用差分进化算法进行驰豫信号反演。通过上述两 种方法，可以减小测量误差为驰豫

17、信号反演结果带来的不良影响。六、模型的评价与推广6.1模型的评价6.1.1模型的优点本文在传统的5反演算法的基础之上，提出了改进的SVD反演算法，更容易实 现非约束，计算速度快，优化项的提出，更好的解决了T组分离散且分布较宽的问题。本文采用的非负最小二乘法的特点是数字稳定性好、运算速度快、容易实现；基 于差分进化的优化算法不依赖于初值选取，计算稳定，适用程度好。本文将二者进行结 合，得到了误差较小的结果。本文采用多种方法对同一问题求解，再对不同的模型进行比较，综合得出在符合 题意下的最优解，使得模型更具有推广性。6.1.2模型的缺点非负最小二乘法在信噪比较小的情况下，结果偏差小，在信噪比较大的

18、情况下有 其局限性，本文采用差分进化算法对其进行优化，得到了较好的结果。S0算法适合信噪比较高（SNR 80）的数据的反演，当数据信噪比较低时，反演 结果的分辨率比较低，有可能造成勺谱的畸形，解会出现不规则的跳动，并且在实现非 负约束时，可能出现谱线不连续的问题。6.2模型的推广本文建立的模型可用于解决多指数反演问题，采用基于DE的优化算法不依赖于初 值的选取，计算稳定，本文创新性的将非负最小二乘法和DE的优化算法相结合，不仅 对核磁共振勺谱反演问题有效，还可应用于其他非指数多项式非线性拟合问题中，也可 用于非负约束条件的优化问题，如：神经网络优化、阵列天线方向图综合等方面。参考文献王鹤，李鲠

19、颖。反演与拟合相结合处理核磁共振弛豫数据的方法。物理学报，第54 卷第3期2005年3月。王为民，李培，叶朝辉。核磁共振弛豫信号的多指数反演。中国科学院武汉物理与 数学研究所，2001年8月。王鹤，李鲠颖。非负最小二乘法在弛豫时间谱中的应用。中国新疆乌鲁木齐，2004。4附录附件1:运用Matlab软件随机取点：对数随机取点及求解2的幕指数随机取点及求解clcclear allload data.txtt=data(:,1);y=data(:,2);h,l=size(t);m=30;d=logspace(log10(0.0001),log10(10),m);A=;for i=1:hfor j=

20、1:mA(i,j)=exp(-t(i)/d)；endendf=A*(A*A)A(-1)*y;plot(t,f,k)xlabel(驰豫时间/s);ylabel(幅度 f);grid onfigureplot(1:30,d)xlabel(取点 /个);ylabel(T2 时间/s);clcclear allload data.txtt=data(:,1);y=data(:,2);h,l=size(t);m=30;d=linspace(log2(0.0001),log2(10),m);A=;for i=1:hfor j=1:mA(i,j)=exp(-t(i)/d(j);endendf=A*(A*A)

21、A(-1)*y;plot(t,f,k)xlabel(驰豫时间/s);ylabel(幅度 f);grid onfigureplot(1:30,d)xlabel (幕指数取点/个);ylabel(T2 时间/s);线性随机取点及求解clcclear allload data.txtt=data(:,1);y=data(:,2);h,l=size(t);m=30;d=linspace(0.0001,10,m);A=;for i=1:hfor j=1:mA(i,j)=exp(-t(i)/d(j);endendb=y;f=A*(A*A)A(-1)*y; plot(t,f,k)xlabel(驰豫时间/s); ylabel(幅度 f);grid on figure plo