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

（武漢大學 衛(wèi)星導航定位技術(shù)研究中心，武漢 430079）

1 引言

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

2 數(shù)據(jù)仿真

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

2.1 地震仿真模型

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

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

其中

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

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

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

f（t）可根據(jù)下式確定

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

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

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

2.2 儀器相關誤差參數(shù)設置

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

表1 GNSS噪聲參數(shù)

表2 IMU噪聲參數(shù)

2.3 仿真信號

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

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

圖1 仿真IMU加速度

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

3 數(shù)據(jù)處理

3.1 卡爾曼濾波

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

狀態(tài)方程為

量測方程為

圖3 仿真GNSS位置誤差

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

卡爾曼濾波方程為

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

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

3.2 算法模型建立

3.2.1 卡爾曼濾波系統(tǒng)方程

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

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

3.2.2 卡爾曼濾波量測方程

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

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

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

3.3 固定區(qū)間最優(yōu)平滑算法

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

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

3.4 數(shù)據(jù)處理及結(jié)果分析

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

圖6 濾波后位置誤差

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

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

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

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

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

表3 位置誤差對比統(tǒng)計／m

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

4 總結(jié)

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

5 致謝

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

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

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

［3］ 黃世敏，王亞勇，丁潔民，等.GB50011－2010，建筑抗震設計規(guī)范［S］.北京：中國建筑工業(yè)出版社，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］ 李睿佳，李榮冰，劉建業(yè)，等.衛(wèi)星／慣性組合導航事后高精度融合算法研究［J］.系統(tǒng)仿真學報，2010（增刊1）：75－78.