王樂洋， 李志強

1 東華理工大學(xué)測繪工程學(xué)院， 南昌 330013 2 山東建筑大學(xué)測繪地理信息學(xué)院， 濟南 250101

0 引言

地震是地殼巖層受力突破彈性極限(抗剪切強度)后快速破裂錯動引起的地表振動或破壞，是對人類文明威脅最大的自然災(zāi)害之一.震源參數(shù)反演是研究地震斷層錯動特征、震源機制、觸發(fā)關(guān)系等的重要方式，也是地震危險性評價及地震預(yù)警的基礎(chǔ).地表形變是地震發(fā)生時最直接的表現(xiàn)，是研究地震斷層構(gòu)造、破壞程度及地震機理等的重要資料(王樂洋等，2017).在獲取地表形變數(shù)據(jù)的過程中，由于各種噪聲的干擾導(dǎo)致形變數(shù)據(jù)的準確性和可靠性受到影響，進而導(dǎo)致震源參數(shù)的精度損失.因此，獲取震源參數(shù)的精度信息、誤差水平及可信程度顯得尤為重要，對震源參數(shù)進行精度評定研究成為了提升成果質(zhì)量的重要內(nèi)容，對準確、深刻的理解地震也具有一定的實際意義.

由于依據(jù)彈性半空間矩形位錯模型(Okada，1985)建立的關(guān)于地表形變與震源參數(shù)之間的關(guān)系為復(fù)雜多維的非線性函數(shù)，因此需要依靠搜索算法進行非線性反演.常用的搜索算法有模擬退火法(Kirkpatrick et al.，1983)、粒子群算法(Kennedy and Eberhart，1995；馮萬鵬和李振洪，2010)、神經(jīng)網(wǎng)絡(luò)法(Barradas and Vieira，2000)、遺傳算法(Genetic Algorithm，GA)(Holland，1975)等，該類方法無需獲取參數(shù)的解析表達式，也無需求導(dǎo)計算，因此常用于非線性反演獲取參數(shù)最優(yōu)值.遺傳算法不必明確地描述優(yōu)化問題的全部特征，且不存在連續(xù)性、可導(dǎo)或單峰等特殊要求，因此具有通用性.此外，該搜索算法全局性性強，不易陷入局部最優(yōu)，并且適合并行計算(王新洲，2002)，因此較適用于數(shù)據(jù)量大的大地測量震源參數(shù)非線性反演問題研究(Nolte and Frazer，1994；萬永革等，1997；王文萍和王慶良，1999).

將震源參數(shù)反演過程視為一個非線性函數(shù)，目前已有的非線性精度評定方法有近似函數(shù)法(徐培亮，1986；Xue et al.，2015；Wang and Zhao，2019)和近似非線性函數(shù)的概率密度分布方法(王樂洋和趙英文，2019；Wang and Yu，2019；Wang and Zou，2019).近似函數(shù)法需要利用泰勒級數(shù)公式對非線性模型展開，根據(jù)誤差傳播定律獲取參數(shù)估值的方差或協(xié)方差.但由于該方法無法避免求導(dǎo)計算，因此在震源參數(shù)精度評定研究中適用性較差.通過采樣來代替復(fù)雜求導(dǎo)計算的非線性精度評定方法包括蒙特卡羅(Monte Carlo, MC)方法、自適應(yīng)蒙特卡羅(Adaptive Monte Carlo, AMC)方法、Jackknife方法、Sterling插值方法等.該類方法一般通過重復(fù)采樣獲取隨機樣本來近似非線性函數(shù)的概率密度分布，進而輸出待統(tǒng)計量的精度信息，這些采樣方法在一般的非線性模型精度評定問題中能夠取得較好的效果.在震源參數(shù)反演精度評定方面，Monte Carlo方法得到了較為廣泛的應(yīng)用(Wright et al.，1999；Funning et al.，2005；Dawson and Tregoning，2007；Xu et al.，2017)，亦有學(xué)者將Jackknife方法、SUT(Scaled Unscented Transformation)方法、Sterling插值方法用于地震大地測量反演精度評估(Ichinose et al.，2003；Prieto et al.，2007；Wang and Zou，2021；Wang and Ding，2020).

Bootstrap方法(Johnson，2001；Chen，2018)也稱為自助法，是一種用于估計和修正統(tǒng)計估計值的偏差或方差信息的一種統(tǒng)計推斷方法.Efron(1979)首次提出了通過重復(fù)采樣得到的自助世界的經(jīng)驗分布來逼近總體真實分布的思想，并給出了其基本采樣策略和相關(guān)理論證明.由于Bootstrap方法通過重采樣代替求導(dǎo)計算，只需選擇采樣次數(shù)并結(jié)合有放回抽樣策略獲取自助樣本，最后解算每個樣本便可求得未知統(tǒng)計量的均值和協(xié)方差陣，因此歷年來有較多學(xué)者將該方法用于偏差分析和方差計算(Sahinler and Topuz，2007；Kitagawa and Konishi，2010；O′Hagan et al.，2019).目前，關(guān)于Bootstrap方法的研究成果較為豐富，主要集中于該方法在統(tǒng)計學(xué)上的性質(zhì)(Martínez-Camblor and Corral，2012；Pan and Politis，2016).目前，已有部分學(xué)者將Bootstrap方法用于地球科學(xué)反演精度評估問題，Wéber(2005)利用貝葉斯方法和Bootstrap重采樣對震源參數(shù)進行了誤差分析；Cesca等(2014)通過Bootstrap采樣來評估發(fā)震斷層位置的不確定性；Jia等(2017)采用Bootstrap方法來估計質(zhì)心深度的不確定性.但Bootstrap采樣方法在大地測量震源參數(shù)反演精度評定方面的性能未得到較為全面的研究，尚未給出將該方法用于獲取震源參數(shù)精度信息的詳細采樣步驟及具體實現(xiàn)流程，并且在震源參數(shù)中誤差計算及區(qū)間估計方面的對比評價也未得到重視.

針對以上問題，本文從為了獲取參數(shù)估值更為合理的精度信息出發(fā)，將Bootstrap方法用于震源參數(shù)非線性反演的精度評定研究中，通過對地表形變數(shù)據(jù)實施自助重采樣并結(jié)合遺傳算法，進而設(shè)計得到震源參數(shù)反演精度評定的Bootstrap方法計算流程.最后通過6個模擬地震、Amatrice地震和Visso地震，對本文精度評定方法的有效性和可靠性進行驗證.

1 震源參數(shù)反演與遺傳算法

1.1 震源參數(shù)反演

均勻彈性半空間矩形位錯模型是反演震源參數(shù)的常用有效模型(Okada，1985)，本文采用Okada矩形位錯模型并結(jié)合GPS臺站地表形變觀測資料進行震源參數(shù)反演.地表形變量觀測值與震源參數(shù)之間為高度非線性關(guān)系，其函數(shù)模型可表示為

d=G(ξ)+e,

(1)

式中，d為GPS觀測到的地表形變量，它包含dx、dy及dz三個分量，分別表示東西、南北和垂直三個不同方向的地表形變.G(·)為利用Okada彈性位錯理論模型建立的關(guān)于GPS地表形變量與震源參數(shù)的復(fù)雜非線性關(guān)系.e為GPS地表形變量觀測誤差.ξ為震源參數(shù)，包括發(fā)震斷層位置(經(jīng)度、緯度、深度)，斷層尺寸(長度、寬度)，斷層幾何(走向角、傾角、滑動角)以及滑動量；本文以頂深、底深、走向角、傾角、長度、滑動角、滑動量以及斷層中心位置為參數(shù)進行反演和精度評定，其中通過頂深和底深可計算獲得深度，長度可聯(lián)立頂深、底深及傾角獲取.

1.2 遺傳算法在震源參數(shù)反演中的應(yīng)用

由于震源參數(shù)未知向量與地表形變觀測值向量之間是復(fù)雜的多維非線性關(guān)系，本文將遺傳算法運用于震源參數(shù)反演問題.遺傳算法最早由美國教授Holland(1975)提出，該算法是一種啟發(fā)式高級全局搜索算法，它所依據(jù)的是自然界生物優(yōu)勝劣汰的生存法則，是基于達爾文的進化論和孟德爾的遺傳學(xué)說建立的.遺傳算法僅通過目標函數(shù)而不需要其他附加條件對種群(將待求參數(shù)稱之為種群，每個參數(shù)稱之為個體)進行全局搜索.對提高算法的搜索效率并保證反演得到的震源參數(shù)具有物理意義，反演前通常需要給定震源參數(shù)的搜索區(qū)間作為約束.本文在使用遺傳算法進行震源參數(shù)反演時，首先在全參數(shù)空間內(nèi)進行全局搜索，當獲取初步的反演結(jié)果后再縮減搜索空間并通過反復(fù)迭代計算最后求得穩(wěn)定的最優(yōu)解.

將遺傳算法用于震源參數(shù)反演的大致步驟是首先使用二進制對區(qū)間內(nèi)的數(shù)進行編碼；其次，隨機生成一組任意排列的字符串，即產(chǎn)生初始種群；然后根據(jù)提前設(shè)置的適應(yīng)度函數(shù)判斷個體的優(yōu)劣性，淘汰適應(yīng)度低的個體，選擇復(fù)制適應(yīng)度較高的個體作為父代并且進行交叉和變異操作進而生成下一代個體；最后通過不斷繁殖、選擇、交叉及變異，便可朝最優(yōu)解的方向進化獲得最優(yōu)的震源參數(shù).遺傳算法實現(xiàn)非線性反演的具體過程總結(jié)為

1)選擇初始群體大小pop及最大迭代次數(shù)maxgen.

2)編碼：編碼是將震源參數(shù)解的搜索范圍轉(zhuǎn)換到遺傳算法所能認知的搜索空間的一種轉(zhuǎn)換方式，通過編碼將具體的非線性反演問題與自然界生物進化建立聯(lián)系.

3)生成初始種群：在編碼區(qū)間內(nèi)隨機生成數(shù)量為pop的初始參數(shù)，作為遺傳算法求解優(yōu)化問題時的初始解的集合.

4)個體適應(yīng)度函數(shù)值的計算：適應(yīng)度函數(shù)是用來評價個體優(yōu)劣的指標.震源參數(shù)反演中常將殘差的均方根即式(2)作為目標函數(shù)，計算每個個體的適應(yīng)度函數(shù)值，區(qū)分出每個個體的優(yōu)劣性.

5)選擇復(fù)制：將步驟4)所計算出的個體適應(yīng)度函數(shù)值作為評價標準，選擇出適應(yīng)度較好的個體作為父代來繁衍后代種群.

6)交叉：交叉是個體兩兩配對的過程即基因重組過程，具體實現(xiàn)過程為從父代中選取部分基因，并以某個交叉概率交換部分染色體進而生成新的個體.該過程的目的是為了生成新一代個體以提升算法搜索能力.

7)變異：基因遺傳算法中的變異過程是指等位基因替換引起的基因結(jié)構(gòu)的變化，是自然界生物進化的必然要求.變異的目的是保證種群多樣性，防止陷入局部最優(yōu).

8)更新：將步驟7)獲取到的子代作為更適應(yīng)環(huán)境的父代種群，重復(fù)步驟4)～步驟7)，重新計算個體適應(yīng)度函數(shù)值、選擇、交叉、變異，直到迭代次數(shù)大于設(shè)定的最大迭代次數(shù)即Itei>maxgen則跳出迭代.那么，最新一代中適應(yīng)度函數(shù)值最小所對應(yīng)的個體即為震源參數(shù)的最優(yōu)參數(shù)估值.

在使用基因遺傳算法反演震源參數(shù)的過程中，本文采用地表形變觀測值與Okada模型預(yù)測值差值的均方根(Root Mean Square, RMS)誤差作為目標函數(shù)來評估反演結(jié)果，其表達式為

(2)

式中，n為地表GPS觀測點個數(shù)，Pi表示每個觀測值的權(quán)值，di為形變實際觀測值，ci為通過反演模型得到的地表形變預(yù)測值.

2 震源參數(shù)反演及精度評定的Bootstrap方法

Bootstrap方法將原始觀測值樣本作為總體，利用有放回隨機抽樣法從總體分布函數(shù)中得到一系列獨立樣本(稱為自助樣本)，并通過計算每個自助樣本來獲取未知統(tǒng)計量抽樣分布的經(jīng)驗估計(Efron，1979；Sahinler and Topuz，2007).該方法的優(yōu)勢性主要表現(xiàn)在它不需要對未知模型及分布做任何假設(shè)，也無需推導(dǎo)估計量的精確表達式，只需通過檢驗樣本內(nèi)統(tǒng)計量的變化便能夠估計未知參數(shù)的均值和方差.Bootstrap方法通過重采樣獲取自助樣本的方式來代替復(fù)雜的求導(dǎo)運算，因此本文將Bootstrap方法引入到震源參數(shù)精度評定問題，試圖通過在避免求導(dǎo)計算的情形下依然能夠獲取合理的精度評定結(jié)果.

2.1 精度評定的Bootstrap方法

(3)

震源參數(shù)估值的方差-協(xié)方差陣計算公式為

(5)

=(1-α),

(6)

式中，pr表示概率，α為顯著性水平，Zα/2為對應(yīng)的標準分數(shù).

2.2 精度評定計算步驟

為充分利用總體包含在樣本中的先驗信息和數(shù)據(jù)性質(zhì)，Bootstrap方法通過有放回隨機抽樣將GPS觀測點三個方向形變量數(shù)據(jù)采樣成多個自助樣本，原始樣本集中的每個樣本點被采樣到的概率相同，每個自助樣本的樣本容量與原始樣本相同.顧及GPS站點形變量含有誤差，故對三個方向形變量數(shù)據(jù)按照其對應(yīng)的觀測精度進行加權(quán)處理，通常將三個方向形變量所對應(yīng)的觀測中誤差作為定權(quán)的依據(jù).在重采樣獲得自助樣本的過程中，需保持自助樣本中的形變觀測值與其權(quán)值Pi相對應(yīng)，具體表現(xiàn)為重采樣三個方向形變量數(shù)據(jù)的中誤差.形變值與對應(yīng)中誤差兩者的采樣方式一致.根據(jù)以上分析，總結(jié)得到震源參數(shù)精度評定的Bootstrap方法的完整計算步驟：

1)假設(shè)n個GPS觀測點的三個方向形變量為隨機項樣本N=(ux,uy,uz)，其對應(yīng)的觀測中誤差為ε=(σx,σy,σz)；

2)將1/n的概率分別設(shè)置在N中的每個樣本值上；

3)產(chǎn)生服從均勻分布的n個隨機數(shù)(i1,i2,…in)，其中1≤ij≤n,j=1,2,3,…,n；

根據(jù)以上計算步驟，可得Bootstrap方法評定震源參數(shù)精度的計算流程圖，如圖1所示.

圖1 震源參數(shù)精度評定的Bootstrap方法迭代流程圖Fig.1 The iteration flow chart of the Bootstrap resampling method for precision estimation of fault source parameters

由以上的震源參數(shù)精度評定的Bootstrap方法計算流程可知，本文精度評定方法獲得的所有自助樣本中的數(shù)據(jù)均來源于原始固定的GPS觀測點三個方向形變量數(shù)據(jù)，并未根據(jù)更多的觀測形變信息進行計算.Bootstrap方法的采樣過程不會產(chǎn)生額外的模型誤差，也不會改變模型態(tài)性.盡管每個自助樣本的樣本容量與原始樣本相同，但不是改變GPS觀測點三個方向形變量數(shù)據(jù)的先后排列順序，有放回隨機抽樣過程使得每個自助樣本中可能存在重復(fù)的數(shù)據(jù)點，而另外一些樣本點則沒有出現(xiàn).因此，每個自助樣本將隨機地異于原始樣本，導(dǎo)致每個自助樣本反演得到的震源參數(shù)估值存在細微差異.

3 實驗與分析

為驗證Bootstrap方法在震源參數(shù)非線性反演精度評定中的可靠性和有效性以及探討該方法在實際震例中的性能，本文共設(shè)計了6個模擬實驗，同時也將本文精度評定方法用于獲取Amatrice地震和Visso地震的震源參數(shù)精度信息，并與Jackknife方法和Monte Carlo方法進行對比分析.當使用到遺傳算法時，為同時保證較高的反演精度和反演效率，初始群體大小取值為450，最大迭代次數(shù)選擇200次，交叉概率為0.85，變異概率選擇0.1.Bootstrap采樣次數(shù)的選擇取決于統(tǒng)計量的具體研究問題并且依賴于需要做的檢驗(例如方差或置信區(qū)間)(Efron，1981)，諸多統(tǒng)計學(xué)家給出了抽樣次數(shù)的經(jīng)驗取值，若要獲取未知參數(shù)較為合理的方差或中誤差，采樣次數(shù)應(yīng)為50～200之間，當采樣次數(shù)超過200時，方差近似值的改進很小(Efron，1981；Efron and Tibshirani，1986；Léger et al.，1992).經(jīng)實驗，通過增加自助采樣次數(shù)在一定程度上能夠提升精度評定結(jié)果質(zhì)量，但隨著采樣次數(shù)的增加，中誤差結(jié)果趨于較為穩(wěn)定的值，而其計算量將成倍遞增.本文依據(jù)已有研究的經(jīng)驗取值并顧及復(fù)雜多維的地震非線性反演實際問題，在保證較高的反演精度條件下，同時考慮較高的反演效率，因此本文將Bootstrap方法的采樣次數(shù)設(shè)置為300.

圖2 仿真設(shè)置1的斷層平面位置及觀測點位分布Fig.2 The fault plane position and observation point positions of the first simulation setting

3.1 模擬地震

通過仿真獲取某一地震區(qū)域GPS觀測點三個方向的形變量資料(仿真設(shè)置1)：斷層面的幾何中心位于坐標系原點，37個GPS觀測點分布在斷層面幾何中心的40 km范圍內(nèi)，模擬斷層平面位置及GPS觀測點分布如圖2所示(其中橫坐標E為東方向，縱坐標N為北方向，矩形框表示斷層投影平面)；斷層面幾何參數(shù)的仿真設(shè)置為：斷層頂深為2 km、底深為18 km、走向角為45°、傾角為60°、斷層長度為40 km，滑動角為70°、滑動量為1.5 m；將GPS觀測點坐標及參數(shù)模擬真值代入Okada模型正演獲取形變量真值；最后在形變值真值上施加服從均值為0、單位權(quán)中誤差為3 mm的隨機正態(tài)分布噪聲進而獲取模擬形變量.仿真設(shè)置1的GPS同震形變場在水平方向和垂直方向的模擬形變值如圖3所示.

在獲取GPS形變量數(shù)據(jù)后，分別采用Jackknife方法和Bootstrap方法對地表形變值進行采樣并獲取震源參數(shù)均值和中誤差；為推斷總體均值，利用Bootstrap方法計算震源參數(shù)的置信區(qū)間，并與Jackknife方法的結(jié)果進行對照.各方法獲取仿真設(shè)置1的震源參數(shù)結(jié)果、置信上下限及區(qū)間間隔列于表1.

由表1中的震源參數(shù)估值可以看出，使用遺傳算法搜索到的震源參數(shù)估值與參數(shù)模擬真值基本一致，表明本文采用遺傳算法進行非線性反演是正確且可靠的.Jackknife方法和Bootstrap方法得到的震源參數(shù)在符號和量級上與遺傳算法的結(jié)果相同，且Bootstrap方法的結(jié)果與真值的差值二范數(shù)小于遺傳搜索算法和Jackknife方法的結(jié)果；根據(jù)遺傳算法、Jackknife方法和Bootstrap方法獲取的震源參數(shù)得到的矩震級均為MW6.97，與模擬的矩震級MW6.98的差異僅在小數(shù)點后第2位.這表明基于重采樣的Bootstrap方法能夠適用于震源參數(shù)非線性反演，也說明該方法通過對形變量數(shù)據(jù)實施自助采樣在一定程度上能夠降低反演估值的偏差，因此其震源參數(shù)反演值比遺傳算法和Jackknife方法的結(jié)果更加靠近模擬真值.

圖3 仿真設(shè)置1的GPS同震形變場(a) 水平位移場； (b) 為垂直位移場.Fig.3 The coseismic displacement fields of the first simulation setting obtained by GPS(a) The horizontal displacement field； (b) The vertical displacement field.

表1 仿真設(shè)置1的震源參數(shù)搜索區(qū)間、參數(shù)估值及置信水平為95%情況下的置信區(qū)間計算結(jié)果Table 1 Search ranges, inversion values, and confidence intervals with 95% confidence level of source parameters in the first simulation setting

根據(jù)表1中各參數(shù)的置信上限、置信下限及區(qū)間間隔可以看出，在95%置信水平下，Jackknife方法獲取的震源參數(shù)置信區(qū)間均比Bootstrap方法獲取的置信區(qū)間更寬，說明Jackknife方法反演得到的震源參數(shù)的可靠性較弱，因此Bootstrap方法的反演精度較高.這表明在震源參數(shù)非線性反演問題中，Bootstrap方法獲取的參數(shù)估值比Jackknife方法的結(jié)果更具可信度.由于置信區(qū)間的計算公式取決于所用的統(tǒng)計量，Jackknife方法的樣本量為總觀測個數(shù)，而Bootstrap方法的樣本量可根據(jù)模型復(fù)雜度或樣本容量進行調(diào)整.因此，與Jackknife方法相比，Bootstrap方法在震源參數(shù)反演中更具可靠性和適用性.

圖4 仿真設(shè)置6的GPS同震形變場(a) 水平位移場； (b) 垂直位移場.Fig.4 The coseismic displacement fields of the sixth simulation setting obtained by GPS(a) The horizontal displacement field； (b) The vertical displacement field.

表2 六種實驗方案及其說明Table 2 The six experimental schemes and their explanations

在獲取以上6個仿真設(shè)置的形變量數(shù)據(jù)后，分別利用Bootstrap方法、Jackknife方法和Monte Carlo方法獲取震源參數(shù)的中誤差，3種精度評定方法的計算結(jié)果列于表3，并將6個不同仿真設(shè)置的中誤差可視化繪于圖5((a—f)所繪中誤差分別表示仿真設(shè)置1～6的結(jié)果).

通過表3中的參數(shù)中誤差結(jié)果可以看出，Bootstrap方法在6個不同的仿真設(shè)置下所獲得的震源參數(shù)中誤差與Monte Carlo方法的結(jié)果在符號和量級上相同、在數(shù)值上相近，表明基于重采樣的Bootstrap方法能夠?qū)Ψ囱莸玫秸鹪磪?shù)進行有效的精度評定.從計算結(jié)果上看，相比于Monte Carlo方法所獲得的震源參數(shù)中誤差，Jackknife方法得到的結(jié)果偏離較大，是由于Jackknife法易受原始樣本容量大小的影響且其重采樣過程獲取的刀切樣本較少所致.表明Jackknife方法的結(jié)果可能會低估震源參數(shù)反演估值的精度，因此無法反映參數(shù)估值的真實精度.Bootstrap方法獲取的震源參數(shù)中誤差小于Jackknife法的結(jié)果且與Monte Carlo方法的結(jié)果較為接近.從圖5也可直觀看出，在不同仿真設(shè)置情況下，Jackknife方法獲取的參數(shù)中誤差總是表現(xiàn)出較多深色，而Bootstrap方法的結(jié)果與Monte Carlo方法所得中誤差顏色較為相近且淺色居多，這表明將Bootstrap方法用于解決震源參數(shù)精度評定問題是可行的，也表明Bootstrap方法在獲取震源參數(shù)精度信息方面較Jackknife方法更加精確有效.

表3 六種不同仿真設(shè)置下震源參數(shù)的中誤差計算結(jié)果Table 3 The standard errors of source parameters under the six different simulation settings

圖5 不同仿真設(shè)置情況下的震源參數(shù)中誤差圖圖中，“TD”表示頂深(Top Depth)、“BD”表示底深(Bottom Depth)、“Strike”表示走向角(Strike)、“Dip”表示傾角(Dip)、“Length”表示長度(Length)、“Rake”表示滑動角(Rake)、“Slip”表示滑動量(Slip)Fig.5 The standard error of source parameters under the six different simulation settingsLet “TD” denote the Top Depth, and “BD” indicate the Bottom Depth. “Strike”, “Dip”, “Length”, “Rake” and “Slip” corresponding to the source parameters Strike, Dip, Length, Rake and Slip, respectively.

通過表3中仿真設(shè)置1和仿真設(shè)置2的結(jié)果對比、以及仿真設(shè)置3和仿真設(shè)置4的結(jié)果對比可以看出，在GPS觀測點個數(shù)相同的情況下，正態(tài)分布隨機噪聲越大，通過執(zhí)行不同的采樣方法獲取的震源參數(shù)的中誤差也均相應(yīng)增大.對比圖5a與5b、圖5c與5d也可以直觀看出，單位權(quán)中誤差為5 mm的隨機噪聲情況下得到的震源參數(shù)中誤差比3 mm情況下誤差顏色更深，表明不同誤差水平對震源參數(shù)精度的影響不同：隨機誤差越大，震源參數(shù)的精度越低.

對比表3中仿真設(shè)置1、3、5的結(jié)果可以看出，當GPS觀測點三個方向形變量數(shù)據(jù)的隨機噪聲為同分布情況下，將仿真設(shè)置1中的GPS觀測點個數(shù)增加一倍，并未對仿真設(shè)置5的參數(shù)中誤差結(jié)果造成較大的影響.并且隨著觀測點個數(shù)的增加，參數(shù)中誤差并未呈現(xiàn)出較明顯的遞減趨勢.通過對比仿真設(shè)置2、4的結(jié)果，圖5a、5c、5e，以及圖5b與5d的中誤差也可發(fā)現(xiàn)，當誤差水平相同條件下，觀測點的個數(shù)對震源參數(shù)的中誤差影響較小.因此，若試圖通過增加GPS觀測點數(shù)量來提升震源參數(shù)的精度可能并不十分顯著.

從表3中仿真設(shè)置6的結(jié)果可以看出，在GPS站點隨機分布情況下，Bootstrap方法依然能夠估出合理的震源參數(shù)中誤差，而Jackknife方法的結(jié)果與Monte Carlo方法結(jié)果偏離較大.這表明Bootstrap方法同樣適用于站點分布不均勻的普遍情況，并且能夠獲取比Jackknife方法更加合理的震源參數(shù)精度信息.對比仿真設(shè)置1與仿真設(shè)置6的結(jié)果可以看出，在相同的隨機誤差水平及GPS觀測點個數(shù)的條件下，通過執(zhí)行不同采樣方法獲取的不均勻分布站點情況下的震源參數(shù)中誤差較均勻分布情況下的結(jié)果更大，這可能是因為不均勻分布站點對斷層面的約束力不及均勻分布情況，進而導(dǎo)致更大的震源參數(shù)中誤差.

3.2 Amatrice地震

2016年8月24日，意大利中部阿馬特里斯地區(qū)發(fā)生一次MW6.2地震，造成約300人傷亡且震中附近大量房屋被毀.美國地質(zhì)調(diào)查局(United States Geological Survey, USGS)發(fā)布該地震的震中位置為(13.188°E，42.723°N)，震源深度為4.4 km，屬于淺層地震.震后，許多研究機構(gòu)和學(xué)者快速對該地震進行了研究，分析了地震的震源機制、觸發(fā)關(guān)系以及庫侖應(yīng)力(Wang et al., 2018)；亦有學(xué)者對該地震進行了震源參數(shù)反演和精度評估(Xu et al.，2017)，但針對不同的重采樣精度評定方法的性能評估對比、以及Jackknife方法和Bootstrap方法區(qū)間估計的對比分析未得到深入的研究.

本文獲取了Amatrice地震的106個GPS觀測臺站的形變量數(shù)據(jù)(https:∥ring.gm.ingv.it/)，將Amatrice地震GPS同震形變場在水平方向和垂直方向的實際形變值、利用Okada模型正演得到的預(yù)測值及其對應(yīng)殘差分別繪于圖6(a、b、c).GPS形變整體擬合較好，水平方向的擬合殘差RMS為0.11 cm，垂直方向的擬合殘差RMS為0.36 cm，略小于Wang等(2018)得到的結(jié)果.水平方向的擬合精度整體上高于垂直方向，僅在主要形變區(qū)域的南側(cè)及西北側(cè)存在少量較大垂直殘差.

在獲取形變量數(shù)據(jù)后，經(jīng)過坐標轉(zhuǎn)換將原始數(shù)據(jù)中的經(jīng)緯度坐標轉(zhuǎn)換為以USGS發(fā)布的震中位置為原點的大地坐標；在均勻位錯模型和遺傳算法的基礎(chǔ)上，利用Bootstrap方法對該地震的形變量數(shù)據(jù)進行重采樣，通過遺傳算法搜索震源參數(shù)并獲取斷層參數(shù)的置信區(qū)間；再采用Bootstrap方法獲取參數(shù)的中誤差并與Jackknife方法和Monte Carlo方法的計算結(jié)果進行對比評估.參數(shù)的搜索區(qū)間、部分已有文獻研究結(jié)果、Monte Carlo方法獲取的中誤差，以及Jackknife方法和Bootstrap方法的反演估值、置信上下限、區(qū)間間隔和參數(shù)中誤差均列于表4.

表4 Amatrice地震震源參數(shù)搜索區(qū)間、反演估值、95%置信區(qū)間及參數(shù)中誤差計算結(jié)果Table 4 Search ranges, inversion values, and confidence intervals with 95% confidence level of source parameters for the Amatrice earthquake

圖6 Amatrice地震和Visso地震GPS同震形變場及殘差(a) Amatrice地震水平位移場； (b) Amatrice地震垂直位移場； (c) Amatrice地震的殘差圖； (d) Visso地震水平位移場； (e) Visso地震垂直位移場； (f) Visso地震的殘差圖；“Dis.”為實際形變值(Displacement)，“Pre.”為模型預(yù)測值(Prediction)，“Hor.”為水平方向殘差，“Ver.”為垂直方向殘差.Fig.6 The coseismic displacement fields and the corresponding residuals of Amatrice earthquake and Visso earthquake obtained by GPS(a) The horizontal displacement field of Amatrice earthquake； (b) The vertical displacement field of Amatrice earthquake； (c) The corresponding residuals of Amatrice earthquake; (d) The horizontal displacement field of Visso earthquake； (e) The vertical displacement field of Visso earthquake； (f) The corresponding residuals of Visso earthquake; “Dis.” is the actual displacement, “Pre.” is the model prediction; “Hor.” is the horizontal residual, and “Ver.” is the vertical residual.

根據(jù)表4中的震源參數(shù)估值可以看出，遺傳算法、Jackknife方法及Bootstrap方法獲取的Amatrice地震震源參數(shù)較為接近.由于各研究者采用的反演方法和數(shù)據(jù)的不同，這3種方法得到的震源參數(shù)與已有文獻的研究成果有所差異.通過對比各研究所得的矩震級可以發(fā)現(xiàn)，Bootstrap方法獲取的震源參數(shù)得到的矩震級與已有研究成果較為接近，之間的差異在僅在小數(shù)點后第2位.這證明了將Bootstrap重采樣與遺傳算法相結(jié)合、并用于解決Amatrice地震的震源參數(shù)反演問題是可行的.

由表4中Jackknife獲取的震源參數(shù)置信區(qū)間的上下限可以看出，頂深、底深及滑動量的Jackknife置信區(qū)間的下限均為負值.其中滑動量的Jackknife置信區(qū)間為(-4.15551,5.86072)即存在95%的概率相信該區(qū)間包含滑動量的真值，因此該結(jié)果說明Jackknife方法求取的頂深、底深及滑動量的置信區(qū)間包含的均值信息是模糊的或不精確的.通過對比Jackknife方法和Bootstrap方法獲取到的震源參數(shù)置信區(qū)間的區(qū)間間隔可知，在置信水平同為95%的情況下，采用Bootstrap方法獲取的區(qū)間間隔與Jackknife方法獲取的結(jié)果相差較大，Bootstrap方法獲取的9個震源參數(shù)置信區(qū)間均比Jackknife方法獲取的區(qū)間更窄.這表明與Bootstrap方法相比，Jackknife方法反演得到的震源參數(shù)的精確度較低，而Bootstrap方法獲取的參數(shù)估值更加可靠.

通過表4中的震源參數(shù)中誤差結(jié)果可以看出，Jackknife方法獲取的震源參數(shù)中誤差與Monte Carlo方法所得結(jié)果偏離較大，是因為Jackknife方法的重采樣過程僅僅能夠得到等于總觀測數(shù)的少量刀切樣本量所致，因此Jackknife方法會低估Amatrice地震的震源參數(shù)精度，無法提供反演估值的有效精度信息.Bootstrap方法獲取的震源參數(shù)中誤差小于Jackknife方法的結(jié)果，盡管Bootstrap方法獲取的震源參數(shù)中誤差比Monte Carlo方法的結(jié)果略大，但這兩種方法得到的結(jié)果仍在符號和數(shù)量級上均保持一致、在數(shù)值上也較為接近.造成Bootstrap方法的中誤差結(jié)果略大的原因可能是方法采樣的隨機性以及非線性反演算法自身無法避免的不穩(wěn)定性.實驗結(jié)果說明了本文震源參數(shù)精度評定的Bootstrap方法的有效性，也表明Bootstrap方法在獲取震源參數(shù)精度信息方面較Jackknife方法更加精確.從9個參數(shù)的中誤差數(shù)值結(jié)果可以看出，對于Amatrice地震，走向角、傾角、滑動角的中誤差比其他6個震源參數(shù)更大，這可能是由于本文所采用的106個GPS觀測臺站的形變量數(shù)據(jù)中的隨機噪聲對角度量的反演產(chǎn)生較大影響所致，該現(xiàn)象也從側(cè)面反映出本文給出的震源參數(shù)精度評定的Bootstrap方法能夠合理反映Amatrice地震各個震源參數(shù)的精度信息.

3.3 Visso地震

Amatrice地震發(fā)生的兩個月后，震中30km范圍內(nèi)于2016年10月26日又發(fā)生一次6.1級地震.根據(jù)USGS消息，此次地震的震中位于13.067°E、42.956°N的維索西部約3km地區(qū).地震觸發(fā)后，許多學(xué)者和機構(gòu)采用不同的數(shù)據(jù)和方法對該地震展開深入研究，分析了該地震的觸發(fā)機制和周邊斷層的應(yīng)力變化(Wang et al., 2018)，反演了該地震的震源參數(shù)并獲取了參數(shù)精度(Xu et al.，2017).

本文獲取了127個GPS臺站形變數(shù)據(jù)(https:∥ring.gm.ingv.it/)生成Visso地震同震形變場，GPS水平方向的擬合效果較好，殘差RMS為0.17 cm，垂直方向殘差RMS為0.44 cm，與Wang等(2018)的計算結(jié)果較為一致.殘差的主要來源是GPS臺站的觀測誤差，而垂直方向殘差偏大是由于垂直觀測精度較水平方向精度低造成的.各臺站觀測得到的水平方向和垂直方向的實際形變值、模型正演得到的預(yù)測值及其對應(yīng)殘差繪于圖6(d、e、f).

依據(jù)Okada彈性半空間位錯模型理論及遺傳算法，利用本文設(shè)計的精度評定方法對該地震進行震源參數(shù)搜索、置信區(qū)間計算及震源參數(shù)精度信息獲??；并與Jackknife方法得到的參數(shù)置信區(qū)間和震源參數(shù)中誤差、Monte Carlo方法得到的中誤差進行對比分析，從而驗證Bootstrap方法在實際地震的震源參數(shù)精度評定問題中的廣泛適用性.參數(shù)搜索區(qū)間和部分已有文獻研究結(jié)果，Jackknife方法和Bootstrap方法的反演估值、置信上下限、區(qū)間間隔和參數(shù)中誤差，以及Monte Carlo方法獲取的中誤差均列于表5.

對比表5中的遺傳算法、Jackknife方法及Bootstrap方法獲取的Visso地震矩震級可以發(fā)現(xiàn)，Bootstrap方法的結(jié)果為MW6.04，略小于USGS提供的結(jié)果(MW6.10)及Xu等(2017)的結(jié)果(MW6.08)，且比遺傳算法和Jackknife方法更加接近已有的文獻成果.由表5中的震源參數(shù)95%置信區(qū)間結(jié)果可以發(fā)現(xiàn)，USGS提供的震源參數(shù)及Xu等(2017)的研究結(jié)果均在本文Bootstrap方法求取的置信區(qū)間內(nèi).由Jackknife方法計算得到的傾角置信下限為-31.34522、置信上限為107.78006.而Okada彈性半空間位錯模型中的傾角是斷層面與水平面的夾角，其范圍僅為(0,90].說明Jackknife方法獲取的傾角置信區(qū)間內(nèi)包含了一些無物理意義的信息，因此可信程度較低.由表中數(shù)據(jù)也可看出，Jackknife方法獲取的滑動量置信區(qū)間為(-8.68948, 11.53444)，區(qū)間長度達20.22392，說明Jackknife方法求取的滑動量置信區(qū)間過寬因此可靠性較弱.對比Jackknife方法和Bootstrap方法獲取的Visso地震震源參數(shù)置信區(qū)間能夠發(fā)現(xiàn)，Jackknife方法得到的9個參數(shù)置信區(qū)間寬度均比Bootstrap方法的結(jié)果偏大較多.以上分析說明了Bootstrap方法獲取的震源參數(shù)可信程度更高.

表5 Visso地震震源參數(shù)搜索區(qū)間、反演估值、95%置信區(qū)間及參數(shù)中誤差計算結(jié)果Table 5 Search ranges, inversion values, and confidence intervals with 95% confidence level of source parameters for the Visso earthquake

表5中的震源參數(shù)中誤差顯示，Jackknife方法獲取的震源參數(shù)中誤差與Monte Carlo方法所得結(jié)果差異較大，甚至出現(xiàn)量級不相等的情況.導(dǎo)致該現(xiàn)象是因為Jackknife重采樣容易受原始樣本影響，對于Visso地震，通過Jackknife方法重采樣僅能夠獲取127個少量刀切樣本，使得Jackknife方法獲取的震源參數(shù)中誤差因受樣本量的影響進而偏離參數(shù)的真實精度信息，因此Jackknife方法在實際震例的震源參數(shù)精度評定問題中可能不適用.Bootstrap方法獲取的震源參數(shù)中誤差小于Jackknife方法的結(jié)果，并且與Monte Carlo方法所獲結(jié)果在量級上一致、在數(shù)值上相近，說明Bootstrap方法能夠有效地起到對震源參數(shù)精度評定的作用，表現(xiàn)出比Jackknife方法更強的優(yōu)勢.實驗結(jié)果證明了本文方法用于Visso地震震源參數(shù)精度評定的可行性和有效性，也證明了本文方法具有廣泛的適用性.從參數(shù)的中誤差數(shù)值結(jié)果可以看出，在Visso地震精度評定研究中，分別使用Bootstrap方法和Monte Carlo方法獲取的走向角、傾角、滑動角的中誤差比其它6個參數(shù)的中誤差更大，可能是因為本文采用的127個GPS臺站形變量數(shù)據(jù)中的隨機噪聲對Visso地震的角度量參數(shù)反演產(chǎn)生較大影響所致.該現(xiàn)象也再次從側(cè)面反映了Bootstrap方法的確能夠合理反映實際震例的各個震源參數(shù)的精度信息.

4 結(jié)論

為獲取震源參數(shù)在統(tǒng)計意義上更加合理的精度信息，若采用傳統(tǒng)的泰勒級數(shù)展開方法在計算過程中無法獲取復(fù)雜多維非線性函數(shù)的偏導(dǎo)數(shù).本文在Okada彈性半空間矩形位錯模型及非線性反演遺傳算法的基礎(chǔ)上，引入了基于重采樣的Bootstrap方法，定義了一種通過實施自助采樣并解算自助樣本便可獲取震源參數(shù)中誤差的精度評定方法，并給出了具體的采樣步驟和迭代流程圖.最后將震源參數(shù)精度評定的Bootstrap方法用于模擬地震、Amatrice地震和Visso地震實驗中，通過對震源參數(shù)進行非線性反演、置信區(qū)間計算以及參數(shù)中誤差獲取，并與Jackknife方法計算得到的置信區(qū)間和中誤差、Monte Carlo方法獲取的中誤差進行對比，對本文方法進行了驗證.

模擬實驗結(jié)果顯示，采用本文方法獲取的震源參數(shù)均值與模擬真值更為接近；同時，在不同的隨機噪聲水平、觀測點個數(shù)以及GPS站點分布情況下，本文方法均能得到合理的精度信息.這表明將基于重采樣的Bootstrap方法用于震源參數(shù)非線性反演及精度評定問題是可行且有效的.在Amatrice地震和Visso地震實驗中，本文對比的已有研究結(jié)果均在Bootstrap方法求取的置信區(qū)間內(nèi)，且Bootstrap方法獲取的震源參數(shù)95%置信區(qū)間均比Jackknife方法得到的區(qū)間長度更窄.與Monte Carlo方法獲取的震源參數(shù)中誤差相比，Jackknife方法得到的參數(shù)中誤差與之差異較大.盡管Bootstrap方法獲取的震源參數(shù)中誤差略大于Monte Carlo方法的結(jié)果，但這兩種方法得到的結(jié)果仍在符號和數(shù)量級上均保持一致、在數(shù)值上也較為接近.實際地震結(jié)果驗證了本文震源參數(shù)精度評定的Bootstrap方法的有效性和可靠性.Monte Carlo方法的結(jié)果依賴于三個方向形變量的平差值以及模擬隨機誤差時的協(xié)方差陣，并且該方法必須對總體模型做出具體的假設(shè)，例如需要假設(shè)擾動項的概率分布為高斯分布，而所得的結(jié)論只適用于該特定的數(shù)據(jù)生成過程.此時，不需要對未知模型及分布做任何假設(shè)的Bootstrap方法則表現(xiàn)出更強的適用性.隨著Bootstrap方法引入震源參數(shù)精度評定問題，將拓展該采樣方法在大地測量地震反演領(lǐng)域中的應(yīng)用，對認知震源參數(shù)反演精度也具有一定的實際意義.

致謝感謝審稿專家對本文提出的寶貴意見.文中部分圖件是采用開源軟件GMT繪制.