NLLoc 定位方法在福建地區(qū)的適用性探究

張麗娜， 張紅才， 巫立華， 段 剛， 戴麗金， 廖詩榮

（福建省地震局， 福州 350003）

0 引言

地震定位是地震學中基礎問題之一， 也是地震臺網(wǎng)的基本任務之一。 地震定位結果對于震源幾何構造的研究、 地殼應力場分析等具有重要的意義。 快速準確的地震速報并產出地震目錄， 也可為地震應急救援、 震后減災和救災及震后地震趨勢預測等提供基礎數(shù)據(jù)[1]。 地震定位包括確定震源位置和發(fā)震時刻， 其精度受到定位方法、 地殼速度模型等諸多因素的影響。 地震定位方法研究及地震定位精度的提高， 一直是地震學研究的基本課題。

2000 年， Lomax 等[2]研究學者提出一種非線性搜索定位方法 NLLoc， 并開發(fā)了地震定位程序NonLinLoc（http：//alomax.free.fr/nlloc/）。 Stephan[3]利用NLLoc 定位方法對中國臺灣黃石國家公園地區(qū)的25 267 個地震事件重新定位， 與原始地震位置相比， 重新定位的位置條帶狀分布特征更明顯。Yavor 等[4]采用NLLoc 定位方法對加利福尼亞南部地震事件進行定位， 進而對該地區(qū)活動斷層特性等進行探索。 Theunissen[5]等研究指出三維速度模型能夠更好地解決深度定位問題， 尤其是在臺網(wǎng)邊界地區(qū) （如間隙角大于180°和首臺距離大于15 km）， 并且采用等時差 EDT 的 NLLoc 定位方法，給出震源參數(shù)及其誤差分析。 Antonino[6]應用NLLoc方法分析了臺站位置分布對定位震源位置質量的影響， 指出NLLoc 方法是一個快速有效的自動定位方法， 即使使用較少的地震震相， 定位結果仍然可靠有效， 且該算法對到時中的噪聲影響并不敏感。 區(qū)別于標準的三維線性迭代法， NLLoc 方法還給出了位置不確定性分析和解析度信息。 韓雪君[7-8]也簡要介紹了NLLoc 定位方法， 并將八叉樹搜索方法應用于預警的三維實時定位中， 能夠快速給出震源在三維空間的可能位置。 可見， 該方法有望大幅改善地震定位精度。

本文利用 “福建及臺灣海峽深部構造海陸聯(lián)測” 項目實施的22 次人工爆破事件記錄和88 次福建仙游震群序列M≥1.5 級地震事件記錄， 應用NonLinLoc 地震定位程序重新定位， 檢驗NLLoc 方法在福建臺網(wǎng)地震定位中的適用性。 為方便與現(xiàn)有定位算法進行對比， 本文定位中同樣采用華南地區(qū)一維速度模型[9]， 如表1， 并將其網(wǎng)格化為三維速度模型進行使用， 輸入相關震相文件， 進行重定位， 將其重定位結果與臺網(wǎng)常用定位方法（Hyposat、 HYP2000、 LocSAT 和單純型）作對比，并對該方法定位結果的可靠性進行分析。

表1 華南模型Table 1 South China model

1 定位方法簡介

1.1 NLLoc 方法及主要算法

地震定位方法NLLoc 是一種基于三維速度模型的非線性搜索定位方法， Lomax 開發(fā)的定位軟件NonLinLoc 能夠產出欠擬合函數(shù)、 震源位置和邊際后驗概率密度函數(shù)估計（如式 （1））。 該程序基于Tarantola 等[10]的反演方法以及 Tarantola 等[10]、 T.J.Moser 等[11]和 Wittlinger 等[12]的地震定位方法。 假設震相到時觀測誤差（震相拾取誤差）和理論走時誤差符合高斯分布， 在給定的觀測到時和由觀測臺站和空間網(wǎng)格點計算出的理論走時， 這個假設能夠計算發(fā)震時刻的最大概率值。 四維的地震定位問題就轉化為三維空間的網(wǎng)格搜索問題。

式（1）中， k 為歸一化常數(shù)， D 為數(shù)據(jù)空間， d∈D， p（d）為概率密度函數(shù)， p（m）代表震源位置參數(shù)向量 m 的先驗概率密度函數(shù)， F（d，m）為正演問題的概率密度函數(shù)， μ（d，m）通常設為均勻分布。 該式的積分項一般稱作似然函數(shù)， 記為L（m）。 假設觀測數(shù)據(jù)的概率密度函數(shù)p（d）為高斯分布， 均值為d0和協(xié)方差矩陣為Cd， d 和m 的不確定性是可忽略的且是相互獨立的， 則 μ（d，m）=μ（d）·μ（m)，其中μ（d）通常被認為是常數(shù)。 通過簡化， 似然函數(shù)可定義為：

NLLoc 方法中采用等時差（EDT）格式[13-14]的似然函數(shù)， 當出現(xiàn)異常值時， 該格式是穩(wěn)健的。 在等時差情況下， 該似然函數(shù)可簡化為：

式（3）中， X 是m 的空間上分量； 對于兩個觀測臺站，是觀測到時，是理論走時； σa和σb是觀測到時和理論走時的標準差。

1.2 八叉樹搜索算法

NonLinLoc 程序中， 其搜索方法有多種選擇，如八叉樹法、 網(wǎng)格搜索法等， 本文選擇八叉樹搜索法進行定位。 八叉樹搜索法在三維空間中使一種準確、 高效的全局搜索法， 對產生的樣本細胞不斷遞歸細分。 首先， 初始化一個給定的空間，計算每個網(wǎng)格中心點的概率值， 將其概率值放置有序列表LP中， 接著將概率值最大的點剖分成8個子細胞， 計算8 個子細胞的概率值并繼續(xù)放入概率列表LP中， 然后， 從列表中選出概率最大點再次剖分， 再次循環(huán)， 直到滿足給定的終止法則。這是一個快速收斂的遞歸過程， 比網(wǎng)格搜索法快，比模擬退火法和遺傳算法更具有全局性， 但缺點是依賴于初始網(wǎng)格大小和初始點。

圖1 八叉樹搜索采樣單元（http：//alomax.free.fr/nlloc/）Fig.1 Oct-tree searching method sampling cell

每個細胞中震源位置的概率計算近似如下：

式（5）中， Vi是單元體積， xi是細胞中心坐標。

1.3 NonLinLoc 主程序簡介

NonLinLoc 主程序定位流程圖如圖2（此流程圖僅示意NonLinLoc 部分產出結果）， 首先， 輸入一維（1D）速度模型或三維（3D）速度模型生成一個三維網(wǎng)格速度文件， 進而計算三維網(wǎng)格走時文件；走時文件可產出走時圖， 亦可對給定的震源位置Time2EQ 模塊計算預測走時； 輸入事件的震相文應用NLLoc 模塊定位， 即可尋找最優(yōu)的震源位置和發(fā)震時刻， 并產出結果文件和結果展示圖。

圖2 NonLinLoc 主程序流程圖Fig.2 The flow chart of NonLinLoc main program

1.4 福建臺網(wǎng)常用定位方法

福建臺網(wǎng)中心目前采用Jopens 系統(tǒng)MSDP 軟件作為地震定位軟件， 進行臺網(wǎng)地震速報和地震編目等工作。 福建省常用的網(wǎng)內及網(wǎng)緣事件定位方法為： 單純型、 Hyposat、 HYP2000 和 LocSAT等方法。

單純型是利用基礎數(shù)學上的單純形搜索極值，達到終止法則如殘差最小的理論模型作為震源位置， 該方法在極值附近收斂較慢， 對初始值比較敏感[15]， 適用于地方震、 近震和遠震的定位程序。

HYP2000 和 Hyposat 都采用傳統(tǒng) Geiger 法的基本思路， 即觀測方程組降維后直接用奇異值分解最小二乘法方程組， 實際計算中采用多種數(shù)據(jù)加權[16]。 HYP2000 可采用分區(qū)水平分層速度模型，可以為每個地震臺站指定不同的速度模型， 定位起始位置均采用近臺初值， 該法適用于網(wǎng)內震。而Hyposat 不僅可定位地方震和近震， 也可以定位遠震， 其定位效果也相對較好。

LocSAT 采用傳統(tǒng)Geiger 法的基本思路， 應用阻尼最小二乘法， 即將觀測方程組化為正規(guī)方程組后用主元素消去法求解。 該法無法加權， 為了計算初值， 采用水平分層速度模型[17]。 該方法對于地方震、 近震及遠震均可定位。

2 研究使用數(shù)據(jù)情況

2.1 22 次人工爆破事件

2010 年至 2014 年， 福建省地震局 “福建及臺灣海峽深部構造海陸聯(lián)測” 項目實施期間， 共進行22 次人工爆破激發(fā)實驗， 激發(fā)時刻及震中位置信息見表2。

表2 22 次人工爆破信息Table 2 The catalogs of 22 explosions

本文利用這些人工爆破記錄資料， 首先采用MSDP 軟件進行分析。 由于震中距大于150 km 后大部分Pn 震相較弱， 到時位置拾取誤差較大， 因此本文拾取震中距150 km 內的清晰震相 （主要震相為 Pn、 Pg 和少量 Sg）， 并剔除存在時鐘誤差的臺站的震相數(shù)據(jù)。 隨后， 分別采用目前福建臺網(wǎng)常用的 Hyposat、 HYP2000、 LocSAT 和單純型等四種定位程序及NLLoc 定位程序進行地震定位。 福建臺網(wǎng)常用定位方法和NLLoc 定位程序采用的速度模型均為華南模型。 在定位事件時， Hyposat 定位可以設置自定義初始化選擇和深度反演類型；該初始化選項中， 初始化深度設置為10 km， 初始深度誤差設置為10 km， 深度計算設置為 “同時”；LocSAT 設置選項中， 初始深度設置為10 km， 深度計算設置為 “同時”， 迭代次數(shù)為： 40。

2.2 88 次M≥1.5 級仙游地震序列事件

本文選取2010-01-01—2019-09-01 福建仙游震群序列M≥1.5 級天然地震共88 個， 拾取震中距150 km 內的清晰震相， 并剔除存在時鐘誤差、斷記等震相數(shù)據(jù)， 利用NLLoc 方法進行重定位，檢驗其定位天然地震事件的可靠性。

3 定位結果分析

3.1 22 次人工爆破重定位結果分析

3.1.1 定位震中誤差分析

MSDP 自 帶 的 定 位 方 法 Hyposat、 HYP2000、LocSAT、 單純型和NLLoc 方法均使用華南模型，得到的地震定位結果震中分布如圖3。 從表3 可見， 5 種定位方法中， NLLoc 方法定位震中誤差為0.38±0.19 km， 結果最優(yōu)； 使用 Hyp2000 方法震中誤差均值0.72±1.09 km， 定位效果較差。 而單純型和LocSAT 定位震中誤差結果介于二者之間。

3.1.2 發(fā)震時刻誤差分析

如圖4 所示， NLLoc 方法發(fā)震時刻誤差較小，為 0.35±0.13 s。 除了 HYP2000 定位法有 45.45%的事件發(fā)震時刻誤差超過1s， 福建臺網(wǎng)日常定位的其它方法發(fā)震時刻誤差均值均在1s 以內。

3.1.3深度值分析

如圖5 所示， 采用華南模型的NLLoc 方法深度誤差最?。?0.75±1.62 km。 而 MSDP 自帶的定位方法 Hyposat、 HYP2000、 LocSAT 和單純型中， 其深度誤差均值均超過2 km。 NLLoc 方法定位深度誤差優(yōu)于其它方法， 也表明該方法在定位事件深度上具有明顯的優(yōu)勢。

3.2 88 次M≥1.5 級仙游地震序列重定位結果分析

圖3 5 種定位方法的震中位置圖Fig.3 Epicenter location maps of the 5 location methods

表3 震中誤差值 （單位： km）Table 3 Errors of the located epicenters by 5 location methods （unit： km）

圖4 發(fā)震時刻誤差圖Fig.4 Errors map of the origin time

圖5 深度誤差對比圖Fig.5 Comparison of the depth errors

通常天然地震震源深度比人工地震要深得多，天然地震波所通過的路徑也復雜得多， 而人工爆破深度大多數(shù)為幾米至幾百米， 在地表巖層進行，其介質密度低， 爆破源是作用時間很短的點源瞬時膨脹力， 震源體積也相對小很多。 基于這些不同特征， 本文利用NLLoc 方法對88 次M≥1.5 級仙游地震序列重定位結果分析， 檢驗其定位天然地震事件的可靠性。

3.2.1 定位震中誤差分析

如圖 6， 截止至 2019-09-01， 序列活動在時間上存在叢集-平靜的特征。 利用NLLoc 方法重定位， 仙游序列空間分布如圖7， Ⅱ、 Ⅲ和Ⅳ三個時段地震主要呈北西向線性展布， 地震活動空間主體存在自北西向東南方向遷移特征， Ⅴ時段地震向西擴散， 主要分布在北西向線性展布。 該方法定位空間整體分布與袁麗文等[18]、 許振棟等[19]研究結果較一致。

圖6 福建仙游序列M≥1.5 級M-t 圖Fig.6 M-t diagram of Xianyou earthquake sequence with M≥1.5

圖7 仙游震群序列M≥1.5 空間分布Fig.7 Distribution map of Xianyou earthquake swarm sequence with M≥1.5

利用NLLoc 方法重定位結果與采用Hyposat 或單純型定位的中國臺網(wǎng)正式目錄對比， 如圖8， 震中誤差為： 0.66±0.37 km， 有 2 個地震事件誤差超過1.5 km， 分析這兩個事件均為多個事件重疊，清晰震相少且臺數(shù)少， 定位效果不佳。

3.2.2 發(fā)震時刻誤差分析

重定位的發(fā)震時刻誤差為： 0.16±0.15 s， 誤差小。

3.2.3 深度值分析

在震源深度方面如表4， 利用NLLoc 方法重定位結果中， 地震序列活動初期深度較淺5 km 左右， 與李強等[20]利用 CAP 反演 2012 年 04 月 15 日4.1 級地震的震源深度4 km 一致， 但正式目錄中，地震序列活動初期深度較深， 均值15 km 左右，與CAP 反演震源深度差值較大。 2012.11-2013.3逐漸加深趨勢， 深度均值7.4 km， 2013 年8 月4.2級地震發(fā)生后， 地震震源深度分布穩(wěn)定， 集中分布在 9 km 左右（如圖 9）。

圖8 福建仙游震群序列M≥1.5 震中誤差和發(fā)震時刻誤差圖Fig.8 The epicenter error and origin time error of Xianyou earthquake swarm sequence with M≥1.5 in Fujian

圖9 仙游地震序列M≥1.5 深度時序圖Fig.9 The depth sequence diagram of Xianyou earthquake swarm sequence with M≥1.5

表4 深度均值和標準差 （單位： km）Table 4 the mean and standard deviation of depth （unit： km）

4 討論

本文利用三維非線性NLLoc 定位方法對 “福建及臺灣海峽地殼深部構造海陸聯(lián)測” 項目實施期間的 22 次人工爆破事件重新定位， 并與Hyposat、 HYP2000、 LocSAT 和單純型這四種常用定位算法進行對比， 討論了福建臺網(wǎng)應用NLLoc程序測定人工爆破事件參數(shù)的精度問題。 研究結果表明， NLLoc 方法定位結果在發(fā)震時刻誤差、 震中誤差和震源深度誤差等方面都優(yōu)于現(xiàn)有定位方法結果， 尤其是在震源深度確定方面NLLoc 方法優(yōu)勢明顯。 隨后， 選取 2010 年 1 月至 2019 年 09月福建仙游震群序列M≥1.5 級共88 個天然地震，利用NLLoc 方法重定位， 與中國臺網(wǎng)正式目錄對比， 結果顯示： 深度方面， 地震序列活動初期深度較淺， 逐漸加深趨勢， 集中在 9 km 左右， 震中位置呈北西向線性展布； 但地震序列活動正式目錄的深度初期較深， 與CAP 反演震源深度差值較大。

在震源深度方面， 如果某一臺距離震中較近或就在震中區(qū)， 則該臺震源距與震源深度相近[21]，根據(jù)震源距公式 D=Vφ×Δt， 其中 Vφ為虛波速度，Vφ=（VP*Vs）/（VP-Vs）， Δt 為該臺的 S 波與 P 波到時差。 福建仙游序列事件中能記錄到的最近臺站為“福建仙游石蒼臺 （25.62°N， 118.74°E） ”， 震中距大約0.02°， 該臺的震源距可近似為震源深度。 統(tǒng)計該臺的觀測記錄， 獲得到時差約為1.2s， 結合表1， 推算該臺震源距為10.4 km， 則震源深度應略小于 10.4 km， 如圖 9 所示， 利用 NLLoc 方法定位的震源深度略小于理論估值。 因此， 本文研究認為無論是人工爆破還是天然地震， NLLoc 方法有助于更好地確定事件深度， 有望提高臺網(wǎng)定位精度。

此外， 還存在一些需要進一步討論和研究的細節(jié)。 首先， 擁有已知發(fā)震時刻和絕對位置的爆破的數(shù)據(jù)量有限， 并且選取震中距150 km 以內的臺站作為定位事件的界限僅僅是根據(jù)經驗來確定：震中距在150 km 內， 震相較清晰， 參與定位臺站數(shù)較多； 遠臺震相不清晰， 參與定位殘差較大，有可能會影響NLLoc 定位效果。 第二， 對于網(wǎng)緣（網(wǎng)外）地震的定位問題尚未涉及。 本文的22 次爆破和福建仙游震群序列M≥1.5 事件均為陸上網(wǎng)內事件， 震中包圍相對較好， 定位效果較好。 今后將有針對性研究網(wǎng)外事件， 進一步探討該定位程序對該類事件的定位效果。 第三， 該定位方法輸入三維速度模型時， 存在依賴于初始值的局限性，初始值偏離震中位置較遠， 定位效果不佳。 若今后使用三維速度模型定位， 將解決依賴初始值的問題。

5 結語

通過本文研究結果表明， NLLoc 算法無論是人工爆破還是天然地震， 其定位結果在發(fā)震時刻誤差、 震中誤差和震源深度誤差等方面都優(yōu)于現(xiàn)有定位方法結果， 尤其是在震源深度確定方面NLLoc算法優(yōu)勢明顯。 該法可用于臺網(wǎng)日常地震定位，有助于更好地解決震源深度測定問題， 提高地震定位精度。 今后， 將進一步探究網(wǎng)緣（網(wǎng)外）人工爆破和天然事件的定位效果； 該方法為 “福建及臺灣海峽深部構造海陸聯(lián)測” 項目獲取閩臺三維地殼速度模型提供地震定位方法。