GNSS／INS組合用于地震監測的算法實現和仿真研究

閔 翔，牛小驥，李 由，施 闖

（武漢大學 衛星導航定位技術研究中心，武漢 430079）

1 引言

我國是一個自然災害頻發的國家，近年來地震、滑坡、尾礦庫坍塌、地基不均勻沉陷等事件時有發生，造成大量人員傷亡和巨額經濟損失，迫切需要有效的連續位移監測手段對災害進行預警和記錄；同時國內許多新建的大型建筑結構（如跨海大橋、摩天大樓、大壩等）也迫切需要有效的形變監測手段來監控其健康狀況并對異常狀況提供預警。對比其它連續位移監測手段，高精度全球導航衛星系統（GNSS）接收機具有精度穩定、布設方便、相對成本低等獨特優勢；慣性導航系統（INS）具有短期精度高，不受障礙物影響等優勢，兩者具有很強的互補性。由于地震是一種特殊的運動形式，引發地震需要兩個基本元素：一是地殼的運動；二是局部區域的地殼具有較強的剛性，阻礙地殼的運動，使應變能得以積累。而地殼運動一般可以簡單地分為水平運動和垂直運動［1］。所以，造成地震的地殼運動也可以分為水平運動和垂直運動兩種。在這些運動中不存在角運動的情況，因此慣導的慣性測量單元（IMU）可以省掉用于角運動測量的陀螺，只保留加速度計即可。相應的導航算法和組合算法需要做特殊的簡化設計，數據仿真和誤差建模也要有針對性。另外，GNSS用于地震監測時一般采用精密單點定位（PPP）工作模式，以避免精密定位對主站或周邊參考站的依賴。本文首先通過數據仿真給出有代表性的地震運動數據以及GNSS和加速度計的觀測數據，再通過松組合方法進行解算，得出對地震運動的最優估計結果。

2 數據仿真

根據地震的運動特點以及設備的觀測值誤差參數，進行了如下模型設計和數據仿真。

2.1 地震仿真模型

對地震運動的數據仿真根據特征周期、水平地震波影響系數最大值和地震波幅值等初始條件，采用三角級數法［2］生成人工地震波。

對于給定的功率譜密度函數Sx（ω），按照下面的公式可以方便地生成以Sx（ω）為功率譜密度函數、均值為零的高斯平穩過程a（t）。

其中

式中，φk為（0，2π）內均勻分布的隨機相角，ωu，ωl分別為正ω域內的上、下限值，即認為Sx（ω）的有效功率在（ωu，ωl）范圍內，而范圍外的Sx（ω）值可視為零。

為了反映地面運動的非平穩性，采用包絡函數f（t）乘以平穩過程a（t）得

式（3）即為人工地震波模型。

f（t）可根據下式確定

式中，c為衰減系數、通常取值范圍為0.1～1.0（仿真中取0.15），t1、t2和t3根據不同實際情況取值，T為地震波持時。仿真中給出地震時間為2min，三軸上對應的t1、t2和t3分別取不同的時間。

仿真采用文獻 ［3］中的反應譜作為目標譜，通過Kaul提出的平穩過程反應譜與功率譜的近似關系

式中，SaT（ωk）為規范反應譜，ξ為阻尼比，Td為地震動持時，p為反應不超過反應譜值的概率、仿真中取0.85。通過式（3）和式（5）即可生成仿真地震波。

2.2 儀器相關誤差參數設置

IMU的加速度計輸出a是在仿真給出的真實加速度基礎上疊加白噪聲和零位漂移（建模為一階高斯馬爾科夫過程）誤差得出的，其相關參數按照中等精度IMU（戰術級）的指標［4］給出。GNSS定位結果（PPP）在動態情況下的測量誤差是用ALLAN方差分析真實數據后通過模型辨識和參數估計給出的誤差參數［5］。針對定位模式，主要有兩種誤差，即位置噪聲和系統性的位置漂移，前者用白噪聲建模，后者用一階高斯馬爾科夫過程建模。IMU和GNSS相應的誤差參數如表1及表2所示。

表1 GNSS噪聲參數

表2 IMU噪聲參數

2.3 仿真信號

為了讓仿真中加入的一階馬爾科夫過程誤差和組合解算有充分的收斂時間能夠穩定下來，仿真中也包含了地震前15min的仿真數據。結果包括IMU的加速度和GNSS的位置。相關曲線見圖1、圖2及圖3。

仿真地震運動的真實位置和真實速度是通過無誤差的IMU加速度積分得出的，見圖4及圖5。真實位置主要是用來分析組合解算結果的位置誤差參數。

圖1 仿真IMU加速度

圖2 仿真IMU加速度（局部）

3 數據處理

3.1 卡爾曼濾波

卡爾曼濾波［6］是一種最優遞推估計算法，在組合導航以及其他工程項目中有廣泛的應用。其離散化形式的方程如下。

狀態方程為

量測方程為

圖3 仿真GNSS位置誤差

圖4 地震運動仿真的真實位移

卡爾曼濾波方程為

圖5 地震運動仿真的真實位移（局部）

式（6）－式（8）中：下標k表示時標，Φk，k－1為狀態轉移系數矩陣，Gk－1為系統噪聲驅動系數矩陣，Hk為量測系數矩陣，Kk是增益矩陣，是狀態預測向量，是狀態估計向量，Pk／k－1是預測均方誤差陣，Pk是狀態均方誤差陣。在卡爾曼濾波方程中，更新是在預測的基礎之上進行的，因此增益Kk扮演了一個連接此次觀測數據與上一次預測的權重分配的角色，通過不斷調整權重分配來得出最優結果。

3.2 算法模型建立

3.2.1 卡爾曼濾波系統方程

系統狀態方程反映了模型的狀態演化更新，狀態向量由位置誤差、速度誤差、和加速度計零偏組成，一共有9維［7］。方程的建立如下

式中，Δp、Δv、Δb分別表示位置、速度和零偏的誤差量，為對應量的導數，E、N、U分別代表導航坐標系的三個軸（東北天），n、q分別代表加速度計的白噪聲和加速度計零偏模型（一階高斯馬爾科夫過程）的驅動白噪聲，T表示零偏模型的相關時間。此方程是連續時間的形式，在后面的實際計算中需要離散化。加速度計零偏估計主要是用來在線補償加速度計誤差的。

3.2.2 卡爾曼濾波量測方程

量測方程是用INS和GNSS的位置結果的差異作為量測向量，通過位置差異信息來對速度和位置的預測量進行改正。量測方程如式（10）。

其中，p表示位置，下標E、N、U表示導航坐標系軸向，上標i、g分別代表INS和GNSS，n是白噪聲。而在仿真過程中INS的噪聲包括白噪聲和一階馬爾科夫過程，以符合實際情況。而在上式的觀測方程中將一階馬爾科夫過程也當作白噪聲來處理，以簡化模型以適應Kalman濾波算法。

在量測方程中，本質上是通過引入位置量測信息來對各狀態量進行直接或間接的修正。

3.3 固定區間最優平滑算法

平滑算法主要是針對正向濾波后的數據處理，處理后的數據精度在正向濾波的基礎上有一定的提升。本文采用固定區間最優平滑算法（RTS）［8］。在正向濾波中需要保存濾波中的相關變量，如狀態估計量，狀態預測量，狀態估計量的協方差陣Pk，狀態預測量的協方差陣Pk／k－1，狀態轉移系數陣φk，k－1等。其算法公式如下

式中，k表示時間，N是最后一個時刻，下標s表示平滑的相關變量。在反向平滑中k＝N－1，N－2，…0，作為平滑解算的初始值設為正向濾波的結尾值

3.4 數據處理及結果分析

基于上述特殊的GNSS／INS組合算法，對仿真的數據處理之后發現前100s內存在誤差波幅比較大的現象，這主要是因為在仿真中速度和位置的初始值不準確的影響，本文直接給出初始速度和位置都為0，從而致使INS加速度計積分出來的位置偏差比較大，經過一段時間的濾波更新后誤差逐漸變小。由于仿真是由15min的靜止和2min的地震運動組成，因此在分析誤差時，可以去除開始一段收斂過程的結果，而重點考慮算法收斂后的數據。這里將考察200s收斂過程以后的數據，避免初始收斂過程引入的較大誤差。處理后的相關結果圖見圖6－圖9。

圖6 濾波后位置誤差

圖6－圖9給出了濾波后的位置誤差與平滑后的位置誤差對比圖，可以看出平滑后的誤差要明顯比濾波后的小。

圖7 東向位置誤差平滑前后對比

圖8 北向位置誤差平滑前后對比

圖9 天向位置誤差平滑前后對比

由于單次仿真的時間長度偏短，不足以統計其誤差水平。因此本文在保持各項參數不變的情況下，進行了多次仿真和解算從而得到一個位置誤差的統計值，表3是經過10次仿真和解算后的統計結果。

表3 位置誤差對比統計／m

通過這10次的仿真和解算結果，可以清楚的看出濾波之后的三軸方向上的精度要比GNSS的精度分別高出12.9%、15.7%、10.5%，而平滑后的精度則更高，相對于GNSS的三軸精度要高出19.4%、31.5%、18.7%，達到 了水平位移0.6cm和高程位移1.5cm，能夠滿足地震監測的精度要求。從平滑結果相對于濾波結果精度的改善效果來看，并不是很理想。其主要原因是，在仿真過程中給出的GNSS位置噪聲主要是一階馬爾科夫過程的漂移誤差，而不是白噪聲誤差。

4 總結

本文嘗試了采用由加速度計構成的慣導與GNSS精密單點定位（PPP）的組合算法用于地震監測。在數據仿真的基礎之上，驗證和分析了所提出的算法的正確性、精度和可行性，結果表明基于中等精度（戰術級）加速度計和GPS定位模式的GNSS／INS組合系統能夠達到水平位移0.6cm和高程位移1.5cm的精度水平（中值），滿足地震監測精度要求，具有明確的推廣應用潛力。

5 致謝

本文作者衷心感謝武漢大學衛星導航定位技術研究中心的博士生張全和碩士生陳起金提供了組合導航算法和GNSS誤差模型方面的支持。

［1］ 廖永巖.地球科學原理［M］.北京：海洋出版社，2007.

［2］ 胡聿賢.地震工程學［M］.2版.北京：地震出版社，2006.

［3］ 黃世敏，王亞勇，丁潔民，等.GB50011－2010，建筑抗震設計規范［S］.北京：中國建筑工業出版社，2010.

［4］ NOVATEL INC.User＇s Manual of FSAS［EB／OL］.［2013－01－15］.http：／／www.novatel.com／assets／Documents／Papers／FSAS.pdf.

［5］ NIU Xiao－ji，CHEN Qi－jin.Using Allan Variance to Analyze the Error Characteristics of GNSS Positioning［EB／OL］.［2012－12－26］.http：／／www.ion.org／meetings／gnss2012／gnss2012abstracts.cfm.

［6］ 秦永元.慣性導航［M］.北京：科學出版社，2006.

［7］ BOCK Y，MELGAR D，CROWELL B W.Real－time Strong－motion Broadband Displacements from Collocated GPS and Accelerometers［J］.Bulletin of the Seismological Society of America，2011，101（6）：2904－2925.

［8］ 李睿佳，李榮冰，劉建業，等.衛星／慣性組合導航事后高精度融合算法研究［J］.系統仿真學報，2010（增刊1）：75－78.