線性方程組迭代解法在重力學中的應用

1、線性方程組迭代解法在重力學中的應用肖* 200*9 *學院地球物理學摘要：在重力學中，許多利用地表觀測值求內部質量分布和斷層產狀的反演問題通常存在積分解算。這些積分一般都非常復雜，很難用解析法直接運算，但是我們可以通過線性化來求其近似解。本文討論的就是重力學中的積分問題通過線性方法解算的內容。關鍵詞：線性代數(shù)；迭代解法；重力學一、 重力學中的線性問題1、 重力反演問題在地球重力學中，反演問題是指，由實測重力異常及其倒數(shù)的數(shù)值大小、空間分布和變化規(guī)律，定性和定量的推斷地下客觀存在的異常地質結構、構造和地質體的形狀、產狀及剩余密度分布。反演問題的流程可以描繪為：觀測數(shù)據(jù)反演數(shù)學物理模型場源模型參數(shù)

2、一般地，場源函數(shù)與地球表面物理場之間存在著積分關系： （）其中，是場源的作用區(qū)域，稱為積分核函數(shù)（又稱格林函數(shù)），它刻畫“場”與“源”之間的“關系”。若函數(shù)為已知函數(shù)，積分核函數(shù)根據(jù)場與場源的物理實際給出，求解場源的物理問題就是反演問題。為了敘述方便，將寫成算子形式。若存在逆算子使得，則稱m為的（反演問題意義下的、廣義的）解。一般情況下，積分算子與微分算子互為逆算子。2、 重力線性問題的線性化求解積分方程的解析求解或精密求解時比較困難的事情，對此，人們往往采用迭代方法或線性近似方法。先考慮線性問題，設為一組歸一化的正交多項式基，則一維問題可以變?yōu)椋海?.2.1）系數(shù)待定，帶入積分得：(1.2.

3、2)系數(shù)為已知。為線性的，稱用該方法求解（也就是求解系數(shù)）的問題為線性問題。在上述方程中，若已知足夠多的觀測值，則可能球的所有的。在實際中，不可能使用無限多的觀測值去求定無限多個系數(shù)，實際中常常使用有線展開方法做的近似，即設關于基函數(shù)的選擇問題可以這樣處理，由于M個觀測數(shù)值處的坐標已知，則為一組函數(shù)，由該組函數(shù)構成基函數(shù)，記：為常系數(shù)，為正交基函數(shù)，若則為標準正交基函數(shù)。如果在方程中G為一個維向量，則方程演變?yōu)橐粋€簡單的矩陣方程，即:()方程（）也可以看作是方程（）離散化的結果，這樣，反演問題就變成了求解階線性方程組的問題。3、 重力反演存在的問題反演問題的主要內容有三方面，其一是解得適定問題

4、，包括借的存在性、唯一性及穩(wěn)定性。其二是反演問題的求解方法。其三是反演問題解的評價。在重力反演中，解的存在性已經被大量事實所證實。然而反演問題的解具有不唯一性，為說明解得非唯一性，讓我們來考慮零向量和零空間的概念。若線性反演問題有兩個解和，其解非唯一，則有：，兩個方程相減，得因為假定這兩個解是截然不同的，所以它們的差是非零的，而稱矢量為零矢量，而由零矢量組成的空間稱為零空間。由此看出，任何具有零矢量的線性反演問題的解都是非唯一的。如果（這里表示特解）是的一個非零解，則也是一個非零解。其中是一個不為零的任意常數(shù)。若給定線性問題有q個獨立的零矢量，則其一般解為: (1.3.1)總之，反演一組觀測數(shù)

5、據(jù)就是一個，而中的任意一個都可以加到上而仍然擬合觀測數(shù)據(jù)，因而使解變的非唯一。現(xiàn)在回到方程，且設，r為G的秩，顯然有以下四種情況：（1） 當M=N=r，且存在時，方程為適定問題，有唯一解。（2） 當M>N=r時，方程為超定方程，無常規(guī)意義下的解，但有最小方差解。（3） 當N>M=r時，方程是欠定方程，無常規(guī)意義下的解，但有意義下的解。（4） 當min（M,N）>r時，方程的問題是混定問題。這時，只有同時在兩種限制條件之下，方程才有解。下面將討論每一種情況的具體解法。（1） 適定方程對于適定方程，直接使用下文所介紹的線性方程迭代解法求解即可。（2） 超定方程對于超定方程，最簡單

6、、最常用的反演方法是最小方差法。(1.3.2)r為G的秩。設e是觀測數(shù)據(jù)g與理論計算值Gm之誤差向量，則方差（即目標函數(shù)）為 (1.3.3)最小方差解必須滿足：所以：(1.3.4)注：當有零特征值存在時，方程是奇異的，無法求解；當?shù)奶卣髦岛苄r，方程是病態(tài)的，會使解變得極不穩(wěn)定，使反演非常困難，甚至無法收斂。這也是在重力資料反演中經常遇到的問題。（3） 欠定方程從線性代數(shù)理論可以知道，此時有無限多個解能滿足方程，且其誤差均為零。這是因為，雖然數(shù)據(jù)提供了一些確定參數(shù)模型的信息，但其數(shù)量不足以全部確定模型參數(shù)，或為題夠確定模型參數(shù)足夠充分的信息。因此，解是不唯一的，甚至有無限多能擬合觀測數(shù)據(jù)的解。

7、為了求得重力反演問題的一個解，我們必須從無限多個能擬合觀測數(shù)據(jù)的解中，挑選一個我們所需要的特定解。因此，我們必須加上一些它未包含的信息，這種附加給反演問題的信息叫“先驗信息”。在補充先驗信息的時候，我們需要講究“擇缺補充”的原則。第一類是補充待求參數(shù)的物理性質和可能的數(shù)值范圍，如地球密度的非負性，且根據(jù)地球物理常識可以先定在0.016.0。第二類是先驗信息來源于其他一致的地質、重力資料。比如反演地球基底的埋深，油層的厚度，金屬礦的屬性等。第三類，某些參數(shù)比其他參數(shù)對解決重力反演問題更重要，此時可以對模型參數(shù)經行加權，在一定權系數(shù)約束下求解。第四類，假定地球物理模型最簡單。這里的最簡單是指在保留

8、實際地球物理模型基本特征不變的情況下，最地球物理模型的一種簡化。借的長度，比如說借的歐幾里得長度為最小的模型，應該是一種簡單的模型。當然，這里定義的簡單不一定處處非常合理。因此又出來了其他形式的最簡單模型，如該模型的變化為最小的范數(shù)意義下的模型等。設為一切頂問題，此時的目標函數(shù)，在的約束之下有極小，即根據(jù)極致理論，必須引入拉格朗日算子將條件極值問題化為無條件極值問題。因此目標函數(shù)應為：(1.3.5)顯然求上述目標函數(shù)的極小值可以化為求：故則 (1.3.6)注：在欠定問題求解中也存在這“奇異”“病態(tài)”的問題。此處奇異是指的特征值中有為零的問題；病態(tài)問題只是中有小特征值的問題。這都是在解欠定問題時

9、必須認真對待的。（4） 混定方程鑒于混定問題的復雜性，不難想象其目標函數(shù)中應該兼有方差項和模型長度項，即(1.3.7)求解可以得到： (1.3.8)此處，稱為阻尼系數(shù)或加權因子，它決定預測誤差項和模型范數(shù)長度項再寄笑話目標函數(shù)E時的相對重要性，如果阻尼系數(shù)足夠大，則問題以欠定問題為主，反之則以超定問題為主。二、 線性方程組的迭代解法概述科學研究中及生產實踐中，很多問題都可以歸結為線性方程組的求解，高效求解線性方程組成了許多科學計算及實踐的核心之一。解線性方程組的傳統(tǒng)方法是直接法，理論上可以得到方程組的真解，但是很容易出現(xiàn)大數(shù)“吃掉”小數(shù)、不穩(wěn)定的問題。線性方程組的迭代解法是一種間接解法，它利用

10、某種迭代過程去逐步逼近準確解，而求出方程組具有指定精確度的近似解的方法。它可以有效地提高求解的速度。迭代法的另一個突出特點是可以充分利用和保持系數(shù)矩陣的稀疏性，只需知道系數(shù)矩陣與向量乘積的計算法則，而不必知道具體的系數(shù)矩陣，這樣可以節(jié)省內存的開銷?？梢?，用迭代法求解線性方程組，無論在時間上和空間上都有很大的優(yōu)勢。以下，介紹五種常用的迭代方法。1、 雅可比(Jacobi)迭代法對于n階方程組中的每一個方程（）可以改寫為未知向量x的分量的形式：(2.1.1)設方程的準確解為向量，任取一個向量作為的初始近似，將帶入（）式的右端，求出的結果記為，稱為1次近似。一般的，在求出了第m次近似后，在帶入（）式

11、的右端便可得到m+1次近似(2.1.2)給定初值，雅可比迭代公式格式為：(2.1.3)定義1：如果向量序列中的向量每一個分量當時都趨于向量的對應分量，則稱是該向量的極限，記為在一定條件下，對任意初始向量，按迭代公式（）求出的向量序列的極限存在且等于方程組的解，這種用迭代格式（）求線性方程組近似解的方法成為雅可比迭代法。2、 高斯賽德爾（Gauss-Seidel）迭代法在雅可比迭代公式中，每算出的一個分量便將其代替該次迭代中以下分量的計算式中的相應分量，便得到另一種迭代公式：（）在給定初始向量，按迭代公式（）求出的向量序列的極限存在且等于方程組的解，這種用迭代格式（）求線性方程組近似解的方法稱為

12、高斯賽德爾迭代法。3、 超松弛迭代法（SOR方法）若第次分量迭代已經完成，記第次迭代的第個分量為；第次向量迭代的前個分量已經算出，由高斯賽德爾迭代法算出的第個分量記為：(2.3.1)如果在高斯賽德爾迭代法中引入加權因子，將和進行加權組合作為第m+1次迭代的第i個分量，就得到了松弛迭代法的迭代公式：（）用（）式求方程組解得方法稱為帶有松弛因子的松弛迭代法：當時稱為超松弛迭代方法；時成為低松弛迭代方法；時即為高斯賽德爾迭代方法。4、 對稱超松弛迭代法（SSOR方法）如果我們先按自然次序（）用向前的SOR方法逐點計算：(2.4.1)然后再按相反的次序（）用向后的SOR方法逐點計算：(2.4.2)5、

13、 共軛梯度法（CG方法）對于線性方程組的求解問題，定義二次泛函(2.5.1)則可證明求線性方程組的解等價于求二次泛函（）的極小點。由此，給定了初值，按照某一方向去求式（）的極小值點，就得到先一個迭代值，再由出發(fā)求最終逼近精確解。若取求最小值的方向為在（）處的負梯度方向，就是所謂的最速下降法。然而理論和實際計算表明這個方法的收斂速度較慢，共軛梯度法則是在處的梯度方向和這一步的修正方向所構成的二維平面內，尋找使減少最快的方向作為下一步的修正方向，即求極小值的方向。計算公式為：(2.5.2)再逐次計算：（表示x,y的內積）可以證明當時，有=0，=0。從而，形成一組共軛向量組；形成以正交向量組。后者說明若沒有舍入誤差的話，至多n次迭代就可以得到線性方程組的精確解，然而在實際計算中一般都有舍入誤差，所以并不是真正相交，所以也只能得到近似解。6.線性迭代的收斂性討論設使方程組的解，對于給定的初始向量，若由于某種迭代方法產生的向量序列有，則稱該方法收斂，否則方法發(fā)散。對任意初始向量及任意右端向量f，由迭代產生的迭代向量序列收斂的充要條件是譜半徑r(B)<1判別條件1：若|B|<1, 則迭代對任何初始向量都收斂.判別條件2：如果A為嚴格對角占優(yōu)陣，則其 Jacobi迭代和Seidel迭代對任何初始向量都收斂。 判別條件3： 如果A為對稱正定陣