一种基于微波遥感与时空信息的云下像元LST估算方法与流程

本发明涉及MODIS云下像元LST估算方法
技术领域：
，尤其涉及一种基于微波遥感与时空信息的云下像元LST估算方法。
背景技术：
：地表温度(LandSurfaceTemperature,LST)是衡量地球表面水热平衡的关键参数，对于气候、水文、地球物理等科学领域的相关研究具有重要意义。20世纪70年代以来，国内外学者就如何利用热红外遥感获取地表温度开展了大量研究，提出了包括劈窗算法在内的多种反演算法。尽管热红外遥感覆盖范围广、实时性强、能够快速获取地表比辐射率和温度信息，但其易受大气水汽影响，无法透过云层，因此不能进行地表温度的全天候监测。当前被广泛使用的MODIS地表温度产品，部分影像缺值率高达60％，严重影响了产品使用效率，实现对云覆盖等像元地表温度的重建具有重要意义。而被动微波遥感受大气干扰小，可穿透云层获取地表辐射信息，并具有全天候、多极化特点，在地表温度反演中具有独特的优越性，为重构MODIS等地表温度产品提供了可能。探讨地表亮温与地表温度的相关性，建立二者之间的关系模型反演地表温度成为当前研究的热点和重点，目前被广泛应用的微波传感器主要有三种：搭载在Nimbus-7卫星上的SMMR(ScanningMultichannelMicrowaveRadiometer)传感器、搭载在美国国防气象卫星计划(DMSP)的系列极轨卫星上的SSM/I(SpecialSensorMicrowaveImager)传感器以及地球观测系统(EOS)的Terra和Aqua卫星所搭载的AMSR_E(AdvancedMicrowaveScanningRadiometer-EOS)传感器。其中，SMMR传感器在轨运行时间为1978～1987年，SSM/I在轨运行时间为1987～2005年，AMSR_E传感器从2002年服务至今，且其空间分辨率高于SMMR和SSM/I的空间分辨率，该传感器的亮温产品是当前用于反演地表温度的首选数据源；当前基于被动微波遥感反演地表温度的算法主要包括三类，分别是统计模型法、物理模型法和神经网络算法。其中，物理模型法以辐射传输方程为理论基础，具有明确的物理意义，不受时空限制，但模型参数多且复杂，大多采用条件假设的方法进行简化，影响反演精度；而神经网络算法可以并行求解非线性问题，不需要精确地刻画反演方程，避免了对物理参数中不确定因子的模拟，但该算法只是从数学层面上弱化了地表温度病态反演问题，缺乏实际物理意义；统计模型法则是利用亮温与地表温度之间的内在关系、通过大量数据回归总结规律，算法机理简单，便于实现，已成为当前应用最为广泛的基于被动微波遥感反演地表温度的方法。根据参与拟合的亮温通道数量，统计模型法又可以分为单通道模型和多通道模型，其中多通道模型的反演效果远远好于单通道的反演效果，其算法思想是：基于地表类型实现像元分类，按类别建立LST与地表亮温的回归模型，按照无LST值像元的地表属性选择合适的回归模型估算地表温度值。然而，该算法在实际应用中以现有地表类型产品(如MOD12)为基准，反演精度低，平均误差在3K左右；反演结果明显受环境条件和气象条件限制，缺值率越高，反演误差越大；反演区域以像元为单位，不能实现大区域的地表温度重建，算法的实用性低。技术实现要素：基于
背景技术：
存在的技术问题，本发明提出了一种基于微波遥感与时空信息的云下像元LST估算方法。本发明提出的一种基于微波遥感与时空信息的云下像元LST估算方法，包括如下步骤：S1，地表覆盖分类：Landsat/TM遥感数据集经辐射校正、大气校正、重采样等预处理后，利用ENVI波段运算功能计算其NDVI，公式如下：NDVI＝(b4-b3)/(b4+b3)1式中，b3、b4分别为第三通道和第四通道的反射率。根据NDVI值，将地表覆盖分为6大类，分别是：NDVI＜0；0≤NDVI＜0.1；0.1≤NDVI＜0.3；0.3≤NDVI＜0.5；0.5≤NDVI＜0.7；NDVI≥0.7；S2，像元分割：选用18.7、23.8、36.5与89GHz(H/V)8个通道参与地表温度反演，利用ENVI软件对18.7、23.8和36.5GHz(H/V)数据进行像元分割，使其与89GHz像元大小一致，在LST重建过程中，利用坐标约束实现地表亮温与地表温度的空间匹配；S3，MODISLST与AMSR_E亮度温度回归分析：微波辐射计在给定频段下接收到的辐射亮温表示为：Tf＝τfεfTs+Δ2式中，Tf为亮温，Ts为地表温度，εf为地表辐射率，τf为大气透过率，Δ为其他改进项，包括大气上行辐射和下行辐射；上式简化为：Tf＝εfTs3S4，建立地表温度与地表亮温的多元回归模型：MODIS地表温度与AMSR_E地表亮温的多元回归模型如公式3所示：LSTn×1MODIS=Tn×8AMSR_E×A8×1+Bn×1---4]]>LST=lst1lst2...lstn-1lstnA=a1a2...a7a8B=b1b2...bn-1bn]]>T=t118.7Ht118.7Vt123.8Ht123.8Vt136.5Ht136.5Vt189Ht189Vt218.7Ht218.7Vt223.8Ht223.8Vt236.5Ht236.5Vt289Ht289V.........................tn-118.7Htn-118.7Vtn-123.8Htn-123.8Vtn-136.5Htn-136.5Vtn-189Htn-189Vtn18.7Htn18.7Vtn23.8Htn23.8Vtn36.5Htn36.5Vtn89Htn89V]]>式中，LST为MODISLST列向量，T为同期AMSR_E各通道地表亮温均值矩阵，A、B为待估参数，n为MODIS参与反演的像元数量；S5，基于加权移动平均滤波的估计值改正：设自变量x以步长h做等距观测，对应观测结果y，即：xi＝x0+ih5表1观测序列说明对公式5做变换，得到：t=x-xih---6]]>则：表2观测序列变形构造平滑公式如下：yi'+t＝Ao+A1t+A2t2+...+Amtm7平滑多项式的系数由最小二乘原理确定，即：min(δ)=min[Σt(yi+t′-yi+t)2]---8]]>即：min(δ)=minΣt((A0+A1t+A2t2+...+Amtm)-yi+t)2---9]]>式中，t取最靠近i的2n+1个整数值(即平滑点数)，即t＝-n,-n+1,...0,...n-1,n，且应保证m＜2n+1＜N，根据求解各系数Ai，根据最小二乘原理：∂δ∂A0=Σt=-22(yi+t-A0-A1t-A2t2)=0∂δ∂A1=Σt=-22(yi+t-A0-A1t-A2t2)t=0∂δ∂A2=Σt=-22(yi+t-A0-A1t-A2t2)t2=0---10;]]>通过分析，采用五点加权二次平滑滤波法对含有初始估计值的MODIS地表温度时间序列进行重构，根据公式(10)求解权值系数，将其作用到MODISLST时间序列中，得到公式(11)：y(i)=(-3×y0(i-2)+12×y0(i-1)+17×y0(i)+12×y0(i+1)-3×y0(i+2))/35(3≤i≤44)y(1)=(31×y0(1)+9×y0(2)-3×y0(3)-5×y0(4)+3×y0(5))/35y(2)=(9×y0(1)+13×y0(2)+12×y0(3)+6×y0(4)-5×y0(5))/35y(45)=(-5×y0(42)+6×y0(43)+12×y0(44)+13×y0(45)+9×y0(46))/35y(46)=(3×y0(42)-5×y0(43)-3×y0(44)+9×y0(45)+31×y0(46))/35---11]]>式中，y(i)为第i期LST重构值，y0(i)为含有初步估计值的第i期温度值；S6，LST重建精度评定：计算各区域LST重建均误差，公式如下：ME=Σi=1n|Pi-Oi|/n---12]]>式中，ME为均误差，Pi为第i个像元LST重建值，Oi为原始值，n统计像元量。优选地，所述S1中的1式中，b3、b4分别为第三通道和第四通道的反射率。优选地，所述S5中的10式中，采用的数值为，m＝2,2n+1＝5，t＝-2,-1,0,1,2；公式11为针对MODISLST8Day产品一年数据的计算公式。本发明中声明被动微波遥感反演地表温度的算法主要包括三类，分别是统计模型法、物理模型法和神经网络算法。其中，物理模型法以辐射传输方程为理论基础，具有明确的物理意义，不受时空限制，但模型参数多且复杂，大多采用条件假设的方法进行简化，影响反演精度；而神经网络算法可以并行求解非线性问题，不需要精确地刻画反演方程，避免了对物理参数中不确定因子的模拟，但该算法只是从数学层面上弱化了地表温度病态反演问题，缺乏实际物理意义；统计模型法则是利用亮温与地表温度之间的内在关系、通过大量数据回归总结规律，算法机理简单，便于实现，已成为当前应用最为广泛的基于被动微波遥感反演地表温度的方法。统计模型法分为单通道模型和多通道模型，其算法思想是基于地表类型实现像元分类，按类别建立LST与地表亮温的回归模型估算缺失的地表温度值。然而，该算法在实际应用中以现有地表类型产品(如MOD12)为基准，反演精度低，平均误差在3K左右；反演结果明显受环境条件和气象条件限制，缺值率越高，反演误差越大；反演区域以像元为单位，研究区面积越大误差越大，不能实现大区域的地表温度重建，算法的实用性低。本发明引入时间序列滤波，利用前后时段对无LST时段的影响改正统计模型估计值，通过补充时间信息减少反演结果对空间信息的依赖，提高算法本身的空间适用度。当前公开发布的全球地表覆盖卫星产品误差较大，在进行地表覆盖分类时容易将不同属性或属性相差较大地物错分为一类，影响模型的拟合，降低反演精度。本发明提出利用NDVI(归一化植被指数)作为地物分类的基准。植被指数用定量化的手段描述植被覆盖度，不同地物类型的NDVI具有较高的可区分度，已作为一种遥感手段广泛应用于土地利用覆盖探测、植被覆盖密度评价、作物识别和作物预报等方面。用于解释遥感数据的植被指数有多种，如简单植被指数(LCI)、比值植被指数(RVI)、归一化植被指数(NDVI)、绿度植被指数(GVI)和土壤调节植被指数(SAVI)等，其中最常用的是归一化植被指数(NDVI)，本发明以NDVI作为地物分类的参考标准进行地表覆盖分类，建立基于下垫面类型的、刻画地表温度与地表亮温相关关系的多元回归模型，实现对无LST像元的地表温度值的初估计；考虑地表温度的周期性特征，引入时间序列滤波，利用前后时段对无LST时段的影响改正统计模型的估计值，提高反演精度。附图说明图1为本发明提出的一种基于微波遥感与时空信息的云下像元LST估算方法的LST重建流程图；图2为本发明提出的一种基于微波遥感与时空信息的云下像元LST估算方法的LST时间序列多项式拟合示例图；图3为本发明提出的一种基于微波遥感与时空信息的云下像元LST估算方法的基于不同时间窗口二次加权滤波算法的LST重建结果对比图。具体实施方式下面结合具体实施例对本发明作进一步解说。参照图1-3，本发明提出的一种基于微波遥感与时空信息的云下像元LST估算方法，包括如下步骤：S1，Landsat/TM遥感数据集经辐射校正、大气校正、重采样等预处理后，利用ENVI波段运算功能计算其NDVI，公式如下：NDVI＝(b4-b3)/(b4+b3)1式中，b3、b4分别为第三通道和第四通道的反射率。根据NDVI值，本发明将地表覆盖分为6大类，分别是：NDVI＜0；0≤NDVI＜0.1；0.1≤NDVI＜0.3；0.3≤NDVI＜0.5；0.5≤NDVI＜0.7；NDVI≥0.7；S2，像元分割：研究表明37GHz受水汽、云、雨等散射和吸收影响较小，适宜用来反演地表温度；此外，37GHz与22GHz通道组合有利于消除大气水汽影响，37GHz与19GHz通道组合有利于消除土壤水分影响，而89GHz与MODIS地表温度相关性最高。因此，本发明选用18.7、23.8、36.5与89GHz(H/V)8个通道参与地表温度反演。然而，不同频道的AMSR_E地表亮温数据空间分辨率差异较大，为方便计算，利用ENVI软件对18.7、23.8和36.5GHz(H/V)数据进行像元分割，使其与89GHz像元大小一致，在LST重建过程中，利用坐标约束实现地表亮温与地表温度的空间匹配；S3，MODISLST与AMSR_E亮度温度回归分析：根据辐射传输方程,微波辐射计在给定频段下接收到的辐射亮温可以表示为：Tf＝τfεfTs+Δ2式中，Tf为亮温，Ts为地表温度，εf为地表辐射率，τf为大气透过率，Δ为其他改进项，包括大气上行辐射和下行辐射等。忽略大气影响，上式可以简化为：Tf＝εfTs3众所周知，大部分地物的介电常数主要由其所含水分多少决定。利用被动微波亮温数据获取陆表温度就是在地表辐射率未知的情况下,通过不同频率不同极化的亮温组合,建立各通道亮温与地表温度之间的经验关系；S4，地表温度与地表亮温的多元回归模型：MODIS地表温度与AMSR_E地表亮温的多元回归模型如公式3所示：LSTn×1MODIS=Tn×8AMSR_E×A8×1+Bn×1---4]]>LST=lst1lst2...lstn-1lstnA=a1a2...a7a8B=b1b2...bn-1bn]]>T=t118.7Ht118.7Vt123.8Ht123.8Vt136.5Ht136.5Vt189Ht189Vt218.7Ht218.7Vt223.8Ht223.8Vt236.5Ht236.5Vt289Ht289V.........................tn-118.7Htn-118.7Vtn-123.8Htn-123.8Vtn-136.5Htn-136.5Vtn-189Htn-189Vtn18.7Htn18.7Vtn23.8Htn23.8Vtn36.5Htn36.5Vtn89Htn89V]]>式中，LST为MODISLST列向量，T为同期AMSR_E各通道地表亮温均值矩阵，A、B为待估参数，n为MODIS参与反演的像元数量；S5，基于加权移动平均滤波的估计值改正：众所周知，地表温度具有周期性变化特征，引入时间序列滤波，可利用短时间窗口内地表温度变化的一致性对上一步的估计值进行改正，即补充时间信息降低估计误差。一般来说，同一地区，时间越接近，地表温度值越相似。基于这一前提，本发明采用加权移动平均滤波法对含有初步估计值的LST时间序列加以处理；加权的基本思想是平均区间内中心数据的权值最大，越偏离中心处的数据权值越小，这样就在一定程度上减小了对真实信号本身的平滑作用。设自变量x以步长h做等距观测，对应观测结果y，即：xi＝x0+ih5表1观测序列说明对公式5做变换，得到：t=x-xih---6]]>则：表2观测序列变形构造平滑公式如下：yi'+t＝Ao+A1t+A2t2+...+Amtm7平滑多项式的系数由最小二乘原理确定，即：min(δ)=min[Σt(yi+t′-yi+t)2]---8]]>即：min(δ)=minΣt((A0+A1t+A2t2+...+Amtm)-yi+t)2---9]]>式中，t取最靠近i的2n+1个整数值(即平滑点数)，即t＝-n,-n+1,...0,...n-1,n，且应保证m＜2n+1＜N，根据求解各系数Ai，根据最小二乘原理：∂δ∂A0=Σt=-22(yi+t-A0-A1t-A2t2)=0∂δ∂A1=Σt=-22(yi+t-A0-A1t-A2t2)t=0∂δ∂A2=Σt=-22(yi+t-A0-A1t-A2t2)t2=0---10]]>分析地表温度时间序列性质，发现温度值与时间之间具有较高的二阶非线性相关特征，如图a所示(LST时间序列示例)。此外，平滑窗口越大，平滑效果越好，但信号失真也越大。根据分析，本文最终采用五点平滑窗口对含有初始估计值的MODIS地表温度时间序列进行滤波重构。根据公式10求解权值系数，将其作用到MODISLST时间序列中，以一年数据为例，得到公式(11)。y(i)=(-3×y0(i-2)+12×y0(i-1)+17×y0(i)+12×y0(i+1)-3×y0(i+2))/35(3≤i≤44)y(1)=(31×y0(1)+9×y0(2)-3×y0(3)-5×y0(4)+3×y0(5))/35y(2)=(9×y0(1)+13×y0(2)+12×y0(3)+6×y0(4)-5×y0(5))/35y(45)=(-5×y0(42)+6×y0(43)+12×y0(44)+13×y0(45)+9×y0(46))/35y(46)=(3×y0(42)-5×y0(43)-3×y0(44)+9×y0(45)+31×y0(46))/35---11]]>式中，y(i)为第i期LST重构值，y0(i)为含有初步估计值的第i期温度值。S6，LST重建精度评定：依次将每景影像的有值像元视作云下像元，采用本专利提出的重建算法对其地表温度值进行重建，与原始值对比，分别统计绝对误差在0～1K、0～2K和0～3K的像元比例，计算各区域LST重建均误差，公式如下。ME=Σi=1n|Pi-Oi|/n---12]]>式中，ME为均误差，Pi为第i个像元LST重建值，Oi为原始值，n统计像元量。由于该算法增加了时间维度信息，降低了反演结果对空间信息的依赖程度，重建过程受环境和气象条件的约束明显减少，重建效率高，重建结果满足地表温度的空间一致性。以MOD11A2地表温度产品为研究对象，以新疆、青海、四川、云南、河南、安徽、湖北、湖南、江西等共9个省份作为研究区，对这些区域的MODIS云下像元进行重建研究，并对比基于MODIS地表覆盖产品MOD12Q1的空间域多通道统计模型的反演结果，重建精度明显提高：基于MOD12Q1分类的统计模型平均反演精度为5.14K，而基于时空域联合算法的平均重建精度为1.24K，提高了76个百分点；因此，相较于以往的研究，该算法具有较高区域普适性和重建精度，能很好地实现大区域的云下像元LST重建，实用性强。以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本
技术领域：
的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。当前第1页1 2 3