基于數(shù)據(jù)的前兆儀器分辨力檢測方法研究

成萬里，王秀英

(1.河南省地震局，河南 鄭州 450016；2.中國地震局地殼應(yīng)力研究所，北京 100085)

0 引言

在地震科學(xué)領(lǐng)域，特別是“十五”數(shù)字地震觀測網(wǎng)絡(luò)項目以來，地震前兆觀測數(shù)據(jù)在時間積累、采樣頻率、觀測布點、學(xué)科種類四個維度上激增，各學(xué)科產(chǎn)出的海量數(shù)據(jù)中蘊含著大量價值信息,同時也包含大量的噪聲信息，要從噪聲數(shù)據(jù)中分辨出價值信息，必須正確認(rèn)識儀器分辨力。劉洋君等認(rèn)為地震儀器的信噪比大小決定測量信號的真實度,因而各種噪聲之和的大小決定測量的精度、靈敏度及檢測下限。要降低檢測下限,提高檢測精度,首先應(yīng)該設(shè)法降低各種噪聲的水平[1]。目前關(guān)于信號降噪的研究很多，包括小波變換、經(jīng)驗?zāi)B(tài)分解、濾波技術(shù)等[2-3]，利用小波方法,將信號分解為不同頻段，實現(xiàn)信噪分離。但是,如果對信號的特征頻段不夠了解，則難以獲得滿意的結(jié)果。經(jīng)驗?zāi)B(tài)分解是基于信號本身的時間尺度特征，無需選擇基函數(shù)就可把復(fù)雜信號由精細(xì)尺度到粗大尺度分解為若干本征模態(tài)分量，適合對非平穩(wěn)、非線性信號進(jìn)行平穩(wěn)化處理。對于由偶然因素引起的脈沖干擾噪聲很容易通過限幅濾波法、中值濾波法和算術(shù)平均法等方法去除?？梢钥闯觯鲜鼋翟敕椒ǘ际菑脑夹盘栔腥コ@些加噪聲信息的盲源分離方法。通常情況下，特定周期的環(huán)境動態(tài)變化與噪聲信息的周期或尺度無顯著差異時，盲源分離方法很難真正實現(xiàn)信噪分離[4]。因此，噪聲對數(shù)據(jù)應(yīng)用的影響大小與環(huán)境動態(tài)變化特征有關(guān)，如分鐘值尺度的噪聲信號可能與小時周期動態(tài)變化尺度相當(dāng)，兩者難以分辨，但其可能對日周期環(huán)境動態(tài)變化無明顯影響。據(jù)此特點，文章以觀測數(shù)據(jù)短期變化能否有效區(qū)分長時間環(huán)境變化，來描述噪聲信號與環(huán)境動態(tài)變化速率的差異性，并以此來度量儀器的相對分辨力。這與通常在實驗室中測試完成的儀器分辨力有一定的差別。事實上，實驗室的測試環(huán)境與實際觀測環(huán)境就有很大差別。

收集2006年至 2017年全國595個測點近1TB的地下水溫數(shù)據(jù)，以觀測數(shù)據(jù)短期變化能否有效區(qū)分長時間環(huán)境變化為度量依據(jù)，嘗試用大數(shù)據(jù)統(tǒng)計分析方法從這巨量數(shù)據(jù)中估算3類地?zé)嵊^測儀器的分辨力。一方面探索從數(shù)據(jù)中解決實際監(jiān)測應(yīng)用問題的新方法；同時，用全部實測數(shù)據(jù)參與計算，避免偶然因素對觀測數(shù)據(jù)和計算結(jié)果的影響。

1 方法與數(shù)據(jù)

1.1 方法原理分析

儀器觀測數(shù)據(jù)攜帶的信息IT，可分解為儀器自身影響的信息ID和環(huán)境動態(tài)變化影響信息IE。其中，環(huán)境影響包含地震、地球物理場及各種干擾等，即：

IT=ID+IE

。

(1)

大多數(shù)情況，環(huán)境信息的變化需要一定的持續(xù)時間，認(rèn)為在較短時間內(nèi)觀測環(huán)境IE不變或變化極小可以忽略不計，則較短時間的變化主要由儀器自身因素引起，式(1)表示為IT≈ID，表明較短時間內(nèi)數(shù)據(jù)的變化僅反映儀器自身性能的變化。

較長時間的變化反映出各種因素疊加的環(huán)境變化，如果環(huán)境變化遠(yuǎn)大于儀器性能(分辨力)變化，則式(1)表示為IT≈IE，說明儀器可以滿足觀測環(huán)境變化的要求。如果較長時間的環(huán)境變化與較短時間的儀器性能變化接近，則儀器無法滿足觀測環(huán)境的要求。

按照上述方法原理，將短時間內(nèi)觀測數(shù)據(jù)的波動看作主要是由儀器自身分辨力造成。如果在此較短時間觀測環(huán)境的波動與儀器分辨力接近，得到的儀器分辨力仍可反映儀器自身的特性。

利用該方法計算分辨力時需解決兩個問題。第一是確定觀測數(shù)據(jù)滿足IT≈ID的時間窗長度；第二是如何避免計算時偶然因素的影響。針對問題一，采用不同的時間窗計算數(shù)據(jù)波動情況，綜合大量數(shù)據(jù)選擇一個合理的計算時間窗；針對問題二，全部數(shù)據(jù)參與計算，抵消短時偶然因素對數(shù)據(jù)的影響。如，利用同一測點全部數(shù)據(jù)參與計算，降低局部可能存在的干擾或其他影響；利用同類觀測儀器、不同觀測點的全部數(shù)據(jù)，避免單一臺站或測點可能出現(xiàn)的偏差或偶發(fā)情況。因此，對于兩個問題，最終轉(zhuǎn)化為大量數(shù)據(jù)計算問題，通過全部數(shù)據(jù)參與計算，挖掘出數(shù)據(jù)的動態(tài)變化特征，對儀器技術(shù)指標(biāo)進(jìn)行合理評價。

計算儀器分辨力、環(huán)境變化等反映數(shù)據(jù)動態(tài)變化的指標(biāo)，稱為特征量，可以選取基于統(tǒng)計學(xué)的均值、標(biāo)準(zhǔn)差、中位數(shù)、極差等量，也可以根據(jù)具體的研究目標(biāo)特別設(shè)計，視具體任務(wù)而定[5]。

1.2 研究數(shù)據(jù)

深井溫度短期變化非常微小，波動通常不大于0.000 1 ℃，年漂移變化幅度遠(yuǎn)小于0.1 ℃～ 0.2 ℃[6-7]。地下水溫的這種變化特點，使式(1)中短期瞬時變化與長期環(huán)境變化的關(guān)系更清楚。選取地下水溫數(shù)據(jù)作為實驗數(shù)據(jù)，取“十五”前兆數(shù)據(jù)庫中2006年至2017年所有水溫觀測數(shù)據(jù)，共計595套儀器約128萬條日數(shù)據(jù)記錄，每條日數(shù)據(jù)記錄長度為1440的觀測時間序列。主要包含Szw-1a、Szw-1和ZKGD3000三種型號的觀測儀器，按表層(4311)、中層(4312)、深層(4313)三段水溫劃分測項代碼。因部分測點停測或數(shù)據(jù)不可用，整理出有效參與計算的數(shù)據(jù)信息如表1所示。

由于選取觀測時間跨度較長，儀器觀測起止時間不統(tǒng)一，缺測情況較多，因此研究采取小時段分段的方法選取數(shù)據(jù)，不考慮數(shù)據(jù)的完整率和連續(xù)率，按照大數(shù)據(jù)分析的思想，只要有數(shù)據(jù)片段均可參與計算。

2 計算時間窗確定

據(jù)式(1)，確定一個合理的時間窗，使得在此計算時段內(nèi)，環(huán)境變化盡量足夠小。原理上這個時間越短越好，但時間太短，分鐘采樣數(shù)據(jù)參與運算的太少，會由于偶然數(shù)據(jù)問題對計算結(jié)果產(chǎn)生較大影響。

表1 全國地下水溫觀測數(shù)據(jù)基本信息Table 1 Basic information of underground water temperature observation data in China

因此，選擇的時間窗雖盡量短，但確保有足夠的數(shù)據(jù)參與計算。

為確定實際儀器分辨力評估計算時的時間窗，以沂南臺的水溫觀測數(shù)據(jù)進(jìn)行計算選取。沂南臺的觀測環(huán)境較好，日變幅、月變幅、年變幅分別不超過0.004 ℃、0.01 ℃和0.2 ℃，對于水溫觀測環(huán)境的要求具有代表性。

確定時間窗需要描述數(shù)據(jù)動態(tài)變化的指標(biāo)。用一小段數(shù)據(jù)X的標(biāo)準(zhǔn)差SD(X)來描述儀器自身影響的變化，用滑動平均法(n=3)平滑后X的最大變幅RANG(X)表示水溫受環(huán)境影響的變化。在數(shù)據(jù)上下均勻波動無趨勢變化的理想狀態(tài)下， RANG(X)/2略小于SD(X)。根據(jù)這一數(shù)學(xué)特性，當(dāng)RANG(X)/2明顯大于SD(X)時，表明此時環(huán)境影響與儀器自身影響的信息可以區(qū)分開。通過分別計算不同時間窗(即不同長度的X)下的SD(X)與RANG(X)/2，找到剛好RANG(X)/2大于SD(X)的時間窗。如圖1所示，設(shè)計沂南臺水溫數(shù)據(jù)從10 min～120 min12個Δt的時間窗進(jìn)行對比。

圖1 沂南臺SD(X)與RANG(X)/2的變化關(guān)系Fig.1 The relationship between SD (X) and RANG (X)/2 in Yinan Station

由圖1看出，20 min之前，RANG(X)/2足夠小，很難與SD(X)區(qū)分開。由于環(huán)境變化與儀器自身影響變化速率不同，20 min后，RANG(X)/2明顯大于SD(X)。說明沂南臺水溫數(shù)據(jù)Δt 取值20 min附近能夠較好地區(qū)分出儀器自噪聲和環(huán)境變化。需要說明的是，在計算水溫動態(tài)變化特征時，時間窗選取采用5年以上的數(shù)據(jù)全部參與計算，通過聚類方法排除偶然因素后的計算結(jié)果。如計算長度為15 min的SD(X)與RANG(X)/2，將全部數(shù)據(jù)中滿足15分鐘長度的數(shù)據(jù)段全部計算出結(jié)果，通過聚類方法，將密度最大、個數(shù)最多的一類視為正常，將該類結(jié)果的均值作為最后結(jié)果。圖1表明，環(huán)境變化和儀器自身影響變化都是一種穩(wěn)定規(guī)律的變化，說明計算結(jié)果可信。

3 計算結(jié)果

按照如上方法原理，時間窗確定后，還需解決計算結(jié)果如何避免偶然因素影響的問題?？臻g維度上，有些測點的數(shù)據(jù)變化不正常；時間維度上，一個測點有些時段數(shù)據(jù)變化不正常。研究通過兩步聚類分析來分別篩選掉這些異常數(shù)據(jù)。聚類分析是用機器學(xué)習(xí)的方法將數(shù)據(jù)集中在某些方面相似的數(shù)據(jù)成員進(jìn)行自動分類組織的過程。因此無論在空間維度還是時間維度，都可以將絕大部分正常時段的數(shù)據(jù)聚為一類，將少量異常時段的數(shù)據(jù)聚為一類。具體過程與結(jié)果如下。

以Szw-1a型儀器、測項代碼為4313的中層水溫209個測點573 250條記錄為例(每條記錄包含24 h的數(shù)據(jù))。分別計算3個相隔長度為6 h的均值差來描述環(huán)境變化情況。在每個測點內(nèi)按照均值差進(jìn)行均值聚類，記錄個數(shù)最多的一類作為代表該測點正常環(huán)境變化的記錄，并求均值作為該測點的環(huán)境變化特征值。表2是其中臺站代碼為16002測點1的3 503條記錄的聚類中心。顯然，聚類結(jié)果最多的第1類，其聚類中心的3個值差異最小，樣本個數(shù)最多，數(shù)據(jù)比較符合邏輯，視為正常記錄，將該類的聚類中心3個值作為該測點的環(huán)境變化特征。

表2 臺站代碼為16002測點1的均值差聚類中心Table 2 The mean difference clustering center of the measurement points 1 in the station with code 16002

分別篩選出209個測點的環(huán)境變化正常時段樣本，計算各測點環(huán)境變化特征。從空間維度根據(jù)測點環(huán)境變化特征對209個測點進(jìn)行二次聚類，其聚類結(jié)果如圖2所示(圖中Y軸儀器相對分辨力反映的是209個測點儀器自身因素影響的變化整體情況，X軸反映的是209個測點6小時內(nèi)受環(huán)境影響變化的整體情況)。二次聚類中1類的113個測點均值差聚類中心為0.000 55，說明該類測點受觀測環(huán)境影響較小，認(rèn)為這類測點短時間內(nèi)環(huán)境溫度相對穩(wěn)定。將描述其短時間變化的標(biāo)準(zhǔn)差作為該類測點的相對分辨力。經(jīng)計算該類113個測點時間窗內(nèi)的標(biāo)準(zhǔn)差均值為0.000 38。

圖2 szw-1a型水溫儀觀測數(shù)據(jù)環(huán)境變化特征二次聚類結(jié)果Fig.2 Quadratic clustering results of the environmental change characteristics of szw-1a type water temperature meter observation data

按照該方法，計算其他類型儀器和測項，結(jié)果如圖3所示。不同型號的儀器及不同的入水深度在相對分辨力上均有一定的差異性。圖中結(jié)果表明該計算結(jié)果與儀器實際性能較吻合。

圖3 三套儀器在不同測項的分辨力Fig.3 Resolution of three sets of instruments in different terms

4 討論分析

在實驗方法設(shè)計時，時間窗的選取是重要步驟，其對實驗結(jié)果有影響，是一個求解分離觀測數(shù)據(jù)的宏觀變化與微觀變化最優(yōu)解的過程。研究采取大量實驗對比方法選取的最優(yōu)值，有待于通過設(shè)計算法完成的最優(yōu)解來驗證其精確度。

分析表明，受觀測環(huán)境影響較小的測點和觀測時段的樣本更適合用于計算儀器的分辨力，這類測點對儀器分辨力要求更高。對于相當(dāng)一部分環(huán)境日變化較小的測點，儀器的分辨力不能分辨出日變化以內(nèi)的環(huán)境變化。從圖2看出，測點的分辨力和等距環(huán)境變化沿著聚類中心呈一定斜率展布，其中部分至原點斜率大于1的測點，無論分辨力高低，因儀器自噪聲水平已經(jīng)大于長期環(huán)境變化，該類測點至少無法滿足記錄環(huán)境日變化的要求。

根據(jù)實驗結(jié)果，發(fā)現(xiàn)測點的數(shù)量對實驗結(jié)果有影響。淺層水溫ZKGD3000型測點有5個，Szw-1型16個；深層水溫ZKGD3000型測點有11個，Szw-1型9個。因這幾類測點數(shù)量太少，在自動聚類時穩(wěn)定性不是很好，其結(jié)果與實際情況有一定的偏差。如圖3所示，盡管結(jié)果能將3類儀器差異明顯地展示出來，但對傳感器入水深度差異的區(qū)分效果不明顯。

考慮到儀器類型、傳感器入水深度、季節(jié)變化、地域影響等因素，文章僅對不同儀器和不同測項做差異性假設(shè)檢驗。實驗表明，同一類儀器在不同測項下計算的分辨力有差異，同一測項的不同類儀器的測點分辨力也有差異。因此在設(shè)計實驗時，按照測項進(jìn)行分類。從計算結(jié)果看，每類儀器在不同測項下的分辨力有小幅度的差異，但不影響整體評價。季節(jié)變化、地域影響等受測點數(shù)量的限制，未做一一分類計算，但這些影響因素的差異性分析也是一個非常有意義的研究方向。

5 結(jié)論

從某一類儀器的全部觀測數(shù)據(jù)中挖掘信息，用來分析估計儀器分辨力，是一種新方法的應(yīng)用嘗試。該方法對研究數(shù)據(jù)質(zhì)量要求低，抗干擾性強，避免缺數(shù)插值、臺階處理等復(fù)雜的數(shù)據(jù)預(yù)處理過程，數(shù)據(jù)的統(tǒng)計模型簡單，物理意義明確。相對于實驗室對比觀測的分辨力研究方法，更容易操作。實驗結(jié)果表明，得出各類儀器的相對分辨力與其實際性能都較吻合，說明該研究方法有效、可行。