地震同震滑動(dòng)分布反演平滑因子的確定

王樂洋，趙 雄

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

在地震機(jī)制研究領(lǐng)域中，利用大地測量資料研究地震機(jī)制是地學(xué)領(lǐng)域的前沿研究熱點(diǎn)[1]。大地測量數(shù)據(jù)(如GPS、InSAR數(shù)據(jù)等)因觀測范圍廣、精度高，反演出的滑動(dòng)分布結(jié)果更為精細(xì)、全面，近年來常被國內(nèi)外學(xué)者用于進(jìn)行地震同震滑動(dòng)分布反演研究[2-8]。在震源機(jī)制研究中，位錯(cuò)模型是反演震源參數(shù)的有效模型[9-11]，一般采用彈性半空間矩形位錯(cuò)模型進(jìn)行地震同震滑動(dòng)分布反演[12-13]。

在地震同震滑動(dòng)分布反演過程中，地表形變位移與滑動(dòng)分布參數(shù)之間為線性關(guān)系?；瑒?dòng)分布參數(shù)與地表形變位移的線性模型依賴于格林函數(shù)矩陣的構(gòu)造，而格林函數(shù)矩陣元素通常是確定斷層的幾何參數(shù)后，將斷層破裂面進(jìn)一步離散成若干斷層子塊，由每個(gè)子塊的單位走滑與傾滑引起的地表形變位移構(gòu)成。另外，格林函數(shù)矩陣在地震滑動(dòng)分布反演過程中通常是病態(tài)的，一般采用Tikhonov正則化方法[14]解決。

在正則化過程中，平滑因子的確定是關(guān)鍵問題之一，不同的平滑因子會(huì)使反演結(jié)果大不相同。目前主要方法有廣義交叉核實(shí)法(generalized cross validation,GCV)[15]、嶺跡法[16-17]、L曲線法[18-20]。GCV法理論嚴(yán)密，但有時(shí)GCV方法函數(shù)繪制的曲線過于平緩，很難定位到最優(yōu)平滑因子[21]；嶺跡法計(jì)算簡單，但是在確定平滑因子的時(shí)候具有一定的主觀性；在地震同震滑動(dòng)分布反演中，目前最常用的平滑因子確定方法為L曲線法，但是L曲線法在確定平滑因子過程中具有以下幾個(gè)缺點(diǎn)：①無法從L曲線圖中根據(jù)橫縱坐標(biāo)軸直觀讀取平滑因子的大?。虎贚曲線法過于依賴數(shù)據(jù)擬合精度[22]，求解過程可能不收斂[23]；③計(jì)算時(shí)間過長。

針對以上問題，本文提出一種確定地震同震滑動(dòng)分布反演平滑因子的新方法——折中相交曲線法(eclectic intersection curve method，EI)(簡稱EI曲線法)。通過模擬試驗(yàn)以及拉奎拉與臺灣美濃實(shí)際震例反演驗(yàn)證EI曲線法確定平滑因子的優(yōu)勢與可行性。

1 EI曲線法確定平滑因子的基本原理

1.1 同震滑動(dòng)分布基本方程

在同震滑動(dòng)分布反演中，地表同震位移與斷層滑動(dòng)量之間為線性關(guān)系，如式(1)所示

d=Gm+ε

(1)

式中，d表示地表形變數(shù)據(jù)；G表示格林函數(shù)矩陣；m為滑動(dòng)參數(shù)；ε為隨機(jī)觀測誤差。

在同震滑動(dòng)分布反演中，通常將斷層面劃分成n個(gè)矩形單元，基于矩形位錯(cuò)模型計(jì)算矩形單元的單位走滑量和單位傾滑量引起地表k個(gè)觀測點(diǎn)的形變值，構(gòu)成了k×2n個(gè)格林函數(shù)矩陣G。將斷層面離散化處理，直接用式(1)求解會(huì)出現(xiàn)病態(tài)問題，從而使滑動(dòng)分布解不穩(wěn)定，因?yàn)榇藭r(shí)格林矩陣G是嚴(yán)重秩虧的。為了避免這種情況，通常對斷層面施加一定的平滑約束條件，一方面增加虛擬觀測數(shù)據(jù)解決系數(shù)矩陣G的秩虧問題，另一方面約束了相鄰斷層單元間出現(xiàn)大的梯度變化。

本文采用目前較成功用于斷層面平滑約束的拉普拉斯二階平滑矩陣對斷層單元進(jìn)行約束[24]。根據(jù)Tikhonov正則化原理，相應(yīng)的準(zhǔn)則為

(2)

由式(1)、式(2)可得地震滑動(dòng)分布解為

m=(GTP1G+αTTP2T)-1GTP1d

(3)

1.2 L曲線法的基本原理

圖1 L曲線圖Fig.1 L curve

1.3 EI曲線法基本原理

針對以上問題，EI曲線法首先利用赫爾默特方差分量估計(jì)法[1,25]統(tǒng)一反演過程中各類數(shù)據(jù)單位權(quán)方差，在統(tǒng)一單位權(quán)方差過程中需要加入權(quán)修正因子[26]a、b進(jìn)行修正，使σ01=σ02，此時(shí)相應(yīng)的目標(biāo)函數(shù)可以表達(dá)如下

(4)

式中，a、b為實(shí)際觀測數(shù)據(jù)與虛擬觀測數(shù)據(jù)內(nèi)部權(quán)陣P1、P2的權(quán)修正因子。

圖2 EI曲線圖Fig.2 EI curve

綜上，利用EI曲線法確定地震滑動(dòng)分布平滑因子步驟可總結(jié)如下：

(1) 利用方差分量估計(jì)法統(tǒng)一實(shí)際觀測數(shù)據(jù)與虛擬觀測數(shù)據(jù)單位權(quán)方差：即使σ01=σ02，并求得a、b。

2 模擬試驗(yàn)

為了驗(yàn)證本文方法的可行性，本文作了以下模擬試驗(yàn)。模擬試驗(yàn)中，模擬斷層面幾何參數(shù)如下：斷層面幾何中心位置X=0 km，Y=0 km，斷層頂深0 km，長度30 km，寬度30 km，走向角60°，傾角45°，滑動(dòng)角45°。模擬GPS 3方向觀測數(shù)據(jù)如圖3所示，根據(jù)點(diǎn)位與震中距離設(shè)置形變點(diǎn)密度，并給形變點(diǎn)施加觀測誤差N(0，32mm2)。

在模擬試驗(yàn)中，本文分別利用L曲線法、EI曲線法確定平滑因子，模擬試驗(yàn)中平滑因子α均以0.01為初值、0.001為步長、1為終值繪制L曲線圖、EI曲線圖(其中，L曲線法確定平滑因子大小為0.179，EI曲線法確定的平滑因子大小為0.095。在利用EI曲線法確定平滑因子過程中，利用方差分量估計(jì)法確定實(shí)際觀測數(shù)據(jù)與虛擬觀測數(shù)據(jù)內(nèi)部權(quán)陣修正因子大小分別為a=1，b=0.003 7)，并將兩種方法確定的平滑因子進(jìn)行地震滑動(dòng)分布反演。利用L曲線法、EI曲線法確定平滑因子如圖4(a)、(d)所示；兩種方法反演地震滑動(dòng)分布結(jié)果以及兩種方法反演地震滑動(dòng)分布與模擬值的差值分布如圖4(b)、(e)、(c)、(f)所示。

模擬試驗(yàn)預(yù)設(shè)參數(shù)以及利用L曲線法、EI曲線法反演滑動(dòng)分布各參數(shù)結(jié)果見表1。兩種方法確定平滑因子反演地震滑動(dòng)分布部分參數(shù)結(jié)果與模擬值的差異，以及EI曲線法反演地震滑動(dòng)分布部分參數(shù)結(jié)果，較L曲線法精度改進(jìn)程度見表2。

表1 模擬試驗(yàn)反演地震滑動(dòng)分布結(jié)果

表2模擬試驗(yàn)兩種方法反演各參數(shù)結(jié)果與模擬值的差異

Tab.2Comparisontheresultsofeachparametersandsimulationvaluesoftwomethodsinsimulationexperiments

與模擬值的差值最大滑動(dòng)量/m平均滑動(dòng)量/m平均滑動(dòng)角/(°)L曲線法0.16720.00640.9012EI曲線法0.09990.00340.8826L-EI?5.23%0.95%0.04%

注：L-EI*表示利用EI曲線法確定平滑因子反演地震滑動(dòng)分布各參數(shù)結(jié)果較L曲線法改進(jìn)程度。

表1給出了利用L曲線法和EI曲線法確定平滑因子反演地震同震滑動(dòng)分布各參數(shù)的結(jié)果及兩種方法確定平滑因子所用時(shí)間；表2給出了兩種方法反演地震同震滑動(dòng)分布各參數(shù)結(jié)果與模擬值的差異。在各參數(shù)反演精度上，從表1、表2可以看出，利用EI曲線法確定的平滑因子反演地震滑動(dòng)分布結(jié)果在最大滑動(dòng)量、平均滑動(dòng)量、矩震級等參數(shù)反演結(jié)果均要優(yōu)于L曲線法反演結(jié)果；在計(jì)算時(shí)間上，在相同的初值、步長、終值情況下，EI曲線法比L曲線法在計(jì)算時(shí)間上減少了107.784 3 s，大大地提高了計(jì)算效率；在均方根誤差上，EI曲線法較L曲線法減小0.1 mm，進(jìn)一步說明EI曲線法確定平滑因子反演地震滑動(dòng)分布結(jié)果精度要高于L曲線法。

3 實(shí)際震例反演

3.1 拉奎拉地震滑動(dòng)分布反演

2009年4月6日，意大利中部拉奎拉地區(qū)發(fā)生Mw6.3級地震。據(jù)美國地質(zhì)調(diào)查局(USGS)網(wǎng)站發(fā)布，地震震中位置為(13.334°E，42.334°N)，震源深度為8.8 km。地震造成1000余人的傷亡，大量的房屋被損毀或破壞[27]。本文從文獻(xiàn)[6，28]中獲取拉奎拉地震InSAR降軌數(shù)據(jù)位移1254個(gè)，升軌位移數(shù)據(jù)1282個(gè)，并利用文獻(xiàn)[29]中通過非線性反演方法反演得到的拉奎拉地震斷層的走向角、傾角及斷層中心位置等參數(shù)來固定斷層破裂面。在此基礎(chǔ)上將斷層破裂面沿?cái)鄬幼呦?、傾向把斷層長度、寬度擴(kuò)展至30 km、27 km，并將破裂面延伸至地表，將斷層面均勻剖分成1.5×1.5 km大小的矩形單元，共得到了360個(gè)矩形斷層單元。分別采用EI曲線法與L曲線法確定平滑因子進(jìn)行地震滑動(dòng)分布反演。

在拉奎拉實(shí)際震例中，平滑因子α均以0.05為初值、0.001為步長、1為終值繪制L曲線圖、EI曲線圖(其中，L曲線法確定平滑因子大小為0.588與文獻(xiàn)[28]中用L曲線法確定的平滑因子大小0.6，兩者大小相近；EI曲線法確定的平滑因子大小為0.159。在利用EI曲線法確定平滑因子過程中，利用方差分量估計(jì)法確定實(shí)際觀測數(shù)據(jù)與虛擬觀測數(shù)據(jù)內(nèi)部權(quán)陣修正因子，大小分別為a=1，b=0.002 4)。利用EI曲線法和L曲線法確定平滑因子如圖5(a)、(c)所示；兩種方法反演拉奎拉地震成果及兩種方法反演滑動(dòng)分布結(jié)果差值分布如圖5(b)、(d)、(e)所示。

利用L曲線法、EI曲線法確定的平滑因子反演拉奎拉地震滑動(dòng)分布各參數(shù)結(jié)果見表3；近年來國內(nèi)外部分其他學(xué)者反演拉奎拉地震滑動(dòng)分布各參數(shù)結(jié)果見表4。

表3 用L曲線法、EI曲線法確定平滑因子反演拉奎拉地震滑動(dòng)分布結(jié)果

表4 拉奎拉地震斷層滑動(dòng)參數(shù)

由圖5(b)、(d)可知，利用L曲線法、EI曲線法確定的平滑因子進(jìn)行拉奎拉地震反演時(shí)，其反演結(jié)果均顯示拉奎拉同震滑動(dòng)分布區(qū)域主要分布在地下4～15 km，該結(jié)果同文獻(xiàn)[28]研究成果一致；對比表3、表4可以看出，本文兩種方法確定的平滑因子反演拉奎拉地震滑動(dòng)分布結(jié)果均在其他學(xué)者研究范圍內(nèi)。對于最大滑動(dòng)量，利用L曲線法和EI曲線法反演最大滑動(dòng)量大小分別為1.01、1.02 m，文獻(xiàn)[6]聯(lián)合了不同波長、不同入射角的Enviast和ALOS衛(wèi)星差分干涉數(shù)據(jù)，并聯(lián)立了震區(qū)部分GPS三維形變數(shù)據(jù)，通過彈性三角位錯(cuò)模型，反演得到滑動(dòng)分布的最大滑動(dòng)量為1.07 m，對于其最大滑動(dòng)量略大于本文反演結(jié)果，這可能與約束地表的數(shù)據(jù)源有關(guān)。文獻(xiàn)[28]利用最小二乘和總體最小二乘反演最大滑動(dòng)量均為0.95 m，略小于本文兩種方法反演結(jié)果；對于平均滑動(dòng)角，本文兩種方法反演的結(jié)果分別為-101.3°、-103.9°，該研究結(jié)果均在其他學(xué)者研究范圍內(nèi)，并且與文獻(xiàn)[6]反演結(jié)果相近；對于矩震級參數(shù)，本文兩種方法反演結(jié)果相近，且在其他學(xué)者研究范圍內(nèi)。另外，文獻(xiàn)[27]聯(lián)合ASAR數(shù)據(jù)反演得到的斷層滑動(dòng)量為0.66 m，給出的矩張量為2.80×1018N·m；文獻(xiàn)[31]利用拉奎拉震中附近的GPS觀測站數(shù)據(jù)，研究了拉奎拉地震的同震滑動(dòng)分布和震后的余震滑動(dòng)分布，得到的最大滑動(dòng)量為1.1 m，矩張量為3.90×1018N·m。

3.2 臺灣美濃地震滑動(dòng)分布反演

2016年2月6日3時(shí)57分(北京時(shí)間)我國臺灣省高雄市美濃區(qū)發(fā)生了Mw6.4(global centroid moment tensor,GCMT)地震(下稱美濃地震)。此次地震給臺灣南部地區(qū)帶來巨大的人員傷亡及財(cái)產(chǎn)損失。本文采用文獻(xiàn)[32]獲取美濃GPS站位移數(shù)據(jù)225個(gè)參與計(jì)算美濃地震滑動(dòng)分布反演。利用文獻(xiàn)[32]中均勻滑動(dòng)分布反演的斷層走向、傾向、斷層深度等參數(shù)，作為本文通過多峰值顆粒群算法[29]進(jìn)行非線性斷層參數(shù)反演的初始參數(shù)大小。通過非線性反演最終確定適合本文數(shù)據(jù)的美濃地震斷層參數(shù)走向角為strike=303.6°，傾角為dip=16.2°，震源中心深度depth=12.3 km。將斷層長和寬擴(kuò)展至44、54 km，并將其按照2 km×2 km矩形劃分為594個(gè)單元。為了避免少量滑動(dòng)單元出現(xiàn)方向矛盾問題，采用非負(fù)約束最小二乘法進(jìn)行線性滑動(dòng)分布求解[33]。分別采用EI曲線法與L曲線法確定平滑因子進(jìn)行地震滑動(dòng)分布反演。

在美濃實(shí)際震例反演中，平滑因子α均以0.05為初值、0.002為步長、1為終值繪制L曲線圖、EI曲線圖(其中，L曲線法確定的平滑因子大小為0.445，EI曲線法確定的平滑因子大小為0.328。在利用EI曲線法確定平滑因子過程中，利用方差分量估計(jì)確定實(shí)際觀測數(shù)據(jù)與虛擬觀測數(shù)據(jù)內(nèi)部權(quán)陣修正因子，大小分別為a=1，b=0.831 1)。利用EI曲線法和L曲線法確定平滑因子如圖6(a)、(c)所示；兩種方法反演美濃地震結(jié)果及兩種方法反演結(jié)果差值分布如圖6(b)、(d)、(e)所示。

圖3 GPS 3方向模擬觀測點(diǎn)Fig.3 Distribution of three directions simulated points of GPS

圖4 模擬試驗(yàn)兩種方法反演地震滑動(dòng)分布結(jié)果及解的殘差Fig.4 Slip distribution results of simulation experiment and the residuals of two methods

圖6 兩種方法反演美濃地震滑動(dòng)分布結(jié)果及差值分布Fig.6 The slip distribution inversion results of Meinong earthquake and the residuals of two methods

利用L曲線法、EI曲線法確定的平滑因子反演美濃地震滑動(dòng)分布各參數(shù)結(jié)果見表5；近年來國內(nèi)外部分其他學(xué)者反演美濃地震滑動(dòng)分布各參數(shù)結(jié)果見表6。

表5 用L曲線法、EI曲線法確定平滑因子反演美濃地震滑動(dòng)分布結(jié)果

表6 美濃地震斷層滑動(dòng)參數(shù)

從圖6(b)、(d)可以看出，利用L曲線法、EI曲線法確定的平滑因子進(jìn)行美濃地震反演時(shí)，其反演結(jié)果均顯示美濃地震滑動(dòng)分布區(qū)域主要分布在地下9～14 km，該結(jié)果同文獻(xiàn)[32]反演結(jié)果一致。表6列舉了美濃地震的部分研究成果，表5、表6顯示對于最大滑動(dòng)量，本文兩種方法反演結(jié)果分別為0.453、0.479 m。文獻(xiàn)[32]聯(lián)合InSAR和GPS數(shù)據(jù)反演美濃地震結(jié)果顯示最大傾滑量為0.517 m、最大走滑量為0.553 m，略大于本文兩種方法反演結(jié)果；對于平均滑動(dòng)角，本文兩種方法反演的結(jié)果分別為45.32°、43.45°，其中，利用L曲線法反演結(jié)果與文獻(xiàn)[34]反演結(jié)果相近。而EI曲線法反演結(jié)果在文獻(xiàn)[34]與文獻(xiàn)[32]反演結(jié)果之間；對于矩震級參數(shù)，本文兩種方法反演結(jié)果相近，且在其他學(xué)者研究范圍內(nèi)；文獻(xiàn)[35]反演的矩震級結(jié)果為Mw6.52，對于其反演結(jié)果較其他學(xué)者反演結(jié)果偏大的原因可能與沒有較好地?cái)M合GPS數(shù)據(jù)有關(guān)[32]。

對于本文所做的模擬試驗(yàn)以及拉奎拉與臺灣美濃實(shí)際震例試驗(yàn)，作出以下幾點(diǎn)分析：

(1) 對于兩種方法確定平滑因子反演地震滑動(dòng)分布各參數(shù)精度，由模擬試驗(yàn)可知，利用EI曲線法確定平滑因子反演滑動(dòng)分布的各參數(shù)結(jié)果與模擬值的差值均小于L曲線法。模擬試驗(yàn)及兩個(gè)實(shí)際震例反演的均方根誤差結(jié)果表明EI曲線法反演的均方根誤差均小于L曲線法，從而進(jìn)一步說明EI曲線法確定的平滑因子反演地震滑動(dòng)分布精度要高于L曲線法。

(2) 對于計(jì)算效率，在模擬試驗(yàn)、拉奎拉及美濃實(shí)際震例中EI曲線法確定平滑因子時(shí)間較L曲線法分別減少107.784 3、125.664 1、74.729 4 s。總體上，EI曲線法較L曲線法效率提高2～7倍，故利用EI曲線法確定平滑因子較L曲線法可以有效地提高計(jì)算效率。

4 結(jié) 論

本文提出的EI曲線法為地震滑動(dòng)分布反演中平滑因子的確定增加了一條新的途徑。本文通過模擬試驗(yàn)及拉奎拉與美濃實(shí)際地震反演表明，EI曲線法較L曲線法確定平滑因子具有定位準(zhǔn)確、反演精度高、計(jì)算效率高，以及EI曲線圖橫坐標(biāo)直觀反映平滑因子大小、無需依賴于數(shù)據(jù)擬合度等優(yōu)勢。本文提出的EI曲線法可以應(yīng)用在地震滑動(dòng)分布反演中，在其他領(lǐng)域的應(yīng)用有待于進(jìn)一步的研究與探索。

致謝：感謝武漢大學(xué)溫?fù)P茂博士提供意大利拉奎拉地震的InSAR數(shù)據(jù)。