CN106404178B - 基于小波变换的高光谱热红外地表温度与发射率分离方法 - Google Patents

Abstract

本发明公开了一种基于小波变换的高光谱热红外地表温度与发射率分离方法，由以下步骤组成：1)获取离地辐亮度R_surf；2)估计温度初值T₀；3)将T₀代入公式，计算得到地表发射率ε_s；4)对ε_s进行小波变换，获取重建后的地表发射率估值ε′_s；5)构建代价函数；6)将ε′_s与T₀带入公式计算此时的代价函数C；7)结合牛顿迭代法和退火算法迭代计算地表温度，并将其作为下一次迭代的温度初值T₀，重复步骤(3)至(7)，获取最终的地表温度和发射率。本发明克服了对发射率光谱形状假设的依赖，合理减少了反演中未知数个数，消除了遥感反演的病态性，实现了高光谱热红外地表温度和发射率的高精度遥感分离和准确获取。

Description

基于小波变换的高光谱热红外地表温度与发射率分离方法

技术领域

本发明涉及一种分离方法，尤其涉及基于小波变换的高光谱热红外地表温度与发射率分离方法。

背景技术

地表发射率和地表温度作为定量遥感领域的两个关键参数，已被广泛的应用于水文，生态，环境和生物地球化学等领域。近年来，越来越多的学者开始关注地表温度和发射率在各领域中的关键作用，相关的研究也在不断的开展。而这其中，关于如何通过遥感反演的方式精确的获取高光谱热红外地表温度与发射率一直是人们关注的研究热点之一。

目前，利用遥感方式进行热红外地表温度与发射率分离的最大的难点在于辐射传输方程病态性问题的求解。即便是在已获取精准大气信息的前提下，对于给定的N个通道的星上辐亮度，其可建立的方程组仍有N+1个未知数(N个通道发射率及1个地表温度)，不可解。针对这一问题，现有的高光谱热红外地表温度与发射率的分离算法从消除反演病态的出发点上来看，大致可分为有代表性的两类：一类是引入约束增加方程组个数的反演算法，如光谱平滑算法；另一类是减少未知数个数的反演算法，如分段线性算法。

光谱平滑算法是在地表发射率具有平滑光谱的前提下，引入了平滑度函数的概念，并以此建立约束方程，消除辐射传输方程的病态性，实现地表温度与发射率的分离。光谱平滑算法在一定程度上依赖于其平滑度函数，当实际发射率光谱无法满足平滑条件时，尤其是对于一些特定的人造材料，其算法会导致一定误差。此外，当大气等效温度与地表温度接近时，采用光谱平滑算法往往会存在奇异值的问题。

分段线性地表温度与发射率分离算法则是通过合理的设置波长区间，将区间内的地表发射率作为波长的函数，并用以线性的方式进行表示。通过这样的方式，实现了减少辐射传输方程中未知数个数的目的，使得方程可解。通常以分段线性的方式对地表发射率进行表示所带来的误差可以忽略，其精度与光谱平滑算法相当。但是，受限于分段区间设定的不唯一性，对于地表发射率波动明显的光谱曲线，过大的分段区间必然会导致光谱信息的丢失，而过小的分段区间则会损失线性子函数对发射率曲线的代表性，同时也会增加算法运算的复杂度。此外，分段线性算法将发射率光谱曲线进行了类似于离散化的处理(波长子区间)，这也导致了其无法保证各子区间的连续性，与实际发射率光谱连续的客观事实相违背。

发明内容

为了解决上述技术问题中的不足之处，本发明提供了一种基于小波变换的高光谱热红外地表温度与发射率分离方法。

根据维恩位移定律可以得知，地表的辐射能量主要集中在热红外波段，并且此波段内的太阳辐射能量很小，可忽略不计。因此，晴空条件下，在满足局地热平衡时的热红外波段的辐射传输方程可以表示为：

R_sensor＝(ε·B(T_s)+(1-ε)·R_down)·τ+R_up (1)

其中，R_sensor为卫星接收的星上辐亮度；ε为地表发射率；T_s为地表温度；B(T_s)为温度T_s的普朗克函数；R_down大气下行辐亮度；R_up为大气上行辐亮度；τ为大气透过率。

从式(1)可知，在已完成精准大气校正的条件下，当估算的地表温度T_s与真实地表温度存在差异时，其反演得到的地表发射率光谱曲线将会含有大气特征，尤其是在大气吸收线的位置，会出现明显的峰谷现象。从获取真实地表发射率的角度来看，此时在地表发射率光谱中的大气特征可以看做是一种大气噪声的表现。另一方面，由于卫星数据中往往会存在仪器的噪声影响，而这种影响也会直接体现在地表发射率的光谱曲线之上。也就是说，获取真实地表发射率可以看做是从含有噪声的观测信号中分离出地表温度和发射率信息的过程。基于此，本专利提出了一种基于小波变换的高光谱热红外地表温度与发射率分离算法。

为了解决上述技术问题，本发明采用的技术方案是：一种基于小波变换的高光谱热红外地表温度与发射率分离方法，主要由以下步骤组成：

1)获取离地辐亮度R_surf

通过辐射传输模型，在准确的配套大气温湿度剖面廓线数据的辅助下，估算大气上/下行辐亮度(R_up/R_down)和大气透过率(τ)；随后，利用辐射传输方程计算离地辐亮度R_surf，即

R_surf＝(R_sensor–R_up)/τ (2)

其中，R_surf为离地辐亮度；R_sensor为星上观测辐亮度；R_up为大气上行辐亮度；τ为大气透过率。

2)估计温度初值T₀

(a)首先，分别设定两个地表发射率的初值，ε₁和ε₂；

(b)然后，利用辐射传输方程分别计算两个地表发射率初值所对应的各波长位置的地表温度，并将其地表温度的通道温度最大值的均值作为初始温度T₀，即：

T₀＝(max(B^-1((R_surf-(1-ε₁)·R_down)/ε₁))+max(B^-1((R_surf-(1-ε₂)·R_down)/ε₂)))/2(3)

其中，T₀为温度初值；R_surf为离地辐亮度；R_down为大气下行辐亮度；ε₁和ε₂分别为地表发射率初值；B^-1表示普朗克反函数。

3)将T₀代入下式，计算得到地表发射率ε_s

ε_s＝(R_surf-R_down)/(B(T₀)-R_down) (4)

其中，R_surf为离地辐亮度；R_down为大气下行辐亮度；B表示普朗克函数；T₀为温度初值。

4)对ε_s进行小波变换，获取重建后的地表发射率估值ε′_s

(a)首先，对ε_s进行小波变换，获取相应的高频和低频小波系数：

[C_A,C_D]＝WT(ε_s,db) (5)

其中，C_A,C_D分别表示低频小波系数与高频小波系数；ε_s为地表发射率；db表示所采用的小波基；WT表示小波变换。

(b)然后，进行逆小波变换，将高频小波系数所对应的部分全部视为噪声予以排除，仅利用低频小波系数重建信号，获取重建后的地表发射率ε′_s：

ε′_s＝WT^-1(C_A,db) (6)

其中，WT^-1表示逆小波变换；C_A表示低频小波系数；db表示小波基。

5)构建代价函数

在已知地表温度T₀与重构后的发射率ε′_s后，我们可利用辐射传输方程估算其对应的星上辐亮度值R′_sensor:

R′_sensor＝(ε′_s·B(T₀)+(1-ε′_s)·R_down)·τ+R_up (7)

其中，R′_sensor为估算的星上辐亮度；ε′_s为重构的地表发射率；B(T₀)为温度T₀下的普朗克函数；R_down，R_up分别为大气下行及上行辐亮度；τ为大气透过率。此时，若地表温度与发射率接近真值，则其估算的星上辐亮度应与测量值近似相等。因此，可建立如下的代价函数：

C＝∑((R′_sensor-R_sensor)/mean(R_sensor))² (8)

其中，C为代价函数；R_sensor为星上观测辐亮度；R′_sensor为估算的星上辐亮度；当代价函数C达到最小时，其对应的地表温度T₀与发射率ε′_s便最为接近真实值。

6)将ε′_s与T₀带入式(8)，计算此时的代价函数C。

7)利用牛顿迭代法，并结合退火算法对地表温度进行扰动，迭代计算地表温度，并将其作为下一次迭代的温度初值T₀，重复步骤3)至步骤7)，直至满足代价函数达到最小值，结束迭代，获取最终的地表温度和发射率。

其中，步骤2)中ε₁和ε₂的取值分别为0.9和1.0。

本发明克服了对发射率光谱形状假设的依赖，通过小波变换将地表发射率从谱域信号转为频域信号，并通过获取的地表发射率的低频信息重构原始发射率光谱，这一过程合理减少了反演中未知数个数，消除了遥感反演的病态性，实现了高光谱热红外地表温度和发射率的高精度遥感分离和准确获取。

附图说明

下面结合附图和具体实施方式对本发明作进一步的详细说明。

图1是本发明方法的整体流程图。

具体实施方式

如图1所示，本发明主要由以下步骤组成：

1)获取离地辐亮度R_surf

通过辐射传输模型，在准确的配套大气温湿度剖面廓线数据的辅助下，估算大气上/下行辐亮度(R_up/R_down)和大气透过率(τ)；随后，利用辐射传输方程计算离地辐亮度R_surf，即

R_surf＝(R_sensor–R_up)/τ (2)

其中，R_surf为离地辐亮度；R_sensor为星上观测辐亮度；R_up为大气上行辐亮度；τ为大气透过率。

2)估计温度初值T₀

(a)首先，分别设定两个地表发射率初值，如ε₁＝0.9，ε₂＝1.0；

(b)然后，利用辐射传输方程分别计算两个地表发射率初值所对应的各波长位置的地表温度，并将其对应的波长位置的地表温度的通道温度最大值的均值作为初始温度T₀，即：

T₀＝(max(B^-1((R_surf-(1-ε₁)·R_down)/ε₁))+max(B^-1((R_surf-(1-ε₂)·R_down)/ε₂)))/2(3)

其中，T₀为温度初值；R_surf为离地辐亮度；R_down为大气下行辐亮度；ε₁＝0.9,ε₂＝1.0为地表发射率初值(对于典型地类，地表发射率在热红外谱段的取值范围通常介于0.9到1之间，经过多次试验表明，当其初值分别取值为0.9和1.0时，算法可具有较好的执行效率)；B^-1表示普朗克反函数。

3)将T₀代入下式，计算得到地表发射率ε_s

ε_s＝(R_surf-R_down)/(B(T₀)-R_down) (4)

其中，R_surf为离地辐亮度；R_down为大气下行辐亮度；B表示普朗克函数；T₀为温度初值。

4)对ε_s进行小波变换，获取重建后的地表发射率估值ε′_s

(a)首先，对ε_s进行小波变换，获取相应的高频和低频小波系数：

[C_A,C_D]＝WT(ε_s,db) (5)

其中，C_A，C_D分别表示低频小波系数与高频小波系数；ε_s为地表发射率；db表示所采用的小波基；WT表示小波变换。

由于小波变换的特点决定其不具有单一性，理论上存在无限多种的小波基，这使得本发明提出的方法适用性较强。因此，可根据不同的需求来选择不同的小波基，如更关注算法的整体执行效率，则可选择Haar小波基，而若更关注算法的反演精度，则可选用Daubechies或Symlets小波基等。

(b)然后，进行逆小波变换，将高频小波系数所对应的部分全部视为噪声予以排除，仅利用低频小波系数重建信号，获取重建后的地表发射率ε′_s：

ε′_s＝WT^-1(C_A,db) (6)

其中，WT^-1表示逆小波变换；C_A表示低频小波系数；db表示小波基。

5)构建代价函数

在已知地表温度T₀与重构后的发射率ε′_s后，我们可利用辐射传输方程估算其对应的星上辐亮度值R′_sensor:

R′_sensor＝(ε′_s·B(T₀)+(1-ε′_s)·R_down)·τ+R_up (7)

其中，R′_sensor为估算的星上辐亮度；ε′_s为重构的地表发射率；B(T₀)为温度T₀下的普朗克函数；R_down，R_up分别为大气下行及上行辐亮度；τ为大气透过率。此时，若地表温度与发射率接近真值，则其估算的星上辐亮度应与测量值近似相等。因此，可建立如下的代价函数：

C＝∑((R′_sensor-R_sensor)/mean(R_sensor))² (8)

其中，C为代价函数；R_sensor为星上观测辐亮度；R′_sensor为估算的星上辐亮度。基于非线性最优化理论，当代价函数C达到最小时，即为非线性方程的最优解，因此可以认为此时对应的地表温度T₀与发射率ε′_s便最为接近真实值。

6)将ε′_s与T₀带入式(8)，计算此时的代价函数C。

7)利用牛顿迭代法，并结合退火算法对地表温度进行扰动，迭代计算地表温度，并将其作为下一次迭代的温度初值T₀，重复步骤3)至步骤7)，直至满足代价函数达到最小值(即前后两次迭代过程中的代价函数差值小于一定阈值，而在综合考虑了算法精度和执行效率的基础上，认为10^-6是一个较为合理的数值)，结束迭代，获取最终的地表温度和发射率。

本发明充分考虑了地表发射率光谱曲线具有连续性的特点，无需发射率波谱形状等先验知识，将小波变化引入地表温度与发射率的同步分离过程，充分利用了小波变化在处理高频细节(噪声)信息上的优势，将地表发射率光谱转换至频域进行去噪处理，即仅利用低频信息实现发射率光谱的重构，在不损失反演精度的同时，最大程度的消除了各类噪声并恢复发射率光谱信息，更加简单易行，具有更广泛的适用性。

研究表明，在小波域，一个空间上(或时间域)具有一定连续性的有效信号，其小波系数模值往往较大，而对于具有随机性的噪声来说，其模值则很小。也就是说，当对含有噪声的信号进行小波分解时，含噪声部分主要包含在高频小波系数中，而信号信息则主要包含在低频小波系数中。另一方面，从地表发射率的光谱曲线上我们可以清晰的发现，地表发射率是一个随着波长连续变化的曲线，其正好符合小波变换在消除噪声中的应用特点，具有可行性，也为本发明的提出奠定了理论基础。

小波去噪方法包括三个基本的步骤：对含噪声信号进行小波变换；对变换得到的小波系数进行某种处理，以去除其中包含的噪声；对处理后的小波系数进行小波逆变换，得到去噪后的信号。因此，本发明主要在此基础上，克服了现有算法中对发射率波谱形状等先验知识的依赖，将地表温度与发射率的分离看做是消除观测噪声的过程，并利用发射率是连续性信号的特点，将小波变换引入地表温度和发射率的分离过程，通过在频域对发射率光谱进行分解与重构，成功减少了未知数个数，消除了辐射传输方程的病态性问题，并以循环迭代的方式逐步消除观测噪声，最终实现高光谱热红外地表温度与发射率的高精度分离。

本发明不依赖于任何先验知识，例如现有算法中均存在对光谱发射率曲线形状的先验约束，如光谱平滑或分段的线性表示等，而本发明是将发射率光谱看做是一个连续性的谱信号，通过小波变换的方式将其分解为高频和低频两部分信息，并基于小波变换的特点，认为高频部分对应的信息主要为噪声，可仅通过低频部分来实现原始发射率光谱的恢复，即通过分解与重构的方式，减少辐射传输方程的未知数个数，实现了地表温度与发射率的分离。算法本身不依赖于任何的先验假设或知识，简单易行，具有更好的适用性。此外，本发明受奇异值影响小。当大气等效温度与地表温度近似相等时，现有的光谱平滑算法和分段线性算法都会出现奇异值的问题，而本发明，通过谱域与频域的转换，以分解和重构的方式，仅利用发射率光谱的主要信息(低频)实现其光谱的重构，对波段范围内的噪声进行了综合消除，在一定程度上消弱了奇异值的影响，最大限度的保证了光谱形状的合理性。

上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。

Claims (2)

1.一种基于小波变换的高光谱热红外地表温度与发射率分离方法，由以下步骤组成：

1)获取离地辐亮度R_surf

通过辐射传输模型，在准确的配套大气温湿度剖面廓线数据的辅助下，估算大气上行/下行辐亮度和大气透过率；随后，利用辐射传输方程计算离地辐亮度R_surf，即

R_surf＝(R_sensor–R_up)/τ (2)

其中，R_surf为离地辐亮度；R_sensor为星上观测辐亮度；R_up为大气上行辐亮度；τ为大气透过率；

2)估计温度初值T₀

(a)首先，分别设定地表发射率的两个初值作为定值；

(b)然后，利用辐射传输方程分别计算两个发射率初值所对应的各波长位置的地表温度，并将其地表温度的通道温度最大值的均值作为初始温度T₀，即：

T₀＝(max(B^-1((R_surf-(1-ε₁)·R_down)/ε₁))+max(B^-1((R_surf-(1-ε₂)·R_down)/ε₂)))/2 (3)

其中，T₀为温度初值；R_surf为离地辐亮度；R_down为大气下行辐亮度；ε₁和ε₂分别为地表发射率初值；B^-1表示普朗克反函数；

3)将T₀代入下式，计算得到地表发射率ε_s

ε_s＝(R_surf-R_down)/(B(T₀)-R_down) (4)

其中，R_surf为离地辐亮度；R_down为大气下行辐亮度；B表示普朗克函数；T₀为温度初值；

4)对ε_s进行小波变换，获取重建后的地表发射率估值ε′_s

(a)首先，对ε_s进行小波变换，获取相应的高频和低频小波系数：

[C_A,C_D]＝WT(ε_s,db) (5)

其中，C_A，C_D分别表示低频小波系数与高频小波系数；ε_s为地表发射率；db表示所采用的小波基；WT表示小波变换；

(b)然后，进行逆小波变换，将高频小波系数所对应的部分全部视为噪声予以排除，仅利用低频小波系数重建信号，获取重建后的地表发射率ε′_s：

ε′_s＝WT^-1(C_A,db) (6)

其中，WT^-1表示逆小波变换；C_A表示低频小波系数；db表示小波基；

5)构建代价函数

在已知地表温度T₀与重构后的发射率ε′_s后，利用辐射传输方程估算其对应的星上辐亮度值R′_sensor：

R′_sensor＝(ε′_s·B(T₀)+(1-ε′_s)·R_down)·τ+R_up (7)

其中，R′_sensor为估算的星上辐亮度；ε′_s为重构的地表发射率；B(T₀)为温度T₀下的普朗克函数；R_up，R_down分别为大气上/下行辐亮度；τ为大气透过率；

此时，若地表温度与发射率接近真值，则其估算的星上辐亮度应与测量值近似相等；因此，可建立如下的代价函数：

C＝∑((R′_sensor-R_sensor)/mean(R_sensor))² (8)

其中，C为代价函数；R_sensor为星上观测辐亮度；R′_sensor为估算的星上辐亮度；当代价函数C达到最小时，其对应的地表温度T₀与发射率ε′_s便最为接近真实值；

6)将ε′_s与T₀带入式(8)，计算此时的代价函数C；

7)利用牛顿迭代法，并结合退火算法对地表温度进行扰动，迭代计算地表温度，并将其作为下一次迭代的温度初值T₀，重复步骤3)至步骤7)，直至满足代价函数达到最小值，结束迭代，获取最终的地表温度和发射率。

2.根据权利要求1所述的基于小波变换的高光谱热红外地表温度与发射率分离方法，其特征在于：所述步骤2)中的地表发射率初值ε₁和ε₂的取值分别为0.9和1.0。