土壤水分含量是重要的农业气象要素之一, 它反映土壤的湿润状况, 是农田旱情预报的主要依据。准确监测土壤水分含量具有重要意义。由于影响土壤水分含量的因素很多, 目前大范围实时监测土壤水分含量仍是一个难题。现有的常规土壤水分含量监测方法主要有取土法、中子仪法和时域反射法等^[1-3]。这些监测方法能测定的样点数量有限, 难于实现大范围地区土壤水分含量的实时动态监测。而遥感技术的发展, 为实现大范围地区土壤水分含量的实时或准实时动态监测提供了条件^[4-6]。

应用MODIS数据, 进行裸土土壤湿度的监测, 目前较为成熟的技术方法是热惯量法^[7-11]。虽然, 热惯量法可以获得较为满意的监测结果, 其缺点是必须要有白天和夜间两时相的遥感数据, 这为实际业务应用带来了困难。尤其在青藏高原, 春季TERRA卫星夜间过境时, 一般为多云天气, 很难获得夜间晴空数据。因此, 应用白天单时相数据监测土壤湿度在青海高原地区有一定的应用前景。郭广猛等人研究认为MODIS第7波段的值和土壤表层水分有很好的相关关系^[12], 并依此建立了土壤湿度遥感监测模型。但是, 经试验, 在青海高原地区, 应用此模型监测的土壤表层湿度经检验精度较低, 不能满足实际应用需求。本文沿袭了遥感地物分类中K-L变换的思想, 采用主成分分析法, 应用MODIS数据的多个波段数据, 建立了青海省多维特征空间土壤湿度监测模型^[13-15]。

在遥感地物分类方法中, 最常见的就是应用特征空间进行分类。通过构建特征空间图可以识别地物类型, 梅安新^[15]分析了训练样本在TM第3, 4波段上的不同作物二维特征空间分布, 认为尽管二维特征空间中不同作物域之间有一定的交集, 但不同域的主体仍有明显差别, 如果空间维数增加, 域和域之间的区别将更加明显, 有助于地物类型的识别。

图 1是2005年3月8日对应于地面土壤湿度测点MODIS资料第1、第2波段值。可以看出, 两者有一定的相关关系, 但并不是完全线性相关, 空间点的分布遵循一定规律, 然而分布较为散乱, 从统计学意义分析, 这种二维空间特征分布方差较大, 两个波段既有相关性, 又蕴涵了较多的地物信息。选取和土壤水分相关系数较大的多个波段, 通过主成分分析, 建立多维特征空间, 是可以反演地面土壤水分的^[16-17]。

图 1

图 1. 利用主成分分析方法特征空间变换 Fig 1. Spacial characteristic transformation by principal component analysis

图 1中, 两个波段的亮度值之间有一定的相关性, 同时, 这两个波段又包含不同的信息。两个波段的平均值分别为u₁和u₂, 其数据的分散程度显示两个波段的相关程度, 分散程度越大, 相关性越小, 反之亦然。如果两个波段的二维点状分布图中所有点都集中于一个很小的区域, 那么这些数据很难提供大量有用的信息, 主成分分析的目的就是通过平移和旋转坐标轴, 使原数据信息重新分布到新轴上称第1主成分(PC1), 条件为新轴方向上资料方差最大, 与PC1正交的轴称第2主成分(PC2)。

取MODIS第7、第21、第29、第31波段作为4维特征空间。选取以上波段的依据是4个波段值和地面土壤湿度有较好的线性关系(相关系数绝对值大于0.6, 且通过α=0.001的显著性检验)。另外, 第7波段位于近红外范围(1.0~2.5 μm), 第21波段在中红外(2.5~6 μm)范围之内, 第29和第31波段处于远红外范围(>6 μm)之内, 以保证尽可能覆盖在MODIS红外波段范围。选取4个波段后, 求取对应35个地面测点的4个波段值, 并对其做标准化, 输入矩阵为标准化后的数据。采用式(1)标准化。X_it是标准化后的波段值, X_i是波段的平均值, σ_i是均方差, i为MODIS的第7、第21、第29、第31中的某一波段。这样做的目的是使各波段之间的数据有可比性, 并且做了一次相对辐射订正。

(1)

应用SAS统计软件分析, 得到第1主成分的累积贡献率为0.8757, PC2的累积贡献率为0.9829。用PC1已经可以足够反映原始数据的信息。PC1对应的特征向量为

(2)

特征向量是各波段值标准化后因子权重值, PC1已经可以解释原始数据信息, 因此, 本研究中仅对其进行分析。

完成以上处理之后, 要建立土壤湿度监测模型, 必须使波段值和土壤湿度之间发生联系。通过分析, 建立了PC1和土壤湿度回归模型, 如图 2所示, 横轴为K-L变换的PC1, 纵轴为10 cm土壤湿度。

图 2

图 2. 第1主成分与土壤湿度回归结果 Fig 2. Regression result of principal component factors and soil moisture

(3)

最终, 求得多维特征空间值和10 cm土壤湿度之间的模型为式(4), 其中, M_S为反演所得的10 cm土壤湿度, b_zn代表n波段标准化后的值。

(4)

上述特征空间建模中, 应用的数据是1000 m分辨率的MODIS数据, 分析各波段值和土壤湿度的相关关系发现, 第1、第2波段其实和土壤湿度也有很好的相关关系。因此, 可以通过资料融合技术, 来改变运算后图像的几何精度, 提高模型的空间分辨能力。

图像融合需解决的是数据的空间配准、融合模型的建立与优化以及融合方法的选择等问题。研究中采用的低分辨率图像是多维特征空间模型运行后的土壤湿度监测结果, 高分辨率图像是同一时相MODIS数据的第1或第2波段的数据。由于两种数据都来自同一遥感器、同一时相, 融合过程中不存在空间配准的问题。

分析波段值、K-L变换后的PC1, 测点土壤湿度三者相关关系(表 1), b2为MODIS第2波段值, PC1是第1主成分, M_S10为10 cm土壤湿度。b2和PC1的相关系数为0.67, 和10 cm土壤湿度的相关系数为-0.45, 样本数是30, 通过0.01显著性检验。

表 1

表 1 图像融合各要素相关系数 Table 1 Correlation coefficient of factors in image infusing

表 1 图像融合各要素相关系数 Table 1 Correlation coefficient of factors in image infusing

建立融合模型的方法很多, 有彩色技术、图像变换、数学运算。在这里, 依据b₂, PC1, M_S10的相关关系, 构建融合模型。

(5)

式(5)中, D_f为融合后的图像值, D_g, D_d分别是高、低分辨率的图像值。f₁和f₂分别是权重系数, 式(5)给出了求算方法, R_bS是MODIS第2波段和土壤湿度的相关系数, R_PS是K-L变换后的PC1和土壤湿度的相关系数。

利用已建立的多维特征空间模型, 对2004年2月18日的土壤湿度监测数据进行了检验。青海省中部大部分地区有云覆盖, 但东部地区为晴空区, 并且, 该区域部分低海拔测点有农业气象土壤测湿数据, 可以用于模型检验。

青海省气象灾害标准中利用土壤湿度将干旱程度分为5级, 即偏湿、无旱、轻旱、中旱、重旱, 本研究遵循这种分级方法。采用MODIS第2波段值和空间变换结果实施融合。按照前文所述, 求得f₁=41.3, f₂=58.7。图 3是经过融合处理后青海省土壤湿度分布情况, 其结果使空间分辨率提高了4倍, 可以监测250 m的地理单元。至于融合后土壤水分的监测精度, 取决于辐射校正和第2、第7、第21、第29、第31波段对地面土壤水分的响应程度。2004年青海省东部农业区部分地区有重旱发生, 海西的唐古拉山乡, 玉树的部分地区草原旱情也较为严重。从图 3可以看出, 模型结果与实际情况相符。

图 3

图 3. 2004年2月18日干旱监测结果 Fig 3. Result of monitoring drought on Feb 18, 2004

表 2是实际监测值和模型值的对比, 共有35个测点的对照值。检验预报结果中, 如果允许模型预测结果的误差为正负一个级别。则模型监测的准确度为80%。

表 2

表 2 遥感监测值和实测值对照(单位: %) Table 2 Comparison between RS monitoring value data and observing data(unit : %)

表 2 遥感监测值和实测值对照(单位: %) Table 2 Comparison between RS monitoring value data and observing data(unit : %)

研究结果表明, 建立的青海高原土壤湿度多维特征空间监测模型可以在较高精度水平下大面积监测土壤湿度, 而且所需遥感数据资料为单时相, 较热惯量模型需白天、黑夜两时相资料的要求, 可以更为方便地获取遥感资料, 模型更易于在高原地区使用。

研究过程中尚存在如下问题: ①地面土壤测湿是单点数据, 模型建立的过程中存在以点代面的问题。遥感监测最小单元是1 km×1 km, 得到的反射率值是混合像元的整体反射率, 若需提高模型精度、混合像元的分解问题^[18]须进行深入研究。②研究中的资料处理仅做了简单的太阳高度角订正, 相当于将地面看作郎伯体, 但事实并非如此。另外, 对于辐射校正, 仅仅是做了相对辐射校正, 按照理论, 遥感反演地面温度, 须做绝对大气辐射校正。③利用多维特征空间建立的土壤监测模型, 仅仅进行了一次检验, 模型的运用仍处于试验阶段, 在投入业务应用之前, 仍有大量的工作。如主成分分析波段的选择, 以及采用的波段数量对最终结果的影响等等。