基于背景噪聲互相關(guān)函數(shù)的地震觀測密集臺(tái)陣時(shí)鐘偏差分析方法

田原 王偉濤 李麗 于常青 張海明

基于背景噪聲互相關(guān)函數(shù)的地震觀測密集臺(tái)陣時(shí)鐘偏差分析方法

田原1,2王偉濤2,?李麗2于常青3張海明4

1.中國電子科學(xué)研究院, 北京 100041; 2.中國地震局地球物理研究所, 北京 100081; 3.中國地質(zhì)科學(xué)院地質(zhì)研究所, 北京 100037; 4.北京大學(xué)地球與空間科學(xué)學(xué)院, 北京 100871; ?通信作者, E-mail: wangwt@cea-igp.ac.cn

利用地震觀測臺(tái)陣中相鄰臺(tái)站間背景噪聲互相關(guān)函數(shù)隨時(shí)間的偏移的特性, 考慮參考互相關(guān)函數(shù)、噪聲源變化、反演方法和數(shù)據(jù)誤差的影響, 給出基于背景噪聲互相關(guān)函數(shù)的密集臺(tái)陣時(shí)鐘偏差分析及其誤差評價(jià)方法。對鹽源短周期地震觀測密集臺(tái)陣的單臺(tái)時(shí)鐘偏差分析結(jié)果表明, 基于背景噪聲互相關(guān)函數(shù)的時(shí)鐘偏差分析方法可以給出連續(xù)的單臺(tái)時(shí)鐘偏差, 并能夠篩選出有明顯時(shí)鐘偏差的臺(tái)站及偏差出現(xiàn)的時(shí)段。該臺(tái)陣的 209 個(gè)臺(tái)站中, 有 17 個(gè)臺(tái)站在觀測期內(nèi)出現(xiàn)大于 1s 的時(shí)鐘偏差, 可能與儀器數(shù)據(jù)采集裝置的硬件或軟件故障有關(guān)。

地震觀測密集臺(tái)陣; 背景噪聲互相關(guān)函數(shù); 時(shí)鐘偏差分析; 噪聲源分布

地震波形的數(shù)據(jù)質(zhì)量評估是使用地震數(shù)據(jù)進(jìn)行科學(xué)研究的必要前提。隨著地震觀測密集臺(tái)陣(簡稱密集臺(tái)陣)數(shù)據(jù)處理和應(yīng)用技術(shù)(簡稱臺(tái)陣技術(shù))的快速發(fā)展, 已有大量密集臺(tái)陣數(shù)據(jù)應(yīng)用于高精度地下結(jié)構(gòu)成像和地震定位等對波形記錄的時(shí)間準(zhǔn)確性要求很高的研究中。具備準(zhǔn)確和統(tǒng)一的時(shí)間軸, 是應(yīng)用大規(guī)模密集臺(tái)陣數(shù)據(jù)進(jìn)行深度科學(xué)研究的前提之一。

現(xiàn)今的地震儀大多使用數(shù)據(jù)采集器自動(dòng)連接全球?qū)Ш较到y(tǒng)(GPS)進(jìn)行授時(shí), GPS 信號斷開時(shí)則使用傳感器內(nèi)部時(shí)鐘進(jìn)行授時(shí)。在儀器運(yùn)轉(zhuǎn)正常且 GPS信號連接穩(wěn)定時(shí), 系統(tǒng)的時(shí)間精度可以達(dá)到 50ns, 能夠保證數(shù)據(jù)授時(shí)的精度。授時(shí)系統(tǒng)在兩種情況下可能產(chǎn)生較大的偏差: 1)GPS 信號連接失敗, 導(dǎo)致數(shù)據(jù)采集系統(tǒng)只能使用傳感器內(nèi)部時(shí)鐘, 此時(shí)記錄的波形可能產(chǎn)生連續(xù)的時(shí)間漂移, 直至 GPS 信號重新連接成功[1?2]; 2)由于采集或授時(shí)系統(tǒng)的硬件或軟件出現(xiàn)故障, 數(shù)據(jù)采集系統(tǒng)的時(shí)間可能會(huì)出現(xiàn)較大的變化[3?4]。無論哪種情況產(chǎn)生的時(shí)間偏差, 數(shù)據(jù)使用者都無法獲得準(zhǔn)確的時(shí)間校正信息。

在無法獲知準(zhǔn)確時(shí)間信息的情況下, 可以通過記錄的波形數(shù)據(jù)對臺(tái)站時(shí)鐘進(jìn)行簡單的分析。利用波形的時(shí)鐘偏差分析方法主要有兩類。1)通過地震震相在各臺(tái)站的到時(shí), 估計(jì)臺(tái)站的時(shí)鐘偏差。例如, Anchieta 等[5]利用海底地震儀與臨近海島上的陸地寬頻帶地震儀的 P 波到時(shí)差, 估計(jì)海底地震儀的時(shí)鐘偏差。此類方法依賴于地震震相的清晰度和出現(xiàn)頻率, 因此僅適用于地震頻發(fā)區(qū)域, 對于地震稀少的地區(qū), 不能保證時(shí)鐘的不間斷校正。2)利用臺(tái)站間的相干噪聲信號, 檢測時(shí)鐘偏差。此類方法可用于檢測小孔徑臺(tái)陣、海底臺(tái)陣及區(qū)域臺(tái)網(wǎng)中的臺(tái)站時(shí)鐘偏差。例如, Stehly 等[3]利用地球微震頻帶噪聲信號(5~10s 和 10~20s), 對固定臺(tái)站的時(shí)間偏差進(jìn)行連續(xù)檢測; Sens-Sch?nfelder[1]利用兩個(gè)小孔徑火山監(jiān)測臺(tái)陣的每日互相關(guān)函數(shù)相對于全時(shí)段疊加互相關(guān)函數(shù)的偏移量, 給出每個(gè)檢波器的單臺(tái)時(shí)間日變化, 并指出該臺(tái)陣的時(shí)間漂移主要由 GPS 授時(shí)失敗期間的內(nèi)部時(shí)鐘漂移產(chǎn)生; Xia 等[4]利用幾內(nèi)亞灣的 26s 固定噪聲源信號, 分析全球范圍的臺(tái)網(wǎng)時(shí)鐘偏差; 王俊等[6]對江蘇省的固定臺(tái)網(wǎng)進(jìn)行連續(xù)的時(shí)鐘分析, 并進(jìn)行時(shí)鐘偏差的誤差分析。

隨著臺(tái)陣技術(shù)的進(jìn)一步發(fā)展, 人們開始嘗試布設(shè)密集的短時(shí)觀測臺(tái)陣(臺(tái)間距數(shù)十米至數(shù)公里, 觀測周期十幾天至數(shù)月), 利用噪聲中的高頻成分, 提取高頻(5s 以下)面波乃至體波數(shù)據(jù), 以期研究局部小尺度淺層地殼乃至近地表的精細(xì)結(jié)構(gòu)[7?10]。相對于中長周期臺(tái)陣, 短周期密集臺(tái)陣具有臺(tái)間距小、總臺(tái)站數(shù)大和總觀測時(shí)間短等特點(diǎn), 數(shù)據(jù)對時(shí)間精度要求較高, 因此需要一種能夠同時(shí)檢測觀測時(shí)段內(nèi)所有臺(tái)站時(shí)間精度的快速檢測方法。

本文采用類似 Sens-Sch?nfelder[1]提出的基于臺(tái)陣間噪聲互相關(guān)函數(shù)的單臺(tái)時(shí)間偏差分析方法, 針對短周期密集觀測臺(tái)陣的特點(diǎn), 給出適用于它的臺(tái)陣時(shí)間偏差分析方法, 并對鹽源短周期地震觀測密集臺(tái)陣(簡稱鹽源臺(tái)陣)進(jìn)行時(shí)間偏差分析, 探討參考互相關(guān)函數(shù)計(jì)算誤差、矩陣優(yōu)化方案、數(shù)據(jù)讀取誤差及噪聲源方向改變等因素對分析結(jié)果產(chǎn)生的影響, 最終給出鹽源臺(tái)陣所有臺(tái)站的時(shí)鐘偏差以及偏差的誤差估計(jì)。

1 方法原理

1.1 基本原理

在理想狀態(tài)下, 噪聲源均勻地分布, 兩個(gè)臺(tái)站之間的背景噪聲互相關(guān)函數(shù)代表兩臺(tái)站間的經(jīng)驗(yàn)格林函數(shù), 與兩臺(tái)站間的路徑結(jié)構(gòu)有關(guān), 不隨時(shí)間變化[11?12]。以下因素可能會(huì)導(dǎo)致互相關(guān)函數(shù)的變化: 1)臺(tái)站間路徑結(jié)構(gòu)發(fā)生變化; 2)某個(gè)臺(tái)站的時(shí)鐘出現(xiàn)偏差, 或數(shù)據(jù)采集器出現(xiàn)相位偏差; 3)噪聲源分布發(fā)生變化。

可將某一時(shí)段內(nèi)臺(tái)站間互相關(guān)函數(shù)相對于參考互相關(guān)函數(shù)的時(shí)間偏差用如下式[3]表達(dá):

=()+()+(), (1)

其中,()代表路徑結(jié)構(gòu)變化導(dǎo)致的時(shí)間偏差,()代表由臺(tái)站時(shí)鐘偏差或相位誤差導(dǎo)致的時(shí)間偏差,()表示由噪聲源分布變化產(chǎn)生的時(shí)間偏差。短周期密集臺(tái)陣的觀測時(shí)間通常小于 3 個(gè)月, 可以認(rèn)為觀測期內(nèi)臺(tái)站間的路徑結(jié)構(gòu)幾乎不隨時(shí)間變化, 因此相對于,()可以忽略不計(jì)。此時(shí), 臺(tái)陣中任意兩臺(tái)站與間的時(shí)鐘(或相位)偏差可以表示為

D() = ?() – ?() =δt–ε(), (2)

其中, ?()表示第個(gè)臺(tái)站的單臺(tái)時(shí)鐘偏差,()表征由于噪聲源分布的變化在兩臺(tái)站間產(chǎn)生的時(shí)間偏移。若噪聲源在空間分布均勻, 或空間分布雖不均勻但隨時(shí)間變化不大時(shí), 相對于,()的影響可忽略, 此時(shí)幾乎全部由兩個(gè)臺(tái)站各自的單臺(tái)時(shí)間偏差貢獻(xiàn)。若噪聲源的優(yōu)勢方向在一段時(shí)間內(nèi)發(fā)生明顯的改變, 則()的影響不可忽略, 需要想辦法予以修正。使用較長觀測時(shí)間(1 個(gè)月以上)的互相關(guān)函數(shù)疊加, 可以改善噪聲源的分布情況, 進(jìn)而減少()的影響, 這種做法多應(yīng)用于對長時(shí)間觀測的固定臺(tái)站進(jìn)行的時(shí)鐘偏差分析中[2?3,6], 但會(huì)在一定程度上降低計(jì)算結(jié)果的時(shí)間分辨率。本文利用密集觀測臺(tái)陣的特性, 擬合噪聲源變化在不同方向的臺(tái)站間產(chǎn)生的時(shí)間變化, 并加以去除, 盡可能地減小噪聲源分布變化的影響, 同時(shí)保證計(jì)算結(jié)果的時(shí)間分辨率。

對于整個(gè)臺(tái)陣中的共個(gè)臺(tái)站, 通過臺(tái)陣中所有個(gè)臺(tái)站對互相關(guān)函數(shù)的時(shí)間偏移量(式(2)), 構(gòu)建求解單臺(tái)時(shí)鐘偏差的線性方程組:

式(3)中, 線性方程組的系數(shù)矩陣為不滿秩矩陣, 需要至少選擇臺(tái)陣中的某一個(gè)臺(tái)站, 認(rèn)為它的時(shí)鐘相對準(zhǔn)確, 即?()=0, 此時(shí)式(3)化為典型的超定方程組求解問題, 可以通過最小二乘法求得數(shù)據(jù)誤差最小化的穩(wěn)定解。

1.2 臺(tái)站對選取及互相關(guān)頻帶選擇

臺(tái)站對的選擇取決于密集臺(tái)陣的布設(shè)方式。距離越近的臺(tái)站, 互相關(guān)函數(shù)受介質(zhì)及噪聲源分布變化的影響越小, 也就越能精確地計(jì)算臺(tái)站時(shí)間偏差對互相關(guān)函數(shù)的影響[1]。每個(gè)臺(tái)站應(yīng)擁有足夠的能與之組成臺(tái)站對的臺(tái)站, 同時(shí), 為了更精確地去除噪聲源變化對時(shí)間偏差估計(jì)的影響, 臺(tái)站對應(yīng)在各個(gè)方向上較均勻地分布, 或者有較明顯的優(yōu)勢方向。

利用背景噪聲分析時(shí)鐘偏差時(shí), 計(jì)算互相關(guān)函數(shù)的頻率范圍取決于臺(tái)站間的距離及臺(tái)陣布設(shè)區(qū)域背景噪聲場的強(qiáng)度和穩(wěn)定性。本文通過計(jì)算所有臺(tái)站的單臺(tái)噪聲功率譜密度(power spectrum density, PSD)[13]來考察臺(tái)陣的背景噪聲強(qiáng)度在不同頻帶的穩(wěn)定性, 選擇強(qiáng)度大且變化幅度小的頻帶進(jìn)行互相關(guān)計(jì)算。除單臺(tái)背景噪聲強(qiáng)度外, 還可根據(jù)單日互相關(guān)函數(shù)與觀測期內(nèi)疊加互相關(guān)函數(shù)之間的波形相似度來選擇互相關(guān)頻帶, 相似度高說明噪聲源穩(wěn)定, 相似度過低則說明單日的噪聲源與觀測期內(nèi)的疊加噪聲源有較大差別, 不適宜進(jìn)行時(shí)鐘偏差分析。實(shí)際操作時(shí), 選擇的互相關(guān)頻帶應(yīng)滿足如下條件: 利用該頻帶進(jìn)行計(jì)算時(shí), 臺(tái)陣中每個(gè)臺(tái)站對都有 90%以上的單日互相關(guān)函數(shù)與疊加互相關(guān)數(shù)的波形互相關(guān)系數(shù)大于0.9。

1.3 互相關(guān)函數(shù)相對時(shí)間偏差的計(jì)算

首先按照 Sens-Sch?nfelder[1]的方法, 計(jì)算臺(tái)陣中每個(gè)臺(tái)站對每天的互相關(guān)函數(shù)相對于參考互相關(guān)函數(shù)的時(shí)間偏移量, 操作步驟如下。

3)利用步驟 2 求得的時(shí)間偏移量, 對每日互相關(guān)函數(shù)進(jìn)行時(shí)間修正, 并利用修正后的每日互相關(guān)函數(shù), 重復(fù)步驟 1, 得到修正的參考互相關(guān)函數(shù)。

4)重復(fù)步驟 2, 得到相對于修正參考互相關(guān)函數(shù)的時(shí)間偏移量。

1.4 臺(tái)站對參考互相關(guān)函數(shù)的時(shí)間偏差修正

對于短時(shí)觀測的臺(tái)陣, 利用觀測期內(nèi)互相關(guān)函數(shù)疊加計(jì)算求得的參考互相關(guān)函數(shù)有時(shí)不能代表真正的“正確時(shí)間”[2]。若存在貫穿觀測期的明顯時(shí)鐘或相位偏差, 則簡單疊加后的參考互相關(guān)函數(shù)不能代表臺(tái)站的“正確”時(shí)鐘, 此時(shí)需要對這類情況進(jìn)行修正。

引入絕對時(shí)間偏移量后, 式(2)改寫為

1.5 噪聲源分布改變的影響修正

1.6 單臺(tái)時(shí)鐘偏差計(jì)算及誤差分析

單臺(tái)時(shí)鐘偏差分析結(jié)果的誤差由兩部分組成: 1)線性方程組求解過程中引入的誤差; 2)數(shù)據(jù)誤差造成的結(jié)果誤差, 即臺(tái)站對時(shí)鐘偏差D()的誤差對單臺(tái)時(shí)鐘偏差的影響。在線性方程組的求解過程中, 參考臺(tái)站的選取可能對最終的反演結(jié)果造成影響。為確定此影響的程度, 我們放寬臺(tái)站選擇的第3 個(gè)條件, 選取滿足(D())

本文中反演數(shù)據(jù)(即各個(gè)臺(tái)站對的時(shí)鐘偏差D())的誤差估計(jì)由 3 個(gè)部分組成:

err=errdata+errref+errsrc, (7)

本文使用自舉法估計(jì)數(shù)據(jù)誤差對單臺(tái)時(shí)間偏差反演結(jié)果的影響。在每個(gè)臺(tái)站對的時(shí)鐘偏差D()中加入隨機(jī)誤差(0,), 隨機(jī)誤差的標(biāo)準(zhǔn)差取每個(gè)臺(tái)站對數(shù)據(jù)誤差的一半。我們進(jìn)行100次加入隨機(jī)誤差的線性反演, 并且計(jì)算所有反演結(jié)果的標(biāo)準(zhǔn)差err。最終, 單臺(tái)時(shí)間偏差反演結(jié)果的誤差由數(shù)據(jù)誤差和反演過程誤差組成, 即

ERR?=2(sta+err)。 (8)

2 數(shù)據(jù)與應(yīng)用

作為本文應(yīng)用示例的短周期密集臺(tái)陣數(shù)據(jù)來自由中國地質(zhì)科學(xué)院和中國地震局地球物理勘探中心布設(shè)于中國西南部鹽源盆地及其周邊山區(qū)的短周期密集臺(tái)陣, 臺(tái)站總數(shù)為209個(gè), 以約2km的臺(tái)間距密集地布設(shè)于約40km×60km的菱形區(qū)域內(nèi)(圖1)。臺(tái)陣使用的儀器全部為EPS-2型短周期三分量地震儀, 頻帶寬度為0.2~150Hz, 采樣頻率為200 Hz。臺(tái)陣觀測時(shí)間為2017年6月4日至2017年7月10日, 單臺(tái)最長觀測時(shí)間為36天, 最短18天。

本文采用Bensen等[15]的背景噪聲處理方法, 對背景噪聲記錄數(shù)據(jù)進(jìn)行預(yù)處理。首先, 將臺(tái)站原始記錄的采樣頻率降至50Hz, 并對每天記錄完整度大于80%的數(shù)據(jù)進(jìn)行去均值和去趨勢處理, 再進(jìn)行去除儀器響應(yīng)的操作。然后, 對每個(gè)長度為1天的數(shù)據(jù), 取每段長20分鐘、移動(dòng)步長10分鐘的小時(shí)窗, 在每個(gè)時(shí)窗中使用剪切閾值法進(jìn)行時(shí)間域歸一化處理, 并在0.05~10Hz頻率范圍內(nèi)進(jìn)行譜白化處理。最后, 對臺(tái)站間所有長度為20分鐘的時(shí)窗數(shù)據(jù)進(jìn)行頻率域互相關(guān)計(jì)算, 并疊加成為兩臺(tái)站間的單日噪聲互相關(guān)函數(shù)。

鹽源臺(tái)陣最小臺(tái)站間距為2km, 本文選取臺(tái)間距5km以內(nèi)的681個(gè)臺(tái)站對(方位分布見圖2)參與臺(tái)站時(shí)鐘偏差分析的計(jì)算。除臺(tái)陣邊緣的臺(tái)站以及少量分布不均勻的臺(tái)站外, 大部分臺(tái)站可以與距離最近的8個(gè)臺(tái)站進(jìn)行互相關(guān)計(jì)算。統(tǒng)計(jì)時(shí), 按照臺(tái)陣分布特點(diǎn), 臺(tái)站對的反方位角BAZ∈[0°, 170°]時(shí), 臺(tái)站對的方向=BAZ; BAZ∈[170°, 350°]時(shí), 互相關(guān)波形相對于時(shí)間零點(diǎn)對稱分布后,=BAZ? 180°; BAZ∈[350°, 360°]時(shí),=BAZ?360°。所有臺(tái)站對共有4個(gè)明顯的優(yōu)勢連線方向。

鹽源臺(tái)陣的單臺(tái)背景噪聲強(qiáng)度如圖3(a)所示?？梢钥闯? 臺(tái)陣的單日PSD在頻率大于1Hz時(shí)變化幅度較大, 可達(dá)40dB, 而頻率小于1Hz時(shí)整體強(qiáng)度較穩(wěn)定。因此, 從能量的角度考慮, 應(yīng)選擇1 Hz以下的低頻帶背景噪聲來進(jìn)行互相關(guān)計(jì)算。圖3(b)給出在3個(gè)不同的頻帶范圍(1~2, 0.5~1和0.1~ 0.5Hz), 鹽源臺(tái)陣中684個(gè)臺(tái)間距5km以內(nèi)的臺(tái)站對的每日互相關(guān)函數(shù)與觀測期內(nèi)疊加互相關(guān)函數(shù)之間的相關(guān)系數(shù)?？梢钥吹? 在0.1~0.5Hz頻帶內(nèi), 互相關(guān)函數(shù)相似度非常高, 大部分互相關(guān)系數(shù)均大于0.9; 0.5~1Hz頻帶內(nèi)的互相關(guān)系數(shù)則在0.3~0.9之間, 峰值互相關(guān)系數(shù)為0.75; 1~2Hz頻帶內(nèi)的互相關(guān)函數(shù)相似度很低, 相關(guān)系數(shù)峰值僅為 0.3左右。上述結(jié)果說明, 鹽源臺(tái)陣的相干噪聲場在0.5Hz以上的高頻帶非常不穩(wěn)定。為了在最大程度上保證通過背景噪聲互相關(guān)函數(shù)計(jì)算得出的時(shí)間偏差的準(zhǔn)確性, 我們選取0.1~0.5Hz作為背景噪聲互相關(guān)計(jì)算的頻帶范圍。

圖4展示兩個(gè)不同臺(tái)站對原始記錄的每日互相關(guān)函數(shù)、初始參考互相關(guān)函數(shù)、修正每日互相關(guān)函數(shù)、修正參考互相關(guān)函數(shù)及其對應(yīng)的互相關(guān)系數(shù)?？梢钥吹? 時(shí)鐘可能出現(xiàn)偏移的臺(tái)站對1329-1361, 修正后的參考互相關(guān)函數(shù)形態(tài)發(fā)生改變, 對應(yīng)的每日互相關(guān)函數(shù)與參考互相關(guān)函數(shù)的相關(guān)系數(shù)也大幅度提升。另一個(gè)臺(tái)站對1501-1533, 第26天起數(shù)據(jù)出現(xiàn)明顯錯(cuò)誤, 互相關(guān)系數(shù)大幅度下降。在后續(xù)計(jì)算過程中, 我們只保留互相關(guān)系數(shù)大于0.5的數(shù)據(jù)。

本文對鹽源臺(tái)陣的臺(tái)站對參考互相關(guān)函數(shù)進(jìn)行時(shí)間偏差修正。圖5(a)展示該臺(tái)陣中一個(gè)臺(tái)站1377與其周圍8個(gè)臺(tái)站的參考互相關(guān)函數(shù)。8個(gè)臺(tái)站對中, 有6個(gè)臺(tái)站對的參考互相關(guān)函數(shù)的振幅極大值位置在時(shí)間軸(橫軸)的零點(diǎn)位置附近, 偏離零點(diǎn)的最大值僅為0.26s, 符合我們對較小臺(tái)間距(5km以內(nèi))及較低頻帶范圍(0.1~0.5Hz)條件下參考互相關(guān)函數(shù)形態(tài)的預(yù)期。但是, 存在兩個(gè)異常的臺(tái)站對(1377-1411, 1377-1345), 其參考互相關(guān)函數(shù)的極值位置偏離零點(diǎn) 2s 以上, 遠(yuǎn)大于其余臺(tái)站對, 可以斷定這兩個(gè)臺(tái)站對的參考互相關(guān)函數(shù)出現(xiàn)明顯偏差, 臺(tái)站對中的某個(gè)臺(tái)站可能在整個(gè)觀測期內(nèi)都存在明顯的時(shí)鐘或相位偏差, 此時(shí)參考互相關(guān)函數(shù)的時(shí)間軸已不能代表臺(tái)站對的“正確”時(shí)鐘。

圖5還展示臺(tái)陣中681個(gè)臺(tái)站對的相對時(shí)間偏移量(圖5(b)和(c))以及絕對時(shí)間偏移量(圖 5(d))。為了繪圖方便, 圖5(b)中統(tǒng)一將臺(tái)站對中編號較小的臺(tái)站作為源臺(tái)站, 用起點(diǎn)在源臺(tái)站處的射線代表臺(tái)站對的極值偏移量, 射線長度代表該臺(tái)站對極值偏移量的絕對值, 當(dāng)臺(tái)站間極值偏移量大于0時(shí), 射線方向與臺(tái)站對連線方向一致, 反之則與臺(tái)站對連線方向相反。在鹽源臺(tái)陣中, 共有107個(gè)臺(tái)站對存在絕對時(shí)間偏移量較大的情況, 說明在鹽源臺(tái)陣的時(shí)鐘偏差分析過程中, 參考互相關(guān)函數(shù)的偏差修正是必要的。

圖6展示2017年6月12日相對時(shí)間偏移量隨臺(tái)站對方位的分布、去除噪聲源響應(yīng)的相對時(shí)間偏移量以及修正噪聲源前后單臺(tái)時(shí)鐘偏差的反演結(jié)果, 可以看到噪聲源分布的優(yōu)勢方向發(fā)生明顯的改變, 東南方向的臺(tái)站對產(chǎn)生最大的時(shí)間偏移, 偏移量超過0.3s。這類噪聲源響應(yīng)會(huì)使由式(3)得出的單臺(tái)時(shí)間偏差結(jié)果產(chǎn)生非常明顯的系統(tǒng)性偏差。經(jīng)過噪聲源響應(yīng)修正后, 大部分臺(tái)站對的相對時(shí)間偏移量大幅度減小, 反演結(jié)果中的系統(tǒng)性偏差也得到修正。位于東南角的部分西北?東南方向臺(tái)站對, 經(jīng)過噪聲源響應(yīng)修正后仍殘留少量偏差, 說明此處可能存在一些較明顯的局部噪聲變化。

鹽源臺(tái)陣的時(shí)鐘偏差反演結(jié)果的誤差如圖 7 所示。由式(8)可知, 反演誤差由反演求解造成的誤差和數(shù)據(jù)造成的誤差兩個(gè)部分組成。反演求解造成的誤差主要來源于參考臺(tái)站的選擇。本研究選取滿足(D())<0.02 (即臺(tái)站對時(shí)鐘偏差的標(biāo)準(zhǔn)差小于采樣間隔)的臺(tái)站作為參考臺(tái)站, 鹽源臺(tái)陣中滿足此要求的參考臺(tái)站共有 23 個(gè)。我們選取這 23 個(gè)臺(tái)站中任意一個(gè)作為參考臺(tái)站, 獲得不同的反演結(jié)果, 并計(jì)算結(jié)果的標(biāo)準(zhǔn)差sta(圖 7(a))。可以看出, 除第一天因參與計(jì)算的臺(tái)站較少而使得反演結(jié)果的標(biāo)準(zhǔn)差較大外, 其余時(shí)間改變參考臺(tái)站后的反演結(jié)果標(biāo)準(zhǔn)差均小于 0.03s, 說明變換參考臺(tái)站對反演結(jié)果的影響很小。

按照 1.6 節(jié)的方法, 計(jì)算數(shù)據(jù)誤差造成的反演結(jié)果標(biāo)準(zhǔn)差err, 結(jié)果如圖 7(b)所示?？梢钥闯? 由數(shù)據(jù)誤差造成的反演結(jié)果標(biāo)準(zhǔn)差整體上小于 0.1s, 僅部分臺(tái)站和部分日期達(dá)到 0.2s 左右。臺(tái)陣邊緣的臺(tái)站、參考互相關(guān)函數(shù)有異常的臺(tái)站以及噪聲源異常的臺(tái)站數(shù)據(jù)誤差較大, 由其引起的反演標(biāo)準(zhǔn)差也較大。臺(tái)陣的單臺(tái)時(shí)鐘偏差反演的最終誤差如圖7(c)~(d)所示?？梢钥闯? 反演結(jié)果的誤差在大部分臺(tái)站的大部分時(shí)間內(nèi)都小于 0.1s, 臺(tái)陣邊緣臺(tái)站的總體誤差較大, 單臺(tái)誤差平均值達(dá)到 0.2s 以上。在觀測期內(nèi), 6 月 5 日和 7 月 3 日所有臺(tái)站的反演誤差均較大, 可能是由觀測開始階段和結(jié)束階段參與計(jì)算的臺(tái)站數(shù)量較少導(dǎo)致。臺(tái)陣右下角處臺(tái)站的總體誤差較大, 可能與局部噪音源變化強(qiáng)烈(圖 6(b))有關(guān)。綜上所述, 時(shí)間偏差結(jié)果的準(zhǔn)確度受臺(tái)站在臺(tái)陣中的位置、參與計(jì)算的臺(tái)站數(shù)量以及背景噪聲分布穩(wěn)定性等因素的綜合影響。

3 結(jié)果及討論

鹽源臺(tái)陣的單臺(tái)時(shí)間偏差反演結(jié)果及其誤差分析結(jié)果如圖 8 所示。鹽源臺(tái)陣的 209 個(gè)臺(tái)站中, 有17 個(gè)臺(tái)站出現(xiàn)比較明顯的時(shí)鐘偏差, 其中有 15 個(gè)臺(tái)站的時(shí)間始終比正常時(shí)間提前 2s 以上, 有 1 個(gè)臺(tái)站(1349)的時(shí)間始終比正常時(shí)間提前 1s 以上, 還有一個(gè)臺(tái)站(1329)在 2017 年 6 月 5—19 日的時(shí)間比正常提前 2 s 左右, 自 6 月 20 日起, 時(shí)間回到零點(diǎn)附近的正常范圍。所有時(shí)鐘偏差均為躍變式, 可能與數(shù)據(jù)采集器內(nèi)的軟件或硬件故障有關(guān), 具體原因有待對儀器進(jìn)行檢測后核實(shí)。鹽源臺(tái)陣中未出現(xiàn)類似 Sens-Sch?nfelder[1]和 Gouédard 等[2]的研究中出現(xiàn)單臺(tái)時(shí)鐘偏差隨時(shí)間線性增加或減少的情況, 可能說明鹽源臺(tái)陣的臺(tái)站未出現(xiàn)數(shù)據(jù)采集器在一段時(shí)間內(nèi)無法連接 GPS 的情況, 也可能表明臺(tái)陣使用的儀器內(nèi)部時(shí)鐘比較穩(wěn)定, 沒有明顯的漂移。

我們利用臺(tái)陣觀測期內(nèi)一個(gè)較大的極遠(yuǎn)震的 P波數(shù)據(jù)測試本文方法的準(zhǔn)確性。該地震的震中位于危地馬拉西岸(13.717°N, 90.972°W), 震源深度為28km, 震級為 Mw6.8, 距離臺(tái)陣中心的震中距為135°, 臺(tái)陣中記錄到清晰的 PKP 震相。PKP 震相是近似垂直地入射臺(tái)陣, 可以在最大程度上減少路徑差異對地震信號到時(shí)的影響。我們利用鹽源臺(tái)陣以及臺(tái)陣附近一個(gè)固定臺(tái)站(鹽源臺(tái))接收的該地震波形, 利用臨近臺(tái)站間的波形互相關(guān)函數(shù), 求出臺(tái)站間的到時(shí)差, 并設(shè)定固定臺(tái)的到時(shí)為 0, 求出臺(tái)陣中各個(gè)臺(tái)站的到時(shí)差。本文對該地震的 PKP 波到時(shí)數(shù)據(jù)進(jìn)行方位角校正后, 計(jì)算鹽源臺(tái)陣中各臺(tái)站與鹽源固定臺(tái)之間的 PKP 波到時(shí)差, 估算臺(tái)陣中各臺(tái)站的時(shí)鐘偏差, 結(jié)果如圖9 所示。

從圖 9 可以看到, 由地震信號求得的時(shí)鐘偏差與噪聲互相關(guān)方法求得的時(shí)鐘偏差的差別在±0.2s之間, 大部分臺(tái)站的偏差小于 0.1s, 說明本文提出的利用噪聲進(jìn)行單臺(tái)時(shí)鐘偏差分析的方法是有效的。兩種方法產(chǎn)生差異的可能原因主要有以下 3 個(gè)方面: 1)地震數(shù)據(jù)方法中, 由路徑結(jié)構(gòu)造成的干擾, 不同臺(tái)站下方速度結(jié)構(gòu)的差異會(huì)造成到時(shí)的小幅差異; 2)噪聲互相關(guān)方法中, 反演算法及噪聲源分布造成的干擾; 3)兩種方法代表的意義有所差別, 地震數(shù)據(jù)求得的時(shí)鐘偏差是瞬時(shí)偏差, 噪聲互相關(guān)方法求得的時(shí)間偏差是一天時(shí)鐘偏差的平均值。

我們注意到, 2017 年 6 月 12 日的單臺(tái)時(shí)鐘偏差分析結(jié)果中仍存在一定程度的系統(tǒng)性偏移(圖 6(d)), 臺(tái)陣東南角的臺(tái)站時(shí)鐘的整體偏差大于臺(tái)陣內(nèi)其他區(qū)域的臺(tái)站, 說明在這一天, 即使是方位角相近的臺(tái)站對, 時(shí)間偏移量仍有一定程度的差別, 僅使用不同方位角的平均值, 不能完全去除由噪聲源變化引起的互相關(guān)函數(shù)時(shí)間偏移。這種現(xiàn)象出現(xiàn)的一個(gè)可能原因是, 這一天的主要噪聲源來自離臺(tái)陣較近的地方, 甚至臺(tái)陣內(nèi)部, 因此其對臺(tái)陣中不同位置臺(tái)站對的影響不盡相同。我們查閱當(dāng)?shù)氐男〉卣鹩涗? 未發(fā)現(xiàn)臺(tái)陣東南角周邊在 6 月 12 日有任何地震記錄, 此處臺(tái)站的原始波形記錄中的確存在一些疑似事件記錄, 推測可能與當(dāng)?shù)氐牡V井采集等人為事件有關(guān)。同時(shí), 這一天結(jié)果的整體誤差也比其他時(shí)間偏大, 其原因主要有兩方面: 1)這一天的噪聲源變化較劇烈, 根據(jù)式(6)計(jì)算得出的數(shù)據(jù)誤差中噪聲源對應(yīng)的數(shù)據(jù)誤差 errsrc較大; 2)因受臺(tái)陣布局形態(tài)影響, 臺(tái)陣的東南角的臺(tái)站分布不均, 此處參與反演計(jì)算的臺(tái)站對數(shù)量不足, 臺(tái)站連線的方位角分布也與臺(tái)陣的其他位置有差別, 在進(jìn)行單臺(tái)時(shí)間偏差反演計(jì)算時(shí), 較少的臺(tái)站對會(huì)放大數(shù)據(jù)誤差對結(jié)果的影響, 我們需要后續(xù)更深入的研究來確認(rèn)這一天噪聲源的具體情況及其對臺(tái)陣噪聲互相關(guān)函數(shù)的影響方式。

介質(zhì)波速的變化也會(huì)引起噪聲互相關(guān)函數(shù)的時(shí)間延遲, 但這種影響極其微弱。例如, 地震引起的同震效應(yīng)會(huì)引起介質(zhì)波速的變化, 由噪聲測得的波速相對變化量一般小于 0.1%[16]。Pei 等[17]通過體波成像方法, 發(fā)現(xiàn)在較小的斷裂帶范圍內(nèi), 大地震后地震波的波速可比震前下降 4%。固體潮、降水和大氣壓的變化也會(huì)影響介質(zhì)的波速, 但波速的變化量在 0.1%左右, 遠(yuǎn)小于地震的同震變化[18]。并且, 與波速變化監(jiān)測常用的尾波信號相比, 直達(dá)波對介質(zhì)的波速變化敏感度不高。在鹽源臺(tái)陣觀測期間, 觀測區(qū)域沒有大震發(fā)生, 介質(zhì)的波速發(fā)生較強(qiáng)變化的可能性小。在我們選用的 0.1~0.5Hz 頻帶內(nèi), 研究區(qū)的面波速度約為 2.5km/s。以波速的相對變化為 1%來計(jì)算, 在 5km 臺(tái)站間距上能夠引起的走時(shí)變化約為 0.02s, 不足以造成本文檢測到的時(shí)鐘偏差。實(shí)際介質(zhì)的波速變化遠(yuǎn)達(dá)不到 1%, 因此在短時(shí)間內(nèi)檢測儀器鐘差時(shí), 介質(zhì)波速的微弱變化可以忽略。

4 結(jié)論

本文提出一種利用短周期地震觀測密集臺(tái)陣中相鄰臺(tái)站間的背景噪聲互相關(guān)函數(shù)進(jìn)行臺(tái)陣時(shí)鐘偏差分析的方法, 針對短周期密集臺(tái)陣的特性, 修正了參考互相關(guān)函數(shù)誤差及噪聲源變化對時(shí)鐘偏差分析結(jié)果的影響, 討論了數(shù)據(jù)拾取誤差、反演優(yōu)化方案、參考互相關(guān)函數(shù)修正以及噪聲源變化對結(jié)果的影響。通過在鹽源臺(tái)陣的應(yīng)用, 證明本文提出的方法可以較快速、準(zhǔn)確地給出密集臺(tái)陣全觀測期內(nèi)的連續(xù)時(shí)鐘偏差, 為密集臺(tái)陣數(shù)據(jù)的后續(xù)應(yīng)用提供可靠的依據(jù)。

對鹽源臺(tái)陣的時(shí)鐘偏差分析結(jié)果顯示, 209 個(gè)臺(tái)站中, 有 17 個(gè)臺(tái)站在觀測期內(nèi)出現(xiàn)較明顯的躍變式鐘差, 可能與儀器數(shù)據(jù)采集裝置的軟件或硬件故障有關(guān)。在利用臺(tái)陣數(shù)據(jù)進(jìn)行涉及時(shí)間的研究(如震相到時(shí)拾取、速度結(jié)構(gòu)反演等)時(shí), 應(yīng)特別注意去除時(shí)鐘偏差明顯的臺(tái)站數(shù)據(jù), 或利用時(shí)鐘偏差分析結(jié)果進(jìn)行相應(yīng)的修正。

[1] Sens-Sch?nfelder C. Synchronizing seismic networks with ambient noise. Geophys J Int, 2008, 174: 966–970

[2] Gouédard P, Seher T, McGuire J J, et al. Correction of ocean-bottom seismometer instrumental clock errors using ambient seismic noise. Bull Seismol Soc Am, 2014, 104(3): 1276–1288

[3] Stehly L, Campillo M, Shapiro N M. Traveltime measurements from noise correlations: stability and detection of instrumental time-shifts. Geophys J Int, 2007, 171: 223–230

[4] Xia Y J, Ni S D, Zeng X F, et al. Synchronizing intercontinental seismic networks using the 26 s per-sistent localized microseismic source. Bull Seismol Soc Am, 2015, 105(4): 2101–2108

[5] Anchieta M C, Wolfe C J, Pavlis G L, et al. Seis-micity around the Hawaiian islands recorded by the PLUME seismometer networks: insight into faulting near Maui, Molokai, and Oahu. Bull Seismol Soc Am, 2011, 101: 1742–1758

[6] 王俊, 鄭定昌, 詹小艷, 等. 基于背景噪聲互相關(guān)格林函數(shù)的單臺(tái)時(shí)間誤差估計(jì). 地震學(xué)報(bào), 2013, 35(6): 888–901

[7] Huang Y, Yao H J, Huang B, et al. Phase velocity variation at periods of 0.5–3 seconds in the Taipei basin of Taiwan from correlation of ambient seismic noise. Bull Seismol Soc Am, 2010, 100: 2250–2263

[8] Li C, Yao H J, Fang H J, et al. 3D near-surface shear-wave velocity structure from ambient-noise tomogra-phy and borehole data in the Hefei urban area, China. Seismological Research Letters, 2016, 87(4): 882–892

[9] Lin F C, Li D, Clayton R W, et al. High-resolution 3D shallow crustal structure in Long Beach, California: application of ambient noise tomography on a dense seismic array. Geophysics, 2013, 78: Q45–Q56

[10] Picozzi M, Parolai S, Bindi D, et al. Characterization of shallow geology by high-frequency seismic noise tomography. Geophys J Int, 2009, 176(1): 164–174

[11] Weaver R L, Lobkis O I. On the emergence of the green’s function in the correlations of a diffuse field. J Acoust Soc Am, 2001, 109: 435–439

[12] Shapiro N M, Campillo M. Emergence of broadband Rayleigh waves from correlations of the ambient noise. Geophys Res Lett, 2004, 31: L07614

[13] McNamara D E, Buland R P. Ambient noise levels in the continental united states. Bull Seismol Soc Am, 2004, 94(4): 1517–1527

[14] Paige C C, Saunders M A. LSQR: an algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software, 1982, 8(1): 43–71

[15] Bensen G D, Ritzowoller M H, Barmin M P, et al. Processing seismic ambient noise data to obtain re-liable broad-band surface wave dispersion measure-ments. Geophys J Int, 2007, 169: 1239–1260

[16] Brenguier F, Campillo M, Hadziioannou C, et al. Postseismic relaxation along the San Andreas fault at parkfield from continuous seismological observations. Science, 2008, 321: 1478–1481

[17] Pei S P, Niu F L, Ben-Zion Y, et al. Seismic velocity reduction and accelerated recovery due to earthquakes on the Longmenshan fault. Nature Geoscience, 2019, 12(4): 387–392

[18] Wang B S, Zhu P, Chen Y, et al. Continuous sub-surface velocity measurement with coda wave inter-ferometry. J Geophys Res, 2008, 113: B12313

A Synchronizing Method for Dense Seismic Array Based on Ambient Noise Correlation Function

TIAN Yuan1,2, WANG Weitao2,?, LI Li2, YU Changqing3, ZHANG Haiming4

1. China Academy of Electronics and Information Technology, Beijing 100041; 2. Institute of Geophysics, China Earthquake Administration, Beijing 100081; 3. Institute of Geology, Chinese Academy of Geological Science, Beijing 100037; 4. School of Earth and Space Sciences, Peking University, Beijing 100871; ? Corresponding author, E-mail: wangwt@cea-igp.ac.cn

This paper develops a NCF-based array synchronizing and error-estimation method employing the daily time shift of NCFs between nearby stations to estimate the stability of the timing system of a short-term but densely deployed seismic array. The error related to the reference NCF, the temporal change of ambient noise field, the data error and the optimization method in inversion are all estimated and correspondingly corrected based on the characteristics of a dense seismic array. The proposed method could effectively calculate the daily clock shift of a seismic array and pick out the time period of stations with obvious clock error. As an example, among the 209 stations of Yanyuan Seismic Array, it is found that 17 stations have obvious clock error with time shift over 1 s. These clock errors might be caused by the hardware or software problems related to the data acquisition system.

dense seismic array; ambient noise cross-correlation function; clock error analysis; noise source distribution

10.13209/j.0479-8023.2020.043

國家重點(diǎn)研發(fā)計(jì)劃(2018YFC1503200)和中國地震局地球物理研究所中央級公益性科研院所基本科研業(yè)務(wù)費(fèi)專項(xiàng)(DQJB18B26)資助

2019?06?24;

2019?09?23