何勤业,卢涵宇*,袁咏仪,卢天健
(1.贵州大学 大数据与信息工程学院, 贵州 贵阳 550025;2.贵州六盘水三力达科技有限公司, 贵州 六盘水 532001)
土壤水分是气候系统中的一个关键状态变量,对理解陆地-大气能量、水分和碳循环通量具有重要意义[1]。土壤水分在洪水和干旱预报、气候变化、天气预报、水资源管理和农业生产等方面有着广泛的应用,对人类社会和环境有着重要的影响[2],因此,精确度和高空间分辨率的土壤水分全球测绘对于满足这些应用的需求至关重要。
目前最广泛应用的测量地表土壤水分的方法是被动微波遥感,它能全天候全天时观测数据。由于微波遥感的1.41 GHz(L波段)具有更好的穿透大气和植被覆盖能力及受射频干扰(radio grequency interference,RFI)影响较小特点,因此在大多数与土壤水分相关的研究以该频段为主要研究对象。第一个搭载L波段微波辐射计的卫星是在2009年11月2日由欧空局发射的土壤水分和海洋盐度(soil moisture and ocean salinity,SMOS)卫星。该卫星可以提供地表的双极化和大范围入射角(5°~60°)的亮温数据,同时这也是第一颗监测全球土壤水分和海水盐度的卫星[3];但是L波段辐射计的空间分辨率只有40 km,并不足以满足水文气象、生态学、水资源管理等应用超过空间分辨率为10 km的使用需求,所以美国国家航空航天局(national aeronautics and space administration,NASA)于2015年1月31日在加利福尼亚州发射土壤水分主动-被动(soil moisture active passive,SMAP)任务[4]。该卫星搭载着空间分辨率为3 km的L波段雷达与空间分辨率为40 km的L波段辐射计,可以通过几种主被动结合的降尺度方法,提供空间分辨率为9 km的全球土壤水分[4],但是,在2015年7月7日SMAP雷达故障,无法提供数据。而为了能继续提供高分辨率数据,SMAP使用了Backus-Gilbert(BG)optimal interpolation对原始的SMAP Level 1B亮温产品(L1B_TB)进行插值处理从而得到空间分辨率为9 km的亮温数据[5]。
目前土壤水分的降尺度方法可以大致分为基于经验模型的方法、基于数理模型的方法和基于微波遥感数据融合的方法。其中基于经验模型的方法要先获取的高分辨率的地表参数,再利用获取的地表参数建立与微波土壤水分的模型,并将模型参数用于微波土壤水分的降尺度处理;但降尺度结果的准确程度主要取决于此类经验模型的精确性,而且为了能更好的构建经验模型还需要获取大量的实地数据,因此限制了经验模型在较大区域内的应用[6]。而基于数理模型的方法主要是通过数理统计模型(如谱分析、多重分形和地学统计等空间插值方法)对微波数据进行降尺度处理,其缺陷在于这些数理模型往往只是对单一的微波数据进行空间上的处理,缺少额外的辅助数据作为支撑,并不能确保降尺度的精确程度[7]。本文使用的是基于微波遥感数据融合的方法,这种方法是通过将2种不同分辨率的微波数据进行融合以达到降尺度的目的。其中文献[8]利用傅里叶变换的方法将主被动微波亮温数据进行结合,但是SMAP的主动微波雷达已经损坏,而且傅里叶变换在分析微波信号这类非平稳信号时存在局限性,会导致降尺度结果出现较大的误差。文献[9]使用时间序列观测回归分析降尺度法研究被动微波中L波段与S波段亮温的线性关系,将S波段作为辅助数据对L波段的亮温进行降尺度处理,但是目前还没有搭载S波段的卫星,而且S波段的分辨率较低远远不及10 km,无法满足土壤水分的应用需求。本文在上述2种方法的基础上进行了改进,首先是在选用微波遥感数据时选用了被动微波的X波段作为L波段亮温降尺度的辅助波段,这是因为X波段的分辨率为10 km,可以满足土壤水分的应用需求,而且X波段的亮温数据可以直接通过AMSR2得到;其次在对X波段与L波段进行融合时使用了二维离散小波变换(2-D DWT)的方法,这种方法在分析非平稳信号时没有傅里叶分析的局限,可以更加精准的对微波信号进行分析与变换。具体步骤是通过2-D DWT分辨率为9 km的X波段亮温与分辨率为36 km的L波段亮温进行分解与融合,从而得到降尺度后的L波段亮温数据。实验结果表明降尺度后的L波段亮温可以满足反演土壤水分的需求,且反演的土壤水分的精度与SMAP 9 km亮温反演的土壤水分近似,具有极高的适用性。
本文主要通过结合L波段与X波段亮温的特性来实现亮温降尺度的研究,并且将小波变换技术引入到了亮温降尺度和土壤水分反演这2个领域,从而拓展了亮温降尺度和土壤水分反演的研究方向。
本研究使用的L波段的亮温数据采用了SMAP-3级辐射计(SPL3SMP版本7)中的亮温数据,SMAP数据于“美国国家冰雪数据中心”(https://nsidc.org/data)进行下载。X波段的亮温数据则使用了AMSR2中10.65 GHz的L3亮温数据,其数据来源于“日本宇宙航空研究开发机构”(https://gportal.jaxa.jp/gpr/),并且选择了位于帕里站点获取的2016年1月1日到12月31日的土壤水分数据作为验证数据,该数据来源于“国家青藏高原科学数据中心”(http://data.tpdc.ac.cn)。
本次研究区域选择为帕里区域。该区域的土地覆盖以草地、裸露地面为主,地形平坦,平均年降水量小于400 mm,且降雨大多集中在5—10月。在帕里区域内一共部署了25个观测站点,在本研究中只用到了其中的8个观测站点,并分为了A、B共2个9 km网格,其中黑色三角点表示土壤水分的观测站点[10]。帕里研究区土壤水分采样点分布如图1所示。
图1 帕里研究区土壤水分采样点分布图
土壤水分主动被动(SMAP)卫星于2015年1月31日发射,运行周期为2~3 d,上升轨道时间为18:00,下降轨道时间为上午06:00[4]。SMAP卫星分别搭载了L波段(1.41 GHz)辐射计和L波段(1.26 GHz)雷达,入射角度为40°[4],可以提供空间分辨率为9、36 km的亮温数据。而空间分辨率为9 km的亮温数据在L波段雷达故障后就无法获取,因此基于对BG 最佳插值的应用,将SMAP的原始观测数据进行了修正从而获得了空间分辨率为9 km的亮温数据。
在本研究中使用了SMAP卫星在帕里实验区域的分辨率36、9 km的亮温数据,以及用亮温反演土壤水分需要的辅助数据,为了方便小波变换的应用选取的范围为16×16的9 km网格,具体的经纬度为27.259 5°~28.534 2°、88.879 7°~90.373 4°。
GCOM-W由Japan Aerospace Exploration Agency(JAXA)在2012年5月18日发射,轨道位于地球上方约700 km处,上升轨道时间为13:30,下降轨道时间为01:30[11]。GCOM-W卫星携带的无源微波辐射计Advanced Microwave Scanning Radiometer 2(AMSR2)是由AMSR-E改进而来,AMSR2的天线旋转周期为1.5 s[11],可以获取1 450 km范围内的数据,使得AMSR2能以2 d为一周期获取全球99%以上的数据[11]。AMSR2的入射角为55°,有6个不同频率的微波波段(6.9、10.7、18.7、23.8、36.5、89 GHz),可以提供空间分辨率为10、25 km的双极化亮度温度。本研究中使用了10.65 GHz(X波段)采集的空间分辨率为10 km亮度温度作为降尺度的辅助数据。
土壤水分是评价亮温降尺度结果的一个重要标准,因此本研究使用了帕里区域中按照布设的8个观测站点(如图1黑色三角点所示)获取的土壤水分。同时为了方便进行土壤水分的对比验证还对这8个观测站点根据其所在的网格进行了平均计算,从而得到这2个9 km网格的土壤水分。
虽然L波段对于土壤水分含量的探测更加精准,但是SMAP的L波段的空间分辨率比AMSR2的X波段低,因此为了提高L波段的空间分辨率,使得L波段更加泛用,本文通过二维离散小波变换(2-D DWT)将L波段与X波段的亮温数据融合,从而得到的同时拥有L波段亮温数据特征和X波段空间分辨率的亮温数据。同时为了对这种降尺度方法的结果进行对比还加入了SMAP 9 km亮温数据(通过Backus-Gilbert(BG)optimal interpolation获取)。为了方便将降尺度后的亮温数据与SMAP的分辨率为9 km的亮温数据进行对比,还需要对AMSR2的X波段亮温数据进行插值,得到X波段分辨率为9 km的亮温数据。
Backus-Gilbert(BG)最佳插值基于Backus-Gilbert理论,是一种对成像微波辐射计数据进行插值的技术,因为保持与微波辐射计相关的天线增益函数的空间分辨率,所以其插值过程称为最优的。SMAP 9 km的亮温数据是通过对SMAP 1B级亮温产品(L1B_TB)中的天线温度(TA)的测量值进行Backus-Gilbert(BG)最佳插值得到的,其插值生成的TA数据在经过校准后,就能在9 km网格上生成SMAP 1C级亮温产品(L1C_TB_E)。
小波变换作为一种分析工具,可以对不同空间分辨率下的亮温数据、土壤水分数据等地球物理分布变量数据进行特性变化的分析。
连续小波变换(CWT)可描述为:假设存在一个能量有限的信号空间为L2(R),且f(t)∈L2(R),可以用小波基将f(t)展开[12],得到f(t)的连续小波变换Wf(a,τ),
(1)
式中:a是伸缩因子;τ是平移因子;函数ψa,τ(t)为小波基函数,其逆变换为
(2)
离散小波变换(DWT)是通过对CWT进行二进离散化实现的,可以将信号分解低频与高频2个部分到离散小波基的方法一般是将平移因子a和伸缩因子τ变成幂级数结构[13],即
(3)
式中:a0≠1;τ0是常数;j,k∈Z对应的离散小波基为
(4)
最终f(t)的离散小波变换为
(5)
2-D DWT就是将二维数据分别变成行数据与列数据后再进行一维离散小波变换。二维离散小波变换主要通过Mallat前向金字塔算法的扩展实现,二维小波变换原理如图2所示。每进行一次一维小波变换,原始数据就会分成高频与低频2个分量,因此二维数据在进行一次二维小波变换后,原始数据将会被分成4个分量,分别是原始二维数据的低频分量LL和高频分量HL、LH、HH。
图2 二维小波变换原理图
其中HL包含表示原始数据在水平方向上的高频分量,LH表示原始数据在垂直方向上的高频分量,HH表示原始数据在对角处的高频分量,这些高频分量表示的是原始信号的边缘信息,而LL表示原始数据的低频分量,是原始数据的中心部分,代表了原始数据的主要信息。而为了能够获得更加精确的小波系数,可以将第一次二维小波变换后获得的近似值(低频分量)LL作为下一次二维小波变换的输入数据,最多可进行二维小波变换p=log2N次,N为原始数据的长度。在实际应用中,分解的级数也与小波基的选取有关。
本研究使用了db5小波作为小波基,并进行了2次2-D DWT,最后又使用了小波逆变换重构出了新的亮温数据。之所以使用db5小波作为小波基是因为db5小波具有较好的正则性,且消失矩较大,比较光滑,具有较好的频域局化能力,同时还不会使得时域的紧支撑性过低从而增大计算量使得实时性变差。而进行2次二维小波变换的原因是为了在不失真的情况下获得更准确的小波系数。其中二维小波的重构过程为:先对二维小波变换后获得的高低频分量的每一列进行一维小波逆变换,再对逆变换后的结果的每一行进行一维小波逆变换,即可得到重构的数据。
本研究的技术路线如图3所示。
图3 技术路线
由于没有实地观测的分辨率为9 km亮温数据,无法对降尺度后的L波段的亮温数据进行直接有效的验证。因此用单通道检索算法(single-channel algorithm,SCA)和地表辅助数据[14]。该算法主要是利用对土壤水分变换敏感的亮温和辅助数据(地表温度、地表土壤质地、地表覆盖情况、植被含水量)来实现对土壤水分的反演。
SCA算法主要是基于零阶近似模型ω-τ构建的,其模型为
(6)
其实现过程具体分为以下5个步骤:
(7)
② 通过地表覆盖分类图确定对应的植被类型和参数b。
③ 根据地表覆盖分类图确定的植被类型获取植被单次反照率ω和植被光学厚度τ,通过公式(8)、(9)计算出植被透过率γP和植被含水量VWC,用来校正植被影响
(8)
τ=b·VWC。
(9)
④ 在将植被的影响去除后,使用Choudhury模型[15]对土壤粗糙度进行校正,获取光滑表面的发射率。
⑤ 最后用菲涅尔方程计算土壤介电常数后,使用Dobson模型[16]估算出土壤水分含量。
该算法因为使用了大量的辅助数据来去除土壤温度、大气温度、植被和地表粗糙度的影响,所以反演得到的土壤水分数据精度很高,同时也是目前SMAP卫星的官方的土壤水分反演算法。在本研究中的土壤水分反演用到的地表辅助数据均由SMAP数据集得到,并根据实验区域内的植被分布选取单次散射反照率ω和植被参数b,且b=0.13,ω=0.05。在将亮温数据反演得到土壤水分数据后,就可以与实地探测土壤水分数据进行对比,根据土壤水分数据的对比的结果就可以确定降尺度后亮温数据的精度。
本次实验为了在获取分辨率为9 km的L波段亮温数据的同时证明加入辅助数据对L波段亮温数据降尺度处理的重要性,以及验证使用2-D DWT对亮温降尺度的有效性,因此计划借助AMSR2的X波段亮温数据作为辅助数据,将SMAP 36 km的亮温数据降尺度到9 km,并与SMAP 9 km的亮温数据进行比较。由于AMSR2的X波段只有10 km的亮温数据,因此需要使用线性插值的方法将其改变为9 km的亮温数据,以此满足小波变换降尺度方法的要求。
本次研究选取的时间段为2016-01-01—12-31日,共计366 d;但是,由于SMAP和AMSR2的运行周期不一样,因此为了保证在实验区域内SMAP和AMSR2都能提供亮温数据,在上午时间段选取了55 d作为研究样本。
本研究中使用了二维离散小波变换对L波段亮温进行了降尺度,其降尺度的亮温图像如图4所示,显示了在2016-06-22—09-26,期间,以32 d为间隔选取了4 d的L波段亮温在上午时间段V极化和H极化的降尺度结果。其中第一行是SMAP 36 km原始的L波段亮温图像,第二行是AMSR2的9 km X波段亮温图像,第三行是二维离散小波变换降尺度的L波段亮温图像, 第四行是为了验证算法加入的SMAP 9 km的L波段亮温图像。通过比较可以发现,第三行的降尺度图像基本上可以呈现出SMAP 36 km亮温图像的空间分布情况与AMSR2 9 km亮温图像的空间结构,并与SMAP 9 km的亮温图像相比,二维离散小波变换的降尺度图像的异质性更高,与原始观测的36 km亮温图像的关系更紧密,空间上的细节更多。
(a)V极化
其亮温降尺度的图像分析结果见表1,分别计算了SMAP 9 km的亮温图像与二维离散小波变换的降尺度亮温图像在不同极化下的信息熵、平均梯度、标准差、空间频率。结果表明通过2-D DWT得来的降尺度图像在不同极化的条件计算得到的标准差和空间频率均高于SMAP 9 km的亮温图像,说明通过2-D DWT降尺度得到的亮温图像的有用信息更多,图像空间的总体活跃度更高。而在信息熵和平均梯度的计算上2-D DWT图像和SMAP 9 km图像各有优劣,总体持平,说明这2种亮温图像在平均信息量和图像对比度上有着近似的特征。这些情况表明二维离散小波变换的亮温图像有着更多的细节表达能力,图像信息更多,图像空间的活跃程度更高。
表1 亮温降尺度图像分析结果
为了验证降尺度的结果,需要利用SCA算法与SMAP提供的辅助数据,对2-D DWT降尺度后的9 km亮温数据和SMAP提供的9 km亮温数据进行表面土壤水分反演。
反演的土壤水分和站点土壤水分的长时间序列图与亮温反演土壤水分的散点图如图5所示。其中区域A和区域B分别是选取的2个不同的观测区域。为了确保土壤水分反演的精度,只反演了地表温度大于0 ℃的时期。从长时间序列图中可知这2种不同的亮温反演的土壤水分与站点观测的土壤水分有很高的相似度,特别是V极化下反演的土壤水分相似度更高。从散点图可知,用降尺度后的亮温反演的土壤水分和用SMAP 9 km的亮温反演的土壤水分极为相似,相关性最高为0.984 8,最低为0.970 3,表明两者之间相似度极高,而且偏差值最高仅为0.005 6,表明两者差值极低。其均方根误差最高为0.008 3 cm3/cm3,表明两者间的离散度非常低。以上情况说明用2-D DWT进行降尺度的亮温反演的土壤水分与SMAP 9 km亮温反演的土壤水分有着极高的相似度。
(a)V极化
反演土壤水分与站点土壤水分的对比统计结果见表2。将2种亮温在不同区域、不同极化的情况下反演的土壤水与相对应的站点观测的土壤水进行了对比,计算了反演的土壤水分与站点土壤水分的相关系数(R)、偏差值(BIAS)、均方根误差(RMSE)、无偏均方根误差(ubRMSE)。比较结果表明,反演的土壤水分偏差值均在0.2以内,且相关系数均大于0.68,均方根误差在0.2 cm3/cm3以内,无偏均方根误差最高仅为0.021 9 cm3/cm3,表明反演的土壤水分与站点观测的土壤水分的有着相似的性能,且在V极化情况下能达到最优。
表2 反演土壤水分与站点土壤水分对比统计结果
本文通过二维离散小波变换借助X波段的高分辨率亮温对L波段亮温进行降尺度处理,并与SMAP 9 km的亮温数据进行了对比。其结果表明通过二维离散小波变换得来的亮温图像在整体上要优于SMAP 9 km的亮温图像,且与SMAP 36 km的亮温图像的联系更加紧密,并且在反演的土壤水分的结果对比中两者结果十分近似,且都与站点观测的土壤水分有很高的相关度,无偏均方根误差都在0.02左右。这些数据说明在对L波段亮温进行降尺度时,可以使用二维离散小波变换对L波段亮温进行降尺度处理,而且精度与SMAP 9 km的亮温数据相比略有提高。本文的研究只限于帕里区域,且在反演土壤水分时只在2个像元上进行土壤水分反演结果验证,因此今后的工作需要将该方法改进,争取应用到全球区域的亮温降尺度,为亮温降尺度方法提供一种新的参考思路。