于冬冬，王樂洋

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

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

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

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

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

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

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

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

1.1 垂直位移Mogi模型

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

(1)

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

(2)

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

當(dāng)?shù)貧椴此山橘|(zhì)時，則基于Mogi模型的垂直位移式可等價為[10,20]

(3)

若有n個監(jiān)測點，則有[10]

(4)

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

(5)

1.2 病態(tài)總體最小二乘問題虛擬觀測解法

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

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

4)重復(fù)第(2)，(3)步，當(dāng)

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

1.3 病態(tài)總體最小二乘問題譜修正迭代法

4)重復(fù)第(2)、(3)步，當(dāng)

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

1.4 病態(tài)總體最小二乘問題共軛梯度法

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

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

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

3)構(gòu)造搜索方向dk

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

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

2 算例及分析

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

表1 不同方法的解算結(jié)果

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

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

3 結(jié) 論

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

[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ī)車廠地裂縫強(qiáng)活動深度段位錯模型反演[J].地殼形變與地震,1998,18(3):80-84.

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

[8] 姜友誼,巨小文,胡亞軒,等.借助 Matlab 遺傳算法反演火山區(qū)壓力源參數(shù)[J].西安科技大學(xué)學(xué)報,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] 王樂洋.基于總體最小二乘的大地測量反演理論及應(yīng)用研究[D].武漢:武漢大學(xué)，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] 袁振超,沈云中,周澤波.病態(tài)總體最小二乘模型的正則化算法[J].大地測量與地球動力學(xué),2009,29(2):131-134.

[14] 葛旭明,伍吉倉.病態(tài)總體最小二乘問題的廣義正則化[J].測繪學(xué)報,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] 王樂洋,許才軍.附有病態(tài)約束反演問題的嶺估計解法[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2011,36(5):612-616.

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

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

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

[20] 于冬冬.病態(tài)總體最小二乘問題解算方法及應(yīng)用研究[D].南昌：東華理工大學(xué)，2015.

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

[責(zé)任編輯：李銘娜]

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模型反演的病態(tài)總體最小二乘解算方法[J].測繪工程，2017，26(7)：22-25.

10.19349/j.cnki.issn1006-7949.2017.07.005

2016-06-15

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

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

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

P317

A

1006-7949(2017)07-0022-04