火山形變Mogi模型反演的病態總體最小二乘解算方法

于冬冬，王樂洋

(1.東華理工大學 測繪工程學院，江西 南昌 330013;2.江西建設職業技術學院 土木工程系，江西 南昌 330200;3.流域生態與地理環境監測國家測繪地理信息局重點實驗室，江西 南昌 330013;4.江西省數字國土重點實驗室，江西 南昌 330013)

火山形變Mogi模型反演的病態總體最小二乘解算方法

于冬冬1,2，王樂洋1,3,4

(1.東華理工大學 測繪工程學院，江西 南昌 330013;2.江西建設職業技術學院 土木工程系，江西 南昌 330200;3.流域生態與地理環境監測國家測繪地理信息局重點實驗室，江西 南昌 330013;4.江西省數字國土重點實驗室，江西 南昌 330013)

相對最小二乘方法，總體最小二乘顧及了觀測方程系數矩陣含有誤差的情況，然而，當系統出現病態時，總體最小二乘受病態的影響將更加明顯。因此，針對病態總體最小二乘問題解算方法的研究越來越多受到關注。文中基于總體最小二乘進行火山形變Mogi模型反演，針對反演過程中出現的病態性問題，采用虛擬觀測解法、譜修正迭代解法、共軛梯度解法，通過模擬算例驗證文中方法在抑制病態性方面的有效性。與一般總體最小二乘、正則化總體最小二乘等方法相比存在優勢。

Mogi模型；病態總體最小二乘問題；虛擬觀測法；譜修正迭代法；共軛梯度法

Mogi模型最早是由Kiyoo Mogi自1958年提出[1]，由于其在火山形變反演方面的簡單適用性，提出至今得到廣泛的應用和推廣[2-5]。在利用該模型進行火山地表形變模擬和解釋過程中，往往利用泰勒級數對其進行線性化處理。針對線性化后可能出現的病態或奇異情況，以往的研究當中，主要采用阻尼最小二乘方法[6]、最小二乘嶺估計方法[7]、遺傳算法[8]等方法來解決。上述研究中都是基于系數矩陣不含誤差情況的最小二乘方法的反演，然而，由于線性化后的系數矩陣會在一定程度上失真而含有相應的誤差。此時，利用最小二乘方法得到的結果將不是最優解。目前，基于總體最小二乘方法在火山形變反演方面的研究相對較少。文獻[9]中利用了結構總體最小二乘反演了該火山源的參數，但未顧及反演過程中可能出現的病態情形；文獻[10]利用總體最小二乘嶺估計方法解決了在迭代過程中出現的病態問題。隨著總體最小二乘方法在近幾十年來的發展，關于病態總體最小二乘的解算方法得到了越來越多的探討和研究。如正則化方法[11-14]、嶺估計方法[15-17]、虛擬觀測解法[18]、譜修正迭代法[19]及共軛梯度法[20]等。本文以總體最小二乘方法進行了基于Mogi的火山地表形變反演，利用虛擬觀測解法、譜修正迭代解法及共軛梯度解法來解決反演過程中出現病態性問題。通過模擬算例，對文中應用的三種方法進行比較，通過算例表明本文方法在火山地表垂直位移Mogi模型反演中的適用性和有效性與其他方法相比具有優越性。

1 垂直位移 Mogi模型及病態總體最小二乘問題解算方法

1.1 垂直位移Mogi模型

應用Mogi模型時，在假定條件下可得出由于壓力源引起的地表垂直形變位移量Δh與源的參數表達式為[10,20]

(1)

式中:a為壓力源半徑，ΔP為源內部由于巖漿運動引起的壓力變化量，D為源的中心深度，μ為剪切模量，r為地表徑向距離。若用源的體積變化來表示火山地表垂直位移，其式為[10,20]：

(2)

式中:K為體積彈性模量。

當地殼為泊松介質時，則基于Mogi模型的垂直位移式可等價為[10,20]

(3)

若有n個監測點，則有[10]

(4)

上式利用泰勒級數在初值θ0處展開，并舍去二次項可線性化得[10,20]

(5)

1.2 病態總體最小二乘問題虛擬觀測解法

在病態總體最小二乘虛擬觀測解法中，將實際的觀測作為第一類觀測，把參數之間相互獨立這一信息以觀測方程的形式表達，作為第二類觀測即虛擬觀測，在同精度獨立觀測條件下，將兩類觀測方程聯立進行求解，通過推導可得其迭代求解步驟為[18]

1)采用最小二乘法求得待定參數的初始值

4)重復第(2)，(3)步，當

迭代終止；ε可選較小正值作為閾值。

1.3 病態總體最小二乘問題譜修正迭代法

4)重復第(2)、(3)步，當

迭代終止；ε可選較小正值作為閾值。

1.4 病態總體最小二乘問題共軛梯度法

共軛梯度法的基本原理是在給定的初值、下降方向、步長情況下，通過迭代的方法慢慢逼近參數最佳估值的一種方法[20]。該方法主要用于無約束最優化問題的求解。當正則化矩陣為單位陣時，總體最小二乘正則化問題即可轉化上述無約束最優化問題進行求解。

則基于Armijo非精確線性搜索再開始共軛梯度法來解決上述問題的主要解算步驟為[20]

1)給定初始值X0和迭代精度0≤ε<1，計算一階導數g0=▽f(X0)，k:=0；

3)構造搜索方向dk

4)利用Amijo非精確線搜索方法確定搜索步長λk；

6)令k:=k+1，轉步驟(2)。

2 算例及分析

由理論的Mogi模型正演得到地表垂直位移，[x,y]∈[-10,10]km× [-10,10]km,相鄰垂直形變觀測點間隔為1km，如圖1所示。模型參數的真值：膨脹源的中心坐標為x0=0km，y0=0km，源深D=4km，體積增量ΔV=6×103m3。對正演得到的垂直位移加上均值為0、標準偏差為σ=0.1的獨立、同分布的高斯噪聲，從而得到含有誤差的觀測數據。解算過程中取模型參數x0，y0，D,ΔV的初始值為2km，2km，2km及2×103m3。分別利用不同方法進行反演所得結果見表1[10,20]。

表1 不同方法的解算結果

圖1 由Mogi模型正演得到的地表垂直形變

由表1中結果可以得出，每種方法得到的參數估值與真值之間的差異主要體現在源深及體積增量中。虛擬觀測解法及譜修正迭代解法所得結果均優于正則化方法，且譜修正迭代法在改善反演過程中的病態性效果最好，共軛梯度解法較正則化方法的結果接近。表明了文中所用三種方法在火山形變反演過程中的有效性和優越性。與此同時，在虛擬觀測法中明確了準則子參數的物理意義，即實際觀測方差與虛擬觀測方差之比[18]。譜修正迭代法相比其他方法在編程實現更為容易，且譜修正參數較正則化參數或嶺參數的確定更為簡單。

3 結 論

本文針對病態總體最小二乘問題的虛擬觀測解法、譜修正迭代解法及共軛梯度解法進行了簡單的介紹，并將其用于解決火山形變Mogi模型反演中出現病態性問題。分析了三種方法相比正則化方法的優勢，通過算例表明了在火山形變反演中的有效性。文中算例為模擬算例，且僅模擬了火山的垂直位移形變，關于文中所用方法在實際情況及在水平位移反演及垂直位移同水平位移的聯合反演中的應用還需一步進行研究和探討。

[1] MOGI K.Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them[J].Bulletin of the Earthquake Research Institute,University of Tokyo,1958,36:99-134.

[2] YAMAKAWA N.On the strain produced in a semi-infinite elastic solid by an interior source of stress[J].J.Seism.Soc.Japan,1955,8:84-98.

[3] YOKOYAMA I.A model for the crustal deformations a round volcanoes[J].J Phy s Earth，1971，19(3):199-207.

[4] DECKER R W,KOYANAGI R Y,DVORAK J J,et al.Seismicity and surface deformation of Mauna Loa volcano,Hawaii[J].Eos,Transactions American Geophysical Union,1983,64(37):545-547.

[5] BERRINO G,CORRADO G,LUONGO G,et al.Ground deformation and gravity changes accompanying the 1982 Pozzuoli uplift[J].Bulletin volcanologique,1984,47(2):187-200.

[6] 王慶良,王文萍.大同機車廠地裂縫強活動深度段位錯模型反演[J].地殼形變與地震,1998,18(3):80-84.

[7] 胡亞軒,郝明,王雄,等.L 曲線法在反演火山區壓力源參數中的應用[J].大地測量與地球動力學,2009,29(2):66-70.

[8] 姜友誼,巨小文,胡亞軒,等.借助 Matlab 遺傳算法反演火山區壓力源參數[J].西安科技大學學報,2008,28(1):91-95.

[9] BIFULCO I,RAICONI G,SCARPA R.Computer algebra software for least squares and total least norm inversion of geophysical models[J].Computers & Geosciences,2009,35(7):1427-1438.

[10] 王樂洋.基于總體最小二乘的大地測量反演理論及應用研究[D].武漢:武漢大學，2011.

[11] TIKHONOV A N.Solution of Incorrectly Formulated Problems and the Regularization Method[J].Soviet Mathematics Doklady,1963,4:1305-1308.

[12] GOLUB G H,HANSEN P C,O’LEARY D P.Tikhonov Regularization and Total Least Squares[J].SIAM Journal on Matrix Analysis and Applications,1999,21(1):185-194.

[13] 袁振超,沈云中,周澤波.病態總體最小二乘模型的正則化算法[J].大地測量與地球動力學,2009,29(2):131-134.

[14] 葛旭明,伍吉倉.病態總體最小二乘問題的廣義正則化[J].測繪學報,2012,41(3):372-377.

[15] HOERL A E,KENNARD R W.Ridge regression:Biased Estimation for Non-orthogonal Problem[J].Technometrics,1970,12(1):55-67.

[16] 王樂洋,許才軍.附有病態約束反演問題的嶺估計解法[J].武漢大學學報(信息科學版),2011,36(5):612-616.

[17] 王樂洋,許才軍,魯鐵定.病態加權總體最小二乘平差的嶺估計解法[J].武漢大學學報(信息科學版),2010,35(11):1346-1350.

[18] 王樂洋,于冬冬.病態總體最小二乘問題的虛擬觀測解法[J].測繪學報,2014,43(6):575-581.

[19] 于冬冬，王樂洋.病態總體最小二乘問題的譜修正迭代法[J].大地測量與地球動力學,2015，35(4)：702-706.

[20] 于冬冬.病態總體最小二乘問題解算方法及應用研究[D].南昌：東華理工大學，2015.

[21] 唐利民.非線性最小二乘的不適定性及算法研究[D].長沙：中南大學，2011.

[責任編輯：李銘娜]

The ill-posed total least squares algorithm in volcano inversion of Mogi model

YU Dongdong1,2,WANG Leyang1,3,4

(1.School of Geomatics,East China Institute of Technology,Nanchang 330013,China;2.Department of Civil Engineering,Jiangxi College of Construction,Nanchang 330200,China;3.Key Laboratory of Watershed Ecology and Geographical Environment Monitoring,NASG，Nanchang 330013,China;4.Jiangxi Province Key Lab for Digital Land,Nanchang 330013,China)

As to the least squares method,the total least squares correlates with the coefficient matrix errors.The result of TLS will be more susceptible in ill-conditioned system.The ill-posed TLS problem leads to more and more research.Due to the ill-posed problem existing in volcano inversion of Mogi model,the methods of ill-posed total least squares are applied.Comparing with general total least squares method and regularized total least squares method,the applicability and effectiveness of the algorithms applied are discussed and analyzed in the simulation example.

Mogi model;ill-posed total least squares problem;virtual observation method;iteration method by correcting characteristic values;conjugate gradient method

著錄：于冬冬，王樂洋.火山形變Mogi模型反演的病態總體最小二乘解算方法[J].測繪工程，2017，26(7)：22-25.

10.19349/j.cnki.issn1006-7949.2017.07.005

2016-06-15

國家自然科學基金資助項目(41204003)；國家重點研發計劃(2016YFB0501400)；測繪地理信息公益性行業科研專項(201512026)；流域生態與地理環境監測國家測繪地理信息局重點實驗室開放基金項目(WE2015005)；江西省教育廳科技資助項目(GJJ150595)；對地觀測技術國家測繪地理信息局重點實驗室開放基金項目(K201502)；東華理工大學博士科研啟動基金資助項目(DHBK201113)；江西省杰出人才資助計劃項目.

于冬冬(1992-)，男，碩士研究生.

王樂洋(1983-)，男，博士，副教授，碩士生導師.

P317

A

1006-7949(2017)07-0022-04