陳家樑,張寶劍
(福建省地震局,福建福州,350003)
福建地震臺網(wǎng)中心從1971年開始建設(shè),從最初3個測震臺站的模擬觀測及人工震中定位,發(fā)展到2020年底88個測震臺的數(shù)字觀測臺網(wǎng),得益于國家的發(fā)展壯大。福建地震臺網(wǎng)包含了測震臺、強(qiáng)震臺、烈度臺及GNSS臺網(wǎng),臺站密度走在了全國前列。其中測震臺網(wǎng)中部署有寬頻帶地震計的臺站占據(jù)84%,為組建虛擬臺網(wǎng)我中心接入了鄰省臺站及16個臺灣臺,福建本省測震臺站的平均臺間距約37千米,震級的觀測范圍為0-6級。臺站數(shù)量的增加提高了地震定位的精度,而地震數(shù)據(jù)的可靠性決定了地位精度的準(zhǔn)確性,因此有必要對地震數(shù)據(jù)質(zhì)量進(jìn)行多方面的監(jiān)測。方位角應(yīng)用于包括SKS剪切波分裂,接收函數(shù)研究等多項科學(xué)研究,如果方位角存在偏差將會影響相關(guān)研究結(jié)果的可靠性。以往都是在相關(guān)研究出現(xiàn)異常時才會發(fā)現(xiàn)臺站有故障,這樣的故障維護(hù)就會出現(xiàn)延遲,且會造成該時間段內(nèi)的數(shù)據(jù)不可用。因此需要有專門軟件對臺站方位角進(jìn)行實時監(jiān)測,發(fā)現(xiàn)問題及時處理,這樣不僅可以提高故障修復(fù)的時效性,且可以讓儀器維護(hù)人員的工作更具有針對性。
圖1 福建臺網(wǎng)及臺灣交換臺站分布
方位角的檢測需要用到震中距30度以上的遠(yuǎn)震,所以該程序為準(zhǔn)實時系統(tǒng)。根據(jù)具體的操作步驟,程序分為地震事件選取模塊、數(shù)據(jù)預(yù)處理模塊、臺站屬性讀取模塊、最小能量計算模塊、相關(guān)系數(shù)計算模塊及虛擬臺站構(gòu)建模塊。為簡化代碼數(shù)量,減少程序運行中所占用內(nèi)存,本程序使用Python中的Obspy這個開源的地震數(shù)據(jù)處理框架。
本程序中計算方位角偏差的方法為利用P波質(zhì)點的運動特性法,需要震中距在30度以上的遠(yuǎn)震,因此在準(zhǔn)備地震事件時需要計算震中距是否達(dá)到30度。福建本省面積較小,對于遠(yuǎn)震而言省內(nèi)不同臺站對于同一地震事件的震中不會有太大差別,因此我們選取位于福建較為中心的永安燕西作為所有臺站的參照臺,只要計算該臺與地震事件的震中距即可,這樣可以減少運算量又不影響數(shù)據(jù)的可靠性。計算震中距我們引入Obspy.geodetics中的gps2dist_azimuth及kilometers2degrees兩個方法,前者用于計算震中距,其返回值是一個數(shù)組,震中距為第一個元素且單位為千米,所以需要使用kilometers2degrees將單位從千米轉(zhuǎn)化為度。程序?qū)逊险鹬芯嘁蟮牡卣鹗录丶罢鹬芯噍敵鲋林付窂较碌腡XT文件中。
圖2 程序架構(gòu)
因?qū)掝l帶地震計中有60秒與120秒甚至更高的頻帶,為保證帶寬一致所以需要進(jìn)行帶通濾波。臺站的采樣率為100HZ,為減少計算量程序采取降采樣處理。利用Obspy.Stream類循環(huán)讀取地震離線波形數(shù)據(jù),全部讀取后利用該類中的resample(sampling_rate)方法即可進(jìn)行降采樣處理。Obspy.Stream.filter可以對數(shù)據(jù)進(jìn)行帶通濾波處理,用該方法離線波形數(shù)據(jù)進(jìn)行帶通濾波處理的具體參數(shù)設(shè)置如下所示 filter(‘bandpass’,freqmin=0.02,freqmax=0.2),參數(shù)表示使用帶通濾波處理數(shù)據(jù)且設(shè)置了低頻與高頻的具體數(shù)值。
臺站屬性指的是臺站名、經(jīng)緯度及高程,這些數(shù)據(jù)均存儲于MySQL數(shù)據(jù)庫中。對于省級地震臺網(wǎng)中心,有固定的數(shù)據(jù)庫用于配置各臺站經(jīng)緯度、高程及靈敏度等參數(shù),所以我們用xml文件配置好數(shù)據(jù)庫的地址等,程序就可以通過xml文件進(jìn)入數(shù)據(jù)庫進(jìn)行臺站屬性的讀取。Xml文件的內(nèi)容形式如下所示:
該模塊中我們創(chuàng)建一個專門讀取xml文件的類,在Python中引用xml.dom.minidom類包,就可以按照該類讀取xml文件的步驟方法,將xml文件中的各配置內(nèi)容讀取到指定的數(shù)組中并返回該數(shù)組。
獲取上述數(shù)據(jù)庫的配置信息后連接MySQL數(shù)據(jù)庫,寬頻帶地震計的通道命名有顯著要求,即開頭字母為B,其三通道的命名為BHN、BHE、BHZ。因此在選取寬頻帶臺站的相關(guān)屬性時對數(shù)據(jù)庫的操作語句可寫作“select Sta_code from Channel_info where Chn_code like ‘%BHE%’”,其 中Sta_code為臺站名,Channel_info為表名稱,Chn_code為通道存儲字段名。通過以上操作可以獲取全部寬頻帶地震計臺站名,而后運用循環(huán)方法獲取所有的臺站屬性后封裝到dictionary中返回。
通過以上以上步驟計算每個臺站對應(yīng)每個地震事件的P波到時,選取到時前后五秒內(nèi)經(jīng)過旋轉(zhuǎn)的能量,再以1度為間隔進(jìn)行180度的旋轉(zhuǎn)計算每一度對應(yīng)的能量,再進(jìn)行對比選出最小能量對應(yīng)的旋轉(zhuǎn)度數(shù),認(rèn)為此時對應(yīng)的角度為與正北方向的角度偏差。再將此時的偏差角與經(jīng)過經(jīng)緯度經(jīng)過計算出的理論偏差值進(jìn)行對比,即可得出此臺對應(yīng)該地震事件時間點上的角度偏差。將每個地震事件的偏差值輸出至String類型的數(shù)組,計算完所有事件的偏差值后,再利用“SNR”法計算所有事件中每一度的能量值后選出最小值對應(yīng)的偏差角,即認(rèn)為是最終的地震計與正北方向的偏差值。
圖3 時間序列圖
為更清晰的顯示每個事件對應(yīng)的方位角偏差,在此步驟中將所得到的數(shù)據(jù)進(jìn)行繪圖顯示。因利用了GMT繪圖軟件所以需要在主機(jī)上安裝GMT5.0以上版本,否則會出現(xiàn)繪圖命令錯誤或不兼容的情況。時間序列圖的縱軸單位為的形式,所以需要用到GMT中的轉(zhuǎn)義字符“@~\152\244\050\260\051@~”。將GMT繪圖命令以String形式寫入bat文件后,執(zhí)行該文件即可進(jìn)行繪圖操作,繪制后的圖片將按xml文件中配置的路徑存儲。
為更加精準(zhǔn)的分析異常方位角偏差臺站,該程序還設(shè)計了異常臺站與正常臺站所構(gòu)建的虛擬臺陣對應(yīng)通道的相關(guān)系數(shù)分析模塊。此模塊以故障臺站為中心,計算所有臺站與此臺的臺間距,并選取臺間距最近的15個正常臺站構(gòu)成虛擬臺陣。線性疊加臺陣內(nèi)所有臺站相同通道的數(shù)據(jù)后,與故障臺相同通道進(jìn)行相關(guān)系數(shù)計算,如果故障臺不出現(xiàn)儀器方向擺放錯誤或極性反轉(zhuǎn)等重大問題,則其相關(guān)系數(shù)會接近1相反則會接近0或負(fù)數(shù)。其相關(guān)系數(shù)的計算引入scipy.stats類中的pearsonr方法,可以迅速計算出兩個相同通道的相關(guān)系數(shù)。為了便于觀察,此模塊同樣具備輸出波形形狀及相關(guān)系數(shù)值的功能。
圖4 上圖為故障臺,下圖為正常臺與虛擬臺陣的相關(guān)系數(shù)波形圖
通過該程序?qū)Ω=ǖ卣鹋_網(wǎng)所屬寬頻帶地震計進(jìn)行計算,地震事件為2015年1月到2019年8月之間的遠(yuǎn)震,得到74個寬頻帶地震計臺站中有72個臺的方位角偏差值在-7°~7°之間,這樣的偏差值在該方法的理論中可以忽略不計,因此福建臺網(wǎng)寬頻帶地震計臺站的方位角總體合格率在97%如圖5所示。其中兩個存在較大偏差值的臺站均為BBVS-60地震計,此類型地震計沒有明顯指北標(biāo)志,所以導(dǎo)致維護(hù)人員在更換過程中無法準(zhǔn)確的進(jìn)行指北操作,經(jīng)過實地校正后已使方位角的偏差值減小到要求值1°以內(nèi)。
圖5 方位角偏差臺站分布
通過對程序模塊化開發(fā),可以理順每個模塊的具體功能,使程序整體研發(fā)更清晰,更易調(diào)試出現(xiàn)問題的部分代碼進(jìn)行小范圍的修復(fù),而不至于出現(xiàn)牽一發(fā)動全身式的更改。因程序為處理地震數(shù)據(jù)而開發(fā),所以使用到Obspy庫,因其幾乎支持地震學(xué)界內(nèi)通常使用的所有波形格式的讀寫,可簡化大量的數(shù)據(jù)格式之間轉(zhuǎn)換的代碼,同時他還集成了大量的地震學(xué)及數(shù)學(xué)所用的專有庫,對于開發(fā)者只需要對其進(jìn)行引入就可使用。程序應(yīng)用在福建臺網(wǎng)中心后,經(jīng)過實地校核發(fā)現(xiàn)該程序的計算結(jié)果可靠。該程序的應(yīng)用彌補(bǔ)了福建臺網(wǎng)中心在方位角檢測方面的不足,且提高了方位角偏差的檢測效率及修復(fù)故障的時效性。