李鵬飛 佟喜峰 游博洋

（東北石油大學(xué)計(jì)算機(jī)與信息技術(shù)學(xué)院 大慶 163318）

1 引言

在地球物理、遙感技術(shù)、工業(yè)控制以及信號(hào)（圖像）處理等眾多科學(xué)技術(shù)領(lǐng)域中，一般都存在反問(wèn)題的求解問(wèn)題。這些反問(wèn)題對(duì)應(yīng)的數(shù)學(xué)模型很多都是病態(tài)線性方程組。

在求解病態(tài)線性方程組時(shí)，其數(shù)值解是不穩(wěn)定的。由于系數(shù)矩陣和右端常數(shù)項(xiàng)量大多來(lái)源于原始觀測(cè)資料，而原始資料中存在的觀測(cè)誤差會(huì)因系數(shù)矩陣的病態(tài)性被放大，導(dǎo)致算法的數(shù)值解偏離真解而失去意義。以Ax=b方程組為例說(shuō)明病態(tài)系數(shù)矩陣對(duì)解的誤差的放大作用。在其系數(shù)矩陣和右端向量存在觀測(cè)誤差時(shí)，則變?yōu)槿缦滦问剑?/p>

式（1）有如下估計(jì)式成立［1］：

其系數(shù)矩陣的譜范數(shù)意義下的條件數(shù)為2503，屬于比較嚴(yán)重的病態(tài)程度。精確解為x1=-100，x2=-200。假設(shè)系數(shù)矩陣存在擾動(dòng)，則方程組變?yōu)?/p>

則得到一個(gè)截然不同的解：x1=40000，x2=79800。

當(dāng)病態(tài)系數(shù)矩陣的條件數(shù)很大時(shí)，在不對(duì)系數(shù)矩陣做出預(yù)先處理的情況下，就直接使用常規(guī)的求解方法得到的數(shù)值解誤差會(huì)非常大甚至無(wú)法獲得數(shù)值解。所以，目前對(duì)于病態(tài)線性方程組求解一般都圍繞先對(duì)病態(tài)系數(shù)矩陣進(jìn)行改良這種方法開(kāi)展工作，即先對(duì)病態(tài)系數(shù)矩陣進(jìn)行預(yù)處理，將其改造成良態(tài)矩陣后，再進(jìn)行求解。文獻(xiàn)［2］中SSOR-ICCG算法，其迭代過(guò)程不影響ICCG算法的收斂性，而且通過(guò)SSOR預(yù)處理來(lái)降低系數(shù)矩陣的條件數(shù)，從而提高算法的穩(wěn)定性，加快收斂速度。該算法要求系數(shù)矩陣必須對(duì)稱正定。文獻(xiàn)［3～4］中均提到正則化方法求解病態(tài)的反問(wèn)題，方法是對(duì)病態(tài)矩陣進(jìn)行正則化改造，使之變?yōu)榱紤B(tài)，然后用這個(gè)正則解作為待求方程的近似解。這類(lèi)方法中如何確定正則化參數(shù)，使正則解更好地逼近原問(wèn)題的解比較困難。根據(jù)原始觀測(cè)數(shù)據(jù)中誤差水平是否已知，有后驗(yàn)策略和先驗(yàn)策略兩種方式來(lái)確定正則參數(shù)。黃小為［5］認(rèn)為譜截?cái)嗟恼齽t化方法甚至比Tikhonov正則化方法更為有效，因此基于截?cái)嗥娈愔捣纸夥ǎ?～10］被越來(lái)越多的科研工作者應(yīng)用在各自領(lǐng)域的反問(wèn)題求解過(guò)程中。此外，人工神經(jīng)網(wǎng)絡(luò)和基于進(jìn)化計(jì)算一類(lèi)的智能優(yōu)化算法［11～15］也被應(yīng)用于病態(tài)的反問(wèn)題求解。

本文將病態(tài)系數(shù)矩陣進(jìn)行改造，在降低其病態(tài)性的同時(shí)，構(gòu)造了一種簡(jiǎn)單的預(yù)處理迭代算法，來(lái)求解嚴(yán)重病態(tài)的線性方程組。

2 預(yù)處理迭代算法設(shè)計(jì)

設(shè)有如下形式線性方程組：

其中，(aij)n×n>0 。，構(gòu)建對(duì)角陣D，并將DX疊加在式（3）兩端，得式（4）：

根據(jù)式（4）得：

要確保迭代公式（6）有效則必須滿足D+DA可逆，以及迭代矩陣(D+DA)-1(D-K)收斂?jī)蓚€(gè)條件。由于矩陣D和DA均為元素大于零的對(duì)角陣，所以D+DA必可逆。下面對(duì)迭代矩陣的收斂性予以證明。

3 迭代矩陣收斂性證明

結(jié)論1：設(shè)式（3）中系數(shù)矩陣A為實(shí)正定方陣，則式（6）中迭代矩陣(D+DA)-1(D-K)必收斂。

證明：設(shè)λ為系數(shù)矩陣A的特征值，λ'為迭代矩陣(D+DA)-1(D-K)的特征值，則對(duì)于這兩個(gè)矩陣有如下等式成立：

對(duì)于系數(shù)矩陣和迭代矩陣的任意一個(gè)特征值λi，，有：

根據(jù)所構(gòu)建的對(duì)角陣D可知，式（6）中迭代矩陣嚴(yán)格行對(duì)角元素占優(yōu)，且主對(duì)角線元素均大于0，依據(jù)文獻(xiàn)［16］中定理3的結(jié)論：矩陣中任一特征值實(shí)部必大于0。則迭代矩陣為正定矩陣，所以>0。又因?yàn)锳正定，所以λ>0。a>0，d>0iiii（已知），所以<1。由此得出：0<<1，因此迭代矩陣收斂。證畢。

4 預(yù)處理迭代算法對(duì)比實(shí)驗(yàn)測(cè)試

4.1 Hilbert病態(tài)線性方程組測(cè)試

以20階Hilbert病態(tài)線性方程組為測(cè)試實(shí)例，其中：

Hilbert矩陣元素：

右端向量：

在實(shí)際測(cè)試過(guò)程中，為了驗(yàn)證算法的穩(wěn)定性，取一組xj=1(j=1，2，…，20)作為真解代入方程組，求得右端向量后，在其基礎(chǔ)上添加SNR=50的高斯白噪聲。表1中對(duì)比了本文中設(shè)計(jì)的預(yù)處理迭代法、Matlab中自帶的非負(fù)最小二乘法優(yōu)化函數(shù)、遺傳算法和雅可比迭代法四種方法求解結(jié)果。預(yù)處理迭代算法的停止條件為當(dāng)0.001退出迭代。遺傳算法采用了Matlab中自帶的遺傳算法工具箱，其除了具有簡(jiǎn)單遺傳算法的選擇、交叉和變異算子之外，還有精英個(gè)體復(fù)制策略、遷移策略和懲罰函數(shù)等算子，用來(lái)提高數(shù)值解的精度。

表1 20階Hilbert病態(tài)線性方程組求解對(duì)比

4.2 核磁共振T2譜模型反演測(cè)試

根據(jù)核磁共振測(cè)井理論可知，氫原子核系統(tǒng)磁化強(qiáng)度矢量的橫向分量是按指數(shù)規(guī)律衰減的，由CPMG脈沖序列測(cè)得的回波串也按指數(shù)規(guī)律衰減。儲(chǔ)層巖石通常存在一個(gè)孔隙尺寸分布，并且常常含有多種流體成分，此時(shí)孔隙中存在多種弛豫組份，即橫向弛豫時(shí)間常數(shù)（T2）不是單值，而是一個(gè)分布，稱之為T(mén)2譜［17］。由CPMG脈沖序列測(cè)量記錄的自旋回波串按多指數(shù)規(guī)律衰減，即各單指數(shù)衰減的疊加?；夭ǚ扰c橫向弛豫分量初始幅度之間滿足式（10）關(guān)系［17～18］：

式中：bi為CPMG脈沖序列測(cè)量記錄的第i個(gè)回波幅度；Ej為第j個(gè)橫向弛豫分量的初始幅度（正數(shù)）；T2j為第j個(gè)橫向弛豫分量衰減的時(shí)間常數(shù)；TE為CPMG脈沖序列測(cè)量記錄的回波間隔；N為橫向弛豫分量的個(gè)數(shù)；M為回波個(gè)數(shù)。

T2譜反演就是根據(jù)預(yù)設(shè)的T2值和儀器測(cè)得的回波串?dāng)?shù)據(jù)，求解出Ej來(lái)擬合測(cè)得的回波串?dāng)?shù)據(jù)。這一過(guò)程稱為解譜或多指數(shù)擬合。其實(shí)質(zhì)是求解M個(gè)方程，N個(gè)未知數(shù)的線性方程組。

在實(shí)驗(yàn)室以實(shí)測(cè)的某井巖心回波數(shù)據(jù)為T(mén)2譜反演模型的右端向量，回波間隔表達(dá)式取TE=322.65+(i-1)×300，橫向弛豫分量的個(gè)數(shù)取N=128，回波個(gè)數(shù)取M=2048。某井巖心實(shí)驗(yàn)參數(shù)表如表2所示。

表2 某井巖心實(shí)驗(yàn)參數(shù)表

該井巖心數(shù)據(jù)構(gòu)成的模型為超定線性方程組，將其左乘系數(shù)矩陣的轉(zhuǎn)置矩陣，則原方程組變?yōu)槠湔齽t方程組。使用預(yù)處理迭代算法對(duì)其進(jìn)行反演解譜，并與實(shí)驗(yàn)室中英國(guó)牛津核磁共振分析儀（MARAN）自帶的解譜軟件包的解譜結(jié)果進(jìn)行了對(duì)比，如圖1所示。

圖1 預(yù)處理迭代法和MARAN解譜對(duì)比圖

5 結(jié)語(yǔ)

在20階Hilbert病態(tài)線性方程組的測(cè)試實(shí)例中，添加了SNR=50的高斯白噪聲。此時(shí)，系數(shù)矩陣數(shù)值是精確的，而右端向量相當(dāng)于有了擾動(dòng)δb。經(jīng)過(guò)Matlab中計(jì)算，其系數(shù)矩陣譜范數(shù)意義下的條件數(shù)為1.53×1018，已經(jīng)嚴(yán)重病態(tài)。觀察表1可以得出：預(yù)處理迭代法的相對(duì)誤差為0.007，精度還是相當(dāng)高的。遺傳算法的某些解偏離精確解比較大，導(dǎo)致相對(duì)誤差也比較大。而非負(fù)最小二乘法的數(shù)值解除已經(jīng)嚴(yán)重失真，失去了意義。對(duì)系數(shù)矩陣未經(jīng)處理而直接用雅可比迭代法進(jìn)行求解，算法失敗，這也印證了文中前面所提：如果不對(duì)病態(tài)矩陣進(jìn)行預(yù)處理而直接使用常規(guī)的解法進(jìn)行求解時(shí)通常會(huì)失敗。求從求解的速度來(lái)看，在雙核i3-2130，主頻3.4GHz，4GB內(nèi)存的臺(tái)式機(jī)上，預(yù)處理迭代法和非負(fù)最小二乘法求解時(shí)間均為0.01s的數(shù)量級(jí)，而遺傳算法（進(jìn)化代數(shù)2000代，停滯代數(shù)50）的求解速度最慢，約為10s數(shù)量級(jí)。

對(duì)于用實(shí)驗(yàn)室?guī)r心數(shù)據(jù)構(gòu)建的T2譜反演模型，回波間隔和回波幅度均為儀器測(cè)定，所以反演模型中系數(shù)矩陣和右端向量分別存在擾動(dòng)δA和δb。根據(jù)儀器測(cè)定的參數(shù)，按照式（10）計(jì)算，系數(shù)矩陣譜范數(shù)意義下的條件數(shù)為1.253×1019，病態(tài)性比測(cè)試實(shí)例1更重。將原方程組轉(zhuǎn)換為正則方程組后，系數(shù)矩陣病態(tài)性會(huì)進(jìn)一步加劇（譜范數(shù)意義下的條件數(shù)為7.113×1020）。此時(shí)，對(duì)于任何算法，求出穩(wěn)定的數(shù)值解難度都非常大。預(yù)處理迭代算法是在對(duì)病態(tài)的系數(shù)矩陣改造的基礎(chǔ)上而構(gòu)成一個(gè)迭代公式，并對(duì)其進(jìn)行求解，所以解的精度和穩(wěn)定性得到了有效的提高。通過(guò)觀察圖1可以得知預(yù)處理迭代算法反演結(jié)果與牛津核磁共振分析儀自帶的解譜軟件的反演結(jié)果基本吻合，譜的形狀與變化趨勢(shì)基本一致。因此，與比較昂貴的進(jìn)口儀器相比，采用預(yù)處理迭代算法進(jìn)行反演，具有更好的性價(jià)比。

以上兩個(gè)測(cè)試實(shí)例具有較好的代表性，既有理論數(shù)據(jù)構(gòu)造的系數(shù)矩陣，也有儀器實(shí)測(cè)的原始數(shù)據(jù)形成的具有擾動(dòng)的系數(shù)矩陣和右端向量，能夠更好地檢驗(yàn)算法的數(shù)值解穩(wěn)定性。從解譜結(jié)果對(duì)比來(lái)看，可以得知預(yù)處理迭代算法是求解嚴(yán)重病態(tài)線性方程組的一種有效算法。