李鵬飛 佟喜峰 李鵬舉

（1.東北石油大學(xué)計算機與信息技術(shù)學(xué)院 大慶 163318）（2.東北石油大學(xué)地球科學(xué)學(xué)院 大慶 163318）

1 引言

在地球物理反演問題中，大多數(shù)建立的數(shù)學(xué)模型都是超定病態(tài)線性方程組［1］，此時方程組為矛盾方程組，往往是無解的。即便方程有精確解，由于方程組的系數(shù)矩陣是病態(tài)的，直接解法和迭代解法等求解方法，所得到的數(shù)值解是不穩(wěn)定的［2］，往往因誤差太大而失去意義。

目前，針對病態(tài)線性方程組的求解方法可以分為如下幾類。截斷奇異值分解法［3～5］，正則化方法［6～7］，條件預(yù)優(yōu)法［8～9］，進化一類的算法，比如遺傳算法求解［10～11］。截斷奇異值分解法具有克服病態(tài)能力強，不放大誤差的特點，但是截斷參數(shù)確定較困難。正則化一類的方法是用一組與原不適定問題相“鄰近”的適定問題的解去逼近原問題的解，其正則參數(shù)的選取對問題解的性態(tài)起關(guān)鍵作用。參數(shù)值太小對系數(shù)矩陣的條件數(shù)改善不明顯，近似解的誤差很大。而參數(shù)值太大，新問題的解穩(wěn)定了，但是與原問題相差太大?？梢杂孟闰炘瓌t和后驗原則來確定正則參數(shù)，但不幸的是，這兩種方式確定正則參數(shù)都比較麻煩。條件預(yù)優(yōu)法通過構(gòu)造預(yù)優(yōu)矩陣左乘原方程組，來降低系數(shù)矩陣的條件數(shù)，這樣也是用一個新問題的解去近似原問題的解。但是這種方法尋找預(yù)優(yōu)矩陣比較復(fù)雜。遺傳算法的優(yōu)點是全局優(yōu)化，在多數(shù)情況下，遺傳算子的參數(shù)設(shè)定需要根據(jù)解的分布情況做動態(tài)改變，否則算法容易早熟而收斂到局部最優(yōu)解。

受迭代Tikhonov正則化思想啟發(fā)，本文設(shè)計一種簡單的正則化迭代算法，來求解超定病態(tài)線性方程組。

2 正則化迭代算法設(shè)計

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

其中：A∈RM×N。x∈RN。b∈RM。M＞N 用 AT左乘式（1），得到如下形式正則方程組：

其中：ATA∈RN×N為N階實對稱方陣。如果矩陣A為病態(tài)矩陣，此時ATA的病態(tài)性會加劇?；赥ik?honov正則化的思想，將式（2）系數(shù)矩陣 ATA的主對角線上元素疊加一個正則參數(shù)α，以此來降低ATA的條件數(shù)。同時為了保證與式（2）同解，可構(gòu)建如下公式：

其中：α為正則系數(shù)，I為與ATA同階的單位陣。由于式（3）等式左右兩端分別含有解向量x，在式（3）兩端分別乘(ATA+αI)-1，則可以構(gòu)成一個一般迭代式（4）：

其中：K=ATA。C=α(K+αI)-1。f=(K+αI)-1ATb

要想確定式（4）對應(yīng)算法是否有效必須滿足兩個條件，一是K+αI是否可逆，二是迭代矩陣c是否收斂。下面通過對正則參數(shù)α取值范圍進行判定，進而回答式（4）對應(yīng)算法何時有效。

3 正則參數(shù)取值范圍判定

結(jié)論1：有實方陣K，如果K的特征值為λ，則矩陣K+αI的特征值為λ+α。

證明：設(shè) λi為矩陣 K=(kij)n×n的第i個特征值為矩陣K+αI的第i個特征值，則對于矩陣K+αI有如下等式成立：

對于矩陣K則有有如下等式成立：

比較上式得：

結(jié)論2：如果K為對稱正定方陣，當α＞0時，式（4）所構(gòu)成的迭代公式必收斂。

證明：根據(jù)式（4）知其迭代矩陣為

因為K為對稱正定方陣，所以K的全部特征值均大于0。根據(jù)結(jié)論1和α＞0（已知）得：矩陣K+αI的全部特征值也均大于0，所以矩陣K+αI必可逆。設(shè)K的任一特征值為λ，則存在非零向量x，使得如下等式成立：

上式兩端同乘以α，則：值均小于1，其譜半徑也小于1，所以式（4）所構(gòu)成的迭代公式必收斂。

證畢。

因此，對于對稱正定矩陣選取正則參數(shù)α＞0，就可以保證矩陣 K+αI可逆，且迭代矩陣α(K+αI)-1收斂。

為了有效地降低系數(shù)矩陣的條件數(shù)，提高數(shù)值解的精度，只有α＞0的結(jié)論是不夠的。對于對稱陣為嚴重病態(tài)時，λmax＞＞λmin。此時取 α=λmax，根此時可將嚴重病態(tài)的系數(shù)矩陣改為良態(tài)矩陣，從而提高了數(shù)值解的精度。

4 核磁共振弛豫反演模型

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

式中：bi為CPMG脈沖序列測量記錄的第i個回波幅度；Ej為第j個橫向弛豫分量的初始幅度（正數(shù)）；T2j為第j個橫向弛豫分量衰減的時間常數(shù)；TE為CPMG脈沖序列測量記錄的回波間隔；N為橫向弛豫分量的個數(shù)；M為回波個數(shù)；T2譜反演就是根據(jù)預(yù)設(shè)的T2值和儀器測得的回波串數(shù)據(jù)，求解出Ei，使其滿足式（5）。其實質(zhì)是求解M個方程，N個未知數(shù)的線性方程組。

5 有噪聲理論數(shù)據(jù)反演

考慮如下核磁共振測井T2譜反演模型：

其中：所有算法均在Matlab上實現(xiàn)。

圖1 不同信噪比正則化迭代算法解譜對比圖

圖2 三種不同算法解譜對比圖（信噪比=25%）

表1 不同信噪比和算法誤差表

6 巖心實驗數(shù)據(jù)反演

在實驗室以實測的某井巖心回波數(shù)據(jù)為T2譜反演模型的右端向量，回波間隔取TE=322.65+(i-1)×300，i的意義同式（6）。使用正則化迭代算法對其進行反演解譜，并與實驗室某國外反演軟件包的解譜結(jié)果進行了對比，如圖3所示。

圖3 某井巖心反演結(jié)果與實驗室反演結(jié)果對比圖

7 結(jié)語

對于有噪聲的理論數(shù)據(jù)反演，對比圖1和表1，可知在α一定的情況下，信噪比snr越大，數(shù)值解的精度越高，這個道理是很淺顯的，因為信噪比snr越大，噪聲對病態(tài)矩陣的影響越小，解的精度越高。對比圖2和表1，可以得出，正則化迭代算法數(shù)值解的精度明顯高于簡單遺傳算法和最小二乘法求解精度。這是因為對于病態(tài)線性方程組，如果不對系數(shù)矩陣改造，使之變?yōu)榱紤B(tài)，在求解時右端向量的噪聲或系數(shù)矩陣的誤差會因系數(shù)矩陣病態(tài)性而被放大，從而導(dǎo)致數(shù)值解的誤差非常大。正則化迭代算法是在對病態(tài)的系數(shù)矩陣改造成良態(tài)的基礎(chǔ)上構(gòu)成的迭代公式進行求解，所以解的精度得到了有效的提高。對于某井巖心實驗數(shù)據(jù)反演結(jié)果，通過觀察可以得知正則化迭代算法反演結(jié)果與實驗室反演結(jié)果基本吻合，譜的形狀與變化趨勢基本一致。因此，正則化迭代算法是求解大型超定病態(tài)線性方程組的一種有效算法。