基于Extrapolation Tikhonov正則化算法的重力數(shù)據(jù)三維約束反演

劉銀萍，王祝文，杜曉娟，劉菁華，許家姝

吉林大學地球探測科學與技術學院，長春 130026

1 引 言

重力數(shù)據(jù)反演是根據(jù)地面或空中觀測數(shù)據(jù)求取地球內部物質分布或者異常源形態(tài)結構的方法.隨著密度差解釋模型的發(fā)展，可從重力數(shù)據(jù)中提取的有效信息量顯著增加，因此重力數(shù)據(jù)反演成為定量解釋中必不可少的重要環(huán)節(jié).重力場只在地球表面是已知數(shù)據(jù)，在地面以下有多種等效的密度分布可產生與地面觀測數(shù)據(jù)相同的已知場.所以如果想通過反演從地球物理數(shù)據(jù)求得太多信息，它的解就是不適定的，這是由基于位場的任何地球物理反演方法固有的不穩(wěn)定性和不唯一性造成的.

引入異常源的先驗信息是將不適定問題轉換為適定問題的有效途徑.在反問題中有兩類不同的反演方法.第一類是當我們有足夠的關于異常源形狀、走向分布、密度差隨深度的變換規(guī)律或者密度差確切值等先驗信息時，可利用這些先驗信息作為相對或者絕對約束反演重力數(shù)據(jù)［1－8］.第二類是如果我們的先驗信息不足，在最小化信息的情況下反演異常源的地下信息成為前人研究的重要分支.比如將先驗信息以一個或者多個近似權函數(shù)的形式加入到算法中，通過最小化由模型目標函數(shù)和擬合數(shù)據(jù)誤差構成的全體目標函數(shù)的方式獲得反演結果，結合磁法數(shù)據(jù)反演并用相似的深度加權函數(shù)約束密度隨深度變化的程序來處理近地面的質量集中［9－10］；也可以在考慮核函數(shù)內部固有約束和外加約束前提下，通過隨機子域方法在信息最小化條件下盡可能的減少計算所需存儲空間，獲得良好的反演結果［11－12］；或者通過計算位場數(shù)據(jù)與地下各點在測區(qū)的位場數(shù)據(jù)的歸一化之間的相關性，顯示地下物性分布特點［13－14］；再者基于非線性共軛梯度算法理論提出加權重力梯度反演方法，約束密度上下限的強制轉換函數(shù)的引入，對于獲得有物理和地球物理意義的解起到重要作用［15］.

重力數(shù)據(jù)三維反演的病態(tài)性主要表現(xiàn)是其敏感矩陣的不穩(wěn)定性，觀測數(shù)據(jù)中的微小誤差可使反演結果嚴重偏離真解，Tikhonov指出可在反演過程加入正則參數(shù)，對初始的反演方程進行微調，從而達到穩(wěn)定反演解的目的［16－18］，因此，Tikhonov正則化方法在改變位場反問題的病態(tài)性方面效果顯著.本文在密度差解釋模型基礎上繼續(xù)進行最小化先驗信息情況下的三維加權反演的研究.Tikhonov正則化方法的正則參數(shù)選擇和使用在解決非自共軛情況下不適定問題解的同時，也給近似解引進了新的誤差［16－19］.Hamarik等提出的Extrapolation Tikhonov正則化方法作為一個提高近似方法精度的工具是將不同的正則參數(shù)值對應的不同近似解通過線性結合的方式獲得最終的近似解，因正則參數(shù)的使用使得解在誤差擴大方面受到限制［20－21］.本文第一部分研究用Extrapolation Tikhonov正則化方法解決反演中解的不唯一性和不穩(wěn)定性問題，詳細分析和討論了如何利用三種后驗的Extrapolation Tikhonov正則化參數(shù)的選擇規(guī)則綜合選擇適合的正則參數(shù)以提高反演結果的精度.第二部分在反演過程中用改進后的深度加權函數(shù)作用于核函數(shù)以有效的抵消隨深度增加核函數(shù)迅速衰減對反演結果產生的負面影響，模型試驗表明，改進后的深度加權函數(shù)與原方法相比，改善了異常體底部密度分布存在的發(fā)散問題，對于埋深較深的異常體密度分布特征具有很好的反演效果.本文的第三部分引用上下限約束，用強制轉換函數(shù)將平滑約束引進到反演中，將超出物理意義和地球物理意義范圍的反演參數(shù)限制在解釋者可接受的范圍內，從而有效的限制了反演結果中物性分布的數(shù)值范圍.

2 基于Extrapolation Tikhonov正則化方法的重力反演理論

2.1 重力數(shù)據(jù)三維反演理論

本文將包含異常源的地下區(qū)域離散為多行列的并列單元立方體構成的網格狀區(qū)域，每個具有特定尺寸的立方體單元網格被視為一個密度單元.這些密度單元是重力數(shù)據(jù)三維反演在右手坐標系中的待估計參數(shù).設地下任一立方體單元網格中心坐標為（x′，y′，z′），x軸和y軸為水平面內方向相互垂直的兩個坐標軸，z軸垂直與水平面指向地下.在觀測點r＝（x，y，z）處，單個立方體網格單元產生的萬有引力積分形式［22］為

式中γ是牛頓萬有引力常量，ρ（r′）是異常體中r′＝（x′，y′，z′）處的單元網格立方體的密度差分布，g（r）是測點距離異常源為r的萬有引力.無限小體積dV′＝dx′dy′dz′的萬有引力垂直分量即重力異常用gz（r）表示，（1）式變成

（2）式的積分形式可以表示為如下的離散形式：

重力異?？梢暈殡x散網格單元密度值ρ與常量核函數(shù)A的乘積，因此，（3）式的緊致形式為

其中A∈RM×N是核函數(shù)矩陣，ρ∈RN是地下物性分布向量，g∈RM是地面觀測重力異常.M和N分別代表測量區(qū)域的觀測異常數(shù)據(jù)數(shù)量和測量區(qū)域的網格單元數(shù)量.

2.2 基于Extrapolation Tikhonov正則化方法的物性反演

公式（4）中核函數(shù)A中有一定數(shù)量的含有不同階次的無限接近零的奇異值，因此核函數(shù)矩陣A的逆不存在.為了解決這個問題，前人曾使用過擬廣義逆［23］、神經網絡［24］等多種優(yōu)化方法.本文利用基于原正則化近似方法的Extrapolation Tikhonov正則化方法反演重力數(shù)據(jù)［20－21］.在正則化方法中假設含有噪聲的重力數(shù)據(jù)gδ可代替精確重力數(shù)據(jù)g，且該噪聲數(shù)據(jù)滿足誤差標準‖gδ－g‖≤δ，其中δ為解釋人員設置的擬合誤差標準，那么反演密度的近似值ρ在非自共軛的情況下可用Tikhonov正則化方法ρα（αI＋A＊A）－1＝A＊gδ求得，其中正則化參數(shù)α＞0，I是與A＊A維數(shù)相同的單位矩陣.為了減少因正則化參數(shù)的加入而擴展的反演誤差，Hamarik等提出用拉格朗日插值方法將用不同的正則化參數(shù)獲得的正則近似解經線性結合獲得想要的近似解，即Extrapolation Tikhonov正則化方法.

我們設ραi，i＝1，2，…，n是根據(jù)不同的正則化參數(shù)αi－1＝αi／q（q＞1）計算獲得的密度正則近似解，其中q是用來產生一系列正則化參數(shù)的比例常數(shù)，其值的大小決定了產生的正則化參數(shù)序列收斂的快慢；與上述各參數(shù)相應的計算密度差的Extrapolation Tikhonov正則化公式ρn，α可表示為：

基于Extrapolation Tikhonov正則化方法的重力三維反演，參數(shù)αi和n的選擇直接影響著病態(tài)線性反演結果的好壞，也直接關系到近似解與精確解之間的近似程度.一般情況下我們先確定αi和n中的一個，另外一個成為要求的正則化參數(shù).本文中先確定一組α1，α2，…，αN，N為已知的常數(shù)，然后根據(jù)擬合誤差水平確定α，進一步確定參數(shù)n（1≤n≤N）.本文選取離散原則、單調誤差原則和平衡原則（見附錄A）來確定參數(shù)序列中滿足先驗誤差水平的正則化參數(shù)α和n，并分析由上述選擇方法獲得的參數(shù)值計算的Tikhonov正則近似解及擬合誤差，同時計算相應的Extrapolation Tikhonov正則近似解及擬合誤差，經數(shù)值分析證明后者在相同的擬合誤差水平下能夠獲得更為精確的反演結果.若設定相同的誤差標準，新方法與原方法相比，根據(jù)所使用的數(shù)值模型不同，擬合精度可提高1到5個數(shù)量級不等.

在Extrapolation Tikhonov正則參數(shù)計算中，本文設αn＝0.95，αi－1＝αi／q（q＞1），n＝500，q分別取試算，結果表明q越大收斂速度越快，所需項數(shù)n值越小，計算速度越快，但是相對而言擬合精度受微小的影響.在計算過程中，為了編程方便，我們可將基于Extrapolation Tikhonov方法的密度分布表示為

（6）式中di同（5）式，此處的n值為外推總項數(shù)，經數(shù)值試算表明n可取2～10的任意常數(shù)，k＝1，2，…，n，此處的k值為附錄A中的n值，確定方法詳見附錄A.

2.3 改進的深度加權函數(shù)

在重力數(shù)據(jù)三維反演中，核函數(shù)隨深度增加而自然衰減的特點使反演獲得的密度分布具有明顯的趨膚效應，這種情況下獲得的密度分布與真實異常源的密度分布相背.為了抵消這種衰減，前人討論了多種適用于重磁數(shù)據(jù)反演的深度加權方法［9－10，15］.在重力數(shù)據(jù)反演及重力梯度反演中，這些方法分別用深度加權函數(shù)直接作用于核函數(shù)或者通過再加權抵消核函數(shù)隨深度自然衰減，消除其近地表權重過大而使反演密度分布不符合真實異常源的情況.本文使用的加權函數(shù)由Commer提出的深度加權函數(shù)發(fā)展而來，Commer的深度加權函數(shù)適用于重力梯度三維反演，其表達式為［15］

式中θ＝0.001（由Commer給出）或者其它很小的值，它的大小直接決定近地表頂層壓制作用的大小，其值越小壓制作用越大，反演結果越集中，反之亦然.d＝Lz，即離散區(qū)域垂直深度值，若模型頂部和底部深度已知，在Commer的函數(shù)中zc＝（z1＋z2）／2是一個設定的常數(shù)，其中z1和z2是異常源頂部和底部埋深值.為了獲得適合重力數(shù)據(jù)三維反演的深度加權函數(shù)，我們做了大量的實驗來研究（7）式中各參數(shù)與ρ之間的關系.通過單個立方體模型的實驗我們獲得一系列關于zc與異常體平均深度的數(shù)值關系，見圖1.圖1中x是模型的平均深度值，y是zc隨x變化能獲得好的反演結果時的值.由圖1可見，這兩個參數(shù)之間存在著如下線性關系，因此，適應于重力數(shù)據(jù)三維反演的深度加權函數(shù)可表示為如下形式：

圖1 參數(shù)zc與平均深度之間的線性關系圖Fig.1 The linear relationship between zc and mean depth

式中α的范圍是0.8～0.9，β的范圍是24～26.上述數(shù)據(jù)是根據(jù)大量模型試驗反演擬合初始模型在獲得好的效果時所得，若核函數(shù)為A，深度加權函數(shù)按照設定的尺寸網格化參數(shù)獲得離散值，再作用于核函數(shù)并參與反演過程即可.當應用于實測重力數(shù)據(jù)時，根據(jù)已知地質或者地球物理信息給出異常體的中心埋深，先將數(shù)據(jù)做位場分離，然后用上述方法反演局部重力異常場數(shù)據(jù)，即可獲得異常體形態(tài)和分布特征等信息.

本文理論模型縱向切面圖如圖2a和2b所示，立方體尺寸為200m×200m×200m，x和y坐標范圍分別是：（400，600）和（400，600），測區(qū)范圍是1000m×1000m.異常體的密度值為1g／cm3，背景場密度值為0g／cm3.圖中黑色代表異常體模型縱切面分布特征.圖2c—2n是圖2a和2b在不同深度加權作用下的反演結果切面圖，圖中黑色線框代表異常體實際所在位置.由圖2可見，對于埋深不同的異常體，在改進后的深度加權函數(shù)作用下，反演結果的橫向和縱向分辨率都比較理想（見圖2（c，f，i，l）），與之相比如果不做深度加權處理，反演密度分布都集中在地表面附近，其縱向分辨率非常差（見圖2（d，g，j，m）），趨膚效應十分嚴重，當Commer的深度加權函數(shù)應用到反演中時（見圖2（e，h，k，n）），對于埋深較淺的異常體，反演結果的橫向分辨率較好，異常體頂部位置的確定較為準確，異常體底部密度分布形態(tài)失真，對于埋深較深的異常體，其反演結果明顯高于實際異常體所在位置，反演結果發(fā)散較嚴重.

2.4 密度分布的上下限約束函數(shù)

圖2 不同深度的單個立方體模型在不同的深度加權作用下的反演結果切面圖（a）頂部埋深100m的立方體模型；（b）頂部埋深200m的立方體模型.（c—n）中黑色邊界范圍代表正演模型所在位置.（c—e）是圖（a）分別在本文改進的深度加權、未加權及Commer的深度加權作用下獲得的反演結果縱切面圖；圖（f—h）是（a）在上述深度加權順序作用下獲得的反演結果橫切面 圖；圖（i—k）是（b）在上述深度加權順序作用下獲得的反演結果縱切面圖；圖（l—n）是（b）在上述深度加權順序作用下獲得的反演結果橫切面圖.Fig.2 The inversion density distribution on different depth weighting slices of the single cube model at different depths（a）Single model with the upper depth 100m.（b）Single model with the upper depth 200m.The black boundaries in（c）—（n）represent the location of the forward model.（c）—（e）are the inversion solutions displayed in vertical slices corresponding to（a）and worked by the improved depth weighting function in this paper，unweighting and Commer′s depth weighting function，respectively.（f）—（h）are the inversion solutions displayed in horizontal slices corresponding to（a）；（i）—（k）are the inversion solutions displayed in vertical slices corresponding to（b）.（l）—（n）are the inversion solutions displayed in horizontal slices corresponding to（b）.

地球物理反問題解的不唯一和不穩(wěn)定性導致反演過程中不符合實際的模型計算的響應與觀測數(shù)據(jù)的擬合差與觀測數(shù)據(jù)和由真正的模型計算的響應的擬合差相一致.因此約束的使用對于減少解的奇異性起到至關重要的作用，我們可使用通過不同方式得到的先驗信息約束反演結果.在重力數(shù)據(jù)反演中，早期通過固定上下限的方法將密度值限制在一定的區(qū)域內，后來發(fā)展了強制密度差為正的方法約束反演結果，本文使用的上下限約束是較小程度的限制解的約束方法之一.反演結果中不具有物理意義的負的或者不真實的過高密度值可以通過上下限約束函數(shù)轉換到合理的值域中.Kim等通過對數(shù)約束實現(xiàn)了將電導率值限制在合理的值域范圍內［25］.那么如果我們獲得了關于密度范圍的信息，在反演中可用上下邊界約束處理反演結果，處理后極少數(shù)不具有物理意義的結果再用強制為正的方法加以約束.Cardarelli等將對數(shù)轉換方法成功地應用于電阻率層析成像中［26］.除了對數(shù)轉換，Commer等引進了反雙曲正切轉換將電導率參數(shù)限制在一個轉換空間的無界區(qū)域［27］.這說明這種轉換函數(shù)在限制電導率可能發(fā)生的無界反演的結構邊界方面很有意義.本文中我們使用的上下限約束函數(shù)［27］表示為：

式中a和b分別為強制邊界約束的上下限，通過公式（9）可將反演密度中不具有地球物理意義的密度值限制在由先驗信息確定的有意義的密度分布范圍內.函數(shù)中參數(shù)p控制著轉換區(qū)域中x＝0處的值的轉換結果，此時方程變?yōu)棣眩剑╝＋b）／2.如圖3a所示，不同的p值，對邊界約束的結果是有差異的，其值越大，則由x向ρ轉換的函數(shù)越陡峭，也就是說在x＝0附近很快達到ρ的上下限；反之，轉換函數(shù)的曲線則較緩，在x＝0周圍更大的數(shù)值區(qū)域內才可達到ρ的上下限.參數(shù)p值相同的情況下，雖然不同的上下限對應的轉換空間的值不同，轉換函數(shù)的陡峭程度也不同，如圖3b所示，但是能在相同的x的值域內達到ρ的上下限，本文中p值?。?.2～1.4之間可使反演結果較合理的轉換到符合地球物理意義的數(shù)值范圍內.在處理實測數(shù)據(jù)時，將反演結果視為（9）式中的x，通過該式實現(xiàn)反演數(shù)據(jù)的上下限約束.

圖3 上下限約束函數(shù)中參數(shù)變化對轉換空間值的影響（a）不同的p值的數(shù)據(jù)轉換特征；（b）同p值時不同上下限的數(shù)據(jù)轉換特征.Fig.3 Influence of the change of the parameter on the transform space in the upper and lower function（a）The characteristic of transform values with different p.（b）The characteristic of transform values of different upper and lower functions with the same p.

3 模型計算及誤差分析

下面使用與前期作者（李耀國，姚長利等）相似的巖脈組合模型驗證基于Extrapolation Tikhonov正則化方法的重力數(shù)據(jù)3D反演方法的有效性［10－11］.本文測區(qū)為一1000m×1000m的正方形區(qū)域，由10×10個數(shù)據(jù)點構成，點線距均為100m.網格化地下由10×10×5個網格單元構成，每個單元的尺寸為100m×100m×100m.

地下正演模型如圖4及圖5所示，模型由兩個高密度體構成，左側為一直立棱柱體，兩個角點的坐標分別為（200，400，100）和（300，600，300），密度差為0.9g／cm3；右側為傾角是45°傾斜巖脈，上頂埋深為100m，底部埋深為400m，上頂兩對角角點坐標為（600，300）和（700，700），底部兩對角角點坐標為（320，300）和（420，700），密度差為1g／cm3.圖4a所示為y＝500m處的垂直切面，上面對應的是由該切面處模型體產生的剖面重力異常曲線，異常值單位為mGal；圖4b所示為由上述組合模型產生的原始異常等值線圖.圖中最大值與圖4a中剖面曲線的最大值一致.組合模型在z＝100m及z＝300m處的橫截面示意圖如圖5a和5b所示.

反演時，在模型原始異常中加入3%的高斯噪聲，反演的密度分布在不同切面處的分布特點如圖6所示，圖6（a，b）是基于Extrapolation Tikhonov正則化方法反演的密度分布圖，圖6（c，d）是基于Tikhonov正則化方法反演的密度分布圖，黑色框所圈定的范圍是正演組合模型體所在位置及分布范圍，圖6（b，d）為y＝500m處的切面圖，由圖中結果可看出，利用本文的正則化方法反演的異常位置及分布區(qū)域與模型中異常源的分布特征基本相符，圖6b較圖6d能夠更加準確的給出正演模型體的位置，圖6d的反演結果顯示，埋深較深的反演密度更發(fā)散.圖6（a，c）為z＝250m處的切面圖，兩圖都可較準確的給出模型體所在位置，但是圖6c兩個相鄰的模型體邊界劃分較圖6a更模糊，因此其反演效果較本文中的方法差.圖6所示的反演結果表明，與原正則化方法相比，新方法能將異常體在不同切面處的密度分布準確的呈現(xiàn)，與圖4及圖5所示的模型異常源的位置，傾向及區(qū)域范圍吻合較好，兩種方法的反演結果比模型體略分散一點，新舊正則化方法反演的密度的最大值分別較模型密度大0.06g／cm3和0.052g／cm3.反演時在不同的誤差水平下由原正則化方法及Extrapolation Tikhonov正則化方法計算出的擬合誤差如表1所示，由此可見Extrapolation Tikhonov正則化方法在與原方法耗時幾乎相同的情況下，顯著減小反演密度分布計算的預測數(shù)據(jù)與觀測數(shù)據(jù)之間的擬合誤差，計算中n取5，q取2，迭代次數(shù)最小為12次，最大為27次.

圖5 組合模型不同橫截面示意圖（a）z＝100m；（b）z＝300m．Fig．5 Synthetic model at different depths（a）z＝100m；（b）z＝300m．

表1 Tikhonov正則化方法與Extrapolation Tikhonov正則化方法在不同誤差下的擬合誤差Table 1 The calculated errors with different error levels given by Tikhonov Regularization method and Extrapolation Tikhonov Regularization method

4 結 論

（1）將基于拉格朗日插值方法的Extrapolation Tikhonov正則化方法應用于重力數(shù)據(jù)的三維反演，深入討論三種選擇兩個正則化參數(shù)的原則，詳細計算分析數(shù)值.依據(jù)這些數(shù)值分析結果再結合三種參數(shù)選擇方法，確定最終反演所需的正則化參數(shù).該方法的引入很大程度上減小了因正則化參數(shù)的引入而在反演結果中介入的新的誤差，進一步提高了反演擬合誤差的精度.

（2）改進后的深度加權函數(shù)更適宜重力數(shù)據(jù)的三維反演.在詳細分析函數(shù)中各參數(shù)對趨膚效應的作用后，給出了參數(shù)合理的取值范圍.該函數(shù)對消除核函數(shù)隨深度增加快速衰減的自然屬性起到有效的控制作用，進一步增強了對近地表的物性分布過度集中的壓制效果，與未加權的反演方法及前人的反演方法相比，在反演異常體底部密度分布方面具有更大的優(yōu)勢.

（3）本文引進了上下限約束函數(shù)，深入討論了函數(shù)中參數(shù)選擇與設定的量化細節(jié)，對反演具有重要意義，且該函數(shù)將不具有地球物理意義的反演結果轉化到合理而有地球物理意義的數(shù)據(jù)范圍內，給進一步的解釋工作帶來方便.

（4）通過多組模型反演計算及結果比對，證明基于Extrapolation Tikhonov正則化方法的反演方法較基于原Tikhonov正則化方法的反演方法能夠更準確的反演出異常體所在位置及密度分布特征，在計算時間幾乎不變的情況下，減小了觀測數(shù)據(jù)與預測數(shù)據(jù)的擬合誤差.

附錄A Extrapolation Tikhonov正則化參數(shù)的選擇方法

附錄詳細論述Extrapolation Tikhonov正則化方法三個參數(shù)選擇原則，以用于選擇適當?shù)恼齽t化參數(shù)［20－21］.

1 分別使用三個原則選擇正則化參數(shù)αD、αME以及αL和nD、nME以及nI的方法

1.1 基于離散原則的參數(shù)選擇方法

解釋人員給定誤差標準δ后，在αi序列中第一個滿足Aραi－g＜δ的αD即為所選第一個參數(shù)值，該參數(shù)所對應的k＝iD為第二個參數(shù)值nD.

1.2 基于單調誤差原則的參數(shù)選擇方法

αi序列中第一個滿足下式的αME即為所求第一個參數(shù)值：

當αME確定后，再進一步確定k值.設ρn，α＝A＊wi，i＝1，2，…，是Aρ＝g的近似解序列.iME是第一個使（A1）式成立的序列號：

若（A1）式成立，則有ρi－ρ≤ρi－1－ρ，i＝2，…，iME，因此，在用單調誤差選擇第一個參數(shù)后，滿足（A1）式的k＝iME為Extrapolation Tikhonov正則近似解ρnME的第二個參數(shù)值nME.

1.3 基于平衡原則的參數(shù)選擇方法

2 以nD和nME為例討論綜合確定第二個參數(shù)的方法則項數(shù)n∈｛nD，nME｝能保證使‖ρn，α－ρ‖趨近于0，即如果（A3）式成立，我們認為ρn，α比ρn－1，α更精確，此時的解為所求近似解，n取nD和nME中較大的值作為第二個參數(shù)值，用同樣的方法可從已選擇出的n值及nI中確定最后的參數(shù)值.

（

）

［1］ Blakely R J.Potential Theory in Gravity and Magnetic Applications.Cambridge：Cambridge University Press，1995.

［2］ Last B J，Kubik K.Compact gravity inversion.Geophysics，1983，48（6）：713－721.

［3］ Guillen A，Menichetti V.Gravity and magnetic inversion with minimization of a specific functional.Geophysics，1984，49（8）：1354－1360.

［4］ Barbosa V C F，Silva J B C.Generalized compact gravity inversion.Geophysics，1994，59（1）：57－68.

［5］ Medeiros W E，Silva J B C.Geophysical inversion using approximate equality constraints.Geophysics，1996，61（6）：1678－1688.

［6］ 陳石，王謙身，盧紅艷等.三維重力反演新方法——調和密度成像理論與模型測試.地球物理學進展，2010，25（6）：1917－1925.

Chen S，Wang Q S，Lu H Y，et al.A new method of gravity inversion——theory and model test of harmonic density imaging.ProgressinGeophysics（in Chinese），2010，25（6）：1917－1925.

［7］ 柯小平，王勇，許厚澤等.青藏高原地殼三維密度結構的重力反演.地球物理學進展，2009，24（2）：448－455.

Ke S P，Wang Y，Xu H Z，et al.The three－dimensional crustal density structure of the Tibetan plateau from gravity inversion.ProgressinGeophysics（in Chinese），2009，24（2）：448－455.

［8］ 柯小平，王勇，許厚澤.基于變密度模型的位場界面反演.地球物理學進展，2006，21（3）：825－829.

Ke S P，Wang Y，Xu H Z.Interface depths inversion of potential field data with variable density model.Progressin Geophysics（in Chinese），2006，21（3）：825－829.

［9］ Li Y G，Oldenburg D W.3－D inversion of magnetic data.Geophysics，1996，61（2）：394－408.

［10］ Li Y G，Oldenburg D W.3－D inversion of gravity data.Geophysics，1998，63（1）：109－119.

［11］ 姚長利，鄭元滿，張聿文.重磁異常三維物性反演隨機子域法方法技術.地球物理學報，2007，50（5）：1576－1583.

Yao C L，Zheng Y M，Zhang Y W.3Dgravity and magnetic inversion for physical properties using stochastic subspaces.ChineseJ.Geophys.（in Chinese），2007，50（5）：1576－1583.

［12］ 姚長利，郝天珧，管志寧等.重磁遺傳算法三維反演中高速計算及有效存儲方法技術.地球物理學報，2003，46（2）：252－258.

Yao C L，Hao T Y，Guan Z N，et al.High－speed computation and efficient storage in 3－D gravity and magnetic inversion based on genetic algorithms.ChineseJ.Geophys.（in Chinese），2003，46（2）：252－258.

［13］ 郭良輝，孟小紅，石磊等.重力和重力梯度數(shù)據(jù)三維相關成像.地球物理學報，2009，52（4）：1092－1106.

Guo L H，Meng X H，Shi L，et al.3－D correlation imaging for gravity and gravity gradiometry data.ChineseJ.Geophys.（in Chinese），2009，52（4）：1092－1106.

［14］ 孟小紅，劉國峰，陳召曦等.基于剩余異常相關成像的重磁物性反演方法.地球物理學報，2012，55（1）：304－309.

Meng X H，Liu G F，Chen Z X，et al.3－D gravity and magnetic inversion for physical properties based on residual anomaly correlation.ChineseJ.Geophys.（in Chinese），2012，55（1）：304－309.

［15］ Commer M.Three－dimensional gravity modelling and focusing inversion using rectangular meshes.GeophysicalProspecting，2011，59（5）：966－979.

［16］ Tikhonov A N，Arsenin V Y.Solution of Ill－Posed Problems.New York：V.H.Winston ＆Sons，1977.

［17］ Tikhonov A N.Regularization of Ill－Posed Problems.Doklady Akad.Nauk，SSSR，1963，153：49－52.

［18］ Zhdanov M S.Geophysical Inverse Theory and Regularization Problems.New York：Elsevier Science，2002.

［19］ Lawson C J，Hanson R J.Solving Least Squares Problems.Englewood Cliffs：Prentice－Hall Inc，1974.

［20］ H?marik U，Palm R，Raus T.Extrapolation of Tikhonov and Lavrentiev regularization methods.JournalofPhysics：ConferenceSeries，2008，135（1）：1－8.

［21］ H?marik U，Palm R，Raus T.Use of extrapolation in regularization methods.JournalofInverseandIll－Posed Problems，2007，15（3）：277－294.

［22］ Zhdanov M S.Integral Transforms in Geophysics.New York：Springer Verlag，1988.

［23］ 劉天佑，管志寧.對二維磁異常利用廣義逆矩陣的選擇法.地球物理學報，1986，29（4）：390－398.

Liu T Z，Guan Z N.Selective method of generalized inverse matrix for the interpretation of two－dimensional magnetic anomalies.ChineseJ.Geophys.（in Chinese），1986，29（4）：390－398.

［24］ 管志寧，侯俊勝，黃臨平等.重磁異常反演的擬BP神經網絡方法及其應用.地球物理學報，1998，41（2）：242－251.

Guan Z N，Hou J S，Huang L P，et al.Inversion of gravity and magnetic anomalies using pseduo－BP neural network method and its application.ChineseJ.Geophys.（in Chinese），1998，41（2）：242－251.

［25］ Kim H J，Song Y，Lee K H.Inequality constraint in leastsquared inversion of geophysical data.EarthPlanetsSpace，1999，51：255－259.

［26］ Cardarelli E，F(xiàn)ischanger F.2Ddata modelling by electrical resistivity tomography for complex subsurface geology.GeophysicalProspecting，2006，54（2）：121－133.

［27］ Commer M，Newman G A.New advances in threedimensional controlled－source electromagnetic inversion.GeophysicalJournalInternational，2008，172（2）：513－535.