王樂洋，李海燕，溫揚茂，許才軍

1.東華理工大學測繪工程學院，江西 南昌 330013； 2.武漢大學測繪學院，湖北 武漢 430079； 3.流域生態(tài)與地理環(huán)境監(jiān)測國家測繪地理信息局重點實驗室，江西 南昌 330013

?

地震同震滑動分布反演的總體最小二乘方法

王樂洋1,3，李海燕1，溫揚茂2，許才軍2

1.東華理工大學測繪工程學院，江西 南昌 330013； 2.武漢大學測繪學院，湖北 武漢 430079； 3.流域生態(tài)與地理環(huán)境監(jiān)測國家測繪地理信息局重點實驗室，江西 南昌 330013

同震滑動分布參數(shù)與地表形變間的線性關(guān)系依賴于格林函數(shù)矩陣的構(gòu)造，格林函數(shù)矩陣元素與破裂面位置、幾何參數(shù)、破裂方式及位錯模型假設(shè)等因素有關(guān)。本文嘗試考慮格林函數(shù)矩陣元素的誤差來補償上述原因在一定程度上對反演參數(shù)的影響，采用同時顧及系數(shù)矩陣(格林函數(shù)矩陣)和觀測向量兩者誤差的總體最小二乘方法反演同震滑動分布。首先確定了系數(shù)矩陣元素和觀測向量的協(xié)因數(shù)矩陣，考慮到格林函數(shù)矩陣的病態(tài)性(秩虧)，借助拉普拉斯二階平滑得到正則化矩陣，采用總體最小二乘正則化法反演同震滑動分布。并對2009年意大利中部拉奎拉(L’Aquila)Mw6.3級地震實例進行同震滑動分布反演研究。結(jié)果表明，拉奎拉地震的走向為144.37°,傾角為59.06°，滑動分布的最大滑動量為0.95 m，平均滑動角為-96.4°，主要滑動深度為4～15 km的范圍，地震矩為3.63×1018N·m，對應(yīng)的矩震級為Mw6.34?？傮w最小二乘與最小二乘法的滑動分布解存在一定差別，但差別的量級在10-4以內(nèi)。

同震滑動分布；總體最小二乘；反演；正則化；拉奎拉地震

大地測量反演主要是通過大地測量觀測手段獲得地表形變值求解地球物理模型參數(shù)的過程，在震源機制研究中，位錯模型是反演震源參數(shù)的有效模型[1]，也是最常用的模型之一。利用彈性半空間均勻矩形位錯模型[2]反演地震震源參數(shù)，通常包括非線性反演斷層幾何參數(shù)和線性反演滑動分布兩個過程，在確定斷層位置等幾何參數(shù)后，滑動參數(shù)與形變觀測值為線性關(guān)系，此時滑動參數(shù)模型的系數(shù)矩陣即為格林函數(shù)矩陣。系數(shù)矩陣的元素通常是由剖分得到的每個子斷層塊的單位滑動引起的地表點的位移構(gòu)成，矩陣中的元素與模型本身、模型的幾何參數(shù)、觀測點位等因素有關(guān)，地殼的復(fù)雜性、模型參數(shù)求解存在的不確定性及形變導(dǎo)致的觀測點的坐標誤差等都會引起系數(shù)矩陣的變化，可以將這部分影響看成系數(shù)矩陣存在誤差，并考慮系數(shù)矩陣誤差對滑動參數(shù)解的影響。目前在求解滑動參數(shù)時，通常都是采用普通最小二乘法，即只考慮了觀測形變量的誤差影響；若采用同時顧及系數(shù)矩陣和觀測向量誤差的總體最小二乘方法[3-6]來求解滑動參數(shù)，顯然更具合理性。而國內(nèi)外未見相關(guān)研究發(fā)表。

基于大地測量觀測資料(如InSAR、GPS數(shù)據(jù)等)反演同震滑動分布是研究震源參數(shù)的重要手段，國內(nèi)外許多學者利用InSAR或GPS等數(shù)據(jù)對同震滑動分布進行了研究[7-12]，研究結(jié)果充分體現(xiàn)了采用InSAR或GPS數(shù)據(jù)進行滑動分布反演的優(yōu)越性。

本文研究中，首先給出滑動分布反演的總體最小二乘函數(shù)模型，考慮到系數(shù)矩陣的病態(tài)性，借助拉普拉斯二階平滑得到了正則化矩陣，通過正則化法解決病態(tài)情形下總體最小二乘問題。模擬滑動分布反演試驗，并利用Envisat衛(wèi)星升降軌數(shù)據(jù)約束同震地表形變場研究2009年意大利拉奎拉(L’Aquila)Mw6.3級地震滑動分布，研究了同時顧及系數(shù)矩陣誤差和觀測值誤差對滑動分布反演結(jié)果的影響。

1 滑動分布反演的總體最小二乘方法

1.1 滑動分布總體最小二乘函數(shù)模型及算法

采用Okada矩形位錯模型反演同震滑動分布，首先需要確定斷層的位置、走向等幾何參數(shù)，這些參數(shù)可以利用非線性最優(yōu)化方法反演獲得。當這些參數(shù)確定后，有

d=GM

(1)

式中，d為同震形變觀測值；M為滑動參數(shù)；G為格林函數(shù)矩陣，定義為[12]

(2)

矩陣G中ss、ds分別表示斷層子塊單位走滑和傾滑引起的地表位移，上標表示斷層的編號；下標表示同震形變觀測值個數(shù)。

同時考慮觀測值d和系數(shù)矩陣G含有誤差的情況，得到

d+e=(G+EG)M

(3)

式中，e為觀測向量的誤差向量；EG為系數(shù)矩陣G的誤差矩陣。

精細滑動分布反演時，將斷層面離散化為均勻的小斷層片，此時矩陣G是嚴重病態(tài)(秩虧)的，即為病態(tài)情形下的總體最小二乘參數(shù)估計，國內(nèi)外許多學者開展了相關(guān)研究，主要方法包括Tikhonov正則化法[13-14]、廣義正則化法[15]、嶺估計法[16]、虛擬觀測法[17]等。本文采用Tikhonov正則化解法[14]來反演滑動分布?；瑒臃植挤囱輹r，為克服系數(shù)矩陣的病態(tài)性，對參數(shù)施加拉普拉斯二階平滑約束，使得

HM=0

(4)

式中，H為拉普拉斯平滑矩陣，而R=HTH視為正則化矩陣。

根據(jù)Tikhonov正則化理論[18]，得到估計準則為

(5)

結(jié)合同震滑動分布參數(shù)反演，采用病態(tài)加權(quán)總體最小二乘正則化解法來迭代求解式(5)，具體如下：

(1) 通過最小二乘法得到參數(shù)最初估計值

(6)

(7)

(8)

(3) 迭代求解第(2)步，滿足條件

(9)

時，計算結(jié)束。

值得注意的是，在滑動分布反演算法中，由于系數(shù)矩陣G是嚴重秩虧的，因此式(6)中的N需使用廣義逆N+。

1.2 權(quán)陣P及協(xié)因數(shù)陣QG的確定

權(quán)陣P可通過觀測數(shù)據(jù)的誤差來確定，本文認為觀測數(shù)據(jù)是獨立且不相關(guān)的；假設(shè)QG=Q0?QM，必須確定Q0和QM，Q0表示系數(shù)矩陣每一列元素只存在數(shù)量因數(shù)不同且互不相關(guān)；QM表示系數(shù)矩陣每一行元素只存在數(shù)量因數(shù)不同且互不相關(guān)[4]。本文假設(shè)Q0為單位陣，即Q0=I2n，對于QM，進行如下分析。

根據(jù)位錯模型，地表某一點的形變量di(i=1,2,…,m)是斷層破裂面上每一個子塊斷層的走滑地表響應(yīng)dsj(j=1,2,…,n)和傾滑地表響應(yīng)ddj(j=1,2,…,n)的總和，即

(10)

故可以假設(shè)系數(shù)矩陣每一行元素的協(xié)因數(shù)為

(11)

即有

QM=diag([q(1),q(2),…,q(m)])

(12)

2 試驗與分析

2.1 模擬滑動分布反演

為驗證總體最小二乘方法反演滑動分布的有效性，本文首先對模擬滑動分布進行反演，分別采用總體最小二乘法和最小二乘法反演最優(yōu)滑動分布。

假定斷層破裂面的幾何參數(shù)為：斷層面中心位置X=0 km，Y=0 km，斷層上底深度MinDep=0 km，長度L=36 km，寬度W=16 km，走向角α=90°，傾角Dip=70°。共模擬了328個GPS地表形變觀測點。將斷層面均勻剖分為2 km×2 km的矩形單元，模擬的斷層面的滑動分布見圖1(a)，其中最大滑動量為1.426 8 m，平均滑動量為0.303 9 m。結(jié)合位錯模型正演得到地表形變位移，分別給形變位移觀測量和系數(shù)矩陣元素施加隨機誤差，水平位移誤差為N(0，32mm2)，垂直位移誤差為N(0，52mm2)，相應(yīng)的系數(shù)矩陣元素誤差分別為N(0，32/(2n)2mm2)及N(0,52/(2n)2mm2)?？傮w最小二乘法、最小二乘法及兩者的殘差結(jié)果分別如圖1(b)、圖1(c)及圖1(d)所示。

模擬不同的觀測數(shù)據(jù)，基于蒙特卡羅誤差傳遞法估計公式[19]給出了參數(shù)的均方根誤差(root mean square error，RMSE)作為估計參數(shù)的精度水平。RMSE的計算公式如下

(13)

反演時總體最小二乘法和最小二乘法所選取的最優(yōu)滑動因子分別為0.01、0.1。從圖1及表1的結(jié)果可以看出，總體最小二乘解與最小二乘解的結(jié)果非常接近，其中總體最小二乘結(jié)果中最大滑動量及平均滑動量分別為1.344 5 m、0.243 3 m，誤差的最大值和平均值分別為0.148 6 m、0.038 1 m，殘差中誤差(RMS)為4.7 mm，最小二乘反演結(jié)果的最大滑動量、平均滑動量及RMS與總體最小二乘結(jié)果一致，誤差最大值和平均值分別為0.148 7 m、0.038 8 m。總體最小二乘與最小二乘滑動分布解的差別最大量級在10-4，說明此時系數(shù)矩陣誤差對滑動參數(shù)解有一定影響，但影響程度有限。

表1 反演統(tǒng)計結(jié)果比較Tab.1 Comparison of the statistical inversion results

2.2 意大利拉奎拉地震震源滑動分布反演

2009年4月6日發(fā)生在意大利中部拉奎拉地區(qū)的Mw6.3級地震，據(jù)美國地質(zhì)調(diào)查局(USGS)網(wǎng)站發(fā)布，地震震中位置為(13.334°E，42.334°N)，震源深度為8.8 km。地震造成一千余人的傷亡，大量的房屋被損毀或破壞[20]。地震發(fā)生后，許多研究機構(gòu)通過不同手段對震源機制進行了較詳細的研究。其中，文獻[21—22]分別通過GPS數(shù)據(jù)反演了地震震源參數(shù)，認為此次地震的主斷層為西南走向的正斷層。文獻[20]通過ASAR數(shù)據(jù)反演斷層的幾何參數(shù)和滑動分布，結(jié)果也表明該斷層為正斷層。文獻[23]利用不同的InSAR數(shù)據(jù)集反演比較后獲得了最優(yōu)震源模型。文獻[11]則結(jié)合拉奎拉地區(qū)的地形特點，聯(lián)合Envisat衛(wèi)星的升降軌數(shù)據(jù)、ALOS衛(wèi)星數(shù)據(jù)及GPS形變觀測數(shù)據(jù)，提供比較完整的震區(qū)三維形變場，反演得到斷層幾何模型參數(shù)及滑動分布。

本文研究該地震的InSAR數(shù)據(jù)來源于文獻[11]：項目組對獲取的ASAR數(shù)據(jù)采用二通法對影像進行干涉差分處理后得到地表同震形變場；使用瑞士的Gamma軟件處理，分別采用了ESA和JAXA提供的DOR精密軌道和精密軌道用于修正軌道誤差，并利用90 m分辨率的SRTM DEM去除地形影響；采用基于能量譜的局部自適應(yīng)濾波方法對干涉圖進行了濾波處理；最后采用枝切法進行相位解纏，得到差分干涉相位。本文研究選用了其中Envisat衛(wèi)星降軌Track097和升軌Track129兩幅影像的數(shù)據(jù)，而沒有采用數(shù)據(jù)誤差較大的ALOS衛(wèi)星數(shù)據(jù)，通過對遠場數(shù)據(jù)進行均勻采樣，對近場數(shù)據(jù)進行四叉樹采樣，獲得了降軌數(shù)據(jù)1254個，升軌數(shù)據(jù)1282個(見圖2)，能有效約束斷層滑動參數(shù)的反演。

首先利用多峰值顆粒群算法(MPSO)[24-25]反演得到拉奎拉地震單一均勻滑動斷層的走向、傾角及斷層矩形中心位置等幾何參數(shù)(見表4)，在此基礎(chǔ)上，將斷層破裂面沿走向、傾向?qū)㈤L度、寬度分別拓展至30 km、26 km，并將破裂面延伸至地表，然后將拓展后的破裂面均勻剖分為2 km×2 km的矩形單元，共得到195個矩形斷層單元。采用上述的總體最小二乘方法進行滑動分布的反演研究。

關(guān)于反演涉及的觀測值和系數(shù)矩陣元素定權(quán)問題，根據(jù)InSAR觀測值的先驗中誤差確定觀測值的權(quán)陣P，由式(11)、式(12)可以確定協(xié)因數(shù)陣QM，Q0仍為單位陣，進而可以確定協(xié)因數(shù)陣QG(此時n=390)。總體最小二乘算法和最小二乘算法最優(yōu)滑動分布解對應(yīng)的最優(yōu)滑動因子分別為0.36和0.6，反演結(jié)果如圖3所示。

從圖3可以看出，震源滑動分布總體最小二乘法反演結(jié)果表明，滑動主要發(fā)生在4～15 km的深度范圍內(nèi)，總體最小二乘法的平均滑量為0.22 m，最大滑動量為0.95 m，滑動角為-96.4°，深度位于6～8 km的范圍，總體最小二乘與最小二乘結(jié)果幾乎一致。反演得到此次地震的地震級3.63×1018N·m，對應(yīng)的矩震級為Mw6.34?；瑒臃植冀獾淖畲笳`差為0.17 m，分布在斷層面下部邊緣，不在主要的滑動區(qū)域，平均誤差為0.05 m，反演模型比較穩(wěn)定。

為驗證滑動分布模型與觀測數(shù)據(jù)的擬合程度，圖4給出了升降軌觀測數(shù)據(jù)與滑動分布模型預(yù)測數(shù)據(jù)的散點分布。觀測數(shù)據(jù)擬合殘差的總體中誤差為1.26 cm，其中降軌、升軌數(shù)據(jù)擬合殘差中誤差分別為1.19 cm、1.33 cm，小于均勻滑動的殘差中誤差(數(shù)據(jù)總體擬合中誤差、降軌及升軌擬合中誤差分別為1.43 cm、1.37 cm、1.48 cm)。從圖4可以看出，總體上升降軌數(shù)據(jù)都擬合得較好，只有少數(shù)遠場數(shù)據(jù)未能很好擬合，這是由于這些遠場數(shù)據(jù)點形變值很小且數(shù)據(jù)受到斷層模型簡化、大氣延遲誤差等影響。

圖1 模擬試驗的震源滑動分布反演結(jié)果Fig.1 Slip distribution inversion results of simulation

圖2 InSAR采樣后的點位分布(五角星為震中位置)Fig.2 The distribution of InSAR samples points(The star is the epicenter)

從總體最小二乘解、最小二乘解及誤差可以看出，差別并不明顯。表2及表3比較了兩者的反演結(jié)果， 最大滑動量、平均量、RMSE最大值及平均值都幾乎一致，總體最小二乘法與最小二乘法的滑動分布平均差值量級為10-4，誤差分布平均差值量級為10-5。

表2 拉奎拉地震的TLS和LS反演結(jié)果Tab.2 Inversion results of TLS and LS in L’ Aquila earthquake

表3 拉奎拉地震TLS法與LS反演的滑動量、誤差差值Tab.3 The slip and its error residuals of TLS and LS in L’Aquila earthquake m

圖3 拉奎拉地震的滑動分布、誤差分布及TLS解與LS解的殘差Fig.3 Slip and its error distribution of L’Aquila earthquake, the residuals of TLS-LS

圖4 拉奎拉地震觀測位移值和滑動分布模型預(yù)測值位移的散點分布Fig.4 Scatter points of the observed and the distribution slip modeled displacement in L’Aquila earthquake

表4列舉了此次地震震源參數(shù)的部分研究成果。可以看出，本文確定的斷層走向、傾角、滑動角等參數(shù)與USGS機構(gòu)得到的機制解(走向122°，傾角53°，滑動角-122°)很相似，只是走向相差了22°，余震分布顯示出該發(fā)震斷層為NW-SE結(jié)構(gòu)及SW傾向的正斷層[11]，參數(shù)反演結(jié)果印證了此結(jié)論。本文方法得到的最大滑動量為0.95 m，平均滑動量為0.22 m，地震矩為3.63×1018N·m，對應(yīng)的矩震級為Mw6.34，而均勻滑動反演的滑動量、矩張量及矩震級分別為0.70 m、2.68×1018N·m、Mw6.34。文獻[11]聯(lián)合了不同波長、不同入射角的Enviast和ALOS衛(wèi)星差分干涉數(shù)據(jù)，并聯(lián)立了震區(qū)部分GPS三維形變數(shù)據(jù)，通過彈性三角位錯模型，反演得到滑動分布的最大滑動量為1.07 m，平均滑動量為0.34 m，給出的矩張量矩震級分別為3.43×1018N·m、Mw6.32；其最大滑動量大于本文的結(jié)果，這可能與約束地表的數(shù)據(jù)源有關(guān)。文獻[20]聯(lián)合ASAR數(shù)據(jù)反演得到的斷層滑動量為0.66 m，給出的矩張量2.80×1018N·m；利用L’Aquila震中附近的GPS觀測站數(shù)據(jù)，文獻[22]研究了此次地震的同震滑動分布和震后的余震滑動分布，得到的最大滑動量為1.1 m，矩張量為3.90×1018N·m。

野外地質(zhì)調(diào)查結(jié)果發(fā)現(xiàn)，拉奎拉地震在地表附近區(qū)域發(fā)生了7～10 cm的破裂，已有的GPS數(shù)據(jù)[21-22]、InSAR數(shù)據(jù)[20,27]反演結(jié)果都表明在斷層接近地表處沒有滑動的跡象，本文研究結(jié)果也驗證了這一結(jié)論，只有部分接近地表處的滑動在10 cm左右，此地表破裂可能是由地下的滑動能量傳遞到地表脆性區(qū)域造成[11]。許多研究者認為，Paganica斷層與相鄰斷層比較，其地形活動表現(xiàn)更弱，而來自遙感觀測結(jié)果卻揭示此次地震就發(fā)生在Paganica斷層帶，Paganica斷層正好位于兩個主要的活動斷層之間的位置[22,28]。此次觀測對于該地震區(qū)域和與之相似地質(zhì)條件的區(qū)域確定潛在的危險斷層顯得更為重要，有必要提高觀測手段對地震震后及震間的應(yīng)變積累進行更為細致和長久的研究。

表4 拉奎拉地震斷層幾何及滑動參數(shù)Tab.4 The various source parameters for L’Aquila earthquake

3 結(jié) 論

地震同震滑動分布反演是分析地震發(fā)震機理，研究活動斷層破裂擴展、震后形變、巖石圈應(yīng)力變化及后期地震危險性評估的基礎(chǔ)。本文研究了顧及格林函數(shù)矩陣的誤差的同震滑動分布反演。研究表明，滑動分布反演時顧及格林函數(shù)矩陣誤差是合理的，但系數(shù)矩陣元素與系數(shù)矩陣元素中誤差的比值(信噪比)較大，導(dǎo)致總體最小二乘與最小二乘反演結(jié)果差別僅在10-4以內(nèi)的量級，即本文同震滑動分布反演中格林函數(shù)矩陣含有的隨機誤差對參數(shù)估計的影響很有限。

本文利用總體最小二乘法研究了2009年意大利拉奎拉(L’Aquila)Mw6.3級地震，反演結(jié)果表明，發(fā)震斷層的走向為144.37°，傾角為59.04°，傾向為SW，是一次以正傾滑為主兼少量右旋滑動事件，滑動分布的最大滑動量為0.95 m，平均滑動量為0.22 m，主要滑動深度范圍為4～15 km。研究表明，此次地震的發(fā)震斷層Paganica斷層正好位于兩個主要的活動斷層之間的位置。因此，應(yīng)該對地震震后及震間的應(yīng)變積累進行更為細致和長久的研究。

[1] 許才軍, 尹智.利用大地測量資料反演構(gòu)造應(yīng)力應(yīng)變場研究進展[J].武漢大學學報(信息科學版), 2014, 39(10): 1135-1146, 1178.XU Caijun, YIN Zhi.Progress in Inversion for Tectonic Stress-strain Fields Using Geodetic Data[J].Geomatics and Information Science of Wuhan University, 2014, 39(10): 1135-1146, 1178.

[2] OKADA Y.Surface Deformation DUE to Shear and Tensile Faults in a Half-space[J].Bulletin of the Seismological Society of America, 1985, 75(4): 1135-1154.

[3] VAN HUFFEL S, VANDEWALLE J.The Total Least Squares Problem: Computational Aspects and Analysis[M].Philadelphia: Society for Industrial and Applied Mathematics, 1991.

[4] SCHAFFRIN B, WIESER A.On Weighted Total Least-squares Adjustment for Linear Regression[J].Journal of Geodesy, 2008, 82(7): 415-421.

[5] 王樂洋.基于總體最小二乘的大地測量反演理論及應(yīng)用研究[J].測繪學報, 2012, 41(4): 629.WANG Leyang.Research on Theory and Application of Total Least Squares in Geodetic Inversion[J].Acta Geodaetica et Cartographica Sinica, 2012, 41(4): 629.

[6] WANG Leyang, XU Guangyu.Variance Component Estimation for Partial Errors-in-Variables Models[J].Studia Geophysica et Geodaetica, 2016, 60(1): 35-55.

[7] FUNNING G J, PARSONS B, WRIGHT T J, et al.Surface Displacements and Source Parameters of the 2003 Bam (Iran) Earthquake from Envisat Advanced Synthetic Aperture Radar Imagery[J].Journal of Geophysical Research: Solid Earth, 2005, 110(B9): B09406.

[8] 許才軍, 劉洋, 溫揚茂.利用GPS資料反演汶川Mw7.9級地震滑動分布[J].測繪學報, 2009, 38(3): 195-201, 215.XU Caijun, LIU Yang, WEN Yangmao.Mw7.9 Wenchuan Earthquake Slip Distribution Inversion from GPS Measurements[J].Acta Geodaetica et Cartographica Sinica, 2009, 38(3): 195-201, 215.

[9] 張國宏, 屈春燕, 宋小剛, 等.基于InSAR同震形變場反演汶川Mw7.9地震斷層滑動分布[J].地球物理學報, 2010, 53(2): 269-279.ZHANG Guohong, QU Chunyan, SONG Xiaogang, et al.Slip Distribution and Source Parameters Inverted from Co-seismic Deformation Derived by InSAR Technology of WenchuanMw7.9 Earthquake[J].Chinese Journal of Geophysics, 2010, 53(2): 269-279.

[10] JIANG Zaisen, WANG Min, WANG Yanzhao, et al.GPS Constrained Coseismic Source and Slip Distribution of the 2013Mw6.6 Lushan, China, Earthquake and Its Tectonic Implications[J].Geophysical Research Letters, 2014, 41(2): 407-413.

[11] 溫揚茂, 何平, 許才軍, 等.聯(lián)合Envisat和ALOS衛(wèi)星影像確定L’Aquila地震震源機制[J].地球物理學報, 2012, 55(1): 53-65.WEN Yangmao, HE Ping, XU Caijun, et al.Source Parameters of the 2009 L’Aquila Earthquake, Italy from Envisat and ALOS Satellite SAR Images[J].Chinese Journal of Geophysics, 2012, 55(1): 53-65.

[12] 李志才, 張鵬, 金雙根, 等.基于GPS觀測數(shù)據(jù)的汶川地震斷層形變反演分析[J].測繪學報, 2009, 38(2): 108-113, 119.LI Zhicai, ZHANG Peng, JIN Shuanggen, et al.Wenchuan Earthquake Deformation Fault Inversion and Analysis Based on GPS Observations[J].Acta Geodaetica et Cartographica Sinica, 2009, 38(2): 108-113, 119.

[13] 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.

[14] 魯鐵定.總體最小二乘平差理論及其在測繪數(shù)據(jù)處理中的應(yīng)用[D].武漢: 武漢大學, 2010.LU Tieding.Research on the Total Least Squares and Its Applications in Surveying Data Processing[D].Wuhan: Wuhan University, 2010.

[15] 葛旭明, 伍吉倉.病態(tài)總體最小二乘問題的廣義正則化[J].測繪學報, 2012, 41(3): 372-377.GE Xuming, WU Jicang.Generalized Regularization to Ⅲ-posed Total Least Squares Problem[J].Acta Geodaetica et Cartographica Sinica, 2012, 41(3): 372-377.

[16] 王樂洋, 許才軍, 魯鐵定.病態(tài)加權(quán)總體最小二乘平差的嶺估計解法[J].武漢大學學報(信息科學版), 2010, 35(11): 1346-1350.WANG Leyang,XU Caijun,LU Tieding.Ridge Estimation Method in Ill-posed Weighted Total Least Squares Adjustment[J].Geomatics and Information Science of Wuhan University, 2010, 35(11): 1346-1350.

[17] 王樂洋, 于冬冬.病態(tài)總體最小二乘問題的虛擬觀測解法[J].測繪學報, 2014, 43(6): 575-581.DOI: 10.13485/j.cnki.11-2089.2014.0091.WANG Leyang, YU Dongdong.Virtual Observation Method to Ill-posed Total Least Squares Problem[J].Acta Geodaetica et Cartographica Sinica, 2014, 43(6): 575-581.DOI: 10.13485/j.cnki.11-2089.2014.0091.

[18] TIKHONOV A N.Regularization of Incorrectly Posed Problems[J].Soviet Mathematics Doklady, 1963, 4(1): 1624-1627.

[19] DAWSON J, TREGONING P.Uncertainty Analysis of Earthquake Source Parameters Determined from InSAR: A Simulation Study[J].Journal of Geophysical Research: Solid Earth, 2007, 112(B9): B09406.

[20] WALTERS R J, ELLIOTT J R, D’AGOSTINO N, et al.The 2009 L’Aquila Earthquake (Central Italy): A Source Mechanism and Implications for Seismic Hazard[J].Geophysical Research Letters, 2009, 36(17): L17312.

[21] ANZIDEI M,BOSCHI E,CANNELLI V, et al.Coseismic Deformation of the Destructive April 6, 2009 L’Aquila Earthquake (Central Italy) from GPS Data[J].Geophysical Research Letters, 2009, 36(17): L17307.

[22] CHELONI D, D’AGOSTINO N, D’ANASTASIO E, et al.Coseismic and Initial Post-seismic Slip of the 2009Mw6.3 L’Aquila Earthquake, Italy, from GPS Measurements[J].Geophysical Journal International, 2010, 181(3): 1539-1546.

[23] 馮萬鵬, 李振洪, 李春來.利用InSAR確定2009年4月6日Mw6.3拉奎拉(Italy)地震最優(yōu)震源模型[J].地球物理學進展, 2010, 25(5): 1550-1559.FENG Wanpeng, LI Zhenhong, LI Chunlai.Optimal Source Parameters of the 6 April 2009Mw6.3 L’Aquila, Italy Earthquake from InSAR Observations[J].Progress in Geophysics, 2010, 25(5): 1550-1559.

[24] 馮萬鵬, 李振洪.InSAR資料約束下震源參數(shù)的PSO混合算法反演策略[J].地球物理學進展, 2010, 25(4): 1189-1196.FENG Wanpeng, LI Zhenhong.A Novel Hybrid PSO/Complex Algorithm for Determining Earthquake Source Parameters Using InSAR Data[J].Progress in Geophysics, 2010, 25(4): 1189-1196.

[25] 溫揚茂, 許才軍, 劉洋, 等.升降軌InSAR數(shù)據(jù)約束下的2007年阿里地震反演分析[J].測繪學報, 2015, 44(6): 649-654.DOI: 10.11947/j.AGCS.2015.20140349.WEN Yangmao, XU Caijun, LIU Yang, et al.The 2007 Ali Earthquake Inversion from Ascending and Descending InSAR Observations[J].Acta Geodaetica et Cartographica Sinica, 2015, 44(6): 649-654.DOI: 10.11947/j.AGCS.2015.20140349.

[26] 曾文憲.系數(shù)矩陣誤差對EIV模型平差結(jié)果的影響研究[D].武漢: 武漢大學, 2013.ZENG Wenxian.Effect of The Random Design Matrix on Adjustment of an EIV Model and Its Reliability Theory[D].Wuhan: Wuhan University, 2013.

[27] ATZORI S, HUNSTAD I, CHINI M, et al.Finite Fault Inversion of DInSAR Coseismic Displacement of the 2009 L’Aquila Earthquake (Central Italy)[J].Geophysical Research Letters, 2009, 36(15): L15305.

[28] BONCIO P, LAVECCHIA G, PACE B.Defining a Model of 3D Seismogenic Sources for Seismic Hazard Assessment Applications: The Case of Central Apennines (Italy)[J].Journal of Seismology, 2004, 8(3): 407-425.

(責任編輯：陳品馨)

Total Least Squares Method Inversion for Coseismic Slip Distribution

WANG Leyang1,3，LI Haiyan1，WEN Yangmao2，XU Caijun2

1.Faculty of Geomatics, East China Institute of Technology, Nanchang 330013, China; 2.School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China; 3.Key Laboratory of Watershed Ecology and Geographical Environment Monitoring, NASG, Nanchang 330013, China

The coefficient matrix (Green matrix) is composed of surface point offset caused by unit slip of sub-fault patches.The elements of the coefficient matrix are related to the location, geometry of rupture surface, assumption of model and other factors.In this paper, we attempted to consider the Green’function matrix (coefficient matrix) errors in order to compensate for the effects of above-mentioned factors to some extent.The total least squares (TLS) method, which both errors of coefficient matrix and observation vector are considered, is proposed for fault slip inversion.So we dealt with the errors in both of coefficient matrix and observation at same time.And by analysis of the relations between observation vector and coefficient matrix elements, we obtained the covariance matrix of coefficient matrix elements and observation vector.Considering the coefficient matrix was ill-posed, we used the second-order Laplace smoothing to constrain the slip parameters each other, then we used the regularized total least squares method to estimate slip distribution.the total least squares (TLS)slip inversion method was applied to simulate oblique fault event and Mw6.3 earthquake occurred in L’Aquila (central Italy) on April 6, 2009, respectively.To L’Aquila earthquake, the results by total least squares method indicate that the inverted geodetic moment is 3.63×1018N·m (Mw6.34).With a maximum slip of 0.95 m, and a average rake of -96.4°, the main slip occurred at depth of 4 km-15 km.The difference of slip distribution solutions between total least squares and least squares method is less than 10-4order.

coseismic slip distribution; total least squares inversion; regularization; L’Aquila earthquake

National Natural Science Foundation of China (Nos.41664001; 41204003; 41574002; 41431069); Support Program for Outstanding Youth Talents in Jiangxi Province (No.20162BCB23050); National Department Public Benefit Research Foundation (Surveying, Mapping and Geoinformation) (No.201512026); National Key Research and Development Program(No.2016YFB0501405); Science and Technology Project of the Education Department of Jiangxi Province (No.GJJ150595)

WANG Leyang (1983—), male, PhD, associate professor, majors in geodetic inversion and geodetic data processing.

王樂洋，李海燕，溫揚茂，等.地震同震滑動分布反演的總體最小二乘方法[J].測繪學報，2017,46(3)：307-315.

10.11947/j.AGCS.2017.20160212.

WANG Leyang，LI Haiyan，WEN Yangmao，et al.Total Least Squares Method Inversion for Coseismic Slip Distribution[J].Acta Geodaetica et Cartographica Sinica,2017,46(3):307-315.DOI:10.11947/j.AGCS.2017.20160212.

P227

A

1001-1595(2017)03-0307-09

國家自然科學基金(41664001; 41204003; 41574002; 41431069)；江西省杰出青年人才資助計劃項目(20162BCB23050)；測繪地理信息公益性行業(yè)科研專項(201512026)；國家重點研發(fā)計劃(2016YFB0501405)；江西省教育廳科技項目(GJJ150595)

2016-05-19

修回日期：2016-10-25

王樂洋(1983—)，男，博士，副教授，主要研究方向為大地測量反演及大地測量數(shù)據(jù)處理。

E-mail：wleyang@163.com