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

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

(桂林理工大學 a.測繪地理信息學院; b.廣西空間信息與測繪重點實驗室, 廣西 桂林 541006)

0 引 言

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

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

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

1.1 GPS反射信號原理

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

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

(1)

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

(2)

式中:ψ表示直射信號和反射信號的相位差。

根據上式, 利用二階多項式擬合作為直射分量, 從SNR中分離出衛星反射信號分量。圖1是美國板塊邊界觀測計劃PBO中P041測站PRN26

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

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

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

RSN,m=10RSN/20。

(3)

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

(4)

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

1.2 基于Kalman濾波的衛星反射信號去噪

設SNR反射信號時間序列觀測值X為

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

(5)

式中:k表示觀測歷元個數,k=1,2,…,i。

Kalman濾波是一種通過遞歸方法解決數據線性濾波的問題,它把最小均方誤差作為最佳估計準則,利用前一時刻的估計值和當前時刻的觀測值來更新對狀態變量的估計,得到當前時刻的最優估計值[12]。主要包括描述狀態向量的過程方程和描述觀測向量的觀測方程,設在歷元時刻狀態GPS反射信號向量估計值為Xk,其對應協方差矩陣為Qk,k+1時刻的GPS反射信號估計值Xk+1可表示為[13-14]

Xk+1=AkXk+Wk;

(6)

Zk=HkXk+Vk。

(7)

式中:Ak代表狀態轉移矩陣;Wk代表過程合成噪聲向量;Zk為動態系統在時間k的觀測向量;Hk代表觀測矩陣；Vk觀測噪聲向量。正態分布的白色噪聲,過程噪聲協方差矩陣為Q, 觀測噪聲協方差矩陣為R, 即

(8)

(9)

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

2 實驗分析

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

圖2 P041測站Fig.2 P041 station

選擇P041測站2013年第132～240天數據進行實驗。利用TEQC軟件解算觀測文件得到SNR

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

由圖3可見, 衛星反射信號與高度角正弦值存在顯著的線性關系, 衛星反射信號在5°～20°下波動明顯。進一步對比發現, 同一顆衛星在上升、下降段波動趨勢并不一致, 這是由于不同時間段內測站周圍地表環境和方位角差異造成的。

圖3 衛星上升、下降段衛星反射信號分離結果Fig.3 Separation results of satellite reflection signals from rising and falling satellites

采用方案2對衛星反射信號去噪后與方案1結果對比分析。限于篇幅, 僅給出PRN10和PRN20衛星采用方案1、2處理后的反射信號與衛星高度角正弦值之間的線性關系, 如圖4所示。

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

圖4 方案1、2衛星反射信號處理結果Fig.4 Satellite reflection signal processing results of Scheme 1 and 2

圖5 P041測站不同衛星的相對延遲相位與土壤濕度反演估算結果關系Fig.5 Relationship between relative delay phase and soil moisture of different satellite at Station P041

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

表1 P041測站土壤濕度反演精度評定統計Table 1 Statistics of soil moisture inversion accuracy assessment of Station P041

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

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

3 結 論

基于GNSS連續運行參考站P041的背景下, 利用GPS-IR進行土壤濕度準確監測對于環境和水循環具有重要的意義和應用前景。本文利用Kalman濾波對衛星反射信號去噪后反演土壤濕度, 經理論分析和實驗表明:

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

(2)利用Kalman濾波能夠有效消除各衛星反射信號中的噪聲, 能夠更加精確地得到衛星反射信號, 同時能夠有效削弱各衛星的異常跳變。實現較高時空分辨率的土壤濕度反演, 彌補了傳統數據處理的缺陷。

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