線性方程組迭代解法在重力學(xué)中的應(yīng)用

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

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

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

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

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

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

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

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

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

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

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

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

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