方位角檢測程序的設計及應用

陳家樑，張寶劍

(福建省地震局,福建福州，350003)

0 引言

福建地震臺網中心從1971年開始建設，從最初3個測震臺站的模擬觀測及人工震中定位，發展到2020年底88個測震臺的數字觀測臺網，得益于國家的發展壯大。福建地震臺網包含了測震臺、強震臺、烈度臺及GNSS臺網，臺站密度走在了全國前列。其中測震臺網中部署有寬頻帶地震計的臺站占據84%，為組建虛擬臺網我中心接入了鄰省臺站及16個臺灣臺，福建本省測震臺站的平均臺間距約37千米，震級的觀測范圍為0-6級。臺站數量的增加提高了地震定位的精度，而地震數據的可靠性決定了地位精度的準確性，因此有必要對地震數據質量進行多方面的監測。方位角應用于包括SKS剪切波分裂，接收函數研究等多項科學研究，如果方位角存在偏差將會影響相關研究結果的可靠性。以往都是在相關研究出現異常時才會發現臺站有故障，這樣的故障維護就會出現延遲，且會造成該時間段內的數據不可用。因此需要有專門軟件對臺站方位角進行實時監測，發現問題及時處理，這樣不僅可以提高故障修復的時效性，且可以讓儀器維護人員的工作更具有針對性。

圖1 福建臺網及臺灣交換臺站分布

1 程序設計

方位角的檢測需要用到震中距30度以上的遠震，所以該程序為準實時系統。根據具體的操作步驟，程序分為地震事件選取模塊、數據預處理模塊、臺站屬性讀取模塊、最小能量計算模塊、相關系數計算模塊及虛擬臺站構建模塊。為簡化代碼數量，減少程序運行中所占用內存，本程序使用Python中的Obspy這個開源的地震數據處理框架。

1.1 地震事件選取模塊

本程序中計算方位角偏差的方法為利用P波質點的運動特性法，需要震中距在30度以上的遠震，因此在準備地震事件時需要計算震中距是否達到30度。福建本省面積較小，對于遠震而言省內不同臺站對于同一地震事件的震中不會有太大差別，因此我們選取位于福建較為中心的永安燕西作為所有臺站的參照臺，只要計算該臺與地震事件的震中距即可，這樣可以減少運算量又不影響數據的可靠性。計算震中距我們引入Obspy.geodetics中的gps2dist_azimuth及kilometers2degrees兩個方法，前者用于計算震中距，其返回值是一個數組，震中距為第一個元素且單位為千米，所以需要使用kilometers2degrees將單位從千米轉化為度。程序將把符合震中距要求的地震事件三要素及震中距輸出至指定路徑下的TXT文件中。

圖2 程序架構

1.2 數據預處理

因寬頻帶地震計中有60秒與120秒甚至更高的頻帶，為保證帶寬一致所以需要進行帶通濾波。臺站的采樣率為100HZ，為減少計算量程序采取降采樣處理。利用Obspy.Stream類循環讀取地震離線波形數據，全部讀取后利用該類中的resample（sampling_rate）方法即可進行降采樣處理。Obspy.Stream.filter可以對數據進行帶通濾波處理，用該方法離線波形數據進行帶通濾波處理的具體參數設置如下所示 filter(‘bandpass’,freqmin=0.02,freqmax=0.2),參數表示使用帶通濾波處理數據且設置了低頻與高頻的具體數值。

1.3 臺站屬性讀取

臺站屬性指的是臺站名、經緯度及高程，這些數據均存儲于MySQL數據庫中。對于省級地震臺網中心，有固定的數據庫用于配置各臺站經緯度、高程及靈敏度等參數，所以我們用xml文件配置好數據庫的地址等，程序就可以通過xml文件進入數據庫進行臺站屬性的讀取。Xml文件的內容形式如下所示：

數據庫密碼

數據庫名

圖片存儲路徑

該模塊中我們創建一個專門讀取xml文件的類，在Python中引用xml.dom.minidom類包，就可以按照該類讀取xml文件的步驟方法，將xml文件中的各配置內容讀取到指定的數組中并返回該數組。

獲取上述數據庫的配置信息后連接MySQL數據庫，寬頻帶地震計的通道命名有顯著要求，即開頭字母為B，其三通道的命名為BHN、BHE、BHZ。因此在選取寬頻帶臺站的相關屬性時對數據庫的操作語句可寫作“select Sta_code from Channel_info where Chn_code like ‘%BHE%’”,其 中Sta_code為臺站名，Channel_info為表名稱，Chn_code為通道存儲字段名。通過以上操作可以獲取全部寬頻帶地震計臺站名，而后運用循環方法獲取所有的臺站屬性后封裝到dictionary中返回。

1.4 最小能量計算

通過以上以上步驟計算每個臺站對應每個地震事件的P波到時，選取到時前后五秒內經過旋轉的能量，再以1度為間隔進行180度的旋轉計算每一度對應的能量，再進行對比選出最小能量對應的旋轉度數，認為此時對應的角度為與正北方向的角度偏差。再將此時的偏差角與經過經緯度經過計算出的理論偏差值進行對比，即可得出此臺對應該地震事件時間點上的角度偏差。將每個地震事件的偏差值輸出至String類型的數組，計算完所有事件的偏差值后，再利用“SNR”法計算所有事件中每一度的能量值后選出最小值對應的偏差角，即認為是最終的地震計與正北方向的偏差值。

圖3 時間序列圖

為更清晰的顯示每個事件對應的方位角偏差，在此步驟中將所得到的數據進行繪圖顯示。因利用了GMT繪圖軟件所以需要在主機上安裝GMT5.0以上版本，否則會出現繪圖命令錯誤或不兼容的情況。時間序列圖的縱軸單位為的形式，所以需要用到GMT中的轉義字符“@～\152\244\050\260\051@～”。將GMT繪圖命令以String形式寫入bat文件后，執行該文件即可進行繪圖操作，繪制后的圖片將按xml文件中配置的路徑存儲。

1.5 虛擬臺陣

為更加精準的分析異常方位角偏差臺站，該程序還設計了異常臺站與正常臺站所構建的虛擬臺陣對應通道的相關系數分析模塊。此模塊以故障臺站為中心，計算所有臺站與此臺的臺間距，并選取臺間距最近的15個正常臺站構成虛擬臺陣。線性疊加臺陣內所有臺站相同通道的數據后，與故障臺相同通道進行相關系數計算，如果故障臺不出現儀器方向擺放錯誤或極性反轉等重大問題，則其相關系數會接近1相反則會接近0或負數。其相關系數的計算引入scipy.stats類中的pearsonr方法，可以迅速計算出兩個相同通道的相關系數。為了便于觀察，此模塊同樣具備輸出波形形狀及相關系數值的功能。

圖4 上圖為故障臺，下圖為正常臺與虛擬臺陣的相關系數波形圖

2 程序在福建臺網中的應用

通過該程序對福建地震臺網所屬寬頻帶地震計進行計算，地震事件為2015年1月到2019年8月之間的遠震，得到74個寬頻帶地震計臺站中有72個臺的方位角偏差值在-7°～7°之間，這樣的偏差值在該方法的理論中可以忽略不計，因此福建臺網寬頻帶地震計臺站的方位角總體合格率在97%如圖5所示。其中兩個存在較大偏差值的臺站均為BBVS-60地震計,此類型地震計沒有明顯指北標志，所以導致維護人員在更換過程中無法準確的進行指北操作，經過實地校正后已使方位角的偏差值減小到要求值1°以內。

圖5 方位角偏差臺站分布

3 結論

通過對程序模塊化開發，可以理順每個模塊的具體功能，使程序整體研發更清晰，更易調試出現問題的部分代碼進行小范圍的修復，而不至于出現牽一發動全身式的更改。因程序為處理地震數據而開發，所以使用到Obspy庫，因其幾乎支持地震學界內通常使用的所有波形格式的讀寫，可簡化大量的數據格式之間轉換的代碼，同時他還集成了大量的地震學及數學所用的專有庫，對于開發者只需要對其進行引入就可使用。程序應用在福建臺網中心后，經過實地校核發現該程序的計算結果可靠。該程序的應用彌補了福建臺網中心在方位角檢測方面的不足，且提高了方位角偏差的檢測效率及修復故障的時效性。