Matlab在天津市強震動臺網(wǎng)烈度速報中的應(yīng)用①

劉雙慶，聶永安，高武平，李雅靜

（天津市地震局，天津 300201）

0 引言

烈度速報是一項非常重要的地震災(zāi)害快速評估工作，快速產(chǎn)出準(zhǔn)確的烈度評估圖可以有效地指導(dǎo)地震救援，從而極大地減輕地震造成的二次災(zāi)害。隨著我國強震動臺網(wǎng)的建設(shè)與完善，烈度速報工作已獲得了長足的進(jìn)展，并在近幾年向地震預(yù)警工作方面進(jìn)行開拓［1］?！熬盼濉逼陂g，中國地震局建立了我國第一個由72個實時數(shù)據(jù)傳輸強震臺和80個電話撥號數(shù)據(jù)傳輸強震臺構(gòu)成的地震動強度（烈度）速報臺網(wǎng)。在“十五”期間，中國地震局除了繼續(xù)對北京市、天津市的地震烈度速報臺網(wǎng)進(jìn)行加密建設(shè)外，還在昆明、蘭州、烏魯木齊等城市新建地震烈度速報臺網(wǎng)，部分省市地震局也正在計劃建設(shè)具有地震烈度速報功能的強震臺網(wǎng)［2］。烈度速報的形式也從過去的宏觀調(diào)查逐步向儀器烈度側(cè)重［3－4］，發(fā)揮儀器烈度在震害快速評估中的作用。烈度等震線也加入了場地因素并對不同烈度段的PGA、PGV、PGD進(jìn)行了融合，這些技術(shù)方法在美國的Shakemap和日本的UrEDAS都有體現(xiàn)［5－7］。澤仁志瑪還就等震線的顯示問題進(jìn)行了討論［8］，劉雙慶［9］則利用蒙特卡羅法分析了強震臺網(wǎng)布局對等震線的繪制精度的影響。

隨著強震動臺網(wǎng)的高密度建設(shè)，不同廠家不同型號的強震儀器得到了使用，從而產(chǎn)生了多種數(shù)據(jù)回收方式及數(shù)據(jù)保存格式。如何快速地處理這些不同型號不同存儲格式的大批量強震數(shù)據(jù)成為了快速烈度速報的阻礙，就目前天津市強震動臺網(wǎng)110個強震臺而言，就涵蓋了四種類型的儀器。其部分回傳的數(shù)據(jù)還存在著壞死而無法用廠家提供的軟件打開的問題，且直接利用廠家提供的軟件無法批量處理強震數(shù)據(jù)并快速歸算烈度值。因此本文擬利用已成為國際控制界的標(biāo)準(zhǔn)計算軟件Matlab來解決上述的主要難題，并對其中的核心處理技巧進(jìn)行說明。

1 Matlab及天津市強震動臺網(wǎng)簡介

1．1 Matlab簡介

20世紀(jì)70年代，美國新墨西哥大學(xué)Cleve Moler為了減輕學(xué)生編程的負(fù)擔(dān)，用FORTRAN編寫了最早的Matlab。1984年由Little、Moler、Steve Bangert合作成立的MathWorks公司正式把Matlab推向市場。到20世紀(jì)90年代，Matlab已成為國際控制界的標(biāo)準(zhǔn)計算軟件，并在2011年推出了Matlab2011b（7．13版本）。從 Matlab6．5開始，其主要代碼已改由Java語言編寫，并與Fortran，C，C＋＋，Java，Vb等語言兼容。由于 MathWorks公司允許世界各地的專家將自己開發(fā)的Matlab子程序上傳到該公司網(wǎng)站，并由MathWorks公司的工程師對程序代碼進(jìn)行復(fù)核，通過全面測試后在新版的Matlab中追加上去，因此目前的Matlab涵蓋了多個模塊的應(yīng)用程序，主要包括：Web服務(wù)、信息通信、財政金融、控制系統(tǒng)設(shè)計、曲線擬合、數(shù)據(jù)庫使用、濾波器設(shè)計、模糊邏輯分析、仿生算法、圖像處理、地圖工具、神經(jīng)網(wǎng)絡(luò)、Mu－分析、優(yōu)化分析、Robust控制、信號處理、樣條工具箱、統(tǒng)計工具箱、符號數(shù)學(xué)、小波分析、偏微分方程分析、線性矩陣不等式控制工具箱等多個模塊，并有良好的界面操作功能和非常優(yōu)秀的使用指導(dǎo)說明。

1．2 天津市強震動臺網(wǎng)簡介

天津市強震臺網(wǎng)是在“九五”實時數(shù)據(jù)傳輸強震臺的基礎(chǔ)上，在“十五”期間重新建設(shè)起來的，總共110個強震臺（工程力學(xué)研究所的SLJ加速度儀暫停用），涵蓋了美國凱尼、瑞士GeoSig、SysCom、港震GSMA－2400IP四種類型的儀器，平均間距14．5 km。其中有31個臺（光纖通訊）具有準(zhǔn)實時波形回傳的條件，另外79臺當(dāng)觸發(fā)條件滿足后能在15分鐘內(nèi)自動回傳數(shù)據(jù)到天津市地震局?jǐn)?shù)據(jù)處理中心數(shù)據(jù)庫。基于Oracle數(shù)據(jù)庫、Java運行環(huán)境及 Web管理的強震動臺網(wǎng)實現(xiàn)了強震動數(shù)據(jù)自動匯集、存儲，準(zhǔn)實時或定時監(jiān)控儀器運行狀態(tài)和臺站遠(yuǎn)程控制與管理等功能。強震臺網(wǎng)運行以來，獲得了2008年5月汶川M8．0大地震較好的加速度記錄和2010年4月9日豐南M4．1地震加速度記錄。但此前烈度速報效率低下的情況比較嚴(yán)重。

2 天津市強震動臺網(wǎng)烈度速報的幾個技術(shù)難點及Matlab的應(yīng)用

通過分析研究，歸納天津市強震動臺網(wǎng)速報的技術(shù)難點如下：

（1）數(shù)據(jù)格式不統(tǒng)一，需要廠家的界面化軟件對事件文件進(jìn)行一個一個人工交互式處理（除個別廠家儀器數(shù)據(jù)能批處理外）。

（2）數(shù)據(jù)文件內(nèi)容存在壞死情況，軟件無法打開，用解碼程序運行則顯示錯誤。

（3）數(shù)據(jù)文件命名不規(guī)范，數(shù)據(jù)處理后存放路徑不統(tǒng)一。

（4）用 Web界面調(diào)用Oracle數(shù)據(jù)，操作不便利，影響下載速度。

（5）數(shù)據(jù)預(yù)處理操作及峰值提取無法批處理。

（6）等值線繪制方法單一且沒有融合理論烈度。

以上過程嚴(yán)重影響了烈度速報的速度和效率。針對以上難點，本文利用Matlab處理了下面幾個問題：

（1）Matlab對Oracle數(shù)據(jù)庫的直接調(diào)用

這個部分主要解決Web調(diào)用的數(shù)據(jù)不能快速下載的缺陷。由于對Oracle數(shù)據(jù)文件命名作了新的設(shè)定（圖1），將觸發(fā)事件和儀器標(biāo)定事件的命名做了區(qū)分，因此可以通過時間段和事件命名的差異直接提取出Oracle某個時間段的地震觸發(fā)事件文件或標(biāo)定文件。其涉及的Matlab代碼有：

當(dāng)在本地機器安裝好Oracle客戶服務(wù)端后，配置Matlab對Oracle的驅(qū)動，通過以上代碼就可以直接從遠(yuǎn)程Oracle數(shù)據(jù)庫批量讀取數(shù)據(jù)文件。

（2）Matlab與exe文件的融合

該部分主要是充分利用廠家提供的exe格式的解碼程序，對從Oracle數(shù)據(jù)庫讀出的數(shù)據(jù)文件進(jìn)行批量轉(zhuǎn)換成ascii碼。由于臺站代碼和命名固定，所以直接通過文件名就可以確定采用哪個廠家的exe文件進(jìn)行數(shù)據(jù)轉(zhuǎn)換。命名方法見圖1，如果是標(biāo)定文件則在上述文件名的基礎(chǔ)上加前綴FT＿，即function test的縮寫。

Matlab與exe的調(diào)用可以通過Matlab提供的system或dos兩個命令來操作．由于存在數(shù)據(jù)文件內(nèi)容壞死而無法解碼的情況，因此可以進(jìn)一步使用Matlab提供的try．．．catch．．．end的命令來嘗試解碼，當(dāng)解碼不成功時直接跳過，從而保證后續(xù)程序的連續(xù)運行。

圖1 強震事件文件命名規(guī)則Fig．1 Naming rule of strong motion record file．

（3）Matlab的批量顯示與通道選擇

這部分可通過Matlab界面化函數(shù)uicontrol（gcf，‘style’，‘slider’，‘unit’，‘normalized’，‘position’），uicontrol（gcf，‘style’，‘listbox’，‘Position’）組合的方式，提供一個左側(cè)的滑動選擇框及右側(cè)的圖形曲線顯示框（圖5）來確認(rèn)哪道數(shù)據(jù)可以使用。

（4）儀器烈度的換算

該部分先對第（3）部分選擇出來的數(shù)據(jù)道保存，取每道觸發(fā)點前20s的數(shù)據(jù)均值作為基線，以去除各道基線漂移。然后利用ButterWorth帶通濾波器壓制低頻和高頻的部分波形能量，考慮強震儀器的寬頻與高通特性以及建筑物反應(yīng)譜設(shè)計曲線［10］，本文通頻帶設(shè)為20Hz～2s。之后程序?qū)⑻崛∶總€臺站兩個水平道上的第20個最大點作為有效峰值，按每套儀器滿量程時對應(yīng)的gal值，計算出每個儀器兩水平道第20個最大點對應(yīng)的gal值的均值。提取絕對值排序后的第20個最大值作為峰值，可以有效避免儀器尖脈沖的影響（單道記錄有時存在多個尖脈沖干擾），同時因為采樣率為100sps，不會嚴(yán)重漏用記錄峰值及其附近區(qū)數(shù)據(jù)。程序經(jīng)多次測試，該選點一般可取第5～25點，換算出的gal值較穩(wěn)定。然后按中國地震烈度表 GB／T 17742－2008［11］中烈度與gal的數(shù)值表，采用線性內(nèi)插的方式（Matlab命令：interp1）獲取gal對應(yīng)的儀器烈度值。當(dāng)然interp1命令還提供了最近鄰、三次多項式、樣條法等多種插值方法。由于烈度表不提供V以下的PGA參考值，因此這些外推的低烈度區(qū)儀器等值線僅作參考，且其本身因烈度低，其災(zāi)評意義亦不明顯（Ⅵ以上才有明顯地震災(zāi)情）。另外，本文亦嘗試?yán)妹绹鳤TC－3規(guī)范中的加速度有效峰值反應(yīng)譜EPA（線性化方法計算程序見附錄）轉(zhuǎn)換儀器烈度，但由于地震樣本量太少，尚未能深入分析本地區(qū)PGA，EPA的具體差異。

（5）Matlab的儀器烈度顯示

對于烈度的顯示，除了利用Delaunay三角剖分三角形顯示之外［9］，本部分還提供了兩種繪制方案．一是利用Matlab的scatter命令來對不同空間點上的烈度值按數(shù)據(jù)大小決定顯示點大小及顏色的離散方式顯示（圖7（a）），這種方式比等值線法更能體現(xiàn)單點的儀器烈度。二是利用華北地區(qū)理論烈度衰減模型（需要人工輸入地震震中的經(jīng)緯度、震級和長軸方向這些參數(shù)，圖4（c））構(gòu)建大網(wǎng)格點（本地區(qū)取0．2°）的理論烈度場，然后將儀器記錄烈度值加密并取代附近的理論烈度，之后利用griddata命令對非規(guī)則的加密后的數(shù)據(jù)插值，生成烈度等值線。這種方法既利用了儀器烈度也使用了理論烈度，是一種綜合效果（圖7（b））。而地理底圖的顯示則直接利用Matlab與arcgis地圖格式兼容的命令shaperead和mapshow進(jìn)行地理圖讀取與顯示。當(dāng)然，低烈度區(qū)（IV以下）的理論計算值亦僅作為畫圖參考，其本身災(zāi)評意義較弱。

（6）可移植性問題

考慮程序的可移植性，程序使用過程中所涉及的程序路徑都是相對路徑，不存在放在不同路徑無法運行的問題。其次，對于臺站名命名、地理底圖、廠家提供的解碼程序等靜態(tài)資料，程序統(tǒng)一從一個cfg文件夾獲取，cfg隨整個程序綁定，只要修改cfg的內(nèi)容并按圖1的命名方法對數(shù)據(jù)文件命名，且替換相應(yīng)的地理底圖，則本文編寫的程序可以移植到其他省臺網(wǎng)使用。同時在最后烈度成圖的操作上本文采用了調(diào)用計算過程中生成的文件內(nèi)容方式，因此當(dāng)有野外實地烈度調(diào)查結(jié)果時可以補充到這個烈度結(jié)果文件中，從而進(jìn)一步突出烈度等值線對災(zāi)情現(xiàn)場烈度的反映。

利用Matlab人工處理天津市強震臺網(wǎng)數(shù)據(jù)與烈度顯示流程及使用的Matlab主要命令如圖2所示。

（7）“事件偵聽”

由于地震是突發(fā)事件，為了保證烈度速報的產(chǎn)出速度，本文的Matlab程序除了提供人機交互處理的方式外，還提供了自動觸發(fā)處理方式——“事件偵聽”（圖4）。這是建立在天津市強震動臺網(wǎng)數(shù)據(jù)上傳的具體特征的基礎(chǔ)上。一般地，地震到來后在5～13分鐘內(nèi)觸發(fā)的臺站基本把數(shù)據(jù)自動上傳到了中心數(shù)據(jù)庫。而無地震的情況下，上傳到數(shù)據(jù)庫的文件在5～13分鐘內(nèi)一般小于2個，因此可以每5分鐘對數(shù)據(jù)庫上傳事件的個數(shù)進(jìn)行統(tǒng)計，如果超過一定的限度，則認(rèn)為有叢集事件記錄上傳，并按如下4個步驟進(jìn)行自動處理（圖3）。

圖2 Matlab處理天津市強震動臺網(wǎng)數(shù)據(jù)流程圖Fig．2 Flow chart of dealing with the data of Tianjin strong motion network using Matlab．

圖3 自動掃描處理示意圖Fig．3 Sketch of auto－scaning and auto－processing．

①如果有地震事件發(fā)生，下次循環(huán)時間由5分鐘改為4．5分鐘并等待8分鐘，預(yù)留0．5分鐘提供程序自動處理前14分鐘的數(shù)據(jù)，并且提取滿足掃描條件時刻的前1分鐘的數(shù)據(jù)，因此總共處理8＋1＋5（其中5＝4．5＋0．5）分鐘內(nèi)上傳到數(shù)據(jù)庫里的數(shù)據(jù)。當(dāng)沒有重復(fù)事件在13（8＋5）分鐘內(nèi)發(fā)生時，數(shù)據(jù)將提取14分鐘的內(nèi)容，與上一個5分鐘的循環(huán)有1分鐘的重疊。目前我們設(shè)定5分鐘內(nèi)上傳數(shù)大于5個則觸發(fā)，這個約束條件將可能丟失4個文件，因為1分鐘內(nèi)最多可以上傳4個文件。提前一分鐘的重疊，至少可以挽回2個文件。

②如果事件發(fā)生后，在13分鐘之后又發(fā)生一次地震事件（即第13～18分鐘再次出現(xiàn)叢集事件上傳，則提取這個后續(xù)事件的前14分鐘的事件文件附加到上一次事件中，此時有1．5分鐘的重疊）?？偣蔡幚?4＋14－1．5＝26．5分鐘內(nèi)上傳的文件數(shù)據(jù)。

③如果在18分鐘后才再次發(fā)生地震事件叢集上傳，則作為第二個地震事件處理。

④如果不發(fā)生其他事件，則正常進(jìn)入下一個5分鐘的循環(huán)。

經(jīng)過一年的實踐運行結(jié)果顯示，未出現(xiàn)誤自動處理現(xiàn)象。當(dāng)“設(shè)定5分鐘內(nèi)上傳數(shù)大于2個”時，偶爾出現(xiàn)自動處理結(jié)果，這些結(jié)果常是由儀器工作不正常，或職工當(dāng)時進(jìn)行儀器維修造成。

3 應(yīng)用實例

本實例采用2010年4月9日豐南M4．1地方震加速度記錄，觸發(fā)臺站36個，可利用的數(shù)據(jù)臺27個（由于當(dāng)時外圍網(wǎng)絡(luò)數(shù)據(jù)傳輸程序有缺陷，SYSCOM公司的儀器記錄文件壞死情況較嚴(yán)重）。地震發(fā)生后，當(dāng)時數(shù)據(jù)處理過程耗用了多個小時。本程序編寫完成后，在3分鐘內(nèi)就可以通過界面化（圖4～6）的鼠標(biāo)操作完成了所有的操作，其中點擊“人工數(shù)據(jù)提取”可以彈出圖4（b）的參數(shù)輸入界面。數(shù)據(jù)繪制的烈度顯示結(jié)果（圖7）也與當(dāng)時各區(qū)縣群眾打電話詢問的頻度很相近，其中寧河縣，塘沽區(qū)和市中心高層住戶有明顯震感。需說明的是，圖5、圖6的縱坐標(biāo)做了如下處理：先對各臺站記錄數(shù)據(jù)歸一化以壓制異常記錄，然后對每臺記錄的基線遞增進(jìn)行顯示以避免曲線堆疊。數(shù)據(jù)提取按內(nèi)存變量值為準(zhǔn)。

圖4 天津強震動臺網(wǎng)速報軟件界面（a），獲取事件選擇框（b）和理論烈度計算參數(shù)輸入框（c）Fig．4 GUI of intensity rapid report of Tianjin strong motion（a），GUI of selecting events（b）and GUI of input parameters to calculate theoretical intensity（c）．

4 認(rèn)識與討論

本文將Matlab軟件提供的調(diào)用數(shù)據(jù)庫命令，調(diào)用儀器廠家的exe命令，數(shù)據(jù)分析處理命令和界面化操作等命令進(jìn)行有機組合，從而比較成功地解決了天津市強震動臺網(wǎng)在烈度速報中遇到的速度與效率問題，大大縮減了儀器烈度發(fā)布的時間，由此得到以下認(rèn)識：

圖5 Uicontrol函數(shù)實現(xiàn)的數(shù)據(jù)顯示及通道選擇操作Fig．5 Data imaging and selecting with uicontrol function．

圖6 Uicontrol函數(shù)實現(xiàn)的數(shù)據(jù)顯示局部放大（由右滑鍵控制）Fig．6Local zooming－in with uicontrol function（Operated by right－side slider）．

（1）Matlab簡明的程序代碼書寫格式和豐富的命令式程序庫，可以有效降低應(yīng)用程序開發(fā)難度，大大縮減程序代碼量。但Matlab軟件再強大，也需要和研究者所研究對象的物理模型、數(shù)學(xué)模型充分融合才能發(fā)揮其作用。本文中涉及的理論烈度計算的非線性迭代方程，數(shù)據(jù)ButterWorth濾波頻段設(shè)計，“事件偵聽”的循環(huán)判斷方案等等都是具體物理模型的數(shù)學(xué)化，是Matlab本身不具備的。

圖7 豐南4．1級地震的儀器烈度散點顯示（a）和等值線顯示（b）Fig．7 Equipment intensity imaging by scattering（a）and contouring method（b）for M4．1earthquake at Fennan county．

（2）儀器烈度與傳統(tǒng)烈度的“多因子綜合判定”、“宏觀分級”、“以結(jié)果表述原因”三個特征尚有一定的銜接困難，本文的PGA直接轉(zhuǎn)化儀器烈度亦有較多值得商榷之處，利用EPA方法尚還在嘗試，由于缺乏足夠的地震樣本未能深入研究其與PGA的差別。日本目前已不再過分強調(diào)儀器烈度與傳統(tǒng)烈度的災(zāi)害評估差異，并減少現(xiàn)場宏觀調(diào)查的工作量，突顯儀器烈度在快速災(zāi)評的指導(dǎo)，他們這種處理方式不失為一種明智之舉［3，12］。

由于在本文中并未對天津地區(qū)場地進(jìn)行烈度校正，且沒有進(jìn)一步討論利用反應(yīng)譜的結(jié)果來刻畫儀器烈度的可行性，因此儀器烈度速報的結(jié)果還是比較粗糙。另外，Matlab的通訊工具箱直接對串口實時數(shù)據(jù)，對數(shù)字采集器直接訪問的研究也未深入展開，在人機交互處理界面上也不太美觀，這些工作都是后續(xù)需要努力的。

致謝：本文程序編寫過程中得到了長安大學(xué)鄧亞虹博士，北京航空航天大學(xué)吳鵬博士的幫助，在此非常感謝！

［1］ 李山有，金星，馬強，等．地震預(yù)警系統(tǒng)與智能應(yīng)急控制系統(tǒng)研究［J］．世界地震工程，2004，20（4）：21－26．

［2］ 金星，廖詩榮，陳緋雯．區(qū)域數(shù)字地震臺網(wǎng)實時速報系統(tǒng)研究［J］．地震地磁觀測與研究，2007，28（1）：64－72．

［3］ 張敏政．地震烈度及其評定［J］．防災(zāi)科技學(xué)院學(xué)報，2010，12（1）：1－6．

［4］ 金星，張紅才，韋永祥．基于地震臺網(wǎng)資料快速發(fā)布的震動烈度標(biāo)準(zhǔn)及其應(yīng)用研究［J］．國際地震動態(tài)，2008，（10）：20－27．

［5］ 李俊，蘇楓，米宏亮，等．ShakeMap及其在地震動快速預(yù)估中的應(yīng)用［J］．中國地震，2010，26（1）：103－111．

［6］ D J Wald，V Quitoriano，T H Heaton，et al．TriNet“ShakeMaps”：Rapid Generation of Peak Ground Motion and Intensity Maps for Earthquakes in Southern California［J］．Earthquake Spectra，1999，15（3）：537－555．

［7］ K T Shabestari，F(xiàn) Yamazaki．A proposal of instrumental seismic intensity scale compatible with MMI evaluated from three－component acceleration records［J］．Earthquake Spectra，2001，17（4）：711－723．

［8］ 澤仁志瑪，陳會忠，何加勇，等．震動圖快速生成系統(tǒng)研究［J］．地球物理學(xué)進(jìn)展，2006，21（3）：809－813．

［9］ 劉雙慶，謝靜，聶永安，等．Delaunay三角剖分算法在地震烈度速報中的應(yīng)用－－以天津市強震臺網(wǎng)為例［J］．地震地磁觀測與研究，2011，32（5）：115－122．

［10］ 國家標(biāo)準(zhǔn)建筑抗震設(shè)計規(guī)范管理組 編．建筑抗震設(shè)計規(guī)范（GB50011－2010）統(tǒng)一培訓(xùn)教材［M］．北京：地震出版社，2010．

［11］ 中國國家標(biāo)準(zhǔn)化管理委員會．中國地震烈度表GB／T 17742－2008［S］．北京：中國標(biāo)準(zhǔn)出版社，2009．

［12］ 顏心儀．利用臺灣寬頻地震網(wǎng)從事強震預(yù)警研究［D］．臺灣?。号_灣大學(xué)地質(zhì)科學(xué)研究所，2007．

附錄：Newmark方法計算地震動反應(yīng)譜子程序（Matlab版）：