第三次模擬賽

1、2012年四川理工學(xué)院暑期培訓(xùn)第三次模擬賽承諾書我們仔細(xì)閱讀了四川理工學(xué)院大學(xué)生數(shù)學(xué)建模競(jìng)賽的競(jìng)賽規(guī)則.我們完全明白，在競(jìng)賽開始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、網(wǎng)上 咨詢等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競(jìng)賽規(guī)則的，如果引用別人的成果或其他公開的 資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考 文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競(jìng)賽規(guī)則，以保證競(jìng)賽的公正、公平性。如有違反競(jìng)賽規(guī) 則的行為，我們將受到嚴(yán)肅處理。我們參賽選擇的題號(hào)是（從A/B/C/D中選擇一項(xiàng)填寫）：A我們的參賽隊(duì)號(hào)為所屬學(xué)

2、校（請(qǐng)?zhí)顚懲暾娜?四川理工學(xué)院參賽隊(duì)員（打印并簽名）:1.郭姝霖徐高飛胡婷婷日期：2012 年 08月 29 日評(píng)閱編號(hào)（由評(píng)委團(tuán)評(píng)閱前進(jìn)行編號(hào)）:2012年四川理工學(xué)院暑期培訓(xùn)第三次模擬賽編號(hào)專用頁(yè)評(píng)閱編號(hào)（由評(píng)委團(tuán)評(píng)閱前進(jìn)行編號(hào)）:核磁共振弛豫信號(hào)反演問題摘要本文針對(duì)核磁共振弛豫信號(hào)反演問題，在S0反演算法的基礎(chǔ)之上，創(chuàng)新性地提出 了變換反演算法，解決了當(dāng)數(shù)據(jù)信噪比較低時(shí)，反演結(jié)果的分辨率比較低，有可能造成 T譜的畸形的情形，建立了非負(fù)最小二乘法和DE的優(yōu)化算法相結(jié)合的模型，很好的解 決了問題。針對(duì)問題一，本文首先選用對(duì)數(shù)、2的幕指數(shù)、線性均勻布點(diǎn)三種方式進(jìn)行布點(diǎn)， 對(duì)T譜建立數(shù)學(xué)

3、模型，再將三種模型相互比較，在此基礎(chǔ)之上，對(duì)模型進(jìn)行修改，得出 ，2 ，一 了很好的結(jié)果，針對(duì)問題二，要求在弛豫時(shí)間未知的情況下對(duì)T譜進(jìn)行求解，2針對(duì)問題三，誤差的主要來源是測(cè)量誤差，而會(huì)造成測(cè)量數(shù)據(jù)誤差的有等待時(shí)間、 回波間隔、輕烴或稠油、順磁物質(zhì)以及采集方法等因素。而這些都是在測(cè)井的時(shí)候就已 經(jīng)存在，無法避免，只能通過改進(jìn)現(xiàn)有的計(jì)算模型以及相關(guān)的算法來使得結(jié)果更加精確。 當(dāng)原始數(shù)據(jù)的信噪比SNR 80時(shí)，可選用SVD反演算法，然而為保證T反演結(jié)果的真 實(shí)性，本文選用30個(gè)T的布點(diǎn)數(shù)來建立改進(jìn)的SVD反演算法，最后將非負(fù)最小二乘法 與差分進(jìn)化的優(yōu)化算法相結(jié)合，很好地排除了不同條件下造成的誤差

4、。本文建立了反演 優(yōu)化模型，并通過差分進(jìn)化算法求解，最終得出在建模和算法上合理的更加精確的解決 方法。最后，本文對(duì)結(jié)果進(jìn)行了分析，對(duì)模型的優(yōu)缺點(diǎn)進(jìn)行了評(píng)價(jià)，并對(duì)模型進(jìn)行了推廣和改 進(jìn)。關(guān)鍵詞：核磁共振SVD反演算法非負(fù)最小二乘法一、問題重述1.1背景核磁共振（NMR）在測(cè)井技術(shù)和巖心分析中已經(jīng)得到廣泛應(yīng)用，正在為油氣資源的勘 探和開發(fā)發(fā)揮重要作用。油氣井中NMR技術(shù)提供的油氣藏流體特性和儲(chǔ)集層參數(shù)，如地 層孔隙度、孔徑分布、束縛水與可動(dòng)流體孔隙體積、滲透率、以及流體的擴(kuò)散系數(shù)和黏 度等信息，都需要經(jīng)過一個(gè)基本的反演處理，即NMR弛豫信號(hào)的多指數(shù)擬合，得到弛豫 時(shí)間的分布。在核磁共振測(cè)井中，一般

5、采用CMPG方法測(cè)量自旋回波串，縱向弛豫時(shí)間廣和橫向 弛豫時(shí)間T都是用來描述自旋回波信號(hào)的弛豫特征。由于T譜能夠提供許多巖石物性和 流體特性的信息，越來越受到人們的關(guān)注。根據(jù)核磁共振理論，單個(gè)空隙中的磁化強(qiáng)度 信號(hào)的衰減滿足單指數(shù)衰減規(guī)律，但由于巖石內(nèi)部是一系列大小不等的孔隙群體組成， 所以在巖石核磁共振中測(cè)得的中磁化強(qiáng)度信號(hào)yt）為一系列單個(gè)孔隙磁化強(qiáng)度信號(hào)的 疊加，因此y（t）可描述為：yG）=Z f - e-T；j, t = n -tj=1 TOC o 1-5 h z 其中f為第j類孔隙在總孔隙中所占的份額，T為第j類孔隙的T弛豫時(shí)間，通常范 j一2 j 一2圍為0.1ms T 為測(cè)量的

6、磁化強(qiáng)度衰減信號(hào)，A =A = e*，1 2 nij nxm |_-nxmf = (f, f，f)為弛豫時(shí)間T對(duì)應(yīng)各點(diǎn)的幅度值。T (j = 1,2,.,m)為預(yù)先指定的T時(shí)1 2 m2 j2 j2間分布系列。本文采用在T . , T )區(qū)間內(nèi)對(duì)數(shù)均勻地選取m個(gè)點(diǎn)，即弛豫時(shí)間布點(diǎn)，同 時(shí)也采用2的幕指數(shù)，線性均勻布點(diǎn)的三種方式。為使算法簡(jiǎn)便化，本文采用變換的反演算法，假設(shè)每個(gè)點(diǎn)的測(cè)量誤差是隨機(jī)獨(dú)立且 標(biāo)準(zhǔn)差都一樣，則可給定如下的目標(biāo)函數(shù)：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)運(yùn)用Matlab軟件求解得C的值，再經(jīng)過線性變換式回代而獲得，選擇合適的人， 以保證矩陣J - At +人L)的可逆，則求解得方程(1)的最小二乘解為：f = At - Ca - At +人I )】 j則分別可得如下如圖4，圖5，圖6結(jié)果，其具體程序見附件1：由圖重結(jié)果可以看出，f的值有正有負(fù)，這在物理意義上來說不符合邏輯，故本文 i對(duì)模型改進(jìn)如下：T譜非負(fù)限制的實(shí)現(xiàn)，為了克服傳統(tǒng)的缺點(diǎn)，不破壞T分布的連續(xù)性，防止T譜的.一.222畸形，本文米用迭代法逐漸消去負(fù)數(shù)的組分，其思想為：通過T域

8、到時(shí)域的變換后，將 第k項(xiàng)相鄰點(diǎn)的T分布信息疊加到負(fù)數(shù)負(fù)數(shù)的T組分中，重新進(jìn)行反演計(jì)算，直到滿22 k足所有的f 0（j = 1,2質(zhì)）為止。具體結(jié)果見下圖：i由模型一的反演結(jié)果分析可知反演的分辨率低的問題，是因?yàn)門在運(yùn)算前就已經(jīng)4.1問題分析2 j選定，所以本模型中將T也看成未知數(shù)。并采用將T和s成為了要求解的未知數(shù)，將2 j2 j jT2和s.間的關(guān)系由線性轉(zhuǎn)為非線性關(guān)系求解。為使求解過程簡(jiǎn)單，本文采用非線性擬 合求解；求解過程采用L M （高斯一牛頓）法。非線性擬合過程為：實(shí)際問題中，匕是非負(fù)的，于是采用非負(fù)最小二乘法進(jìn)行擬合，首先設(shè)置C的初值，最后對(duì)初值進(jìn)行修正，直到X 2相對(duì)于前一次

9、不再增長(zhǎng)后者不再明 顯增長(zhǎng)時(shí)為止。針對(duì)T譜中每一個(gè)峰對(duì)應(yīng)一個(gè)種類孔隙的特征，T譜中峰的個(gè)數(shù)就等于 弛豫時(shí)間T.的種類數(shù)。對(duì)于k的選取，本文在多次擬合的情況下，選用30個(gè)2的布 點(diǎn)數(shù)來建立改進(jìn)的SVD反演算法，將非負(fù)最小二乘法與差分進(jìn)化的優(yōu)化算法相結(jié)合，很 好的排除了不同條件下造成的誤差。4.2模型的建立及求解4.2.1模型的建立進(jìn)行如下改進(jìn)得: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.在運(yùn)算前

10、由Matlab軟件隨機(jī)選定，而反演算法分辨率很低， 為提高本文的分辨率，本文對(duì)模型：其中，k = m / 2，的關(guān)系變?yōu)閥與C之間的非線性關(guān)系，于是本文采用線性擬合求解。i i 的函數(shù)值在點(diǎn)G = 1,2, ., m設(shè)有實(shí)驗(yàn)數(shù)據(jù)(t, y )(i = 1,2,m)，尋找函數(shù)f (t, y)，使得函數(shù)在點(diǎn)x (i = 1,2,m)夕卜 )處的函數(shù)值與觀測(cè)數(shù)據(jù)偏差的平方和達(dá)到最小，即：min (f (t, y )- y= 代iii其中， r Q=f (t, y)- y.。 可得系數(shù)的計(jì)算公式為：i=1i=1 (y - yiR 2 = 1 -i=1 (y- ii=1其中， = 1 lLyi，而R 2

11、越趨于1表明擬合效果越好。.i=1本文采用基于改進(jìn)的L - M (高斯一牛頓法)算法對(duì)其進(jìn)行求解，將其轉(zhuǎn)化為求解誤 差平方和最小。即有：min f (x )= 1 代=12 i 2記：r(x)= r (x), r()., r 1)1, A(x)= (Vr ., Vr (x)。12m1m本文采用改進(jìn)的L-M算法，其形式為：s = I/t Co)j (w)+ pTL J JU)其中，比例系數(shù)p 0為常數(shù)，I是單位矩陣。編程求解的結(jié)果如下:五、問題三5.1問題分析針對(duì)問題三，由于測(cè)量方式和測(cè)量工具的差異，測(cè)量數(shù)據(jù)必然會(huì)存在誤差。前兩問 中考慮到測(cè)量數(shù)據(jù)的誤差，均通過分析去噪信號(hào)，求解測(cè)量信號(hào)的尸頻譜

12、。在實(shí)際檢測(cè) 計(jì)算過程中，需要考慮測(cè)量誤差對(duì)結(jié)果的影響。因此在數(shù)據(jù)處理時(shí)應(yīng)盡可能的減小測(cè)量 誤差的程度。本文將測(cè)量誤差視為測(cè)量過程中各種干擾導(dǎo)致的結(jié)果。5.2模型的建立及求解信噪比SNR :從測(cè)量數(shù)據(jù)中估計(jì)出的信號(hào)強(qiáng)度與噪音之比，定義為第一個(gè)回波的幅 度值除以誤差矢量的標(biāo)準(zhǔn)差。為了衡量測(cè)量誤差的程度，此處引入信噪比（SNR），信噪比（SNR）按如下方式定義： SNR = 20% 產(chǎn)|2其中y表示有效信號(hào)，即去噪信號(hào)，表示原信號(hào)，即未經(jīng)去噪處理的信號(hào)。SNR 數(shù)值越大，表明信號(hào)受干擾越少，即測(cè)量誤差越小。本文采用計(jì)算機(jī)模擬來分析測(cè)量誤差對(duì)于前面兩個(gè)問題帶來的影響。首先給定一個(gè) 具有雙峰特征的T譜

13、分布y（t）= 200/16 + 800/1024 （單位：ms ），由頻譜可計(jì)算出不同 時(shí)刻匕，計(jì)算出不同信號(hào)強(qiáng)度y,再對(duì)y加上一定大小的噪音得到，由數(shù)據(jù)可計(jì)算出 受干擾信號(hào)的信噪比，模擬信號(hào)的噪聲由Matlab軟件隨機(jī)生成。對(duì)于模擬所得的信號(hào)， 在預(yù)先給定弛豫時(shí)間分布T.和未預(yù)先給定弛豫時(shí)間分布兩種情況下，分別采用非負(fù)最 小二乘算法和差分進(jìn)化算法反演頻譜。最終通過對(duì)比得出測(cè)量誤差對(duì)于這兩種情況的頻 譜計(jì)算結(jié)果的影響。由問題一的結(jié)果分析可得，采用對(duì)數(shù)均勻布點(diǎn)方式曲線擬合相似度較高，因此本問 采用對(duì)數(shù)均勻布點(diǎn)方式，取L = 2 j（j = 1,2, ,15）（單位：ms ）。分別在信號(hào) y（t

14、）= 200/16 + 800/1024上加以不同程度的噪聲，計(jì)算模擬信號(hào)的信噪比SNR，在利 用問題一中的NNLS算法對(duì)模擬信號(hào)進(jìn)行反演，得出反演曲線的殘差和對(duì)應(yīng)的f與原信 號(hào)f.的方差作為評(píng)判反演結(jié)果好壞的指標(biāo)。j（1）加為強(qiáng)度為0-1的隨機(jī)噪音利用Matlab計(jì)算結(jié)果如下表1：表1強(qiáng)度為0-1隨機(jī)噪音信號(hào)反演結(jié)果信噪比112.2516殘差34.5041TT 22T 23T 24T 25T 26T 210T211） 時(shí)間（ms ）24163225610242048f.f2ffff6f10f數(shù)值12.78644.6748185.03178.47050.6549794.40954.5702J

15、.方差76.2473（2）加為強(qiáng)度為0-5的隨機(jī)噪音利用Matlab計(jì)算結(jié)果如下表2：表2強(qiáng)度為0-5隨機(jī)噪音信號(hào)反演結(jié)果信噪比86.3802殘差561.2619TTTTTT孕 23242529210時(shí)間ms )416325121024ffffffj345910數(shù)值5709524127.902226.7878224.7314434.0134f j方差38744.9(1)加為強(qiáng)度為0-20的隨機(jī)噪音 利用Matlab計(jì)算結(jié)果如下表3：表3強(qiáng)度為0-20隨機(jī)噪音信號(hào)反演結(jié)果信噪比64.6336殘差8520.4T2jT 24T 29T 210時(shí)間ms )165121024fff9f10數(shù)值154.8

16、43017.7149786.9547f,方差841.0507結(jié)論：通過對(duì)比不同程度噪音干擾下的反演結(jié)果可知，NNLS算法在信噪比較高的情況 下計(jì)算比較準(zhǔn)確，但在測(cè)量誤差比較大的情況下，反演得到的結(jié)果與實(shí)際信號(hào)差距很大， 已經(jīng)不能正確計(jì)算出較好的結(jié)果?？偨Y(jié)：根據(jù)上述分析可知，要減小測(cè)量誤差帶來的不良影響，可從兩方面下手。一方面 在進(jìn)行馳豫信號(hào)反演前，對(duì)測(cè)量數(shù)據(jù)進(jìn)行去噪處理。另一方面，在去噪處理后，估算去 噪信號(hào)的SNR，如果信號(hào)SNR較高，進(jìn)行馳豫信號(hào)反演時(shí)采用NNLS算法或者差分進(jìn)化 算法均可，如果信號(hào)SNR較低，則應(yīng)采用差分進(jìn)化算法進(jìn)行馳豫信號(hào)反演。通過上述兩 種方法，可以減小測(cè)量誤差為馳豫

17、信號(hào)反演結(jié)果帶來的不良影響。六、模型的評(píng)價(jià)與推廣6.1模型的評(píng)價(jià)6.1.1模型的優(yōu)點(diǎn)本文在傳統(tǒng)的5反演算法的基礎(chǔ)之上，提出了改進(jìn)的SVD反演算法，更容易實(shí) 現(xiàn)非約束，計(jì)算速度快，優(yōu)化項(xiàng)的提出，更好的解決了T組分離散且分布較寬的問題。本文采用的非負(fù)最小二乘法的特點(diǎn)是數(shù)字穩(wěn)定性好、運(yùn)算速度快、容易實(shí)現(xiàn)；基 于差分進(jìn)化的優(yōu)化算法不依賴于初值選取，計(jì)算穩(wěn)定，適用程度好。本文將二者進(jìn)行結(jié) 合，得到了誤差較小的結(jié)果。本文采用多種方法對(duì)同一問題求解，再對(duì)不同的模型進(jìn)行比較，綜合得出在符合 題意下的最優(yōu)解，使得模型更具有推廣性。6.1.2模型的缺點(diǎn)非負(fù)最小二乘法在信噪比較小的情況下，結(jié)果偏差小，在信噪比較大的

18、情況下有 其局限性，本文采用差分進(jìn)化算法對(duì)其進(jìn)行優(yōu)化，得到了較好的結(jié)果。S0算法適合信噪比較高（SNR 80）的數(shù)據(jù)的反演，當(dāng)數(shù)據(jù)信噪比較低時(shí)，反演 結(jié)果的分辨率比較低，有可能造成勺譜的畸形，解會(huì)出現(xiàn)不規(guī)則的跳動(dòng)，并且在實(shí)現(xiàn)非 負(fù)約束時(shí)，可能出現(xiàn)譜線不連續(xù)的問題。6.2模型的推廣本文建立的模型可用于解決多指數(shù)反演問題，采用基于DE的優(yōu)化算法不依賴于初 值的選取，計(jì)算穩(wěn)定，本文創(chuàng)新性的將非負(fù)最小二乘法和DE的優(yōu)化算法相結(jié)合，不僅 對(duì)核磁共振勺譜反演問題有效，還可應(yīng)用于其他非指數(shù)多項(xiàng)式非線性擬合問題中，也可 用于非負(fù)約束條件的優(yōu)化問題，如：神經(jīng)網(wǎng)絡(luò)優(yōu)化、陣列天線方向圖綜合等方面。參考文獻(xiàn)王鶴，李鯁

19、穎。反演與擬合相結(jié)合處理核磁共振弛豫數(shù)據(jù)的方法。物理學(xué)報(bào)，第54 卷第3期2005年3月。王為民，李培，葉朝輝。核磁共振弛豫信號(hào)的多指數(shù)反演。中國(guó)科學(xué)院武漢物理與 數(shù)學(xué)研究所，2001年8月。王鶴，李鯁穎。非負(fù)最小二乘法在弛豫時(shí)間譜中的應(yīng)用。中國(guó)新疆烏魯木齊，2004。4附錄附件1:運(yùn)用Matlab軟件隨機(jī)取點(diǎn)：對(duì)數(shù)隨機(jī)取點(diǎn)及求解2的幕指數(shù)隨機(jī)取點(diǎn)及求解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(馳豫時(shí)間/s);ylabel(幅度 f);grid onfigureplot(1:30,d)xlabel(取點(diǎn) /個(gè));ylabel(T2 時(shí)間/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(馳豫時(shí)間/s);ylabel(幅度 f);grid onfigureplot(1:30,d)xlabel (幕指數(shù)取點(diǎn)/個(gè));ylabel(T2 時(shí)間/s);線性隨機(jī)取點(diǎn)及求解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(馳豫時(shí)間/s); ylabel(幅度 f);grid on figure plo