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

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

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

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

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

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

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

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

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

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

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

11、+1次近似(2.1.2)給定初值，雅可比迭代公式格式為：(2.1.3)定義1： 如果向量序列中的向量每一個(gè)分量當(dāng)時(shí)都趨于向量的對(duì)應(yīng)分量，則稱是該向量的極限，記為在一定條件下，對(duì)任意初始向量，按迭代公式（）求出的向量序列的極限存在且等于方程組的解，這種用迭代格式（）求線性方程組近似解的方法成為雅可比迭代法。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、 對(duì)稱超松弛迭代法（SSOR方法）如果我們先按自然次序（）用向前的SOR方法逐點(diǎn)計(jì)算：(2.4.1)然后再按相反的次序（）用向后的SOR方法逐點(diǎn)計(jì)算：(2.4.2)5、 共軛梯度

13、法（CG方法）對(duì)于線性方程組的求解問(wèn)題，定義二次泛函 (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。從而，形成一組共軛向量組；形成以正交向量組。后者說(shuō)明若沒(méi)有舍入誤差的話，至多n次迭代就可以得到線性方程組的精確解，然而在實(shí)際計(jì)算中一般都有舍入誤差，所以并不是真正相交，所以也只能得到近似解。6.線性迭代的收斂性討論設(shè)使方程組的解，對(duì)于給定的初始向量，若由于某種迭代方法產(chǎn)生的向量序列有，則稱該方法收斂，否則方法發(fā)散。對(duì)任意初始向量及任意右端向量f，由迭代產(chǎn)生的迭代向量序列收斂的充要條件是譜半徑r(B)1 判別條件1： 若|B|1, 則迭代對(duì)任何初始向量都收斂. 判別條件2： 如果A為嚴(yán)格對(duì)角占優(yōu)陣，則其 Jacobi迭代和Seidel迭代對(duì)任何初始向量都收斂。 判別條件3： 如果A為對(duì)稱正定陣