水稻株高是水稻本身以及土壤、 水文、 气象等因素的综合反映, 是水稻长势监测的重要指标。 遥感以其宏观性、 时效性、 周期性等特点, 在水稻株高反演中具有不可替代的作用。 光学遥感能够直接获取植被水平分布信息, 利用光谱反射率、 植被指数以及纹理因子等参数, 与水稻株高建立定量关系, 实现大范围快速反演^[1]。 但光学遥感对植被垂直分布信息探测难度较大, 而且水稻生长季内多云雨天气, 又给光学数据获取带来较大的困难, 这些原因大大限制了其在水稻株高反演中的应用^[2]。 合成孔径雷达(synthetic aperture radar, SAR)具有全天时、 全天候、 穿透性的优势, 成为水稻株高反演的重要手段之一。

目前基于SAR的水稻株高反演方法主要有三类: 基于辐射传输理论(radiative transfer theory, RTT)的水稻株高反演方法、 基于SAR干涉测量(interferometry SAR, InSAR)以及基于极化干涉测量(polarimetric SAR interferometry, PolInSAR)的水稻株高反演方法^{[3, 4, 5, 6]}。 其中, 基于PolInSAR的水稻株高反演方法是通过干涉和极化信息的有效组合实现了观测维度的拓展, 将目标竖直方向上不同位置的散射中心分离开来, 估计确定具有最高相干性的散射分量, 然后利用通道差分法^[7]、 模型解算法^{[8, 9]}, 计算水稻株高, 该方法既有InSAR对地表植被散射体的空间分布敏感的特性, 又具有极化SAR对植被散射体的形状和方向敏感的特性^[10]。 近年来, PolInSAR已成为植被高度反演研究热点, 尤其是在森林高度反演方面得到验证和应用^{[8, 11]}。

不同的研究者从不同的角度进行了PolInSAR理论、 方法和应用研究。 目前基于PolInSAR的植被高度反演方法主要有两类: 通道差分法、 模型解算法。 通道差分法是通过直接选择或计算得到代表植被冠层散射和代表地表散射相位并作差, 进而得到水稻株高。 该方法计算简单, 但选择和计算得到的代表植被冠层散射和代表地表散射分量并不准确, 存在较大误差。 模型解算法是基于PolInSAR散射模型开展的反演算法, 有严密的物理模型的支撑, 反演精度较为准确。 目前采用的PolInSAR散射模型主要以双层散射模型最为广泛, 而双层体散射主要有随机体散射体/地表二层模型RVoG模型^{[8, 12]}、 方向体积层模型(oriented volume over ground, OVoG)^[13]。 尽管OVoG模型更贴合真实场景, 水稻层内部结构是垂直分布的, 消光系数极化状态存在依赖性, 但是这种依赖性并不强烈, 且OVoG模型应用需要许多假设, 涉及参数多, 导致对模型的推演和应用较困难。 RVoG模型简单, 涉及的参数少, 反演精度较高, 是森林研究中最常见的模型^[14]。 然而, 水稻作为低矮植被的代表, 具有其独特的生长环境, 下垫面常常为水面, 地表散射容易发生镜面散射, 传统的RVoG模型并不完全适用^{[15, 16]}。

传统的RVoG模型未考虑植被生长环境, 反演结果精度不高, 对水稻株高反演应用有局限性。 因此, 针对水稻植株结构与生长环境等特点, 利用多时相TanDEM-X数据, 对9个时相的TanDEM-X极化干涉SAR数据进行距离向频谱滤波, 求取8个复相干系数, 进行去相干处理, 建立适用于水稻田特性的RVoG模型, 构建基于RVoG模型的极化干涉SAR水稻株高反演方法, 进而评定反演结果。

研究区位于西班牙塞维利亚西南部, 瓜达尔基维尔河的南部, 地中海西北沿岸, 面积约为30 km², 地理位置位于36.85° — 37.38° N, 5.91° — 6.23° W之间。 该区农作物主要为水稻, 水稻常年在五月左右种植, 十月前后收割, 一年一熟。 水稻品种大都是普通栽培稻, 播种方式为撒播种植, 在整个水稻物候期间稻田内有一层水覆盖。 图1为Google earth获取QuickBird卫星数据制作的研究区正射影像。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区地理位置与范围(右图4个红图框为研究区所选的水稻样区)Fig.1 Location and scope of study site and four rice sampling areas of ground measurement campaign coloured in red

由于基于PolInSAR植被株高反演需要较低的时间去相干、 大空间基线数据^[17], 星载卫星TanDEM-X系统, 是双站模式, TanDEM-X卫星发射电磁波, TanDEM-X和TerraSAR-X同时接受回波信号。 因此, TanDEM-X系统可以忽略时间去相干和大气的影响, 另外2015年4月到9月, TanDEM-X系统进行一次科学实验, 提高到2~3 km的空间基线, 是反演农作物的最佳数据。 因此, 在水稻生育期内获取了9个不同时相的TanDEM-X双极化CoSSC数据, 中心波长9.65 GHz, 带宽为150 MHz, 极化方式为HH和VV极化, 该数据为一发双收的双站SAR数据, 包括已配准的单视复数据对, 详细信息如表1所示。 图2给出了TanDEM-X双极化CoSSC数据彩色合成图。

表1

Table 1

表1(Table 1)

表1 TanDEM-XCoSSC信息列表 Table 1 List of TanDEM-X CoSSC image information 获取时间 (年-月-日) 年积日 入射角 /(° ) 分辨率/m (距离向× 方位向) 像元大小 (距离向 m) 像元大小 (方位向m) 垂直波数kz (弧度· m-1) 2015-06-04 155 22.71 1.17× 6.6 0.91 2.18~2.45 2.48 2015-06-15 166 22.71 1.17× 6.6 0.91 2.18~2.45 2.48 2015-06-26 177 22.73 1.17× 6.6 0.91 2.18~2.45 2.48 2015-07-07 188 22.73 1.17× 6.6 0.91 2.18~2.45 2.47 2015-07-18 199 22.73 1.17× 6.6 0.91 2.18~2.45 2.48 2015-07-29 210 22.74 1.17× 6.6 0.91 2.18~2.45 2.48 2015-08-09 221 22.73 1.17× 6.6 0.91 2.18~2.45 2.49 2015-08-20 232 22.73 1.17× 6.6 0.91 2.18~2.45 2.48 2015-08-31 243 22.73 1.17× 6.6 0.91 2.18~2.45 2.48

表1 TanDEM-XCoSSC信息列表 Table 1 List of TanDEM-X CoSSC image information

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 TanDEM-X双极化CoSSC数据彩色合成图(获取时间2015年7月7日, R=VV; G=HH-VV; B=HH) (a): 主影像; (b): 辅影像Fig.2 TanDEM-X color composite images of CoSSC data (time: 7 July 2015, R=VV, G=HH-VV, and B=HH) (a): Master image; (b): Slave image

SAR卫星过境的同时, 获取了详细的地面数据。 9个时相中地面数据的采集大部分为卫星过境当天采集, 个别时相样区地面采集数据与卫星过境时间相差1~2 d。 地面数据主要来自均匀分布于研究区的4块水稻样区, 图1中红色方块为样区。 对每一块样区, 采集数据包括水稻物候和株高、 样区面积(hm²)、 测量/播种/收割日期、 稻株密度(株· m^-2)、 稻穗密度(穗· m^-2)、 产量(kg· hm^-2)、 灌溉方式、 盐碱性、 病虫害等信息。 其中, 利用皮尺测量水稻株高, 其精度为cm级。 此外, 研究区日值气象数据(温度、 湿度、 风速、 风向、 降水等)由西班牙农业气候信息系统(SIAR)提供。

采用基于RVoG模型的极化干涉植被高度反演方法, 构建复相干系数与水稻株高之间的定量关系。 复相干系数与目标场景、 传感器特性等多种因素密切相关, 可以表示为

γ=γtempγgeomγprocγSNRγBQγ~(1)

其中, γ _temp为时间去相干系数; γ _geom为几何去相干系数; γ _proc为干涉处理过程中造成的去相干系数; γ _SNR为信噪比去相干系数; γ _BQ为比特量化去相干系数; 为体去相干。 在极化干涉SAR水稻株高反演之前, 需要剔除水稻场景以外去相干因素对复相干系数的影响。 图3给出了极化干涉SAR水稻株高反演基本流程。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 极化干涉SAR水稻株高反演基本流程Fig.3 Basic process of rice crop height inversion using PolInSAR data

①首先对9景TanDEM-X CoSSC双极化SAR数据进行距离向频谱滤波补偿^[18], 对于TanDEM-X CoSSC数据, 由于其空间基线大, 传感器空间方位差异对复相干系数的影响较大, γ _geom小于0.8, 因此必须考虑γ _geom的影响。 ②对距离向频谱滤波后的TanDEM-X数据进行共轭相乘得到干涉图, 并利用主辅影像斜距、 入射角、 中心波长等获取影像时的轨道参数计算平地相位并去除平地效应。 ③全极化干涉SAR系统中每个散射单元为2× 2的散射矩阵, 但对于TanDEM-X双极化系统, 每个散射单元为2× 1的散射矩阵。 因此将TanDEM-X数据对应的S矩阵转化为T₄矩阵, 并利用Boxcar方法对T₄矩阵进行滤波处理, 窗口大小: 21× 21, 以便减小相干斑点噪声, 提高图像的辐射分辨率。 ④计算8个复相干系数, 分别为主辅影像HH, VV, HH-VV, HH+VV复相干系数、 两个优化的相干性系数^[7](OPT1, OPT2)、 PD相位差分法得到的复相干系数(PDHigh, PDLow)。 然后进行去相干因素影响的剔除。 由于TanDEM-X CoSSC为一发双收的双站模式数据产品, 其主辅图像同时获取, 不受时间去相干的影响, 即γ _temp为1。 γ _proc主要是由主辅图像配准引起的, 而TanDEM-X CoSSC数据产品具有较高的配准精度, 因此γ _proc可以忽略。 对于TanDEM-X and TerraSAR-X数据, 考虑到成像区域主要为农田场景, γ _BQ约为0.965^[15]。 在干涉测量研究中, 信噪比去相干γ _SNR通常被忽略, 因为它只对低散射区造成影响。 但对于稻田场景, 由于其X波段SAR图像上的后向散射系数通常在-25~-5 dB之间, 尤其是幼苗期后向散射系数很低, 因此, 信噪比去相干γ _SNR不能忽略。 基于TanDEM-X头文件数据计算SNR, 采用公式(2)计算信噪比去相干 γSNR[4]。

γSNR=11+SNR-1(2)

⑤利用剔除去相干因素影响惑的8个复相干系数, 构建一种基于RVoG模型的水稻植株高度反演算法, 实现研究区9个时相水稻株高反演。

2.2 基于RVoG模型的水稻植株高度反演算法

RVoG模型包含植被层及其下垫面, 其中植被层厚度近似为h_v且包含随机取向粒子集合。 对于一般植被目标, RVoG主要考虑植被层体散射以及地表面散射。 对于水稻而言, 由于其生长周期内下垫面多为水面, 易发生镜面反射, 因此地表面散射可忽略不计, 此外水稻与下垫面之间的二次散射较强, 需要考虑。 因此, 对于水稻目标, RVoG模型可以表示为^[15]

γ~=eiφ0γV+sinkzhvkzhvmDB(ω)1+mDB(ω)(3)

其中, k_z为垂直波数, φ ₀为地形相位, h_v为植被层厚度也就是我们所要反演的水稻株高, m_DB(ω )为二次散射对应的地体幅度比。

γ _v表示体相干, 取决于趴标地物散射体的消光系数和植被层厚度h_v

γV=∫0hve2σz/cosθe(ikzz)dz∫0hve2σz/cosθdz(4)

其中, σ 为消光系数, θ 为雷达数据入射角。

于是, 式(3)和式(4)中包含了4个未知数: h_v, m_DB, σ , φ ₀。 式(3)可以写成复平面内的线性形式

γ~=eiφ0γV+mDB(ω)1+mDB(ω)(γDB-γV)=eiφ0[γV+L(ω)(γDB-γV)](5)

其中 γDB=sinkzhvkzhv, L(ω)=mDB(ω)1+mDB(ω)∈[0, 1]。由此可以看出, 理想情况下在复平面上不同极化方式所对应的复相干系数都落于以γ _DB eiφ0和γ _V eiφ0为端点的直线段上, 其中φ ₀为植被下的地形相位, 而γ _V eiφ0则为具有完全体散射的复相干系数。 然而在实际获取的全极化SAR影像对上, 地形相位和完全体散射相干系数都无法直接得到。 鉴于此, 利用不同极化方式对应的8个复相干系数进行最小二乘线性拟合, 获得式(5)对应的直线, 然后再计算地形相位φ ₀。 对于传统的RVoG模型, 通过式(5)对应的直线与复平面内单位元的交点来确定地形相位; 然而, 对于水稻而言, 还需要考虑二次散射的去相干影响, 因此, 利用式(5)对应的直线与复平面内半径为γ _DB的圆的交点, 来确定地形相位。

基于上面构建的极化干涉SAR水稻株高反演模型可以看出, 共有5个模型参数。 这里采用模式搜索法^[19]的迭代算法进行水稻株高反演。 由于模式搜索法是直接优化方法, 不需要求解目标函数的导数。 因此, 该方法具有方法简单, 计算方便快速等优点。 总体流程如图4所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 基于RVoG模型的针对水稻特性的水稻株高反演算法的流程图Fig.4 Flow chart of rice crop height inversion algorithm based on RVoG model for rice characteristics

①基于剔除各种去相干因素影响的8个不同极化通道的复相干系数, 利用最小二乘拟合获得式(5)在复平面内对应的直线。 ②由于模式搜索法的搜索步长影响搜索速度和全局搜索能力, 根据经验知识, 确定未知量的搜索步长、 搜索范围。 给定h_v初值、 范围、 步长, h_v∈ (0.05, 1.5), 步长Δ h_v为0.05。 然后计算γ _DB, 进而获得地形相位φ ₀。 ③设定σ , m_DBmax, m_DBmin取值范围、 初始值、 步长。 σ ∈ (0, 10), 步长Δ σ 为1 dB· m^-1; m_DBmax∈ (-10, 10), 步长Δ m_DBmax为0.1 dB; m_DBmin∈ (-10, 10), 步长Δ m_DBmin为0.1 dB。 利用设定的h_v, σ , m_DBmax, m_DBmin根据式(5), 计算模拟每一个像元的复相干系数。 需要注意的是, 消光系数σ 在代入模型中要进行转换, 其转换系数为0.23。 ④确定模式搜索法基本参数。 确定轴向方向E, 加速因子α , 缩减率β ∈ (0, 1)。 由于是四个自变量参数, 因此, 设定轴向方向为E=(0, 0, 0, 1; 0, 0, 0, -1; 0, 0, 1, 0; 0, 0, -1, 0; 0, 1, 0, 0; 0, -1, 0, 0; 1, 0, 0, 0; -1, 0, 0, 0), 加速因子α 设为2, 缩减率β 设为0.5。 ⑤求出每一个像元8个复相干系数中表真实地表散射的最大复相干系数、 最小复相干系数与模拟的复相干系数的距离和的最小值, 其具体公式见式(6)。

minφ0, hv, σ, mDBmax, mDBminγ(kz, w→max)-γsim(kz, φ0, hv, σ, mDBmax)γ(kz, w→min)-γsim(kz, φ0, hv, σ, mDBmin)(6)

其中 w→为极化通道, γ _sim为利用模型模拟的复相干系数。

基于模式搜索方法求出每一个像元8个复相干系数中表真实地表散射的最大复相干系数、 最小复相干系数与模拟的复相干系数的距离的最小值, 如果满足2条件(条件1: 设定迭代次数为2 000次; 条件2: 设定函数阈值为0.001)其中之一, 终止迭代, 提取此时的给定的h_v; 如果不满足, 回到第2步继续运算。 ⑥遍历每景影像中的每一个像元, 从而得到该时相下的水稻株高反演结果。

基于9个时相的TanDEM-X数据, 利用第2节所述方法进行研究区水稻株高反演, 并利用地块水稻样区的地面数据进行验证与精度评价, 如图5所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 样区1, 2, 3, 4水稻株高反演结果与实测数据对比 红色圆点折线表示地面实测水稻株高数据曲线, 绿色圆点表示反演结果, 其误差线表示反演结果标准差Fig.5 Comparison of rice crop heights retrieved and measured for sampling area 1, 2, 3 and 4 The red dotted line represents the measured rice crop height, and the green dots represent the inversion result. What’ s more, the error line represents the standard deviation of the inversion result

对样区1来说, 在水稻生长前期, 第155天(6月4日)、 第166天(6月15日)株高实测数据在0.1~0.3 m之间, 而反演结果大于0.8 m, 明显高于真实值, 而且这两个时相反演结果对应的标准差也较大。 由此, 在水稻生长前期, 当株高较小时(< 0.3 m), 该方法反演结果较差, 株高被过高估计。 第177天(6月26日)株高在0.3 m左右, 反演结果与实测数据吻合较好, 但反演结果标准差较大。 在水稻生长中期、 成熟期, 即第188天(7月7日)、 第199天(7月18日)、 第210天(7月29日)、 第221天(8月9日)、 第232天(8月20日)、 第243天(8月31日), 株高反演与地面实测值非常接近, 绝对差值在0.1 m之内, 并且其标准误差也远低于第155天、 第166天、 第177天的对应值。 由此可以看出, 当株高大于0.4 m时, 该方法反演精度较高, 而且稳定性好。

对样区2来说, 在第155天、 第166天、 第188天、 第199天, 株高实测数据在0.1~0.4 m之间, 而这4个时相的反演结果均比其实测结果高0.2 m以上, 虽然第177天反演结果较好, 但是当株高小于0.4 m时, 该方法稳定性较差。 第210天、 第221天、 第232天、 第243天株高实测数据在0.4~1.0 m之间, 其反演结果与实测株高吻合较好, 尤其是第210, 221, 232和243天, 其反演结果与实测值绝对误差在0.1 m以内, 而且反演结果标准差也较小。 因此, 该方法在株高大于0.4 m时, 反演结果较好, 稳定性较高。

样区3和样区4的反演结果与样区1和样区2类似, 当水稻株高小于0.4 m时, 反演结果过高估计, 当株高高于0.4 m时, 反演精度较高(绝对误差< 0.15 m), 稳定性也较好。

通过上述4个水稻样区反演结果分析可以看出, 当株高低于0.4 m时, 反演结果较差, 稳定性不高; 然而当株高大于0.4 m时, 反演结果精度较高而且稳定性较好, 与实测值的绝对误差在0.1 m以内, 最差不超过0.15 m。 图6(a)给出了四个样区9个时相下的水稻株高反演结果。 由图6(a)看出, 当株高低于0.4 m时, 反演结果基本都分布在直线“ y=x” 上侧, 即反演结果被明显高估, 绝对误差在0.1~0.8 m之间; 当株高大于0.4 m时, 反演结果较好, 决定系数R²为0.86, 均方根误差RMSE为6.79 cm[图6(b)]。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 4个样区9个不同时相水稻株高反演结果与实测水稻株高的对比Fig.6 Comparison of rice crop heights retrieved and measured from four sampling areas under 9 time phases

通过反演结果验证与分析发现, 基于TanDEM-X数据的水稻株高反演算法, 在株高低于0.4 m以下时, 效果不好。 这主要是由于所使用的TanDEM-X数据, 空间基线为2~3 km。 根据Cloude^[20]提出的株高反演最优范围指数k_v=k_zh_v/2(其中k_z为垂直波数, h_v为水稻株高)及其最优范围1≤ k_v≤ 1.5, 结合表1中给出的的k_z值, 可以得出本文所使用的TanDEM-X数据对应的最佳植被高度反演范围为0.81 m≤ h_v≤ 1.21 m。 因此, 当水稻株高低于0.4 m时, 与所使用的TanDEM-X数据对应的最佳植被高度反演范围差距较大, 无法较好地反映出水稻植株的较大体散射量, 进而导致反演结果较差。

对于株高低于0.4 m的水稻, 基于TanDEM-X数据, 利用极化干涉的方法进行高度反演需要考虑两方面因素, 一是要尽可能的增大基线, 提高数据对目标体散射敏感性; 二是要保持较高的相干性, 不能无限制增大基线。 根据频谱偏移超过系统带宽的理论^[20], 得出TanDEM-X数据极限空间基线为3.2 km, 即当空间基线超过3.2 km时, 相干性将为0。 根据这一极限条件, TanDEM-X数据, 对于高度大于0.33 m的水稻株高反演是有效的。 通过k_v最优范围计算得出, 当TanDEM-X数据基线范围在2.625 km< B< 3.2 km之间时, 利用极化干涉的方法能够较好地反演株高在0.33~0.4 m的水稻株高。 对于高度低于0.33 m的水稻, 基于TanDEM-X数据, 利用极化干涉的株高反演方法, 很难得到较高精度的水稻株高反演结果。

利用9个时相的TanDEM-X CoSSC数据, 构建了一种基于RVoG模型的极化干涉水稻株高反演算法, 实现了多时相的水稻株高反演, 结果表明: 当株高高于0.4m时, 反演精度较高, 决定系数R²为0.86, 均方根误差RMSE为6.79 cm; 然而当株高低于0.4 m时, 反演结果被高估而且稳定性差, 这主要是因为本研究所使用的TanDEM-X数据空间基线对株高小于0.4 m的水稻敏感性不足。 当TanDEM-X数据空间基线范围在2.625 km< B< 3.2 km之间时, 利用极化干涉的方法能够较好地反演株高在0.33~0.4 m的水稻株高。

致谢: 感谢德国航空太空中心(The German Aerospace Center, DLR)(项目编号: NTI-POLI6736)提供本文所用的TanDEM-XCoSSC数据。