求解超定病態線性方程組的一種正則化迭代算法

李鵬飛 佟喜峰 李鵬舉

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

1 引言

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

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

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

2 正則化迭代算法設計

設有如下形式的超定線性方程組：

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

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

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

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

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

3 正則參數取值范圍判定

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

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

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

比較上式得：

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

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

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

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

證畢。

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

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

4 核磁共振弛豫反演模型

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

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

5 有噪聲理論數據反演

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

其中：所有算法均在Matlab上實現。

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

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

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

6 巖心實驗數據反演

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

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

7 結語

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