吳華燈1, 2, 3) 謝劍波1, 2, 3)
1)廣東省地震局,廣州 510070
2)中國地震局地震監測與減災技術重點實驗室,廣州 510070
3)廣東省地震預警與重大工程安全診斷重點實驗室,廣州 510070
EDAS-C24型數字測震儀實時數據的解碼及加速度和位移仿真的實現
吳華燈1, 2, 3) 謝劍波1, 2, 3)
1)廣東省地震局,廣州 510070
2)中國地震局地震監測與減災技術重點實驗室,廣州 510070
3)廣東省地震預警與重大工程安全診斷重點實驗室,廣州 510070
本文嘗試使用了一種圖形編程語言進行數據解碼和仿真研究,闡述了用該語言實現EDAS-C24型數字測震儀實時數據流的解碼過程,提出了在G語言下利用數字濾波器逼近模擬積分器及模擬微分器響應實現對解碼數據實時仿真的方法,并通過設定頻帶寬度,比較了設計的補償濾波響應和實際的幅頻響應的一致性。結果表明,在設定的頻帶內,補償濾波響應和實際的幅頻響應是一致的,仿真的精度是理想的,成果已經在廣東省地震科普館的地震互動區部署運行,取得了一定的實效。
LabVIEW G編程語言 測震儀 實時數據 解碼 仿真
Kanamori等(1999)利用美國TriNet強震臺網,曾對實時速度記錄通過單邊差分方法計算位移等問題做過研究,但是其計算精度和相位仍存在問題;金星等(2003)利用單自由度體系地震反應的時域遞歸方法,系統研究了由數字加速度記錄實時仿真地動速度和地動位移的問題,從理論和數值計算上取得了理想的結果;同時金星等(2004)還系統研究了利用數字化寬頻帶速度型記錄仿真地面位移和地面加速度等問題,從時域上提出了一套實時計算公式;謝劍波(2014)也用時域遞歸方法做過不同帶寬的地震記錄的仿真研究。以上的實現都是利用地震臺網匯集的實時數據或事件數據,在傳統文本編程語言的環境下實現數據的仿真。本文將基于LabVIEW圖形編程語言——G編程語言,在數據采集現場,直接解碼數據采集器的實時速度數據,并設計合適的數字濾波器,實現加速度和位移的實時仿真,最終繪制出速度、加速度和位移三種物理量的實時波形。
1.1 數據采集現場設備的架設
地震動數據采集現場如圖1所示,由EDAS-C24型地震數據采集器、FBS-3A型反饋式寬頻帶地震計、GPS授時系統、12V電源供電系統和集成開發電腦平臺等設備架設而成。FBS-3A型反饋式寬頻帶地震計是機電一體化的測震儀器,它在自振周期為2s的機械拾震器上附加電子線路,構成測量頻帶約為0.05—40Hz的寬頻帶測量系統。在儀器結構上,它由一個豎直向和兩個相互正交的水平向組成(劉慶偉等,2001)。該儀器用于拾取地面的震動信號,并把信號轉換成電壓值,供數據采集器數字化。EDAS-C24型地震數據采集器將FBS-3A型地震計輸出的模擬電壓信號轉換成數字信號,并經串口線把數字信號傳輸到電腦,供計算機程序處理。

圖1 數據采集現場設備的架設Fig.1 The set-up of equipment
1.2 數據解碼及仿真環境
在圖1所示的集成開發平臺電腦上,安裝了LabVIEW2010集成開發軟件。LabVIEW是Laboratory Virtual Instrument Engineering Workbench(實驗室虛擬儀器集成環境)的簡稱,是一個功能強大而又靈活的圖形開發環境(喬瑞萍等,2008)。由于EDAS-C24型地震數據采集器是基于串口通信,本文用LabVIEW中的VISA(VISA是Virtual Instruments Software Architecture的縮寫,中文譯為虛擬儀器軟件架構或可視化儀器軟件架構)編程技術,建立與數據采集器串口的通信,實現實時數據的接收,為數據解碼和仿真提供數據源。
2.1 EDAS-C24型地震數據采集器實時數據格式
EDAS-C24型地震數據采集器實時數據有壓縮和非壓縮兩種數據格式,本文主要針對非壓縮數據格式進行解碼。該型號數據采集器實時數據幀是由幀頭、數據區、校驗和組成。幀頭長16個字節,包括幀同步標識碼、臺站編號、幀類型標示字、幀長、秒記數、標志字節、頭段區字節檢查和;數據區包括3個通道的數據,當采樣率是125Hz時,還有一個字節的填充字;幀的尾部是兩個字節的檢查和。數據幀格式見表1,數據幀的具體內容見圖2。從表1中可看出,在EDAS-C24中,字節序是小端,每一個數據幀都是以“74 97 13 BF”同步標識碼開頭,BF是高字節位。圖2中的數據幀內容,是在計算機程序緩沖區中的字節序。因此,程序解碼時,將用小端的字節序來判斷數據幀和提取3通道的數據。

表1 EDAS-C24數據幀格式Table 1 Data format of EDAS-C24

續表

圖2 EDAS-C24數據幀內容Fig.2 Content of EDAS-C24 data frame
2.2 實時數據解碼
解碼是在LabVIEW2010專業版開發系統下,采用狀態機、事件隊列消息處理器和事件結構結合的架構進行編程,用到了While循環、條件結構和事件結構,架構代碼如圖3所示。

圖3 程序架構Fig.3 Program frame
狀態機里定義了Init、No Action、Receive Data和Exit四種事件狀態及對應的事件分支,優先執行Init事件分支,用于初始化串口名稱、串口波特率、采樣率、滿量程、事件觸發模式、權重等配置參數和讀取串口的數據,為數據的解碼提供數據源。通過獲取下一個事件的子VI(圖形編程語言的子程序),依次改變事件狀態為No Action、Receive Data和Exit。在No Action事件分支里,用了事件結構。事件結構的超時端子為50ms,有“超時”、“前面板關閉”、“菜單選擇”、“運行狀態”4個分支,主要用到了“超時”事件分支。當程序運行超過50ms,將執行“超時”事件框中的程序,如果收到串口的數據并且程序處于運行狀態,則驅動Receive Data事件,否則置錯誤信息。Receive Data事件分支的程序是解碼和仿真EDAS-C24型地震數據采集器實時數據的核心。在解碼方面,程序首先由圖標為GET FRAME、名稱為Receive Frame的子VI,從VISA資源里找出數據幀頭。然后,通過Parse Frame Header子VI分析數據幀頭的有效性。如果數據幀頭是有效的,則進一步求解出三個通道的原始速度數據,解碼的具體實現過程如下所述。
2.2.1 數據幀頭解碼的實現
通過隧道和移位寄存器把Init事件分支初始化時的VISA資源輸入到GET FRAME子VI的VISA resource name入口(見圖4),然后利用“VISA讀取”函數,讀取16個字節大小的數據幀頭到header buffer緩沖區,接著將header buffer緩沖區的數據輸入圖5所示的Parse Frame Header子VI,實現數據幀頭Frame Header的輸出及其有效性的判斷。

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

圖5 Parse Frame Header子VIFig.5 Parse Frame Header sub VI
“從字符串還原”函數有類型、二進制字符串、字節序等參量的輸入。類型就是按照“表1 EDAS-C24數據幀格式”以簇方式定義的幀頭;二進制字符串就是header buffer;由于在EDAS-C24中,字節序是小端,所以字節序選擇little endian。該函數的輸出結果是數據幀頭Frame Header。數據幀頭Frame Header有效性的判斷用了3個判斷條件:其一,從數據幀頭Frame Header的簇中求出幀長,判斷幀長是否小于1000;其二,取Frame Header前四個字節,判斷是否與“BF 13 97 74”相等;其三,取Frame Header第七、第八個字節,判斷是否與“AA55”相等,當3個判據同時滿足時,認為Frame Header是有效的,至此實現了數據幀頭解碼。
2.2.2 數據區解碼的實現
在解碼出數據幀頭和確認幀頭有效后,程序進入Case結構。依據EDAS-C24數據幀格式,數據區的字節數byte count等于數據幀長度減8個字節,從VISA resource out讀取byte count個字節的數據到讀寫緩沖區read buffer。當采樣率為100bps時,每個通道每個采樣點有3個字節,3個通道每個采樣點一共9個字節。可以用一個循環,循環的次數為byte count字節/9字節,在循環體里結合“索引數組”函數、24bit子VI和“創建數組”函數,分解出三個通道的數據,如圖6所示。“索引數組”函數的作用是建立每個通道LMH的索引;24bit子VI的作用是將每個通道的LMH重新組合成24位的數據,并轉換成計算機的32位數據,再右移8位,得到每個通道的最終數據;“創建數組”函數的作用是把每個通道的最終數據添加到數組,為繪制波形提供數據源。

圖6 三個通道數據的解碼Fig.6 Decoding of three channels data
3.1 實時數據的仿真方法
實時數據的仿真,一般通過計算單自由度系統的地震反應即可得到(胡聿賢,2006)。單自由度系統地震反應的計算方法可分為二大類:一類是在頻域內進行;另一類是在時域下進行。由于頻域分析方法必須在獲取完整的地震記錄后才能進行,不適合于實時計算的要求,人們一般采用時域方法進行仿真計算。即用時間上有限個離散的動力平衡方程來替換連續形式的動力平衡方程,進而在一定假定下依據離散動力平衡方程來尋找數字濾波器,使得所尋找的數字濾波器的傳遞函數逼近理論傳遞函數。同時,利用數字濾波尤其是遞歸濾波的特點實現實時計算。
3.1.1 常用的時域計算方法
有關單自由度系統地震反應的時域計算方法有很多,如中心差分方法、New-mark方法、變換方法、Duhamel逐步積分法和wilson-0法等。許多學者開展了相關的研究工作,并取得了很好的成果。如金星等(2004)的分析與研究,得到了計算單自由度系統相對位移、相對速度、相對加速度的遞歸公式,并求得由地面速度時程計算位移時程的公式:

式中,b1,b2是遞歸系數;xk-3/2和xk-5/2是前兩個輸出點位移;S0是與傳遞函數的零頻約束條件有關的參數;Δt是離散時間間隔;δ 是與Δt/T0有關的參數;T0是自振周期;Vk是k點的速度。
以及求得由地面速度時程計算相對加速度時程的公式:

利用式(1)和式(2)進行了實時仿真位移與加速度時程的計算研究。
3.1.2 本文作者提出及使用的仿真方法


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

輸出的Sinc信號一組取輸出值,另一組取輸出值的倒數。斜坡信號與兩者捆綁成簇數組后,按照設定的頻帶0.4Hz,提取一定數量的數據,重新組合成簇數組,再由Remez的交換方法建立equi-ripple等紋波濾波器,并輸出實際紋波值。等文波濾波器的濾波輸出經過DFD plot Freq Response子VI,輸出并繪制補償濾波響應,再將filter out輸入到獲取傳遞函數的子VI,取出系統的傳遞函數。接著,由給定的頻率和截止頻率,計算增益、零點、極點參數,通過DFD build filter from zeros-poles-gain子VI設計出新的濾波器。最后,級聯兩個濾波器,再經DFD plot Freq Response子VI得到最終的濾波輸出,實現逼近模擬積分器響應的設計。
微分濾波器的設計,也是首先生成斜坡信號和Sinc信號,再由Remez數字濾波器設計方法設計任意頻響的數字濾波器。任意頻響的數字濾波器經DFD Get Transfer Function子VI后,取出傳遞函數。利用傳遞函數,通過DFD Build Filter From Transfer Function子VI,建立新的濾波器。最后,級聯兩個濾波器,再經DFD plot Freq Response子VI得到最終的幅頻響應和濾波輸出,實現逼近模擬微分器響應的設計。
3.2 實時數據仿真的實現
如圖7所示,微分濾波器和積分濾波器的濾波輸出分別接入DFD Filter子VI,速度數據便可實時仿真成加速度和位移數據,同時將加速度、速度和位移數據在波形圖上顯示出來,見圖8。

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

圖8 加速度、速度和位移數據波形圖Fig.8 Waveform of acceleration, velocity and displacement
在進行濾波器設計過程中,為了觀察仿真的精度,設定了頻帶寬度,將實際的幅頻響應和設計的補償濾波響應繪制了出來(見圖9)。此外,還輸出了實際紋波值和繪制了紋波波形(見圖10),輸出了濾波系數,繪制了最終的幅頻響應。從圖中可看出,補償濾波響應在設定的頻帶內和實際的幅頻響應是一致的,設計的紋波值在0.005左右,從設計的原型上看,仿真的精度是理想的。

圖9 設計的補償濾波響應Fig.9 Designed compensation filter response

圖10 設計的紋波波形輸出Fig.10 Output of designed ripple waveform
為了驗證仿真數據的正確性及精度,在同一臺址架設了珠海泰德企業有限公司生產的強震動數據采集器TDE-324CA和力平衡加速度計TDA-33M,進行仿真加速度記錄與實測加速度記錄對比觀測實驗,對比觀測所使用的儀器信息如表2所示。

從傅立葉幅值譜(圖12)和自相關系數(圖13)的結果可以看出,仿真的加速度數據在頻帶和相位上與實測的加速度數據吻合較好,說明本文設計的濾波器是符合預期的,具有一定的可靠性。
本文在解碼出數據采集器的速度數據后,利用數字濾波器逼近模擬積分器和模擬微分器響應的方法,實現了由速度數據實時仿真位移和加速度數據,數據的仿真精度是理想的。但本文設計的濾波器更多的是考慮加速度數據所關心的頻帶,仿真系統在有效頻帶方面尚應做更多的分析研究。此外,在沒有實測位移數據的情況下,還應考慮利用現有的一些典型地震事件數據進行相關性試算驗證,計算不同阻尼比、不同自振周期等條件的結果,逐步修正和優化濾波器,使其具有實用性,方便應用于地震觀測數據的實時仿真。本文的成果已經在廣東省地震科普館的地震互動區部署運行,取得了一定的實效。
胡聿賢,2006.地震工程學(第二版).北京:地震出版社.
金星,馬強,李山有,2004.利用數字化速度記錄實時仿真位移與加速度時程.地震工程與工程振動,24(6):9—14.
金星,馬強,李山有,2003.四種計算地震反應數值方法的比較研究.地震工程與工程振動,23(1):18—30.
劉慶偉,莊燦濤,劉慧寧,2001.FBS-3A型反饋式寬帶地震計的傳遞函數.地震學報,23(1):79—86
喬瑞萍,林欣等譯,2008.LabVIEW 8實用教程.北京:電子工業出版社.
萬永革,2007.數字信號處理的MATLAB實現.北京:科學出版社.
謝劍波,2014.地震記錄的時間域反褶積仿真及在地震計方位角相對測量中的應用.地球物理學報,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型數字測震儀實時數據的解碼及加速度和位移仿真的實現.震災防御技術,9(3):508—517.
10.11899/zzfy20140318
2013-11-21
吳華燈,男,生于1980年。理學學士,工程師。主要從事地震觀測研究和軟硬件開發工作。E-mail:gdea_whd @aliyun.com