李海燕 王樂洋,2

（1.東華理工大學(xué)測(cè)繪工程學(xué)院 江西南昌 330013；2.流域生態(tài)與地理環(huán)境監(jiān)測(cè)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室 江西南昌 330013）

總體最小二乘方法在地震位錯(cuò)模型中的應(yīng)用探討

李海燕1王樂洋1,2

（1.東華理工大學(xué)測(cè)繪工程學(xué)院 江西南昌 330013；2.流域生態(tài)與地理環(huán)境監(jiān)測(cè)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室 江西南昌 330013）

分析了目前基于位錯(cuò)模型反演地震斷層參數(shù)的研究現(xiàn)狀，嘗試?yán)每傮w最小二乘方法反演地震斷層問(wèn)題，通過(guò)探討總體最小二乘方法在非線性反演震源參數(shù)和線性反演震源滑動(dòng)分布兩個(gè)方面的應(yīng)用，進(jìn)一步促進(jìn)位錯(cuò)模型下地震斷層參數(shù)反演的理論和應(yīng)用的研究。

位錯(cuò)模型；總體最小二乘；震源參數(shù)；滑動(dòng)分布

1 引言

地震是地球上最常見的自然現(xiàn)象之一，也是最嚴(yán)重的自然災(zāi)害之一。研究地震的手段通常有地震學(xué)、地質(zhì)學(xué)、野外地質(zhì)調(diào)查和大地測(cè)量方法等。地震學(xué)方法是通過(guò)地震波資料反演震源參數(shù)，由于地震波輻射場(chǎng)與斷層的滑動(dòng)量、斷層破裂的速度及震源時(shí)間函數(shù)有關(guān)［1］，因此，通過(guò)地震波資料確定震源參數(shù)及地震震級(jí)通常并不能準(zhǔn)確反映地震發(fā)生所釋放的真實(shí)能量。地質(zhì)學(xué)主要通過(guò)地質(zhì)構(gòu)造來(lái)分析斷層的演化過(guò)程，這類方法涉及的時(shí)間尺度一般都是幾十萬(wàn)、幾百萬(wàn)甚至上億年，無(wú)法精確反映短期的地殼活動(dòng)。野外地質(zhì)調(diào)查可以獲得斷層破裂延伸至地表時(shí)的斷層走向方向及斷層近似長(zhǎng)度，而對(duì)于斷層破裂深度或破裂未能延伸至地表時(shí)，該方法將無(wú)法確定震源參數(shù)。大地測(cè)量方法能夠有效彌補(bǔ)前述幾種方法的不足，特別是近二三十年，隨著空間大地測(cè)量技術(shù)（如GPS、INSAR)的發(fā)展，使得獲得精確的同震前后地表位置的變化成為現(xiàn)實(shí)，通過(guò)彈性位錯(cuò)理論，建立了地表形變數(shù)據(jù)和震源幾何學(xué)、運(yùn)動(dòng)學(xué)的聯(lián)系［2-4］，通過(guò)擬合地表形變數(shù)據(jù)，可以比較精確地獲得斷層震源參數(shù)。

本文首先介紹同震矩形位錯(cuò)模型，綜述已有的基于位錯(cuò)模型的斷層參數(shù)研究進(jìn)展，針對(duì)目前位錯(cuò)模型斷層參數(shù)反演大都沒有考慮模型系數(shù)矩陣誤差的問(wèn)題，探討了基于總體最小二乘的地震震源參數(shù)和滑動(dòng)分布反演方法。

2 同震位錯(cuò)模型

斷層位錯(cuò)模型是地球物理模型中描述斷層運(yùn)動(dòng)的主要模型，有矩形位錯(cuò)和三角位錯(cuò)等。位于均勻彈性半空間的矩形位錯(cuò)模型是最常用的位錯(cuò)模型。1958年Steketee［5］首次將位錯(cuò)理論引入地震學(xué)中用來(lái)描述斷層運(yùn)動(dòng)。1985年Okada系統(tǒng)總結(jié)了已有的研究成果，進(jìn)一步推導(dǎo)了均勻彈性半空間介質(zhì)的矩形位錯(cuò)理論，給出了斷層矩形位錯(cuò)引起的地表形變關(guān)系式。位錯(cuò)模型已成為利用大地測(cè)量觀測(cè)數(shù)據(jù)反演地震參數(shù)和震源機(jī)制非常有效的手段。

如圖1，在o-xyz為斷層坐標(biāo)系，斷層面的長(zhǎng)度為L(zhǎng)，寬度為W，斷層深度為d，傾角為δ，斷層走向α，U1，U2，U3分別是斷層面沿走向、傾向和張裂方向的滑動(dòng)量，假設(shè)斷層坐標(biāo)系原點(diǎn)在地面坐標(biāo)系中的位置為（x0，y0），斷層參數(shù)X＝［x0，y0，L，W，d，δ，α，U1，U2，U3］T其中，U1，U2，U3為的滑動(dòng)參數(shù)，其余為幾何參數(shù)。則斷層引起的地表形變量f為

圖1 均勻彈性半空間矩形位錯(cuò)示意圖［6］

Okada矩形位錯(cuò)公式中，斷層幾何參數(shù)與地表形變量是非線性關(guān)系，而斷層的滑動(dòng)參數(shù)與地表形變形變量之間是線性的。所以斷層參數(shù)的反演包含非線性參數(shù)反演和線性滑動(dòng)分布反演。非線性反演就是在單一均勻滑動(dòng)模型下，確定斷層的幾何參數(shù)及均勻滑動(dòng)量；線性滑動(dòng)分布反演是在確定斷層幾何參數(shù)基礎(chǔ)上，得到破裂面上更精細(xì)的滑動(dòng)量分布。通常將斷層面沿?cái)鄬幼呦蚝蛢A向方向分別適當(dāng)延長(zhǎng)，將延長(zhǎng)后的斷層面離散剖分為M個(gè)小斷層片，然后通過(guò)Okada位錯(cuò)公式計(jì)算單位走滑和單位傾滑（1m)在每個(gè)地表觀測(cè)點(diǎn)上產(chǎn)生的地表位移，由這些位移量構(gòu)造N×2M階的格林函數(shù)H。N為地表觀測(cè)點(diǎn)個(gè)數(shù)。得到地表位移與震源參數(shù)的關(guān)系：

其中，d表示地表點(diǎn)形變觀測(cè)值；m （2M×1階）表示斷層分別沿走向和傾向方向的滑動(dòng)量；Hs和Hd分別表示其中某一小斷層走滑和傾滑引起的地表位移。

3 矩形位錯(cuò)斷層反演研究現(xiàn)狀

3.1 非線性斷層參數(shù)反演

由于地表形變值與斷層參數(shù)是非常復(fù)雜的非線性關(guān)系，很難給出其線性化形式，同時(shí)線性方法對(duì)模型參數(shù)的初始值依賴性很大，所以更多的是采用非線性反演方法求解斷層參數(shù)，傳統(tǒng)的非線性反演方法有牛頓法、梯度法，還有啟發(fā)式的搜索算法，如模擬退火法、遺傳算法等［7］。溫?fù)P茂［6］通過(guò)獲得1997年瑪尼Mw7.5級(jí)地震地表點(diǎn)的InSAR形變數(shù)據(jù)，運(yùn)用遺傳算法反演了的該地震的斷層參數(shù)，并用蒙特卡洛誤差傳遞方法評(píng)定反演參數(shù)的精度；得到均勻滑動(dòng)模型下反演的震級(jí)為Mw7.51級(jí)（地震矩Nm），與地震學(xué)方法得到的結(jié)果相近。李爽［8］對(duì)基于位錯(cuò)模型的多類數(shù)據(jù)聯(lián)合反演問(wèn)題進(jìn)行了研究，采用區(qū)間模擬退火算法研究了地震定位問(wèn)題，并利用重力數(shù)據(jù)和GPS數(shù)據(jù)，聯(lián)合反演了川西鮮水河斷裂帶的斷層參數(shù)。有學(xué)者也研究了線性化迭代反演，由于線性化得到的系數(shù)矩陣容易產(chǎn)生結(jié)構(gòu)性擾動(dòng)和受線性化誤差的影響，故Bifulco等［9］利用最小范數(shù)（LN）算法來(lái)求解參數(shù)反演這一結(jié)構(gòu)總體范數(shù) （structured total least norm，STLN）問(wèn)題，同時(shí)為了抵抗觀測(cè)值粗差對(duì)結(jié)果的影響，采用具有抗差能力的1范數(shù)進(jìn)行參數(shù)求解。并通過(guò)模擬數(shù)據(jù)驗(yàn)證了該算法在斷層參數(shù)反演中的有效性。王樂洋等［10］在Bifulco等研究的基礎(chǔ)上，利用InSAR數(shù)據(jù)和LN算法反演了2008年青海大柴旦Mw6.3級(jí)地震的斷層參數(shù)。

3.2 斷層滑動(dòng)分布反演

已知斷層的幾何參數(shù)，那么斷層面上的滑動(dòng)值與觀測(cè)得到的形變場(chǎng)之間則是線性關(guān)系，可以通過(guò)最小二乘方法來(lái)求解斷層滑動(dòng)分布問(wèn)題。李志才等［11］基于GPS觀測(cè)數(shù)據(jù)反演了汶川地震斷層參數(shù)，利用彈性半空間均勻位錯(cuò)理論，采用有界變量最小二乘法反演了汶川地震斷層滑動(dòng)分布。Elliottetal［12］聯(lián)合升降軌Envisat雷達(dá)影像，獲得了2008年大柴旦Mw6.3級(jí)地震的InSAR數(shù)據(jù)同震形變場(chǎng)，采用約束最小二乘法反演了斷層面的滑動(dòng)分布；與Elliottet al對(duì)斷層面進(jìn)行均勻剖分不同的是，溫?fù)P茂等［13］利用自動(dòng)剖分技術(shù)得到斷層面的非均勻子塊，并通過(guò)最小二乘法反演了斷層的精細(xì)滑動(dòng)分布；最小二乘作為最普遍、最常用的數(shù)據(jù)處理方法，已在全球范圍內(nèi)的中強(qiáng)地震滑動(dòng)分布反演中得到了廣泛應(yīng)用。

3.3 斷層參數(shù)及斷層滑動(dòng)分布反演算法的不足

從已有的斷層參數(shù)反演研究可以發(fā)現(xiàn)，無(wú)論是采用非線性方法還是線性迭代方法，都只是考慮了地表形變觀測(cè)值即觀測(cè)向量的誤差，上述的利用LN算法解決了參數(shù)反演的結(jié)構(gòu)總體最小范數(shù)問(wèn)題，在觀測(cè)值含有粗差的情況算法具有較好抵抗粗差的能力，但是當(dāng)觀測(cè)向量不存在明顯粗差值時(shí)，反演效果并不理想。通過(guò)分析可知，由于Okada矩形位錯(cuò)模型本身存在某些假設(shè)條件，這就使得模型不可避免地存在誤差，進(jìn)一步認(rèn)為這些誤差是模型參數(shù)系數(shù)矩陣的誤差，此時(shí)，需要考慮系數(shù)矩陣誤差對(duì)斷層參數(shù)反演的影響，最小二乘方法無(wú)法處理系數(shù)矩陣和觀測(cè)向量同時(shí)含有誤差的情況，這時(shí)可以把最小二乘的擴(kuò)展算法——總體最小二乘 （total least squares，TLS）應(yīng)用于斷層參數(shù)反演。

滑動(dòng)分布反演中，格林函數(shù) 即模型系數(shù)矩陣是通過(guò)斷層幾何參數(shù)計(jì)算的單位位錯(cuò)在地表引起的形變，地球介質(zhì)、形狀、斷層幾何形狀的復(fù)雜性、斷層剖分方式對(duì)系數(shù)矩陣構(gòu)成有很大影響，不同的斷層幾何結(jié)構(gòu)對(duì)應(yīng)不同的系數(shù)矩陣，分析式（2），待求參數(shù)的系數(shù)矩陣 存在誤差最終將會(huì)影響滑動(dòng)模型參數(shù)的反演結(jié)果。而已有滑動(dòng)分布反演研究，都是基于最小二乘或約束最小二乘準(zhǔn)則，顯然，如果把總體最小二乘方法應(yīng)用在滑動(dòng)分布反演中，將會(huì)得到更加更加合理的滑動(dòng)參數(shù)，更好地?cái)M合地表形變位移。

4 位錯(cuò)模型參數(shù)的總體最小二乘反演法

最小二乘法僅在觀測(cè)值存在誤差時(shí)可以得到無(wú)偏最優(yōu)估計(jì)值，而總體最小二乘是20世紀(jì)80年代發(fā)展起來(lái)的同時(shí)顧及系數(shù)矩陣和觀測(cè)值誤差的測(cè)繪數(shù)據(jù)處理方法，其理論與應(yīng)用是目前國(guó)內(nèi)外研究的熱點(diǎn)問(wèn)題［14］。

總體最小二乘反演的函數(shù)模型為

式中，L為觀測(cè)值，eL表示觀測(cè)誤差，系數(shù)矩陣A及誤差矩陣EA，x為待估計(jì)參數(shù)。

當(dāng)觀測(cè)值和系數(shù)矩陣元素等精度、不相關(guān)時(shí)，則

總體最小二乘的平差準(zhǔn)則為

構(gòu)造拉格朗日函數(shù)求條件極值，得到如下的公式［15］

得到可待估計(jì)參數(shù)的解

當(dāng)觀測(cè)值與系數(shù)矩陣元素不等精度且相關(guān)時(shí)［14］，則有

得到待估參數(shù)的加權(quán)總體最小二乘解［14］

因?yàn)榈仁剑?）及等式（10）的右端都與參數(shù)x有關(guān)，在求解時(shí)需要迭代計(jì)算。主要步驟是把基于最小二乘的參數(shù)解作為初始值，反復(fù)計(jì)算待估參數(shù)，直至所求參數(shù)滿足迭代終止條件結(jié)束。

5 結(jié)束語(yǔ)

地震斷層參數(shù)反演、地震滑動(dòng)分布反演等大地測(cè)量反演問(wèn)題，通常是借助最小二乘方法完成。而模型系數(shù)矩陣和觀測(cè)向量值同時(shí)存在誤差時(shí)，研究總體最小二乘方法在地震斷層參數(shù)反演、地震斷層滑動(dòng)分布等問(wèn)題的應(yīng)用，具有十分重要的意義。本文在綜述位錯(cuò)模型下地震斷層參數(shù)反演和斷層滑動(dòng)分布的研究現(xiàn)狀的基礎(chǔ)上，分析了位錯(cuò)模型中觀測(cè)值和模型的誤差，我們認(rèn)為對(duì)于斷層參數(shù)反演和滑動(dòng)分布反演問(wèn)題，都應(yīng)同時(shí)考慮系數(shù)矩陣和觀測(cè)值中的誤差。因此，探討了將總體最小二乘反演法應(yīng)用于在地震位錯(cuò)模型參數(shù)反演，并給出了總體最小二乘反演法的一般算法。借助大地測(cè)量技術(shù)，將總體最小二乘方法應(yīng)用于地震學(xué)反演問(wèn)題，將是大地測(cè)量學(xué)科與地震學(xué)科交叉融合的新的研究方向。

［1］劉洋.顧及模型誤差的震源參數(shù)InSAR反演［D］.博士學(xué)位論文，武漢大學(xué)，2012.

［2］OKADA Y.Surface Deformation due to Shear and Tensile Faults in a Half-space［J］.Bulletin of the Seismologica Socie-ty of America，1985，75(4):1135-1154.

［3］OKADA Y.Internal Deformation duo to Shear and Tensile Faults in a Half-space［J］.Bulletin of the Seismological Society of America，1992，82(2):1018-1140.

［4］Pollitz F.F.Coseism ic Deformation from Earthquake Faulting on A Layered Spherical Earth，Geophysical Journal International，1996，125(1):1-14.

［5］Steketee，J.A.On Volterras Dislocation in A Semi-infinite Elastic Medium，Elastic Medium Can.J.Phys.，1958，36: 192-205.

［6］溫?fù)P茂.利用InSAR資料研究若干強(qiáng)震的同震和震后形變［D］.博士學(xué)位論文，武漢大學(xué)，2009.

［7］王家映.地球物理反演理論［M］.高等教育出版社，2002.

［8］李爽.大地測(cè)量聯(lián)合反演的模式及算法研究［D］.博士學(xué)位論文，武漢大學(xué)，2005.

［9］BIFULCO I，RAICONIG，SCARPA R.Computer Algebra Software for LeastSquaresand Total LeastNorm Inversion of Geophysical Models［J］.Computers&Geosciences，2009，35 (7):1427-1438.

［10］王樂洋，許才軍，溫?fù)P茂.利用STLN和InSAR數(shù)據(jù)反演2008年青海大柴旦Mw6.3級(jí)地震斷層參數(shù)［J］.測(cè)繪學(xué)報(bào)，2013，42(2):168-176.

［11］李志才，張鵬，金雙根等.基于 GPS觀測(cè)數(shù)據(jù)的汶川地震斷層形變反演分析［J］.測(cè)繪學(xué)報(bào)，2009，38(2):108-113

［12］ELLIOTT，J.，B.PARSONS，J.JACKSON，et al.Depth Segmentation of the Seismogenic Continental Crust:the 2008 and 2009 Qaidam Earthquakes［J］.Geophysical Research Letters，2011，38(6):L06305.

［13］溫?fù)P茂，許才軍，劉洋等.利用斷層自動(dòng)剖分計(jì)技術(shù)的2008年青海大柴旦Mw6.3級(jí)地震InSAR［J］.2012，37 (4):458-462.

［14］王樂洋.基于總體最小二乘的大地測(cè)量反演理論及應(yīng)用研究［J］.測(cè)繪學(xué)報(bào)，2012，41(4):629.

［15］魯鐵定.總體最小二乘平差理論及其在測(cè)繪數(shù)據(jù)處理中的應(yīng)用［D］.博士學(xué)位論文，武漢大學(xué)，2010.