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

成萬里，王秀英

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

0 引言

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

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

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

1.1 方法原理分析

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

IT=ID+IE

。

(1)

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

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

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

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

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

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

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

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

2 計(jì)算時(shí)間窗確定

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

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

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

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

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

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

3 計(jì)算結(jié)果

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

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

表2 臺(tái)站代碼為16002測(cè)點(diǎn)1的均值差聚類中心Table 2 The mean difference clustering center of the measurement points 1 in the station with code 16002

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

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

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

圖3 三套儀器在不同測(cè)項(xiàng)的分辨力Fig.3 Resolution of three sets of instruments in different terms

4 討論分析

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

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

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

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

5 結(jié)論

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