劉瓊瓊 張洪勝

（周口市規(guī)劃建筑勘測(cè)設(shè)計(jì)院，河南 周口 466000）

20世紀(jì)60年代，美籍匈牙利數(shù)學(xué)家卡爾曼首先提出了用一個(gè)狀態(tài)方程和一個(gè)量測(cè)方程來(lái)完整描述線性動(dòng)態(tài)方程，即卡爾曼濾波方程，采用時(shí)域法及狀態(tài)方程為其數(shù)學(xué)工具，以多變量控制、最優(yōu)控制與估計(jì)以及自適應(yīng)為主要內(nèi)容。這是一種動(dòng)態(tài)進(jìn)行實(shí)時(shí)數(shù)據(jù)處理的有效方法，無(wú)需存儲(chǔ)不同時(shí)刻的觀測(cè)數(shù)據(jù)。在實(shí)際應(yīng)用中，如何選擇濾波模型，選擇濾波初值、精度評(píng)定，本文將進(jìn)行介紹，并結(jié)合應(yīng)用實(shí)例說(shuō)明。

1 卡爾曼濾波基本原理

卡爾曼濾波屬于軟件濾波方法，是以最小均方誤差為最佳估計(jì)準(zhǔn)則，采用信號(hào)與噪聲的狀態(tài)空間模型，利用前一時(shí)刻的估計(jì)值和當(dāng)前時(shí)刻的觀測(cè)值來(lái)更新對(duì)狀態(tài)變量的估計(jì)，求出當(dāng)前時(shí)刻的估計(jì)值，根據(jù)建立的系統(tǒng)方程和觀測(cè)方程對(duì)需要處理的信號(hào)做出滿(mǎn)足最小均方誤差的估計(jì)[1]。最初的卡爾曼濾波算法適用于解決隨機(jī)線性離散系統(tǒng)的狀態(tài)或參數(shù)估計(jì)問(wèn)題。

在實(shí)際應(yīng)用中，可以將物理系統(tǒng)的運(yùn)行過(guò)程看成一個(gè)狀態(tài)轉(zhuǎn)換過(guò)程，卡爾曼濾波將狀態(tài)空間理論引入物理系統(tǒng)的數(shù)學(xué)建模中，其假設(shè)系統(tǒng)狀態(tài)可以用n維空間的一個(gè)向量X∈Rn來(lái)表示。為了描述方便，可以作以下假設(shè)：物理系統(tǒng)的狀態(tài)轉(zhuǎn)換過(guò)程，描述為一個(gè)離散時(shí)間的隨機(jī)過(guò)程。系統(tǒng)狀態(tài)受控制輸入影響；系統(tǒng)狀態(tài)及觀測(cè)過(guò)程受噪聲影響；對(duì)系統(tǒng)狀態(tài)是非直接可觀測(cè)的[2]。在假設(shè)前提下，定義系統(tǒng)狀態(tài)變量為X∈Rn，系統(tǒng)控制輸入為UK，系統(tǒng)過(guò)程激勵(lì)噪聲為WK，得出系統(tǒng)的狀態(tài)隨機(jī)差分方程：

定義觀測(cè)變量Zk∈Rm，觀測(cè)噪聲為Vk，得到量測(cè)方程：：

假設(shè)Wk、Vk為相互獨(dú)立，正態(tài)分布的白色噪聲，過(guò)程激勵(lì)噪聲協(xié)方差矩陣為Q，觀測(cè)噪聲協(xié)方差矩陣為R。A、B、H統(tǒng)稱(chēng)為狀態(tài)變換矩陣，是狀態(tài)變換過(guò)程中的調(diào)整系數(shù)，是從建立的系統(tǒng)數(shù)學(xué)模型中導(dǎo)出來(lái)的，假設(shè)它們是常數(shù)，得出卡爾曼濾波迭代的五個(gè)公式如下：

時(shí)間更新方程：

狀態(tài)更新方程：

式中， A是作用在Xk-1上面的n×n狀態(tài)變換矩陣；B是作用在控制向量Uk-1上的n×1輸入控制矩陣；H是m×n觀測(cè)模型矩陣，把真實(shí)狀態(tài)空間映射成觀測(cè)空間；為n×n先驗(yàn)估計(jì)誤差協(xié)方差矩陣；Pk為n×n后驗(yàn)估計(jì)誤差協(xié)方差矩陣；Q是n×n過(guò)程噪聲協(xié)方差矩陣；R是m×m過(guò)程噪聲協(xié)方差矩陣；I是n×n階單位矩陣；Kk是n×m階矩陣，稱(chēng)為卡爾曼增益或混合因數(shù)，可使后驗(yàn)估計(jì)誤差協(xié)方差最小。

上述給出的是離散卡爾曼濾波的時(shí)間更新方程和狀態(tài)更新方程（-代表先驗(yàn)，代表估計(jì)）?？柭鼮V波主要分兩步進(jìn)行，先預(yù)測(cè)后修正，對(duì)應(yīng)兩個(gè)方程，預(yù)測(cè)方程（1）和修正方程（2）。預(yù)測(cè)方程功能：由系統(tǒng)上一狀態(tài)最優(yōu)估計(jì)狀態(tài)推算出系統(tǒng)當(dāng)前狀態(tài)，即當(dāng)前狀態(tài)預(yù)測(cè)值＝系統(tǒng)狀態(tài)轉(zhuǎn)移系數(shù)（A）×系統(tǒng)上一狀態(tài)的最優(yōu)估計(jì)值＋狀態(tài)估計(jì)噪聲。修正方程功能：由系統(tǒng)當(dāng)前狀態(tài)的預(yù)測(cè)狀態(tài)和實(shí)際觀測(cè)值求取當(dāng)前狀態(tài)的最優(yōu)估計(jì)狀態(tài)，即系統(tǒng)當(dāng)前狀態(tài)的最優(yōu)值＝當(dāng)前狀態(tài)預(yù)測(cè)值＋修正系數(shù)（K）×（實(shí)際觀測(cè)量－當(dāng)前狀態(tài)預(yù)測(cè)值）。這個(gè)過(guò)程是預(yù)估－校正過(guò)程，對(duì)應(yīng)的這種估計(jì)稱(chēng)為預(yù)估－校正算法。

2 濾波模型的選擇

卡爾曼濾波模型誤差的影響是卡爾曼濾波發(fā)散的重要原因之一。為了避免模型誤差造成的不良影響，測(cè)量工作者要嚴(yán)格按照規(guī)定進(jìn)行觀測(cè)，保證觀測(cè)過(guò)程的合理性和準(zhǔn)確性，盡量給出接近實(shí)際的模型[3]。

在變形測(cè)量實(shí)時(shí)監(jiān)測(cè)項(xiàng)目中，常用的卡爾曼濾波模型為位置、速度模型和位置、速度、加速度模型，這些是將監(jiān)測(cè)點(diǎn)的變形位移看成一個(gè)隨機(jī)過(guò)程。位置速度濾波模型認(rèn)為在變形測(cè)量中，監(jiān)測(cè)點(diǎn)位移速度的均值不變，在濾波中將監(jiān)測(cè)點(diǎn)的位置、移速度作為狀態(tài)量，將位移加速度視作動(dòng)態(tài)噪聲，其狀態(tài)方程為：

位置、速度、加速度濾波模型認(rèn)為在變形測(cè)量中，監(jiān)測(cè)點(diǎn)位移加速度的均值不變，濾波中將監(jiān)測(cè)點(diǎn)的位置、移速度和加速度都作為狀態(tài)量，將加速度變化率視作動(dòng)態(tài)噪聲，，其狀態(tài)方程為：

式中 xk、、分別為監(jiān)測(cè)點(diǎn)的坐標(biāo)、位移速度和加速度，Ωk-1為動(dòng)態(tài)噪聲，t表示時(shí)間。

觀測(cè)噪聲通過(guò)模擬觀測(cè)實(shí)驗(yàn)獲得。取若干期監(jiān)測(cè)數(shù)據(jù)，計(jì)算監(jiān)測(cè)點(diǎn)各期坐標(biāo)及對(duì)應(yīng)的速度、加速度，對(duì)速度和加速度作圖分析。當(dāng)各期速度在其平均速度的兩側(cè)波動(dòng)變化，無(wú)明顯增、減速跡象時(shí)，采用位置、速度模型；當(dāng)各期速度有增、減速跡象，而其加速度無(wú)明顯增、減速跡象時(shí)，采用位置、速度、加速度模型，最后根據(jù)各期加速度的變化率估算狀態(tài)方程的動(dòng)態(tài)噪聲的方差。

3 模型精度的評(píng)定

當(dāng)精確已知觀測(cè)噪聲和動(dòng)態(tài)噪聲的隨機(jī)模型時(shí)，不需要專(zhuān)門(mén)進(jìn)行精度評(píng)定，但實(shí)際應(yīng)用中，所采用的先驗(yàn)系統(tǒng)狀態(tài)誤差協(xié)方差陣Q和量測(cè)誤差協(xié)方差陣R與其實(shí)際值可能不完全一致，這時(shí)Q和R的含義與其說(shuō)是協(xié)方差矩陣不如說(shuō)是權(quán)逆矩陣更適當(dāng)，式（7）中的Pk就是系統(tǒng)狀態(tài)量的權(quán)逆矩陣，因此要計(jì)算單位權(quán)方差的估計(jì)值，以獲得狀態(tài)參數(shù)估計(jì)的協(xié)方差矩陣[4]。

對(duì)于第k期濾波，單位權(quán)方差σ2，可按照下式估計(jì)：

式中，nk為第k期濾波觀測(cè)值的個(gè)數(shù)，為ξk預(yù)測(cè)殘差向量的權(quán)逆矩陣，這時(shí)狀態(tài)參數(shù)向量濾波值得協(xié)方差：

在施工變形測(cè)量中，nk往往較小，這時(shí)按式（10）求得的可能不穩(wěn)定。理想方法是通過(guò)平滑的各期狀態(tài)參數(shù)向量的平滑值，根據(jù)平滑后的殘差向量來(lái)估計(jì)，但平滑計(jì)算較復(fù)雜，計(jì)算量也很大。為了化簡(jiǎn)計(jì)算，同時(shí)消除值的不穩(wěn)定性，可將式（10）計(jì)算的各期值，按權(quán)（各期觀測(cè)值的個(gè)數(shù)）取平均值近似地計(jì)算的估值，若以Nk表示第k期以前（包括第k期）各期觀測(cè)個(gè)數(shù)之和，Wk表示第k期以前（包括第k期）各期之和，經(jīng)推導(dǎo)可得下列公式：

用卡爾曼濾波法處理數(shù)據(jù)的一般步驟是：首先，獲取部分對(duì)象的監(jiān)測(cè)數(shù)據(jù)來(lái)確定系統(tǒng)狀態(tài)量模型。其次，提取建立卡爾曼濾波模型所需的預(yù)測(cè)噪聲量和量測(cè)噪聲量，確定系統(tǒng)狀態(tài)量初值。再次，建立系統(tǒng)卡爾曼濾波模型，處理觀測(cè)數(shù)據(jù)[5]。其中，根據(jù)對(duì)監(jiān)測(cè)數(shù)據(jù)的分析選擇合適模型，提取相應(yīng)的預(yù)測(cè)噪聲量和量測(cè)噪聲量是關(guān)鍵，需要根據(jù)具體工程數(shù)據(jù)分析獲得。

4 應(yīng)用實(shí)例

以濱州大橋運(yùn)行監(jiān)測(cè)項(xiàng)目為例。獲得采用動(dòng)態(tài)RTK模式解算的大橋監(jiān)測(cè)三維坐標(biāo)數(shù)據(jù)，經(jīng)過(guò)統(tǒng)計(jì)分析，此處采用位置、速度、加速度模型，設(shè)t為監(jiān)測(cè)周期，設(shè)某時(shí)刻觀測(cè)點(diǎn)的位置為Xk，其瞬時(shí)速率為，瞬時(shí)加速度為，將加速度的瞬時(shí)變化看作隨機(jī)干擾（動(dòng)態(tài)噪聲）。記此點(diǎn)的狀態(tài)向量為：為9×1的矩陣，I為3×3的單位陣，系統(tǒng)狀態(tài)方程為：

可得系統(tǒng)狀態(tài)方程：

Lk+1為3×1觀測(cè)值矩陣，Bk+1為3×9觀測(cè)模型矩陣，為觀測(cè)噪聲，可得卡爾曼濾波迭代的五個(gè)公式：

三維坐標(biāo)數(shù)據(jù)處理采用卡爾曼濾波處理結(jié)果（如圖1所示），綠色曲線為濾波后結(jié)果，紅色曲線為實(shí)測(cè)數(shù)據(jù)，濾波效果明顯，數(shù)據(jù)波動(dòng)符合實(shí)際監(jiān)控情況。

圖1 觀測(cè)數(shù)據(jù)濾波效果

濾波數(shù)據(jù)單位權(quán)中誤差會(huì)發(fā)生變化，濾波值的中誤差越來(lái)越小，精度逐步提高，趨于穩(wěn)定值。有關(guān)結(jié)果如表1所示（這里僅列出10組數(shù)據(jù)）。

由表1看出，監(jiān)測(cè)值與濾波值之間殘差的最大值為-4.588mm，最小值為0.377mm，殘差有正有負(fù)，有很強(qiáng)的隨機(jī)性，預(yù)測(cè)誤差為1.08cm，對(duì)采用動(dòng)態(tài)RTK模式解算數(shù)據(jù)而言，預(yù)測(cè)效果較好。

表1 變形監(jiān)測(cè)值與相應(yīng)的濾波值

5 結(jié)束語(yǔ)

本文詳細(xì)探討了卡爾曼濾波在監(jiān)測(cè)應(yīng)用中如濾波模型的選擇，模型初值的設(shè)定和精度評(píng)價(jià)方法，并結(jié)合具體的工程實(shí)例，論證了該方法在動(dòng)態(tài)變形監(jiān)測(cè)中的可行性，濾波解算結(jié)果具有較高的精度?？柭鼮V波適合估計(jì)線性隨機(jī)差分方程描述的隨機(jī)過(guò)程的狀態(tài)變量，在變形監(jiān)測(cè)工程領(lǐng)域有廣泛的應(yīng)用前景。