近年来, 在全球气候变化的大背景下, 世界各地水资源严重短缺, 干旱、洪涝等极端事件发生频繁, 水分限制条件下的大范围作物生长动态模拟研究成为作物生长研究的热点之一。传统的作物生长动态模拟主要通过作物生长模拟模型实现。作物生长模拟模型是对作物生长发育及产量形成过程中一系列生理生化过程的数学描述, 是对实际生长情况的模拟。然而, 气候条件的变化以及模型的初始值及作物和土壤参数不确定性的存在, 往往导致模拟结果出现较大偏差^[1]。

Wiegand等^[2]在20世纪70年代指出, 遥感信息可以用来改进作物生长模型的模拟精度。此后, 人们开展了一系列将遥感信息(如植被指数、后向散射系数等)应用于作物生长模拟模型的研究^[2-9](表 1), 但由于逐日高分辨率的土壤水分数据难以获得, 大多数研究都没有考虑土壤水分对作物生长模拟模型的影响。Ines等^[9]将AMSR-E土壤湿度产品以及叶面积指数LAI (leaf area index)同时同化到DSSAT (decision support system for agrotechnology transfer)作物生长模拟模型, 在Iowa州Story郡的实验结果表明, 该方法的相关系数R=0.8, 最终产量模拟的均方根误差为1400 kg/ha。由于使用的AMSR-E土壤湿度产品空间分辨率低, 模拟精度有待于进一步提高。

土壤水分指土壤渗流层的水分含量, 是陆地和大气能量交换过程中的关键因子, 影响植被的蒸散发及光合作用, 在降水、径流、下渗和蒸散发水文过程中发挥着至关重要的作用^[10]。目前, 遥感手段已成为监测大区域范围内土壤水分时空分布和变化的主要方法^[11]。

基于遥感技术对土壤水分时空分布的大范围监测, 按遥感测量手段的不同, 可分为光学遥感、主动微波和被动微波三类^[12]。由于单一手段无法获得逐日高空间分辨率的土壤水分信息, 因此人们开展了许多利用多种传感器进行土壤水分联合反演的研究。张显峰等^[13]用MODIS的光学与热红外波段反演土壤水分的基准值, 利用AMSR-E传感器的X波段反演土壤水分的日变化量, 然后集成二者建立土壤水分协同反演模型, 得到具有较高时空分辨率的土壤水分反演结果。该方法仍然需要地面土壤水分观测数据, 且空间分辨率为1 km。

为获取高时空分辨率的土壤水分数据, 本研究探讨利用AMSR-E (advanced microwave scanning radiometer for EOS)土壤湿度产品和ASAR (advanced synthetic aperture radar)后向散射系数数据, 基于高级积分方程模型AIEM (advanced integral equation model), 在不依赖地面土壤水分观测数据的情况下, 发展适合大区域尺度的土壤水分反演方法。最后, 以陕西省渭北台塬为研究区, 对本文方法进行验证。

1 主、被动微波土壤水分协同反演方法

主动微波数据不受云层的影响, 空间分辨率高, 但由于其重访周期长, 反演结果在时间上不连续。被动微波遥感虽然具有较好的时间连续性, 但其空间分辨率较低, 受地表粗糙度和植被覆盖等的影响较大。针对两者的优点与不足, 本文提出一种主被动微波高时空分辨率土壤水分协同反演方法。本方法包括两步, 技术路线如图 1所示。

第1步, 利用不依赖地面观测数据的土壤水分基准日期值反演方法获得高分辨率土壤水分基准日期值:用补偿后的AMSR-E土壤湿度产品作为“实测土壤水分”, 并将ASAR后向散射系数匹配到与AMSR-E相同的空间分辨率, 最后基于AIEM模型进行数据拟合, 得到土壤水分的表达式, 并应用到日期A的30 m分辨率的ASAR数据中, 获得ASAR一个重访周期内土壤水分的基本情况。

第2步, 依据逐日土壤水协同反演方法获得逐日高分辨率的土壤水分:在第一步获得的基准日期值的基础上, 利用补偿后的AMSR-E土壤湿度产品计算日期B与基准日期A的土壤水分的变化量, 并将其重采样到与ASAR数据相同大小的像元, 最后将基准日期值与日变化量相加即可得到该重访周期内任意日期B的高分辨率土壤水分。

1.1 不依赖地面数据的土壤水分基准日期值反演方法

1.1.1 AMSR-E土壤湿度产品补偿

研究表明, AMSR-E土壤湿度产品虽然能够提供全球粗分辨率的土壤水分数据, 但其值相对地面实测的土壤水分偏小^[14-15]。为使AMSR-E土壤湿度产品更好地与地面实测数据保持一致, 本研究对AMSR-E土壤湿度产品进行补偿。陈洁^[14]对2006-2008年我国西北陕西(包含本文研究区域)、新疆、甘肃等5省148个农气站0~10 cm的土壤水分观测数据与AMSR-E土壤湿度产品进行了对比分析, 结果表明, 2006-2008年AMSR-E反演的西北地区土壤湿度数据整体变化不大, 但相比地面实测数据要小40%左右。本文研究区域内只有一个农业气象站, 无法进行统计分析, 因此采用上述研究结论对AMSR-E土壤水分进行补偿, 用M_v表示补偿后的土壤水分, M_v⁰表示AMSR-E土壤湿度产品的值, 两者的关系式为

$ {M_{\rm{v}}} = 1.67 \times M_{\rm{v}}^0。 $

(1)

1.1.2 ASAR与AMSR-E数据匹配

ASAR数据的像元大小为30 m, 而AMSR-E数据的像元大小为25 km。为得到与AMSR-E土壤湿度产品对应的后向散射系数, 需要将ASAR数据与AMSR-E数据匹配到相同的空间尺度, 本文采取的方法是将ASAR数据分辨率重采样至与AMSR-E相同的分辨率。由于重采样后的一个像元内包含ASAR数据原始像元众多, 地物类型多样, 而AIEM模型对建筑、水体和高大植被区域反演效果较差, 因此在重采样前, 采用该区域同一时期的Landsat真彩色合成影像, 在目视判读的基础上, 对原始ASAR数据进行粗分类, 将ASAR地物分为建筑、裸地农田和水体山林三类, 然后将建筑与水体山林区域进行掩膜, 剔除模型不适用的地物覆盖区域。掩膜处理后, 对ASAR数据进行重采样, 相当于仅对裸地农田区域的后向散射系数取均值作为重采样后ASAR像元的后向散射系数。这样不仅可以保证模型的可靠性和适用性, 而且能大幅度减少参与运算的原始ASAR像元数。

1.1.3 简化AIEM模型参数拟合与反演

在主动微波土壤水分反演中, 积分方程模型(integrated equation model, IEM)能在一个很宽的地表粗糙度范围内再现真实地表后向散射情况, 已广泛应用于微波地表散射、辐射的模拟和分析^[16]。近年来, IEM模型经过不断改进和完善, 发展成为高级积分方程模型AIEM, 使模型模拟结果和精度得到不断提高。已有的研究表明, 裸露地表的后向散射系数主要受地表粗糙度、地表含水量和入射角的影响, AIEM模型可简化^[17]为

$ {\sigma _{{\rm{pq}}}}\left( {{\rm{dB}}} \right) = {A_{{\rm{pq}}}}\left( \theta \right) \times {\log _{10}}\left( {{M_{\rm{v}}}} \right) + {B_{{\rm{pq}}}}\left( \theta \right) \times {\log _{10}}\left( {{Z_{\rm{s}}}} \right) + {C_{{\rm{pq}}}}\left( \theta \right)。 $

(2)

本文采用的数据为HH/VV的双极化组合。渭北台塬地势相对平坦, 地表粗糙度变化不明显, 因此可以利用两种极化组合消去Z_s, 从而得到针对ASAR双极化数据的土壤水分反演模型:

$ {\log _{10}}\left( {{M_{\rm{v}}}} \right) = A\left( \theta \right) \times {\sigma _{{\rm{HH}}}} + B\left( \theta \right) \times {\sigma _{{\rm{vv}}}} + C\left( \theta \right), $

(3)

式中, M_v(%)表示补偿后的AMSR-E土壤体积含水量, σ_HH和σ_vv分别表示重采样到相同空间分辨率后ASAR水平极化和垂直极化的后向散射系数, A(θ), B(θ)和C(θ)是待拟合的参数。对式(3)进行参数拟合, 将拟合得到的模型应用于30 m分辨率的ASAR后向散射系数中, 即可获得高空间分辨率的土壤水分反演基准日期值。

1.2 逐日土壤水协同反演方法

根据张显峰等^[13]提出的陆表土壤水分协同反演策略, 陆地表面土壤水分的变化可以分解为两部分:基准值和变化量。本文借鉴其协同思想, 将陆地表面土壤水分的变化分为基准日期值和日变化量两部分, 表达式如下:

$ {\rm{s}}{{\rm{m}}_{i, j}}\left( t \right) = {m_{i, j}}\left( t \right) + \Delta {m_{i, j}}\left( t \right), $

(4)

式中, m_{i, j}(t)为一定时间尺度T内, 遥感影像的第i列第j像元的土壤水分基准日期值, Δm_{i, j}(t)为时间尺度T内, t时刻第i列第j像元土壤水分与基准日期值的差异变化量, sm_{i, j}(t)为该像元通过遥感反演得到的高分辨率土壤水分。

对于同一区域, 如果时间尺度T较小, 植被和地表粗糙度没有显著变化, 那么估算日期B地表土壤水分对基准日期A土壤水分的变化量, 就只需计算补偿后的日期B与日期A的AMSR-E土壤水分的差值。因此, 本文利用日期A的ASAR数据反演土壤水分的基准日期值m, 利用补偿后的日期B的AMSR-E土壤水分与补偿后的日期A的土壤水分差值作为土壤水分的变化量Δm (图 1), 并将后者重采样至前者的像元大小, 最后利用式(4)即可得到日期B高分辨率的土壤水分。

2 研究区域概况与数据预处理 2.1 研究区概况

渭北台塬地处陕西省中部, 涵盖宝鸡、咸阳、西安等市。本文研究区为渭北台塬西部地区, 主要包括千阳县、宝鸡市辖区、凤翔县、岐山县、眉县和扶风县等六县区(图 2)。渭北台塬西部农作物以冬小麦、玉米等为主, 根据联合国粮农组织(FAO)和维也纳国际应用系统研究所(IIASA)构建的世界土壤数据库HWSD (Harmonized World Soil Data-base version 1.1), 其土壤类型以土娄土和雏形土为主, 在渭河及其支流周围形成小面积的冲击土。渭北台塬属于温带大陆性季风气候区, 年降水量为550~700 mm, 年日照时数为2300~2500小时, 年均温为7.0~12.8 C, 昼夜温差可达10 C, 全年无霜期170天以上^[18]。然而, 降水量在空间, 尤其是时间上存在明显的波动性, 往往是旱涝灾害的多发区, 因此获取较高精度的大范围连续时间的土壤水分分布数据, 对该区水分限制条件下的作物长势动态监测具有重要意义。

2.2 数据获取及预处理

ASAR是ENVISAT上搭载的合成孔径成像雷达系统, 工作在C波段, 过境赤道的时间约为当地时间11:30 (降轨)和23:30 (升轨)。本实验的ENVISAT ASAR数据是HH, VV交替极化模式的ASA_APP_1P数据, 像元大小为30 m, 时间为2005年12月12日和2006年1月16日, 共计两景数据。该数据的产品级别为L1B级, 记录的是DN值, 需要进行辐射定标、几何校正以及滤波等操作后方可使用。本文采用欧洲空间局(European Space Agency, ESA)提供的NEST (Next ESA SAR Tool-box)软件包进行辐射校正, 主要步骤为:原始数据提取, 由幅度图像得到能量图像, 辐射定标得到后向散射系数图像。由于渭北台塬西部地势相对平坦, 因此直接使用NEST软件包中的Geometry/ Terrain Correction进行地形校正。

斑点噪声会降低图像的空间分辨率和信噪比, 严重影响图像的可解译性, 甚至可能导致地物特征的消失^[19]。为减少SAR斑点噪声对图像准确度的影响, 本实验利用Gamma滤波器对图像进行滤波处理, 选取3×3, 7×7, 10×10的滤波窗口。由于滤波后要对ASAR数据进行尺度匹配, 最终不同大小的窗口得到尺度匹配后的ASAR数据像元的后向散射系数差异均在0.02 dB以内, 对简化的AIEM参数拟合的决定系数影响不明显, 因此最终选取处理速度较快的3×3的滤波窗口。

AMSR-E传感器于2002年由美国NASA发射的Aqua卫星搭载升空。AMSR-E的赤道过境时间约为当地时间下午13:30 (升轨)和凌晨01:30 (降轨)。本文中AMSR-E土壤水分值来自美国冰雪数据中NSIDC (National Snow and Ice Data Center)的AMSR-E-L3-DailyLand-V06数据集中的重采样结果。

在地面观测资料方面, 本文从国家气象局气象数据共享平台获取土壤水分数据集中0~10 cm的土壤水分相对湿度, 对反演结果进行验证。由于AMSR-E土壤湿度数据是体积含水量, 国家气象局土壤水分数据集测定的是每旬第8天的相对土壤湿度, 因此在进行验证前, 需要将相对湿度数据转换成体积含水量, 公式为

$ {Q_{\rm{v}}} = {Q_{\rm{r}}} \times \rho \times {f_{\rm{c}}}, $

(5)

式中, Q_v(%)表示土壤体积含水量; Q_r(%)表示土壤相对湿度; ρ(g/cm³)表示土壤容重; f_c(m³/m³)表示田间持水量。其中, 田间持水量和土壤容重与土壤类型有关。

由于本文采用的凤翔县农业气象站数据仅用作验证, 且凤翔县土壤类型单一, 因此, 根据王春娟^[20]关于凤翔县0~10 cm土壤水文特征的数据资料, 土壤容重取1.08 g/cm³, 田间持水量取0.263 m³/m³。

3 结果与讨论 3.1 土壤水分基准日期值反演结果与分析

本文获得渭北台塬西部地区2005年12月12日和2006年1月16日的两景ASAR数据, 数据获取时间均为北京时间23:30左右, AMSR-E数据则分别选取2005年12月13日01:30和2006年1月17日01:30的数据, 对应日期的两种数据分别相差两个小时左右。本文假设土壤水分在这两个小时内基本上无变化。按照本文提出的不依赖地面数据的土壤水分基准日期值反演方法, 将ASAR数据匹配到AMSR-E相同的像元大小。由于覆盖研究区域的AMSR-E数据像元数为9个, 其中2个像元没有有效值, 予以剔除, 因此得到7组2005年12月12日、2006年1月16日ASAR数据匹配后的后向散射系数和补偿后的AMSR-E土壤湿度, 详见表 2。

利用两天的数据, 式(3)的参数拟合结果分别为

$ {\log _{10}}\left( {{M_{\rm{v}}}} \right) = 0.0292 \times {\sigma _{{\rm{HH}}}} + 0.0343 \times {\sigma _{{\rm{vv}}}} + 3.0883, $

(6)

$ {\log _{10}}\left( {{M_{\rm{v}}}} \right) = 0.0397 \times {\sigma _{{\rm{HH}}}} + 0.0158 \times {\sigma _{{\rm{vv}}}} + 2.9915, $

(7)

其中, 式(6)为2005年12月12日的拟合结果, 拟合的决定系数R²=0.82;式(7)为2006年1月16日的拟合结果, 拟合的决定系数为R²=0.81。

将拟合得到的模型进行变换, 得到土壤水分基准日期值的表达式分别为

$ {M_{\rm{v}}} = {10^{0.0292 \times {\sigma _{{\rm{HH}}}} + 0.0343 \times {\sigma _{{\rm{vv}}}} + 3.0883}}, $

(8)

$ {M_{\rm{v}}} = {10^{0.0397 \times {\sigma _{{\rm{HH}}}} + 0.0158 \times {\sigma _{{\rm{vv}}}} + 2.9915}}, $

(9)

其中, 式(8)是2005年12月土壤水分的表达式, 式(9)是2006年1月土壤水分的表达式。

将上述公式分别应用到对应的30 m分辨率的ASAR后向散射系数中, 得到高空间分辨率的土壤水分基准日期值反演结果(图 3)。

对像元进行统计, 发现反演得出的土壤体积含水量在170‰~286‰之间, 约占总数的80%。2005年12月12日的平均土壤水分值为233‰, 方差为0.054; 2006年1月16日的平均土壤水分值为244‰, 方差为0.049。可以看出, 2006年1月16日土壤体积含水量更高、更集中。从空间分布上看, 两者均呈现南高北低的整体态势, 这可能与渭河穿过研究区域南部有关。土壤水分基准日期值反演结果与AMSR-E土壤水分的分布基本上一致。

3.2 逐日土壤水分协同反演与验证

凤翔县农业气象站(简称凤翔站)是研究区域内唯一的农业气象站。2005年12月, 凤翔站有3次(3, 18和28日)地面观测土壤水分值, 2006年1月只有一次观测。在反演出的2005年12月12日土壤水分基准日期值的基础上, 利用前面提到的逐日土壤水分协同反演方法, 反演了2005年12月有地面农业气象站观测资料的3, 18和28日共3天的土壤体积含水量(图 4)。从图 4可以看出, 12月3日的土壤含水量整体上高于其他3天, 渭河流域的土壤水分明显高于其他区域。图 4中绿色三角形表示凤翔站的位置, 其地面观测相对湿度与协同反演的土壤体积含水量见表 3。

利用式(5), 将凤翔站0~10 cm的土壤相对湿度观测值转换成土壤体积含水量, 并与土壤体积含水量反演值进行对比, 结果如图 5所示, 两者的决定系数R²=0.92, 均方根误差为0.025, 说明反演结果与地面观测值较为接近, 具有较高的可信度。

3.3 讨论

本文采用AMSR-E土壤湿度产品代替地面实测土壤水分, 利用ASAR和AMSR-E协同反演的方法, 反演出30 m分辨率的逐日土壤水分, 解决了地面土壤水分同步数据难于获取的问题, 提高了土壤水分反演的空间分辨率, 取得较好的模拟效果。但是, 由于真实土壤水分的不确定性, 难免有一定的偏差。

首先, 鉴于AMSR-E土壤湿度产品数值整体上偏小, 本文采用前人包含该区域同时期的148个农业气象站资料, 对AMSR-E土壤水分标定的结果进行补偿, 以使其更接近真实值, 但结果仍然存在误差。其次, 在反演建模时, 进行粗分类去除模型不适用的区域, 尽量使参与建模的地物适用于模型, 但仍然可能会有模型不适用的地物混杂其中, 增加了建模的不确定性。最后, 在计算土壤水分变化量时, 采用AMSR-E土壤水分变化量的方法估算日期B土壤水分与日期A土壤水分基准日期值的差值, 这也具有一定的不确定性。上述不确定性的存在, 使反演结果与真实值之间存在一定误差, 但经过处理后, 建模和反演结果的误差在可接受的范围内, 实例验证结果也证明了结果的可靠性。

特别地, 2006年1月8-19日的最高气温均高于零度, 12-15日的最高气温高于最低气温的绝对值, 属于冻融交替状态, 此外, 研究区渭河流域之外的区域土壤水分普遍偏低。在土壤含水量较低时, 大部分水分以吸附水的形式存在, 一般不会冻结, 只有少量的自由水会冻结, 从而影响土壤含水量估计。根据曹梅盛等^[21]与张立新等^[22]对冻融土壤的微波辐射的研究, 野外土壤, 尤其是农田表层会有一层很薄的干土, 即使干土下湿土含水量较大, 土壤冻结后总体的发射率、比辐射率和亮温等也不会有很大变化。因此, 在反演土壤水分时, 不考虑土壤冻融对土壤水分的影响, 在一定程度上会增加模型的不确定性, 但是影响较小。根据晋锐等^[23]对冻土土壤水分的研究, 同化SSM/I 19 GHz亮温后可以使土壤水分的均方根误差减小0.018。后续研究中, 可以考虑采用这种同化策略, 进一步减小不确定性。

4 结论

本文利用主被动微波遥感结合的方式, 在缺乏土壤水分地面观测数据的条件下, 估算渭北台塬地区30 m分辨率的逐日土壤水分, 是不依赖地面观测数据反演土壤水分的一种有益尝试。主动微波穿透能力强, 土壤的介电常数和土壤含水量具有很强的相关性, 因此利用主动微波数据反演表层土壤水分具有较高的机理性和可信度。实验结果也表明, 本文提出的主被动微波协同反演土壤水分基准日期值的方法有较高的相关性。同时, 利用被动微波时间分辨率高、主动微波微空间分辨率高的特性, 在充分借鉴前人方法的基础上, 提出一套主被动微波结合反演高空间分辨率日土壤水分的方法, 具有一定的实用性。凤翔县农业气象站的土壤水分观测数据验证结果表明, 土壤水分反演值与地面观测值之间有很高的一致性。

本文提出的方法在一定程度上可以监测逐日高分辨率土壤水分的真实状况, 有利于逐日模拟作物生长模型以及土壤水分耦合的研究, 为实现水分限制条件下大范围作物生长模拟奠定了良好的基础。后续研究将通过土壤水分与作物生长模型的结合, 实现对区域尺度高分辨率作物生长发育情况的实时监测。

致谢:

感谢欧洲太空局提供的ENVISAT影像以及中国气象资料共享中心的地面实测数据支持。