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

基于背景噪聲互相關(guān)函數(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

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

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

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

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

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

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

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

1 方法原理

1.1 基本原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.6 單臺時鐘偏差計算及誤差分析

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

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

err=errdata+errref+errsrc, (7)

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

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

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

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

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

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

鹽源臺陣的單臺背景噪聲強度如圖3(a)所示?？梢钥闯? 臺陣的單日PSD在頻率大于1Hz時變化幅度較大, 可達40dB, 而頻率小于1Hz時整體強度較穩(wěn)定。因此, 從能量的角度考慮, 應(yīng)選擇1 Hz以下的低頻帶背景噪聲來進行互相關(guān)計算。圖3(b)給出在3個不同的頻帶范圍(1~2, 0.5~1和0.1~ 0.5Hz), 鹽源臺陣中684個臺間距5km以內(nè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é)果說明, 鹽源臺陣的相干噪聲場在0.5Hz以上的高頻帶非常不穩(wěn)定。為了在最大程度上保證通過背景噪聲互相關(guān)函數(shù)計算得出的時間偏差的準(zhǔn)確性, 我們選取0.1~0.5Hz作為背景噪聲互相關(guān)計算的頻帶范圍。

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

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

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

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

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

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

3 結(jié)果及討論

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

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

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

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

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

4 結(jié)論

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

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

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

2019?06?24;

2019?09?23