吳華燈1, 2, 3) 謝劍波1, 2, 3)

1）廣東省地震局，廣州 510070

2）中國地震局地震監(jiān)測與減災(zāi)技術(shù)重點實驗室，廣州 510070

3）廣東省地震預(yù)警與重大工程安全診斷重點實驗室，廣州 510070

EDAS-C24型數(shù)字測震儀實時數(shù)據(jù)的解碼及加速度和位移仿真的實現(xiàn)

吳華燈1, 2, 3) 謝劍波1, 2, 3)

1）廣東省地震局，廣州 510070

2）中國地震局地震監(jiān)測與減災(zāi)技術(shù)重點實驗室，廣州 510070

3）廣東省地震預(yù)警與重大工程安全診斷重點實驗室，廣州 510070

本文嘗試使用了一種圖形編程語言進行數(shù)據(jù)解碼和仿真研究，闡述了用該語言實現(xiàn)EDAS-C24型數(shù)字測震儀實時數(shù)據(jù)流的解碼過程，提出了在G語言下利用數(shù)字濾波器逼近模擬積分器及模擬微分器響應(yīng)實現(xiàn)對解碼數(shù)據(jù)實時仿真的方法，并通過設(shè)定頻帶寬度，比較了設(shè)計的補償濾波響應(yīng)和實際的幅頻響應(yīng)的一致性。結(jié)果表明，在設(shè)定的頻帶內(nèi)，補償濾波響應(yīng)和實際的幅頻響應(yīng)是一致的，仿真的精度是理想的，成果已經(jīng)在廣東省地震科普館的地震互動區(qū)部署運行，取得了一定的實效。

LabVIEW G編程語言 測震儀 實時數(shù)據(jù) 解碼 仿真

引言

Kanamori等（1999）利用美國TriNet強震臺網(wǎng)，曾對實時速度記錄通過單邊差分方法計算位移等問題做過研究，但是其計算精度和相位仍存在問題；金星等（2003）利用單自由度體系地震反應(yīng)的時域遞歸方法，系統(tǒng)研究了由數(shù)字加速度記錄實時仿真地動速度和地動位移的問題，從理論和數(shù)值計算上取得了理想的結(jié)果；同時金星等（2004）還系統(tǒng)研究了利用數(shù)字化寬頻帶速度型記錄仿真地面位移和地面加速度等問題，從時域上提出了一套實時計算公式；謝劍波（2014）也用時域遞歸方法做過不同帶寬的地震記錄的仿真研究。以上的實現(xiàn)都是利用地震臺網(wǎng)匯集的實時數(shù)據(jù)或事件數(shù)據(jù)，在傳統(tǒng)文本編程語言的環(huán)境下實現(xiàn)數(shù)據(jù)的仿真。本文將基于LabVIEW圖形編程語言——G編程語言，在數(shù)據(jù)采集現(xiàn)場，直接解碼數(shù)據(jù)采集器的實時速度數(shù)據(jù)，并設(shè)計合適的數(shù)字濾波器，實現(xiàn)加速度和位移的實時仿真，最終繪制出速度、加速度和位移三種物理量的實時波形。

1 實時數(shù)據(jù)解碼及仿真開發(fā)環(huán)境的搭建

1.1 數(shù)據(jù)采集現(xiàn)場設(shè)備的架設(shè)

地震動數(shù)據(jù)采集現(xiàn)場如圖1所示，由EDAS-C24型地震數(shù)據(jù)采集器、FBS-3A型反饋式寬頻帶地震計、GPS授時系統(tǒng)、12V電源供電系統(tǒng)和集成開發(fā)電腦平臺等設(shè)備架設(shè)而成。FBS-3A型反饋式寬頻帶地震計是機電一體化的測震儀器，它在自振周期為2s的機械拾震器上附加電子線路，構(gòu)成測量頻帶約為0.05—40Hz的寬頻帶測量系統(tǒng)。在儀器結(jié)構(gòu)上，它由一個豎直向和兩個相互正交的水平向組成（劉慶偉等，2001）。該儀器用于拾取地面的震動信號，并把信號轉(zhuǎn)換成電壓值，供數(shù)據(jù)采集器數(shù)字化。EDAS-C24型地震數(shù)據(jù)采集器將FBS-3A型地震計輸出的模擬電壓信號轉(zhuǎn)換成數(shù)字信號，并經(jīng)串口線把數(shù)字信號傳輸?shù)诫娔X，供計算機程序處理。

圖1 數(shù)據(jù)采集現(xiàn)場設(shè)備的架設(shè)Fig.1 The set-up of equipment

1.2 數(shù)據(jù)解碼及仿真環(huán)境

在圖1所示的集成開發(fā)平臺電腦上，安裝了LabVIEW2010集成開發(fā)軟件。LabVIEW是Laboratory Virtual Instrument Engineering Workbench（實驗室虛擬儀器集成環(huán)境）的簡稱，是一個功能強大而又靈活的圖形開發(fā)環(huán)境（喬瑞萍等，2008）。由于EDAS-C24型地震數(shù)據(jù)采集器是基于串口通信，本文用LabVIEW中的VISA（VISA是Virtual Instruments Software Architecture的縮寫，中文譯為虛擬儀器軟件架構(gòu)或可視化儀器軟件架構(gòu)）編程技術(shù)，建立與數(shù)據(jù)采集器串口的通信，實現(xiàn)實時數(shù)據(jù)的接收，為數(shù)據(jù)解碼和仿真提供數(shù)據(jù)源。

2 實時數(shù)據(jù)解碼的實現(xiàn)

2.1 EDAS-C24型地震數(shù)據(jù)采集器實時數(shù)據(jù)格式

EDAS-C24型地震數(shù)據(jù)采集器實時數(shù)據(jù)有壓縮和非壓縮兩種數(shù)據(jù)格式，本文主要針對非壓縮數(shù)據(jù)格式進行解碼。該型號數(shù)據(jù)采集器實時數(shù)據(jù)幀是由幀頭、數(shù)據(jù)區(qū)、校驗和組成。幀頭長16個字節(jié)，包括幀同步標(biāo)識碼、臺站編號、幀類型標(biāo)示字、幀長、秒記數(shù)、標(biāo)志字節(jié)、頭段區(qū)字節(jié)檢查和；數(shù)據(jù)區(qū)包括3個通道的數(shù)據(jù)，當(dāng)采樣率是125Hz時，還有一個字節(jié)的填充字；幀的尾部是兩個字節(jié)的檢查和。數(shù)據(jù)幀格式見表1，數(shù)據(jù)幀的具體內(nèi)容見圖2。從表1中可看出，在EDAS-C24中，字節(jié)序是小端，每一個數(shù)據(jù)幀都是以“74 97 13 BF”同步標(biāo)識碼開頭，BF是高字節(jié)位。圖2中的數(shù)據(jù)幀內(nèi)容，是在計算機程序緩沖區(qū)中的字節(jié)序。因此，程序解碼時，將用小端的字節(jié)序來判斷數(shù)據(jù)幀和提取3通道的數(shù)據(jù)。

表1 EDAS-C24數(shù)據(jù)幀格式Table 1 Data format of EDAS-C24

續(xù)表

圖2 EDAS-C24數(shù)據(jù)幀內(nèi)容Fig.2 Content of EDAS-C24 data frame

2.2 實時數(shù)據(jù)解碼

解碼是在LabVIEW2010專業(yè)版開發(fā)系統(tǒng)下，采用狀態(tài)機、事件隊列消息處理器和事件結(jié)構(gòu)結(jié)合的架構(gòu)進行編程，用到了While循環(huán)、條件結(jié)構(gòu)和事件結(jié)構(gòu)，架構(gòu)代碼如圖3所示。

圖3 程序架構(gòu)Fig.3 Program frame

狀態(tài)機里定義了Init、No Action、Receive Data和Exit四種事件狀態(tài)及對應(yīng)的事件分支，優(yōu)先執(zhí)行Init事件分支，用于初始化串口名稱、串口波特率、采樣率、滿量程、事件觸發(fā)模式、權(quán)重等配置參數(shù)和讀取串口的數(shù)據(jù)，為數(shù)據(jù)的解碼提供數(shù)據(jù)源。通過獲取下一個事件的子VI（圖形編程語言的子程序），依次改變事件狀態(tài)為No Action、Receive Data和Exit。在No Action事件分支里，用了事件結(jié)構(gòu)。事件結(jié)構(gòu)的超時端子為50ms，有“超時”、“前面板關(guān)閉”、“菜單選擇”、“運行狀態(tài)”4個分支，主要用到了“超時”事件分支。當(dāng)程序運行超過50ms，將執(zhí)行“超時”事件框中的程序，如果收到串口的數(shù)據(jù)并且程序處于運行狀態(tài)，則驅(qū)動Receive Data事件，否則置錯誤信息。Receive Data事件分支的程序是解碼和仿真EDAS-C24型地震數(shù)據(jù)采集器實時數(shù)據(jù)的核心。在解碼方面，程序首先由圖標(biāo)為GET FRAME、名稱為Receive Frame的子VI，從VISA資源里找出數(shù)據(jù)幀頭。然后，通過Parse Frame Header子VI分析數(shù)據(jù)幀頭的有效性。如果數(shù)據(jù)幀頭是有效的，則進一步求解出三個通道的原始速度數(shù)據(jù)，解碼的具體實現(xiàn)過程如下所述。

2.2.1 數(shù)據(jù)幀頭解碼的實現(xiàn)

通過隧道和移位寄存器把Init事件分支初始化時的VISA資源輸入到GET FRAME子VI的VISA resource name入口（見圖4），然后利用“VISA讀取”函數(shù)，讀取16個字節(jié)大小的數(shù)據(jù)幀頭到header buffer緩沖區(qū)，接著將header buffer緩沖區(qū)的數(shù)據(jù)輸入圖5所示的Parse Frame Header子VI，實現(xiàn)數(shù)據(jù)幀頭Frame Header的輸出及其有效性的判斷。

圖4 GET FRAME子VIFig.4 GET FRAME sub VI

圖5 Parse Frame Header子VIFig.5 Parse Frame Header sub VI

“從字符串還原”函數(shù)有類型、二進制字符串、字節(jié)序等參量的輸入。類型就是按照“表1 EDAS-C24數(shù)據(jù)幀格式”以簇方式定義的幀頭；二進制字符串就是header buffer；由于在EDAS-C24中，字節(jié)序是小端，所以字節(jié)序選擇little endian。該函數(shù)的輸出結(jié)果是數(shù)據(jù)幀頭Frame Header。數(shù)據(jù)幀頭Frame Header有效性的判斷用了3個判斷條件：其一，從數(shù)據(jù)幀頭Frame Header的簇中求出幀長，判斷幀長是否小于1000；其二，取Frame Header前四個字節(jié)，判斷是否與“BF 13 97 74”相等；其三，取Frame Header第七、第八個字節(jié)，判斷是否與“AA55”相等，當(dāng)3個判據(jù)同時滿足時，認(rèn)為Frame Header是有效的，至此實現(xiàn)了數(shù)據(jù)幀頭解碼。

2.2.2 數(shù)據(jù)區(qū)解碼的實現(xiàn)

在解碼出數(shù)據(jù)幀頭和確認(rèn)幀頭有效后，程序進入Case結(jié)構(gòu)。依據(jù)EDAS-C24數(shù)據(jù)幀格式，數(shù)據(jù)區(qū)的字節(jié)數(shù)byte count等于數(shù)據(jù)幀長度減8個字節(jié)，從VISA resource out讀取byte count個字節(jié)的數(shù)據(jù)到讀寫緩沖區(qū)read buffer。當(dāng)采樣率為100bps時，每個通道每個采樣點有3個字節(jié)，3個通道每個采樣點一共9個字節(jié)?？梢杂靡粋€循環(huán)，循環(huán)的次數(shù)為byte count字節(jié)/9字節(jié)，在循環(huán)體里結(jié)合“索引數(shù)組”函數(shù)、24bit子VI和“創(chuàng)建數(shù)組”函數(shù)，分解出三個通道的數(shù)據(jù)，如圖6所示。“索引數(shù)組”函數(shù)的作用是建立每個通道LMH的索引；24bit子VI的作用是將每個通道的LMH重新組合成24位的數(shù)據(jù)，并轉(zhuǎn)換成計算機的32位數(shù)據(jù)，再右移8位，得到每個通道的最終數(shù)據(jù)；“創(chuàng)建數(shù)組”函數(shù)的作用是把每個通道的最終數(shù)據(jù)添加到數(shù)組，為繪制波形提供數(shù)據(jù)源。

圖6 三個通道數(shù)據(jù)的解碼Fig.6 Decoding of three channels data

3 實時數(shù)據(jù)仿真

3.1 實時數(shù)據(jù)的仿真方法

實時數(shù)據(jù)的仿真，一般通過計算單自由度系統(tǒng)的地震反應(yīng)即可得到（胡聿賢，2006）。單自由度系統(tǒng)地震反應(yīng)的計算方法可分為二大類：一類是在頻域內(nèi)進行；另一類是在時域下進行。由于頻域分析方法必須在獲取完整的地震記錄后才能進行，不適合于實時計算的要求，人們一般采用時域方法進行仿真計算。即用時間上有限個離散的動力平衡方程來替換連續(xù)形式的動力平衡方程，進而在一定假定下依據(jù)離散動力平衡方程來尋找數(shù)字濾波器，使得所尋找的數(shù)字濾波器的傳遞函數(shù)逼近理論傳遞函數(shù)。同時，利用數(shù)字濾波尤其是遞歸濾波的特點實現(xiàn)實時計算。

3.1.1 常用的時域計算方法

有關(guān)單自由度系統(tǒng)地震反應(yīng)的時域計算方法有很多，如中心差分方法、New-mark方法、變換方法、Duhamel逐步積分法和wilson-0法等。許多學(xué)者開展了相關(guān)的研究工作，并取得了很好的成果。如金星等（2004）的分析與研究，得到了計算單自由度系統(tǒng)相對位移、相對速度、相對加速度的遞歸公式，并求得由地面速度時程計算位移時程的公式：

式中，b1，b2是遞歸系數(shù)；xk-3/2和xk-5/2是前兩個輸出點位移；S0是與傳遞函數(shù)的零頻約束條件有關(guān)的參數(shù)；Δt是離散時間間隔；δ 是與Δt/T0有關(guān)的參數(shù)；T0是自振周期；Vk是k點的速度。

以及求得由地面速度時程計算相對加速度時程的公式：

利用式（1）和式（2）進行了實時仿真位移與加速度時程的計算研究。

3.1.2 本文作者提出及使用的仿真方法

式中，i=0，1，2…，n-1，n是采樣數(shù)；0x是首端，xΔ=(末端-首端)/m，不包括末端時，m=n，否則，m=n–1。Sinc信號的輸入采樣數(shù)為101，tΔ=0.005，如用序列y表示Sinc信號，可依據(jù)下式生成Sinc信號：

輸出的Sinc信號一組取輸出值，另一組取輸出值的倒數(shù)。斜坡信號與兩者捆綁成簇數(shù)組后，按照設(shè)定的頻帶0.4Hz，提取一定數(shù)量的數(shù)據(jù)，重新組合成簇數(shù)組，再由Remez的交換方法建立equi-ripple等紋波濾波器，并輸出實際紋波值。等文波濾波器的濾波輸出經(jīng)過DFD plot Freq Response子VI，輸出并繪制補償濾波響應(yīng)，再將filter out輸入到獲取傳遞函數(shù)的子VI，取出系統(tǒng)的傳遞函數(shù)。接著，由給定的頻率和截止頻率，計算增益、零點、極點參數(shù)，通過DFD build filter from zeros-poles-gain子VI設(shè)計出新的濾波器。最后，級聯(lián)兩個濾波器，再經(jīng)DFD plot Freq Response子VI得到最終的濾波輸出，實現(xiàn)逼近模擬積分器響應(yīng)的設(shè)計。

微分濾波器的設(shè)計，也是首先生成斜坡信號和Sinc信號，再由Remez數(shù)字濾波器設(shè)計方法設(shè)計任意頻響的數(shù)字濾波器。任意頻響的數(shù)字濾波器經(jīng)DFD Get Transfer Function子VI后，取出傳遞函數(shù)。利用傳遞函數(shù)，通過DFD Build Filter From Transfer Function子VI，建立新的濾波器。最后，級聯(lián)兩個濾波器，再經(jīng)DFD plot Freq Response子VI得到最終的幅頻響應(yīng)和濾波輸出，實現(xiàn)逼近模擬微分器響應(yīng)的設(shè)計。

3.2 實時數(shù)據(jù)仿真的實現(xiàn)

如圖7所示，微分濾波器和積分濾波器的濾波輸出分別接入DFD Filter子VI，速度數(shù)據(jù)便可實時仿真成加速度和位移數(shù)據(jù)，同時將加速度、速度和位移數(shù)據(jù)在波形圖上顯示出來，見圖8。

圖7 仿真信號輸出Fig.7 Output of simulation signal

圖8 加速度、速度和位移數(shù)據(jù)波形圖Fig.8 Waveform of acceleration, velocity and displacement

4 仿真數(shù)據(jù)的驗證

在進行濾波器設(shè)計過程中，為了觀察仿真的精度，設(shè)定了頻帶寬度，將實際的幅頻響應(yīng)和設(shè)計的補償濾波響應(yīng)繪制了出來（見圖9）。此外，還輸出了實際紋波值和繪制了紋波波形（見圖10），輸出了濾波系數(shù)，繪制了最終的幅頻響應(yīng)。從圖中可看出，補償濾波響應(yīng)在設(shè)定的頻帶內(nèi)和實際的幅頻響應(yīng)是一致的，設(shè)計的紋波值在0.005左右，從設(shè)計的原型上看，仿真的精度是理想的。

圖9 設(shè)計的補償濾波響應(yīng)Fig.9 Designed compensation filter response

圖10 設(shè)計的紋波波形輸出Fig.10 Output of designed ripple waveform

為了驗證仿真數(shù)據(jù)的正確性及精度，在同一臺址架設(shè)了珠海泰德企業(yè)有限公司生產(chǎn)的強震動數(shù)據(jù)采集器TDE-324CA和力平衡加速度計TDA-33M，進行仿真加速度記錄與實測加速度記錄對比觀測實驗，對比觀測所使用的儀器信息如表2所示。

從傅立葉幅值譜（圖12）和自相關(guān)系數(shù)（圖13）的結(jié)果可以看出，仿真的加速度數(shù)據(jù)在頻帶和相位上與實測的加速度數(shù)據(jù)吻合較好，說明本文設(shè)計的濾波器是符合預(yù)期的，具有一定的可靠性。

5 結(jié)論

本文在解碼出數(shù)據(jù)采集器的速度數(shù)據(jù)后，利用數(shù)字濾波器逼近模擬積分器和模擬微分器響應(yīng)的方法，實現(xiàn)了由速度數(shù)據(jù)實時仿真位移和加速度數(shù)據(jù)，數(shù)據(jù)的仿真精度是理想的。但本文設(shè)計的濾波器更多的是考慮加速度數(shù)據(jù)所關(guān)心的頻帶，仿真系統(tǒng)在有效頻帶方面尚應(yīng)做更多的分析研究。此外，在沒有實測位移數(shù)據(jù)的情況下，還應(yīng)考慮利用現(xiàn)有的一些典型地震事件數(shù)據(jù)進行相關(guān)性試算驗證，計算不同阻尼比、不同自振周期等條件的結(jié)果，逐步修正和優(yōu)化濾波器，使其具有實用性，方便應(yīng)用于地震觀測數(shù)據(jù)的實時仿真。本文的成果已經(jīng)在廣東省地震科普館的地震互動區(qū)部署運行，取得了一定的實效。

胡聿賢，2006．地震工程學(xué)（第二版）．北京：地震出版社.

金星，馬強，李山有，2004．利用數(shù)字化速度記錄實時仿真位移與加速度時程．地震工程與工程振動，24（6）：9—14.

金星，馬強，李山有，2003．四種計算地震反應(yīng)數(shù)值方法的比較研究．地震工程與工程振動，23（1）：18—30.

劉慶偉，莊燦濤，劉慧寧，2001．FBS-3A型反饋式寬帶地震計的傳遞函數(shù)．地震學(xué)報，23（1）：79—86

喬瑞萍，林欣等譯，2008．LabVIEW 8實用教程．北京：電子工業(yè)出版社.

萬永革，2007．?dāng)?shù)字信號處理的MATLAB實現(xiàn)．北京：科學(xué)出版社.

謝劍波，2014．地震記錄的時間域反褶積仿真及在地震計方位角相對測量中的應(yīng)用．地球物理學(xué)報，57（1），167—178．

Kanamori H.，Maechling P.，Hauksson E.，1999．Continuous monitoring of strong motion parameters. Bull. Seism. Socmer.，89 (1)：311—316.

Implement for EDAS-C24 Digital Seismograph Real-time Data Decoding and Simulation of Ground Displacement and Acceleration

Wu Huadeng1,2,3)and Xie Jianbo1,2,3)

1) Earthquake Administration of Guangdong Province，Guangzhou 510070，China
2) Key Laboratory of Earthquake Monitoring and Disaster Mitigation Technology，Guangzhou 510070，China
3) Key Laboratory of Guangdong Province，Earthquake Early Warning and Safety Diagnosis of Major Projects, Guangzhou 510070，China

This paper outlines a new research on data decoding and simulation using graph programming language, which can realize decoding real-time stream of EDAS-C24 digital seismometer. A new method is proposed to realize real-time simulation for decoded data by using digital filter approach to the response of analog integrator and differentiator with G programming language. Through setting the bandwidth, we compared the designed compensation filter response with the actual amplitude-frequency response. We found that the compensation filter response is consistent with the actual amplitude response with a good simulation precision. This method has been well applied to the interaction display area of earthquake science museum of Guangdong province.

LabVIEW; Graphical programming language; Seismograph; Real-time data; Decoding; Simulation

吳華燈，謝劍波，2014．EDAS-C24型數(shù)字測震儀實時數(shù)據(jù)的解碼及加速度和位移仿真的實現(xiàn)．震災(zāi)防御技術(shù)，9（3）：508—517．

10.11899/zzfy20140318

2013-11-21

吳華燈，男，生于1980年。理學(xué)學(xué)士，工程師。主要從事地震觀測研究和軟硬件開發(fā)工作。E-mail：gdea_whd @aliyun.com