张志刚, 梁月吉, 任 超, 黄仪邦, 潘亚龙

(桂林理工大学 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得到迅速发展, 应当有效结合地基和星载多元化数据。对此, 还需要开展更为深入的研究。