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

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

0 引言

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

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

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

表1 華南模型Table 1 South China model

1 定位方法簡(jiǎn)介

1.1 NLLoc 方法及主要算法

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

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

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

式（3）中， X 是m 的空間上分量； 對(duì)于兩個(gè)觀測(cè)臺(tái)站，是觀測(cè)到時(shí)，是理論走時(shí)； σa和σb是觀測(cè)到時(shí)和理論走時(shí)的標(biāo)準(zhǔn)差。

1.2 八叉樹搜索算法

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

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

每個(gè)細(xì)胞中震源位置的概率計(jì)算近似如下：

式（5）中， Vi是單元體積， xi是細(xì)胞中心坐標(biāo)。

1.3 NonLinLoc 主程序簡(jiǎn)介

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

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

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

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

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

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

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

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

2.1 22 次人工爆破事件

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

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

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

2.2 88 次M≥1.5 級(jí)仙游地震序列事件

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

3 定位結(jié)果分析

3.1 22 次人工爆破重定位結(jié)果分析

3.1.1 定位震中誤差分析

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

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

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

3.1.3深度值分析

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

3.2 88 次M≥1.5 級(jí)仙游地震序列重定位結(jié)果分析

圖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ā)震時(shí)刻誤差圖Fig.4 Errors map of the origin time

圖5 深度誤差對(duì)比圖Fig.5 Comparison of the depth errors

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

3.2.1 定位震中誤差分析

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

圖6 福建仙游序列M≥1.5 級(jí)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 方法重定位結(jié)果與采用Hyposat 或單純型定位的中國(guó)臺(tái)網(wǎng)正式目錄對(duì)比， 如圖8， 震中誤差為： 0.66±0.37 km， 有 2 個(gè)地震事件誤差超過(guò)1.5 km， 分析這兩個(gè)事件均為多個(gè)事件重疊，清晰震相少且臺(tái)數(shù)少， 定位效果不佳。

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

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

3.2.3 深度值分析

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

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

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

表4 深度均值和標(biāo)準(zhǔn)差 （單位： km）Table 4 the mean and standard deviation of depth （unit： km）

4 討論

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

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

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

5 結(jié)語(yǔ)

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