成萬(wàn)里,王秀英
(1.河南省地震局,河南 鄭州 450016;2.中國(guó)地震局地殼應(yīng)力研究所,北京 100085)
在地震科學(xué)領(lǐng)域,特別是“十五”數(shù)字地震觀測(cè)網(wǎng)絡(luò)項(xiàng)目以來(lái),地震前兆觀測(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í)間尺度特征,無(wú)需選擇基函數(shù)就可把復(fù)雜信號(hào)由精細(xì)尺度到粗大尺度分解為若干本征模態(tài)分量,適合對(duì)非平穩(wěn)、非線性信號(hào)進(jìn)行平穩(wěn)化處理。對(duì)于由偶然因素引起的脈沖干擾噪聲很容易通過(guò)限幅濾波法、中值濾波法和算術(shù)平均法等方法去除??梢钥闯觯鲜鼋翟敕椒ǘ际菑脑夹盘?hào)中去除這些加噪聲信息的盲源分離方法。通常情況下,特定周期的環(huán)境動(dòng)態(tài)變化與噪聲信息的周期或尺度無(wú)顯著差異時(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)變化無(wú)明顯影響。據(jù)此特點(diǎn),文章以觀測(cè)數(shù)據(jù)短期變化能否有效區(qū)分長(zhǎng)時(shí)間環(huán)境變化,來(lái)描述噪聲信號(hào)與環(huán)境動(dòng)態(tài)變化速率的差異性,并以此來(lá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)用問(wèn)題的新方法;同時(shí),用全部實(shí)測(cè)數(shù)據(jù)參與計(jì)算,避免偶然因素對(duì)觀測(cè)數(shù)據(jù)和計(jì)算結(jié)果的影響。
儀器觀測(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,說(shuō)明儀器可以滿足觀測(cè)環(huán)境變化的要求。如果較長(zhǎng)時(shí)間的環(huán)境變化與較短時(shí)間的儀器性能變化接近,則儀器無(wú)法滿足觀測(cè)環(huán)境的要求。
按照上述方法原理,將短時(shí)間內(nèi)觀測(cè)數(shù)據(jù)的波動(dòng)看作主要是由儀器自身分辨力造成。如果在此較短時(shí)間觀測(cè)環(huán)境的波動(dòng)與儀器分辨力接近,得到的儀器分辨力仍可反映儀器自身的特性。
利用該方法計(jì)算分辨力時(shí)需解決兩個(gè)問(wèn)題。第一是確定觀測(cè)數(shù)據(jù)滿足IT≈ID的時(shí)間窗長(zhǎng)度;第二是如何避免計(jì)算時(shí)偶然因素的影響。針對(duì)問(wèn)題一,采用不同的時(shí)間窗計(jì)算數(shù)據(jù)波動(dòng)情況,綜合大量數(shù)據(jù)選擇一個(gè)合理的計(jì)算時(shí)間窗;針對(duì)問(wèn)題二,全部數(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è)問(wèn)題,最終轉(zhuǎn)化為大量數(shù)據(jù)計(jì)算問(wèn)題,通過(guò)全部數(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]。
深井溫度短期變化非常微小,波動(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ù)庫(kù)中2006年至2017年所有水溫觀測(cè)數(shù)據(jù),共計(jì)595套儀器約128萬(wàn)條日數(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ì)算。
據(jù)式(1),確定一個(gè)合理的時(shí)間窗,使得在此計(jì)算時(shí)段內(nèi),環(huán)境變化盡量足夠小。原理上這個(gè)時(shí)間越短越好,但時(shí)間太短,分鐘采樣數(shù)據(jù)參與運(yùn)算的太少,會(huì)由于偶然數(shù)據(jù)問(wèn)題對(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)境較好,日變幅、月變幅、年變幅分別不超過(guò)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)來(lái)描述儀器自身影響的變化,用滑動(dòng)平均法(n=3)平滑后X的最大變幅RANG(X)表示水溫受環(huán)境影響的變化。在數(shù)據(jù)上下均勻波動(dòng)無(wú)趨勢(shì)變化的理想狀態(tài)下, RANG(X)/2略小于SD(X)。根據(jù)這一數(shù)學(xué)特性,當(dāng)RANG(X)/2明顯大于SD(X)時(shí),表明此時(shí)環(huán)境影響與儀器自身影響的信息可以區(qū)分開(kāi)。通過(guò)分別計(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ū)分開(kāi)。由于環(huán)境變化與儀器自身影響變化速率不同,20 min后,RANG(X)/2明顯大于SD(X)。說(shuō)明沂南臺(tái)水溫?cái)?shù)據(jù)Δt 取值20 min附近能夠較好地區(qū)分出儀器自噪聲和環(huán)境變化。需要說(shuō)明的是,在計(jì)算水溫動(dòng)態(tài)變化特征時(shí),時(shí)間窗選取采用5年以上的數(shù)據(jù)全部參與計(jì)算,通過(guò)聚類方法排除偶然因素后的計(jì)算結(jié)果。如計(jì)算長(zhǎng)度為15 min的SD(X)與RANG(X)/2,將全部數(shù)據(jù)中滿足15分鐘長(zhǎng)度的數(shù)據(jù)段全部計(jì)算出結(jié)果,通過(guò)聚類方法,將密度最大、個(gè)數(shù)最多的一類視為正常,將該類結(jié)果的均值作為最后結(jié)果。圖1表明,環(huán)境變化和儀器自身影響變化都是一種穩(wěn)定規(guī)律的變化,說(shuō)明計(jì)算結(jié)果可信。
按照如上方法原理,時(shí)間窗確定后,還需解決計(jì)算結(jié)果如何避免偶然因素影響的問(wèn)題??臻g維度上,有些測(cè)點(diǎn)的數(shù)據(jù)變化不正常;時(shí)間維度上,一個(gè)測(cè)點(diǎn)有些時(shí)段數(shù)據(jù)變化不正常。研究通過(guò)兩步聚類分析來(lái)分別篩選掉這些異常數(shù)據(jù)。聚類分析是用機(jī)器學(xué)習(xí)的方法將數(shù)據(jù)集中在某些方面相似的數(shù)據(jù)成員進(jìn)行自動(dòng)分類組織的過(guò)程。因此無(wú)論在空間維度還是時(shí)間維度,都可以將絕大部分正常時(shí)段的數(shù)據(jù)聚為一類,將少量異常時(shí)段的數(shù)據(jù)聚為一類。具體過(guò)程與結(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的均值差來(lái)描述環(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,說(shuō)明該類測(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
在實(shí)驗(yàn)方法設(shè)計(jì)時(shí),時(shí)間窗的選取是重要步驟,其對(duì)實(shí)驗(yàn)結(jié)果有影響,是一個(gè)求解分離觀測(cè)數(shù)據(jù)的宏觀變化與微觀變化最優(yōu)解的過(guò)程。研究采取大量實(shí)驗(yàn)對(duì)比方法選取的最優(yōu)值,有待于通過(guò)設(shè)計(jì)算法完成的最優(yōu)解來(lái)驗(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),無(wú)論分辨力高低,因儀器自噪聲水平已經(jīng)大于長(zhǎng)期環(huán)境變化,該類測(cè)點(diǎn)至少無(wú)法滿足記錄環(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類儀器差異明顯地展示出來(lái),但對(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è)非常有意義的研究方向。
從某一類儀器的全部觀測(cè)數(shù)據(jù)中挖掘信息,用來(lái)分析估計(jì)儀器分辨力,是一種新方法的應(yīng)用嘗試。該方法對(duì)研究數(shù)據(jù)質(zhì)量要求低,抗干擾性強(qiáng),避免缺數(shù)插值、臺(tái)階處理等復(fù)雜的數(shù)據(jù)預(yù)處理過(guò)程,數(shù)據(jù)的統(tǒng)計(jì)模型簡(jiǎn)單,物理意義明確。相對(duì)于實(shí)驗(yàn)室對(duì)比觀測(cè)的分辨力研究方法,更容易操作。實(shí)驗(yàn)結(jié)果表明,得出各類儀器的相對(duì)分辨力與其實(shí)際性能都較吻合,說(shuō)明該研究方法有效、可行。