CN104236716A

CN104236716A - 基于成对hj-1b影像时空信息的地表温度反演方法

Info

Publication number: CN104236716A
Application number: CN201310238628.0A
Authority: CN
Inventors: 陈�峰; 赵小锋; 全元
Original assignee: Institute of Urban Environment of CAS
Current assignee: Institute of Urban Environment of CAS
Priority date: 2013-06-17
Filing date: 2013-06-17
Publication date: 2014-12-24

Abstract

本发明公开了一种基于成对HJ-1B影像时空信息的地表温度反演方法。该方法是针对现有的单窗（或单通道）算法在应用中所存在的参数获取困难而提出的一种新策略。实现该方法的步骤为：S1：获取研究区两个成对时相的HJ-1B影像，并进行相应的预处理；S2：对预处理后的影像以对称的方式进行边界延拓处理；S3：识别相似像元及其数目，确定相应的反演方法；S4：依据S3的判定结果，构建辐射传输方程组，并求解地表辐射值；S5：循环进行S3~S4,得的整个研究区在两个时相的地表辐射值，进行异常值识别和处理；S6：由处理后的地表辐射值反演地表温度。本发明能在相关参数无法准确获取的情况下准确反演地表温度，将有效促进HJ-1B热红外数据的定量应用。

Description

基于成对HJ-1B影像时空信息的地表温度反演方法

技术领域

本发明涉及热红外遥感技术领域，更具体的是在地表比辐射率和大气参数均未知的情况下，利用两个成对时相HJ-1B影像的时空信息，从单一热红外波段的辐射亮度直接定量反演地表温度的方法。

背景技术

地表温度是影响生态系统的关键因子之一，它与许多地表过程有关。热红外卫星遥感技术被认为是现阶段获取区域和全球尺度地表温度信息的一个可行且有效的手段。为了实现从热红外卫星遥感影像准确提取地表温度信息，国内外研究人员提出了一系列的反演方法，如单窗（或单通道）算法、劈窗算法、温度/比辐射率分离算法、白天/昼夜法。与其它方法相比，单窗（或单通道）算法的优势在于它能从单一热红外波段影像反演地表温度。因此，在实际应用中，对于Landsat TM/ETM+, CBERS, HJ-1B等仅具有单一热红外波段的传感器，单窗（或单通道）算法是唯一的选择。由辐射传输方程可知，为了从单一热红外波段影像准确地反演地表温度，需要同时对大气影响和地表比辐射率进行校正。已有的研究通过大量的模拟实验和简化假设提出了相应的经验算法模型及其改进形式（Qin 和Karnieli，2001；Jimenez-Munoz 和Sobrino，2003；张勇等，2006；段四波等，2008；Jimenez-Munoz 等，2009；Cristobal 等，2009）。

然而，这些经验的算法均以地表比辐射率和大气参数（如大气平均温度、大气水分含量）已知为前提。目前，对地表比辐射率的估算主要采用地表覆盖分类赋值法（Synder 等，1998）和植被指数法（Valor和Caselles，1996；覃志豪等，2004；Sobrino 等，2008）。基于地表覆盖分类的方法需要预先知道不同地类的比辐射率，并认为每个类别的比辐射率是固定不变的，这显然是不合理的；尽管基于植被指数的方法考虑了像元的组分差异，但依然需要预先确定代表性组分的比辐射率，如土壤，植被等。另外，受准确大气参数获取条件的限制，很多研究多假定整个研究区域的大气状态是均一分布的，在地表温度反演过程中不考虑大气的影响，如Artis和Carnahan（1982）在早期所提出的一个单窗算法就没有考虑大气的影响，因此利用该算法不能得到真实的地表温度。由于大气状态往往存在一定程度的空间异质性，假定大气状况均质，不进行大气校正，或以某一点代替整个区域必将给遥感地表温度反演带来误差（陈峰等，2010）。理论上，可以从MODIS的大气产品（如MOD 05）提取现有经验算法所需的大气参数，但是该产品自身的不确定性可能会在有些情况下给最终的地表温度结果带来不利影响。可见，由于在准确获取相关参数方面存在的困难，将在一定程度上限制现有的单窗（或单通道）算法实施，从而也将影响HJ-1B单一热红外波段数据的定量应用。

本发明针对HJ-1B星的高回访频率的特点，利用两个成对时相HJ-1B影像的时空信息，基于两个合理的假设，提出一种直接从单一热红外波段影像反演地表温度的方法，而无需过度依赖准确的地表比辐射率和大气参数信息。因此，该新方法可以视为对现有的单窗（或单通道）算法的有益补充。

发明内容

要解决的技术问题 本发明主要解决的技术问题是：在无法准确获取研究区地表比辐射率和大气参数（如大气水分含量，大气平均作用温度）的情况下，即无法有效实施现有的单窗（或单通道）算法时，实现从HJ-1B影像准确反演地表温度。

技术方案 为了解决上述技术问题，本发明提出了一种基于成对HJ-1B 影像时空信息的地表温度反演方法，包括以下步骤：

S1 获取研究区两个时相的HJ-1B影像（含CCD多光谱和IRS热红外影像），结合官方公布资料和影像头文件信息，相应地，计算CCD影像的各波段的星上反射率值和IRS热红外影像辐射亮度值；并对HJ-1B CCD反射率影像和热红外波段辐射影像进行地理校正等预处理；

S2 对预处理后的影像以对称的方式进行边界延拓处理；

S3 依据每个目标像元所在的滑动窗口范围内各像元的光谱相似度指数，识别相似像元及其数目，确定相应的反演方法；

S4 依据S3的结果，相应地提取相似或准相似像元的辐射值，构建一个简化的辐射传输方程组，通过一个带约束的非线性优化模型实现对该方程组的求解，得到目标像元在两个成对时相的地表辐射值；

S5 对延拓后的影像循环进行S3~S4, 得的整个研究区两个时相的地表辐射值，并对其进行异常值识别和处理；

S6 利用普朗克反方程或其简化模型从处理后的地表辐射值反演地表温度，并同时保存反演方法的标识信息。

有益效果 本发明可以在地表比辐射率和大气参数均无法准确获取的情况下，仅利用两个成对时相HJ-1B影像的时空信息实现对地表温度的准确反演，克服了现有的方法对辅助信息过于依赖的问题，进而能在一定程度上促进单一热红外波段影像的定量应用。

附图说明

图1是一种基于成对时相HJ-1B 影像时空信息的地表温度反演方法示意图。

图2是基于本发明的两个时相HJ-1B地表温度反演结果图。

具体实施方式

以下结合附图和具体实施例，对本发明提供的基于成对HJ-1B 影像时空信息的地表温度反演方法进行详细描述。以下实施例仅用于说明本发明，而不是用来限制本发明的范围。如图1所示，本发明基于成对HJ-1B影像时空信息的地表温度反演方法，主要包括以下步骤：

S1 获取研究区影像，并进行辐射校正、地理校正、重采样和裁剪预处理

在本实施例中，获取时相相近的两景HJ-1B IRS影像（获取时相分别为2011年12月18日和2011年12月22（path/row：452/91），后文分别记为时相T18和T22），并以珠江三角洲的主要区域为实验区。由于数据质量的影响，仅有与T18相同时相的一景HJ-1B CCD2多光谱影像（path/row: 454/90）是可用的（均来自http://www.cresda.cn）。因此，本实施例中假定两个时相的地表状态相同，则可用T18时相的CCD影像充当T22时相的影像，而不对最终的结果产生影响。一景Landsat5 TM影像（获取时相为2009年1月2日，path/row:122/44），作为地理校正的参考影像。将HJ-1B CCD影像和IRS热红外波段的DN值转换为辐射值，并进一步将CCD辐射值影像转换为星上反射率值。以Landsat5 TM 为参考影像，在ENVI软件平台上，通过人工目视解译和对比判别，在Landsat5 TM影像、以上各HJ-1B CCD星上反射率和IRS热红外辐射影像选取15个均匀分布的特征同名点，并采用二次多项式模型对HJ-1B CCD星上反射率和IRS热红外辐射影像进行地理校正。由于HJ-1B CCD影像和IRS热红外影像存在空间分辨率差异，这里将HJ-1B CCD星上反射率（30m×30m）结果重采样到300m×300m，并裁剪确定最终的研究区影像；

S2对预处理后的影像进行边界延拓处理

为了得到研究区完整的反演结果，采用对称的方式对S1中预处理后的CCD星上反射率影像和IRS热红外辐射亮度值影像均进行边界延拓处理，在本实施例中延拓宽度取值为3个像元；

S3 识别相似像元及其数目，确定相应的反演方法

依据每个目标像元所在的滑动窗口范围（本实施例中定义为7×7像元的范围）内各像元的光谱相似度指数，识别相似像元及其数目，确定相应的反演方法。本发明中采用RMSD指数（Chen 等，2011）来度量某个像元与目标像元的光谱相似度，RMSD指数通过下式计算：

其中，为滑动窗口内像元 P_j在 t 时相CCD影像 b 波段的光谱星上反射率值；同样地，为滑动窗口内目标像元 P₀在 t 时相CCD影像 b 波段的光谱星上反射率值，n为波段数。通过与阈值进行比较来确定像元P_j是否为目标像元的相似像元，当滑动窗口内像元P_j的RMSD值小于阈值时，该像元被判定为与目标像元是相似的。确定目标像元的相似像元后，统计滑动窗口内目标像元的相似像元数N，并依据N确定相应的求解策略。具体而言，若N＜4，则分别对两个时相的窗口内所有像元的RMSD值进行升序排列，并将两个时相的排序号相加，最终确定前面4个具有最小排序号的为准相似像元，相应的反演方法标记为‘M0’；若N≥4，则所有N个相似像元均参与后面的求解过程，相应的反演方法标记为‘M1’；

S4求解目标像元在两个时相的地表辐射亮度值

提取相似或准相似像元的辐射值，构建一个简化的辐射传输方程组，通过一个带约束的非线性优化模型实现对方程组的求解，得到目标像元在两个时相的地表辐射亮度值。对每个目标像元，在t₁和t₂时相，提取该像元及其N个相似或准相似像元的辐射亮度值，构成如下的简化辐射方程式：

该辐射方程组包含2(N+1)个方程式和N+6个未知参数的方程组，即N+1个地表温度值，2个大气综合参数和1个地表温度变化量。由此，当相似像元数N≥4时，该方程组有解；在实际求解过程中，我们先求解出地表辐射量。本发明中构建一个带约束的非线性优化模型来求解上述辐射传输方程组，得到目标像元在两个时相的地表辐射值。优化模型的目标函数定义如下：

其中，；为像元P_i的权重，为了减小目标像元P₀的模型估计误差，将其权重设置为一个较大的值（1000），其他相似像元的权重则设置为1，即；

S5 对地表辐射估算结果中的异常值识别和修正

依据本发明的假设，滑动窗口内相似像元（或准相似像元）具有相同的地表比辐射率，且受到相同的大气影响，则该空间范围内各相似像元间的辐射观测量差异仅与地表温度的差异有关，从而判定某目标像元地表辐射估计值是否异常。以t₁时相的目标像元P₀为例，目标像元辐射亮度值与其相似像元平均辐射亮度值之差可表述为：

其中，为相似像元平均辐射亮度值；

因此，若不满足式，则判定目标像元的地表辐射值为异常值，此时目标像元P₀的地表辐射值用其相似像元的平均地表辐射值来替代；反之，则原估算值被认为是合理的；

S6 从修正后的地表辐射值反演地表温度

依据普朗克反方程可以从地表辐射值来计算地表温度。而针对HJ-1B IRS热红外波段的响应函数，以准确提取温度为前提，本发明中采用Hu 等（2011）的研究结果。此时，地表温温度可通过如下经验关系由地表辐射量求得：

其中，。

为了说明本发明—基于成对时相HJ-1B 影像时空信息的地表温度反演方法的实用性及其反演地表温的精度，本实施例中采用MODIS标准地表温度/比辐射率产品（MOD11）以及海洋表面温度产品来验证反演结果。由于MODIS Terra星与HJ-1B过境时间接近（Terra星为上午10:30左右，HJ-1B为上午10:45左右），理论上二者获得的陆表/海表温度值应该相同（近）。本实施例中采用3个常用的指标来对地表温度的反演结果进行精度评价，即系统误差（SE）、均方根误差（RMSE）和中误差（MAE）。误差评价结果显示（表1）：与MODIS陆表标准差产品相比，本实验反演得到的陆表温度整体上偏高约0.5K~1K，绝对偏差的中值接近1K；同时，与MODIS海表温度产品相比，我们反演得到的海表温度整体上偏高约0.5K，绝对偏差的中值接近0.5K，且均方根误差可控制在1K以内；另外，两个时相的误差评价结果接近。

表1 总体误差指标评价结果。

对本实施例中两个时相HJ-1B地表温度反演结果的误差分析表明，本发明提出的基于两个成对时相HJ-1B 影像时空信息的地表温度反演方法能在没有地表比辐射率和大气参数信息的情况下准确地反演地表温度，可作为现有的单窗（单通道）算法的有益补充，从而能在一定程度上促进HJ-1B热红外数据的定量应用。本发明同样也适用于其他具有单一热红外波段的传感器，如Landsat5 TM/Landsat7 ETM+、CBERS。

以上实施案例仅用于清晰说明本发明，而并非对本发明的限制。对所属领域的技术人员，在不脱离本发明的精神和范围的情况下，在上述说明的基础上还可以做出其他不同形式的变化或变动，因此所有等同的技术方案或方法也属于本发明的范畴，本发明的专利保护范围应有权利限定。

Claims

1.一种基于成对HJ-1B 影像时空信息的地表温度反演方法，其特征在于包括以下步骤：

S1: 获取研究区两个成对时相的HJ-1B影像（含CCD多光谱和IRS热红外影像），结合官方公布资料和影像头文件信息，相应地，计算CCD影像的各波段的星上反射率值和IRS热红外影像辐射亮度值；并对HJ-1B CCD反射率影像和热红外波段辐射影像进行地理校正等预处理；

S2: 对预处理后的影像以对称的方式进行边界延拓处理；

S3: 依据每个目标像元所在的滑动窗口范围内各像元的光谱相似度指数，识别相似像元及其数目，确定相应的反演方法；

S4: 依据S3的结果，相应地提取相似或准相似像元的辐射值，构建一个简化的辐射传输方程组，通过一个带约束的非线性优化模型实现对方程组的求解，得到目标像元在两个时相的地表辐射估计值；

S5: 对延拓后的影像循环进行S3~S4, 得的整个研究区两个时相的地表辐射值，并对其进行异常值识别和处理；

S6: 利用普朗克反方程或其简化模型从处理后的地表辐射值反演地表温度，并同时保存与反演方法相关的标识信息。

2.如权利要求1所述的基于成对HJ-1B影像时空信息的地表温度反演方法, 其特征在于利用了两个成对HJ-1B影像的时空信息，输入影像包括：HJ-1B CCD多光谱星上反射率影像和HJ-1B IRS热红外辐射亮度值影像，所述步骤S1中需分别将CCD多光谱和IRS热红外的原始DN值影像分别转换为星上反射率和辐射亮度，并进行地理校正、CCD影像重采样和研究区裁剪处理；由于目前HJ-1B卫星在地理定位方面存在的偏差问题，因此采用与HJ-1B影像相近时相的Landsat TM/ ETM+作为参考影像，进一步对HJ-1B CCD星上反射率影像和IRS热红外波段辐射亮度值影像进行地理校正；为了解决HJ-1B CCD影像和IRS热红外影像在空间分辨率方面存在的差异，这里将HJ-1B CCD星上反射率（30m×30m）结果重采样到300m×300m，与热红外影像同空间分辨率，并裁剪影像确定最终的研究区影像。

3.如权利要求1所述的基于成对HJ-1B影像时空信息的地表温度反演方法, 其特征在于，所述步骤S2中，为了得到研究区完整的反演结果，采用对称的方式对S1中预处理后的CCD星上反射率影像和IRS热红外辐射亮度值影像均进行边界延拓处理，延拓宽度为像元大小（为反演方法中所采用滑动窗口的大小, 如3，5，7，9，…；鉴于HJ-1B IRS热红外影像的300m×300m的空间分辨率，结合我们以往的研究，我们建议取7个像元，这样延拓宽度为3个像元）。

4.如权利要求1所述的基于成对时相HJ-1B影像时空信息的地表温度反演方法, 其特征在于，所述步骤S3中，定义了滑动窗口范围内的相似或准相似像元，并认为它们具有相同的光谱和热辐射特性，进而在每个像元所在的滑动窗口内选取合理的像元参与对目标像元地表辐射值的优化求解；具体而言，采用了RMSD指数来度量某个像元与目标像元的相似程度，通过与阈值比较来确定该像元是否为相似像元。

5.如权利要求1所述的基于成对HJ-1B影像时空信息的地表温度反演方法, 其特征在于，所述步骤S3中，统计滑动窗口内目标像元的相似像元数N，并依据N确定相应的求解策略。

6.如权利要求1所述的基于成对时相HJ-1B影像时空信息的地表温度反演方法, 其特征在于，本发明中的两个合理假设：（a）在正常天气状况下，短时间内在同一空间位置，地表比辐射率不变；（b）小空间范围内（本发明中用滑动窗口来表示）的大气条件相同或差异可以忽略；并且某一目标像元，在一定空间范围内存在N个与其类别相同（或相似）的像元（即相似像元）；在具体的实施中，本发明的特征还在于，所述步骤S4中，利用目标像元和所选取的相似或准相似像元的辐射亮度值信息，构建一个简化的辐射传输方程组，通过建立一个带约束的非线性优化模型实现对该方程组的求解，得到目标像元在两个成对时相的地表辐射值。

7.如权利要求1所述的基于成对时相HJ-1B影像时空信息的地表温度反演方法, 其特征在于，所述步骤S5中，依据滑动窗口内相似像元（或准相似像元）具有相同的地表比辐射率，且受到的大气影响相同，则该空间范围内各相似像元间的辐射亮度观测量差异仅与地表温度的差异有关，基于这一原则可判定，对某目标像元的地表辐射的估计是否异常。