GNSS／INS組合用于地震監(jiān)測的算法實(shí)現(xiàn)和仿真研究

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

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

1 引言

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

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

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

2.1 地震仿真模型

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

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

其中

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

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

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

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

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

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

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

2.2 儀器相關(guān)誤差參數(shù)設(shè)置

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

表1 GNSS噪聲參數(shù)

表2 IMU噪聲參數(shù)

2.3 仿真信號(hào)

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

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

圖1 仿真IMU加速度

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

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

3.1 卡爾曼濾波

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

狀態(tài)方程為

量測方程為

圖3 仿真GNSS位置誤差

圖4 地震運(yùn)動(dòng)仿真的真實(shí)位移

卡爾曼濾波方程為

圖5 地震運(yùn)動(dòng)仿真的真實(shí)位移（局部）

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

3.2 算法模型建立

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

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

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

3.2.2 卡爾曼濾波量測方程

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

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

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

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

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

式中，k表示時(shí)間，N是最后一個(gè)時(shí)刻，下標(biāo)s表示平滑的相關(guān)變量。在反向平滑中k＝N－1，N－2，…0，作為平滑解算的初始值設(shè)為正向?yàn)V波的結(jié)尾值

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

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

圖6 濾波后位置誤差

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

圖7 東向位置誤差平滑前后對(duì)比

圖8 北向位置誤差平滑前后對(duì)比

圖9 天向位置誤差平滑前后對(duì)比

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

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

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

4 總結(jié)

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

5 致謝

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

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

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

［3］ 黃世敏，王亞勇，丁潔民，等.GB50011－2010，建筑抗震設(shè)計(jì)規(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］ 秦永元.慣性導(dǎo)航［M］.北京：科學(xué)出版社，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)星／慣性組合導(dǎo)航事后高精度融合算法研究［J］.系統(tǒng)仿真學(xué)報(bào)，2010（增刊1）：75－78.