基于Kalman濾波的GPS-IR反演土壤濕度研究

張志剛, 梁月吉, 任 超, 黃儀邦, 潘亞龍

(桂林理工大學(xué) a.測繪地理信息學(xué)院; b.廣西空間信息與測繪重點(diǎn)實(shí)驗(yàn)室, 廣西 桂林 541006)

0 引 言

土壤濕度是水文、氣象和生態(tài)環(huán)境研究中重要的狀態(tài)參數(shù), 開展土壤濕度估算對(duì)于氣候氣象預(yù)報(bào)、洪水災(zāi)害以及水資源循環(huán)等相關(guān)研究具有重要意義[1]。近年來, 隨著全球?qū)Ш叫l(wèi)星系統(tǒng)的高速發(fā)展, GPS多路徑效應(yīng)被證實(shí)可以有效獲取地表環(huán)境變化參數(shù), 因其具有成本低、高效率、高分辨率等優(yōu)點(diǎn), 成為當(dāng)前國內(nèi)外研究的熱點(diǎn)問題。Masters等[2]利用GPS雙基地雷達(dá)進(jìn)行土壤水分遙感的實(shí)驗(yàn), 初步研究表明該技術(shù)能夠有效反映土壤濕度的時(shí)空變化。Larson等[3]利用連續(xù)運(yùn)行GPS接收機(jī)接收到的多徑反射幅度估算土壤濕度, 證實(shí)了多徑反射幅度與土壤濕度存在一定的相關(guān)性, 提出了全球定位系統(tǒng)干涉反射測量(global navigation satellite signal-interferometer and reflectometry, GPS-IR)反演土壤濕度的遙感技術(shù)方法。GPS-IR是一種新型的遙感測量手段, 近年來已被廣泛應(yīng)用于估算近地表土壤濕度、雪深、海平面高度和植被指數(shù)等環(huán)境參數(shù)及其變化[4-5]。Zavorotny等[6]提出了GPS地面多徑模型, 證明了利用相對(duì)延遲相位反演土壤濕度比振幅更有優(yōu)勢, 反射系數(shù)與距離地表5 cm內(nèi)的土壤介電常數(shù)非常敏感。Chew等[7]建立電動(dòng)正向模型成功地模擬了植被冠層對(duì)GPS SNR的影響, 當(dāng)植被較高時(shí), 在低仰角(5°～15°)下多徑幅度感知植被變化的能力較弱。王迎強(qiáng)等[8]證明機(jī)載GPS反射信號(hào)能夠有效地反演土壤濕度。敖敏思等[9]經(jīng)實(shí)驗(yàn)表明, 利用SNR觀測值通過指數(shù)模型能夠有效反映土壤濕度的變化趨勢, 最大有效監(jiān)測范圍約45 m。梁月吉等[10]利用LS-SVM模型滾動(dòng)式多星融合, 實(shí)現(xiàn)了較長時(shí)間高精度的土壤濕度估算, 改善了部分異常跳變現(xiàn)象。隨著美國旋風(fēng)衛(wèi)星導(dǎo)航系統(tǒng)(cyclone global navigation satellite system, CYGNSS)和英國TDS-1衛(wèi)星(technology demonstrate satellite-1)的發(fā)展。星載GNSS-R地表環(huán)境反演成為當(dāng)前的發(fā)展趨勢, 但是由于數(shù)據(jù)量、配置和軌跡異質(zhì)性等限制, 使用星載數(shù)據(jù)進(jìn)行的研究仍處于初步階段, 地基GPS-IR將會(huì)成為有效驗(yàn)證數(shù)據(jù), 因此, 獲取高精度的地基反演值具有重要意義。由于測站周邊的土壤表面粗糙度、植被覆蓋、土壤濕度等地表因素對(duì)GPS反射信號(hào)的影響均不一樣。以上研究較少考慮GPS直反射信號(hào)分離的精度, 導(dǎo)致土壤濕度估算精度難免受到影響。

基于以上研究, 本文提出一種基于卡爾曼濾波(Kalman filter)輔助的土壤濕度反演模型。Kalman濾波是Kalman和Bucy在1958年提出的一種利用線性系統(tǒng)狀態(tài)方程, 通過對(duì)輸入的觀測數(shù)據(jù)濾波處理, 有效消除信號(hào)中的噪聲和干擾的影響, 已在通信、導(dǎo)航、控制等多領(lǐng)域得到了較好的應(yīng)用[11]。與傳統(tǒng)濾波方法相比, Kalman濾波能夠適應(yīng)信號(hào)和噪聲都是非平穩(wěn)的過程, 具有更強(qiáng)的泛化能力。因此, 本文利用Kalman濾波對(duì)GPS反射信號(hào)去噪, 建立基于Kalman濾波的GPS反射信號(hào)去噪模型。

1 GPS-IR土壤濕度反演原理

1.1 GPS反射信號(hào)原理

GPS接收機(jī)接收到來自測站附近地表環(huán)境反射的衛(wèi)星信號(hào)和衛(wèi)星直射信號(hào)相互疊加干涉形成的信噪比(signal to noise ratio, SNR)觀測數(shù)據(jù)。利用SNR中的物理參數(shù)相對(duì)延遲相位可實(shí)現(xiàn)土壤濕度的反演。SNR可以用直射和反射信號(hào)表示為[10]

RSN=Sd(0)+Sm(φ),

(1)

式中:RSN代表信噪比觀測值;Sd(0)代表直射信號(hào);0代表初始相位;Sm(φ)代表反射信號(hào);φ代表反射信號(hào)相位值。當(dāng)天線高度較低時(shí),GPS直反射信號(hào)的頻率相同。SNR可用直角形式表示為

(2)

式中:ψ表示直射信號(hào)和反射信號(hào)的相位差。

根據(jù)上式, 利用二階多項(xiàng)式擬合作為直射分量, 從SNR中分離出衛(wèi)星反射信號(hào)分量。圖1是美國板塊邊界觀測計(jì)劃PBO中P041測站PRN26

圖1 SNR和直射分量擬合圖Fig.1 SNR and direct component fit map

號(hào)衛(wèi)星2013年第135天的SNR和直射信號(hào)擬合圖。可見, 在低衛(wèi)星高度角下, SNR波動(dòng)顯著, 受多路徑效應(yīng)影響較大, 土壤濕度信息主要包含在多路徑影響成分中。其中GPS反射信號(hào)由于電磁波和介質(zhì)之間非線性相互作用, 有很大的噪聲, 信號(hào)本身容易出現(xiàn)失真。如何有效去除噪聲等干擾是本文研究探討的問題。

Chew等[7]研究表明: SNR觀測值與φ之間存在一種正弦或余弦關(guān)系, 而且GPS土壤濕度僅與多路徑反射信號(hào)相關(guān), 那么, 去除GPS衛(wèi)星直射信號(hào)后的多路徑反射信號(hào)隨sinθ呈現(xiàn)指數(shù)變化, 因此, 要將SNR的原始單位(dB)變換為以功率(W)為單位的SNR時(shí)間序列

RSN,m=10RSN/20。

(3)

以功率為單位的SNR反射信號(hào)與高度角正弦值sinθ存在某一固定頻率的正弦(或余弦)函數(shù)關(guān)系,表示為

(4)

式中:RSN, m代表多路徑反射信號(hào);Am代表反射信號(hào)幅度;h代表天線高;λ代表波長;θ代表衛(wèi)星入射高度角;φ代表相對(duì)延遲相位。根據(jù)上式,利用非線性最小二乘擬合求得相對(duì)延遲相位φ。已有研究表明, 相對(duì)延遲相位與地表土壤濕度之間存在較強(qiáng)的相關(guān)性[11]。因此, 本文基于相對(duì)延遲相位建立土壤濕度反演模型。

1.2 基于Kalman濾波的衛(wèi)星反射信號(hào)去噪

設(shè)SNR反射信號(hào)時(shí)間序列觀測值X為

X=[X1,X2,…,Xk],

(5)

式中:k表示觀測歷元個(gè)數(shù),k=1,2,…,i。

Kalman濾波是一種通過遞歸方法解決數(shù)據(jù)線性濾波的問題,它把最小均方誤差作為最佳估計(jì)準(zhǔn)則,利用前一時(shí)刻的估計(jì)值和當(dāng)前時(shí)刻的觀測值來更新對(duì)狀態(tài)變量的估計(jì),得到當(dāng)前時(shí)刻的最優(yōu)估計(jì)值[12]。主要包括描述狀態(tài)向量的過程方程和描述觀測向量的觀測方程,設(shè)在歷元時(shí)刻狀態(tài)GPS反射信號(hào)向量估計(jì)值為Xk,其對(duì)應(yīng)協(xié)方差矩陣為Qk,k+1時(shí)刻的GPS反射信號(hào)估計(jì)值Xk+1可表示為[13-14]

Xk+1=AkXk+Wk;

(6)

Zk=HkXk+Vk。

(7)

式中:Ak代表狀態(tài)轉(zhuǎn)移矩陣;Wk代表過程合成噪聲向量;Zk為動(dòng)態(tài)系統(tǒng)在時(shí)間k的觀測向量;Hk代表觀測矩陣；Vk觀測噪聲向量。正態(tài)分布的白色噪聲,過程噪聲協(xié)方差矩陣為Q, 觀測噪聲協(xié)方差矩陣為R, 即

(8)

(9)

自適應(yīng)Kalman濾波在計(jì)算過程中會(huì)不斷自我修正過程噪聲和觀測噪聲矩陣。根據(jù)上文的分析, 由于SNR觀測數(shù)據(jù)變化頻率高, 相鄰歷元間SNR反射分量變化緩慢, 并且與sinθ呈線性關(guān)系,所以k+1時(shí)刻的GPS反射分量值與k時(shí)刻的GPS反射分量值可以視為相同, 因此狀態(tài)轉(zhuǎn)移矩陣設(shè)置為A=[1]。

2 實(shí)驗(yàn)分析

本文選取來自于美國板塊邊界觀測計(jì)劃PBO (http: //xenon.colorado.edu/portal)提供的P041測站(105.194 267°W, 39.949 493°N)監(jiān)測數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。該測站能夠連續(xù)記錄采樣率為30 Hz 的L2 SNR觀測值, 較早就開展土壤濕度反演分析, 具有一定的代表性。從P041測站(圖2)可看出, 該測站周邊地形較為平坦、開闊且植被稀少, 主要以草為覆蓋植被, 有利于土壤濕度監(jiān)測。站點(diǎn)均采用鋼制三角支架安置, 接收機(jī)型號(hào)為TRIMBLE NERT9, 采用SCIT的天線罩, 天線型號(hào)為TRM59800.80。

圖2 P041測站Fig.2 P041 station

選擇P041測站2013年第132～240天數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。利用TEQC軟件解算觀測文件得到SNR

和高度角數(shù)據(jù)。為了驗(yàn)證Kalman濾波對(duì)衛(wèi)星反射信號(hào)去噪的可行性, 本文設(shè)置兩種方案進(jìn)行對(duì)比分析: 方案1, 用二次多項(xiàng)式擬合分離各衛(wèi)星的直、反射信號(hào), 采用非線性最小二乘法估算各衛(wèi)星的相對(duì)延遲相位; 方案2, 用二次多項(xiàng)式擬合分離各衛(wèi)星的直、反射信號(hào), 并對(duì)衛(wèi)星反射信號(hào)Kalman濾波去噪, 估算相對(duì)延遲相位。限于篇幅, 僅給出PRN10號(hào)第140天利用二次多項(xiàng)式擬合分離得到的去噪前的衛(wèi)星反射信號(hào)結(jié)果。

由圖3可見, 衛(wèi)星反射信號(hào)與高度角正弦值存在顯著的線性關(guān)系, 衛(wèi)星反射信號(hào)在5°～20°下波動(dòng)明顯。進(jìn)一步對(duì)比發(fā)現(xiàn), 同一顆衛(wèi)星在上升、下降段波動(dòng)趨勢并不一致, 這是由于不同時(shí)間段內(nèi)測站周圍地表環(huán)境和方位角差異造成的。

圖3 衛(wèi)星上升、下降段衛(wèi)星反射信號(hào)分離結(jié)果Fig.3 Separation results of satellite reflection signals from rising and falling satellites

采用方案2對(duì)衛(wèi)星反射信號(hào)去噪后與方案1結(jié)果對(duì)比分析。限于篇幅, 僅給出PRN10和PRN20衛(wèi)星采用方案1、2處理后的反射信號(hào)與衛(wèi)星高度角正弦值之間的線性關(guān)系, 如圖4所示。

可以看出, 方案2能夠很好地估計(jì)衛(wèi)星反射信號(hào)與高度角正弦值(sinθ)的線性變化趨勢, 有效削弱部分異常跳變點(diǎn), 如圖4a中下降段高度角正弦值為0.15和圖4b中上升段高度角正弦值為0.14和下降段高度角正弦值為0.22附近的異常值。為了進(jìn)一步評(píng)估Kalman濾波在GPS-IR土壤濕度反演中的性能, 利用非線性最小二乘法解算各衛(wèi)星的相對(duì)延遲相位。僅給出部分衛(wèi)星估算結(jié)果進(jìn)行分析, 各衛(wèi)星相對(duì)延遲相位與土壤濕度關(guān)系如圖5所示。

圖4 方案1、2衛(wèi)星反射信號(hào)處理結(jié)果Fig.4 Satellite reflection signal processing results of Scheme 1 and 2

圖5 P041測站不同衛(wèi)星的相對(duì)延遲相位與土壤濕度反演估算結(jié)果關(guān)系Fig.5 Relationship between relative delay phase and soil moisture of different satellite at Station P041

可見, 隨著土壤濕度的變化, 各衛(wèi)星的相對(duì)延遲相位均能夠作出響應(yīng)。年積日第168～169天的PRN06號(hào)衛(wèi)星、第151～152天的PRN32號(hào)衛(wèi)星的方案1估算結(jié)果均出現(xiàn)較大異常跳變現(xiàn)象。整個(gè)時(shí)段內(nèi), 測站部分衛(wèi)星的相對(duì)相位延遲與土壤濕度之間的線性回歸方程系數(shù)R2如圖6所示。可見, 不同衛(wèi)星的相對(duì)相位延遲與土壤濕度之間的相關(guān)性均不一樣。進(jìn)一步對(duì)比各衛(wèi)星發(fā)現(xiàn), 方案2估算結(jié)果與土壤濕度的相關(guān)性優(yōu)于方案1的估算結(jié)果。為了進(jìn)一步綜合評(píng)定各方案的性能, 本文采用R2、均方根誤差(RMSE)、平均絕對(duì)誤差(MAE)進(jìn)行評(píng)定。只給出了P041測站部分衛(wèi)星2種方案精度評(píng)定結(jié)果， 見表1。

表1 P041測站土壤濕度反演精度評(píng)定統(tǒng)計(jì)Table 1 Statistics of soil moisture inversion accuracy assessment of Station P041

圖6 土壤濕度與估算結(jié)果線性回歸分析Fig.6 Linear regression analysis of estimation results and soil moisture

可見, 各衛(wèi)星與土壤濕度存在顯著的線性相關(guān)性, 因此, P041測站采用相對(duì)延遲相位可以有效監(jiān)測天線周圍土壤濕度的變化趨勢。進(jìn)一步對(duì)比發(fā)現(xiàn), 方案1平均相關(guān)系數(shù)R2=0.496, 方案2平均相關(guān)系數(shù)R2=0.586, 相比之下, 本文方法相關(guān)性提高18.1%。因此, 利用Kalman濾波去噪后反演土壤濕度具有可行性和有效性。

3 結(jié) 論

基于GNSS連續(xù)運(yùn)行參考站P041的背景下, 利用GPS-IR進(jìn)行土壤濕度準(zhǔn)確監(jiān)測對(duì)于環(huán)境和水循環(huán)具有重要的意義和應(yīng)用前景。本文利用Kalman濾波對(duì)衛(wèi)星反射信號(hào)去噪后反演土壤濕度, 經(jīng)理論分析和實(shí)驗(yàn)表明:

(1)對(duì)于P041測站周圍單一地表環(huán)境條件下, 相對(duì)延遲相位與土壤濕度存在顯著相關(guān)性。

(2)利用Kalman濾波能夠有效消除各衛(wèi)星反射信號(hào)中的噪聲, 能夠更加精確地得到衛(wèi)星反射信號(hào), 同時(shí)能夠有效削弱各衛(wèi)星的異常跳變。實(shí)現(xiàn)較高時(shí)空分辨率的土壤濕度反演, 彌補(bǔ)了傳統(tǒng)數(shù)據(jù)處理的缺陷。

(3)針對(duì)土壤濕度的反演, 僅僅針對(duì)數(shù)據(jù)處理過程優(yōu)化還不夠, 要結(jié)合有效的反演模型才能得到較大的提升。同時(shí), 星載GNSS-R得到迅速發(fā)展, 應(yīng)當(dāng)有效結(jié)合地基和星載多元化數(shù)據(jù)。對(duì)此, 還需要開展更為深入的研究。