本發(fā)明涉及一種分離方法，尤其涉及基于小波變換的高光譜熱紅外地表溫度與發(fā)射率分離方法。

背景技術：

地表發(fā)射率和地表溫度作為定量遙感領域的兩個關鍵參數(shù)，已被廣泛的應用于水文，生態(tài)，環(huán)境和生物地球化學等領域。近年來，越來越多的學者開始關注地表溫度和發(fā)射率在各領域中的關鍵作用，相關的研究也在不斷的開展。而這其中，關于如何通過遙感反演的方式精確的獲取高光譜熱紅外地表溫度與發(fā)射率一直是人們關注的研究熱點之一。

目前，利用遙感方式進行熱紅外地表溫度與發(fā)射率分離的最大的難點在于輻射傳輸方程病態(tài)性問題的求解。即便是在已獲取精準大氣信息的前提下，對于給定的N個通道的星上輻亮度，其可建立的方程組仍有N+1個未知數(shù)(N個通道發(fā)射率及1個地表溫度)，不可解。針對這一問題，現(xiàn)有的高光譜熱紅外地表溫度與發(fā)射率的分離算法從消除反演病態(tài)的出發(fā)點上來看，大致可分為有代表性的兩類：一類是引入約束增加方程組個數(shù)的反演算法，如光譜平滑算法；另一類是減少未知數(shù)個數(shù)的反演算法，如分段線性算法。

光譜平滑算法是在地表發(fā)射率具有平滑光譜的前提下，引入了平滑度函數(shù)的概念，并以此建立約束方程，消除輻射傳輸方程的病態(tài)性，實現(xiàn)地表溫度與發(fā)射率的分離。光譜平滑算法在一定程度上依賴于其平滑度函數(shù)，當實際發(fā)射率光譜無法滿足平滑條件時，尤其是對于一些特定的人造材料，其算法會導致一定誤差。此外，當大氣等效溫度與地表溫度接近時，采用光譜平滑算法往往會存在奇異值的問題。

分段線性地表溫度與發(fā)射率分離算法則是通過合理的設置波長區(qū)間，將區(qū)間內的地表發(fā)射率作為波長的函數(shù)，并用以線性的方式進行表示。通過這樣的方式，實現(xiàn)了減少輻射傳輸方程中未知數(shù)個數(shù)的目的，使得方程可解。通常以分段線性的方式對地表發(fā)射率進行表示所帶來的誤差可以忽略，其精度與光譜平滑算法相當。但是，受限于分段區(qū)間設定的不唯一性，對于地表發(fā)射率波動明顯的光譜曲線，過大的分段區(qū)間必然會導致光譜信息的丟失，而過小的分段區(qū)間則會損失線性子函數(shù)對發(fā)射率曲線的代表性，同時也會增加算法運算的復雜度。此外，分段線性算法將發(fā)射率光譜曲線進行了類似于離散化的處理(波長子區(qū)間)，這也導致了其無法保證各子區(qū)間的連續(xù)性，與實際發(fā)射率光譜連續(xù)的客觀事實相違背。

技術實現(xiàn)要素：

為了解決上述技術問題中的不足之處，本發(fā)明提供了一種基于小波變換的高光譜熱紅外地表溫度與發(fā)射率分離方法。

根據(jù)維恩位移定律可以得知，地表的輻射能量主要集中在熱紅外波段，并且此波段內的太陽輻射能量很小，可忽略不計。因此，晴空條件下，在滿足局地熱平衡時的熱紅外波段的輻射傳輸方程可以表示為：

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

其中，R_sensor為衛(wèi)星接收的星上輻亮度；ε為地表發(fā)射率；T_s為地表溫度；B(T_s)為溫度T_s的普朗克函數(shù)；R_down大氣下行輻亮度；R_up為大氣上行輻亮度；τ為大氣透過率。

從式(1)可知，在已完成精準大氣校正的條件下，當估算的地表溫度T_s與真實地表溫度存在差異時，其反演得到的地表發(fā)射率光譜曲線將會含有大氣特征，尤其是在大氣吸收線的位置，會出現(xiàn)明顯的峰谷現(xiàn)象。從獲取真實地表發(fā)射率的角度來看，此時在地表發(fā)射率光譜中的大氣特征可以看做是一種大氣噪聲的表現(xiàn)。另一方面，由于衛(wèi)星數(shù)據(jù)中往往會存在儀器的噪聲影響，而這種影響也會直接體現(xiàn)在地表發(fā)射率的光譜曲線之上。也就是說，獲取真實地表發(fā)射率可以看做是從含有噪聲的觀測信號中分離出地表溫度和發(fā)射率信息的過程?；诖?，本專利提出了一種基于小波變換的高光譜熱紅外地表溫度與發(fā)射率分離算法。

為了解決上述技術問題，本發(fā)明采用的技術方案是：一種基于小波變換的高光譜熱紅外地表溫度與發(fā)射率分離方法，主要由以下步驟組成：

1)獲取離地輻亮度R_surf

通過輻射傳輸模型，在準確的配套大氣溫濕度剖面廓線數(shù)據(jù)的輔助下，估算大氣上/下行輻亮度(R_up/R_down)和大氣透過率(τ)；隨后，利用輻射傳輸方程計算離地輻亮度R_surf，即

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

其中，R_surf為離地輻亮度；R_sensor為星上觀測輻亮度；R_up為大氣上行輻亮度；τ為大氣透過率。

2)估計溫度初值T₀

(a)首先，分別設定兩個地表發(fā)射率的初值，ε₁和ε₂；

(b)然后，利用輻射傳輸方程分別計算兩個地表發(fā)射率初值所對應的各波長位置的地表溫度，并將其地表溫度的通道溫度最大值的均值作為初始溫度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為大氣下行輻亮度；ε₁和ε₂分別為地表發(fā)射率初值；B^-1表示普朗克反函數(shù)。

3)將T₀代入下式，計算得到地表發(fā)射率ε_s

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

其中，R_surf為離地輻亮度；R_down為大氣下行輻亮度；B表示普朗克函數(shù)；T₀為溫度初值。

4)對ε_s進行小波變換，獲取重建后的地表發(fā)射率估值ε′_s

(a)首先，對ε_s進行小波變換，獲取相應的高頻和低頻小波系數(shù)：

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

其中，C_A,C_D分別表示低頻小波系數(shù)與高頻小波系數(shù)；ε_s為地表發(fā)射率；db表示所采用的小波基；WT表示小波變換。

(b)然后，進行逆小波變換，將高頻小波系數(shù)所對應的部分全部視為噪聲予以排除，僅利用低頻小波系數(shù)重建信號，獲取重建后的地表發(fā)射率ε′_s：

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

其中，WT^-1表示逆小波變換；C_A表示低頻小波系數(shù)；db表示小波基。

5)構建代價函數(shù)

在已知地表溫度T₀與重構后的發(fā)射率ε′_s后，我們可利用輻射傳輸方程估算其對應的星上輻亮度值R′_sensor:

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

其中，R′_sensor為估算的星上輻亮度；ε′_s為重構的地表發(fā)射率；B(T₀)為溫度T₀下的普朗克函數(shù)；R_down，R_up分別為大氣下行及上行輻亮度；τ為大氣透過率。此時，若地表溫度與發(fā)射率接近真值，則其估算的星上輻亮度應與測量值近似相等。因此，可建立如下的代價函數(shù)：

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

其中，C為代價函數(shù)；R_sensor為星上觀測輻亮度；R′_sensor為估算的星上輻亮度；當代價函數(shù)C達到最小時，其對應的地表溫度T₀與發(fā)射率ε′_s便最為接近真實值。

6)將ε′_s與T₀帶入式(8)，計算此時的代價函數(shù)C。

7)利用牛頓迭代法，并結合退火算法對地表溫度進行擾動，迭代計算地表溫度，并將其作為下一次迭代的溫度初值T₀，重復步驟(3)至(7)，直至滿足代價函數(shù)達到最小值，結束迭代，獲取最終的地表溫度和發(fā)射率。

其中，步驟2)中ε₁和ε₂的取值分別為0.9和1.0。

本發(fā)明克服了對發(fā)射率光譜形狀假設的依賴，通過小波變換將地表發(fā)射率從譜域信號轉為頻域信號，并通過獲取的地表發(fā)射率的低頻信息重構原始發(fā)射率光譜，這一過程合理減少了反演中未知數(shù)個數(shù)，消除了遙感反演的病態(tài)性，實現(xiàn)了高光譜熱紅外地表溫度和發(fā)射率的高精度遙感分離和準確獲取。

附圖說明

下面結合附圖和具體實施方式對本發(fā)明作進一步的詳細說明。

圖1是本發(fā)明方法的整體流程圖。

具體實施方式

如圖1所示，本發(fā)明主要由以下步驟組成：

1)獲取離地輻亮度R_surf

通過輻射傳輸模型，在準確的配套大氣溫濕度剖面廓線數(shù)據(jù)的輔助下，估算大氣上/下行輻亮度(R_up/R_down)和大氣透過率(τ)；隨后，利用輻射傳輸方程計算離地輻亮度R_surf，即

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

其中，R_surf為離地輻亮度；R_sensor為星上觀測輻亮度；R_up為大氣上行輻亮度；τ為大氣透過率。

2)估計溫度初值T₀

(a)首先，分別設定兩個地表發(fā)射率初值，如ε₁＝0.9，ε₂＝1.0；

(b)然后，利用輻射傳輸方程分別計算兩個地表發(fā)射率初值所對應的各波長位置的地表溫度，并將其對應的波長位置的地表溫度的通道溫度最大值的均值作為初始溫度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為地表發(fā)射率初值(對于典型地類，地表發(fā)射率在熱紅外譜段的取值范圍通常介于0.9到1之間，經(jīng)過多次試驗表明，當其初值分別取值為0.9和1.0時，算法可具有較好的執(zhí)行效率)；B^-1表示普朗克反函數(shù)。

3)將T0代入下式，計算得到地表發(fā)射率ε_s

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

其中，R_surf為離地輻亮度；R_down為大氣下行輻亮度；B表示普朗克函數(shù)；T₀為溫度初值。

4)對ε_s進行小波變換，獲取重建后的地表發(fā)射率估值ε′_s

(a)首先，對ε_s進行小波變換，獲取相應的高頻和低頻小波系數(shù)：

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

其中，C_A，C_D分別表示低頻小波系數(shù)與高頻小波系數(shù)；ε_s為地表發(fā)射率；db表示所采用的小波基；WT表示小波變換。

由于小波變換的特點決定其不具有單一性，理論上存在無限多種的小波基，這使得本發(fā)明提出的方法適用性較強。因此，可根據(jù)不同的需求來選擇不同的小波基，如更關注算法的整體執(zhí)行效率，則可選擇Haar小波基，而若更關注算法的反演精度，則可選用Daubechies或Symlets小波基等。

(b)然后，進行逆小波變換，將高頻小波系數(shù)所對應的部分全部視為噪聲予以排除，僅利用低頻小波系數(shù)重建信號，獲取重建后的地表發(fā)射率ε′_s：

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

其中，WT^-1表示逆小波變換；C_A表示低頻小波系數(shù)；db表示小波基。

5)構建代價函數(shù)

在已知地表溫度T₀與重構后的發(fā)射率ε′_s后，我們可利用輻射傳輸方程估算其對應的星上輻亮度值R′_sensor:

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

其中，R′_sensor為估算的星上輻亮度；ε′_s為重構的地表發(fā)射率；B(T₀)為溫度T₀下的普朗克函數(shù)；R_down，R_up分別為大氣下行及上行輻亮度；τ為大氣透過率。此時，若地表溫度與發(fā)射率接近真值，則其估算的星上輻亮度應與測量值近似相等。因此，可建立如下的代價函數(shù)：

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

其中，C為代價函數(shù)；R_sensor為星上觀測輻亮度；R′_sensor為估算的星上輻亮度?；诜蔷€性最優(yōu)化理論，當代價函數(shù)C達到最小時，即為非線性方程的最優(yōu)解，因此可以認為此時對應的地表溫度T₀與發(fā)射率ε′_s便最為接近真實值。

6)將ε′_s與T₀帶入式(8)，計算此時的代價函數(shù)C。

7)利用牛頓迭代法，并結合退火算法對地表溫度進行擾動，迭代計算地表溫度，并將其作為下一次迭代的溫度初值T₀，重復步驟(3)至(7)，直至滿足代價函數(shù)達到最小值(即前后兩次迭代過程中的代價函數(shù)差值小于一定閾值，而在綜合考慮了算法精度和執(zhí)行效率的基礎上，認為10^-6是一個較為合理的數(shù)值)，結束迭代，獲取最終的地表溫度和發(fā)射率。

本發(fā)明充分考慮了地表發(fā)射率光譜曲線具有連續(xù)性的特點，無需發(fā)射率波譜形狀等先驗知識，將小波變化引入地表溫度與發(fā)射率的同步分離過程，充分利用了小波變化在處理高頻細節(jié)(噪聲)信息上的優(yōu)勢，將地表發(fā)射率光譜轉換至頻域進行去噪處理，即僅利用低頻信息實現(xiàn)發(fā)射率光譜的重構，在不損失反演精度的同時，最大程度的消除了各類噪聲并恢復發(fā)射率光譜信息，更加簡單易行，具有更廣泛的適用性。

研究表明，在小波域，一個空間上(或時間域)具有一定連續(xù)性的有效信號，其小波系數(shù)模值往往較大，而對于具有隨機性的噪聲來說，其模值則很小。也就是說，當對含有噪聲的信號進行小波分解時，含噪聲部分主要包含在高頻小波系數(shù)中，而信號信息則主要包含在低頻小波系數(shù)中。另一方面，從地表發(fā)射率的光譜曲線上我們可以清晰的發(fā)現(xiàn)，地表發(fā)射率是一個隨著波長連續(xù)變化的曲線，其正好符合小波變換在消除噪聲中的應用特點，具有可行性，也為本發(fā)明的提出奠定了理論基礎。

小波去噪方法包括三個基本的步驟：對含噪聲信號進行小波變換；對變換得到的小波系數(shù)進行某種處理，以去除其中包含的噪聲；對處理后的小波系數(shù)進行小波逆變換，得到去噪后的信號。因此，本發(fā)明主要在此基礎上，克服了現(xiàn)有算法中對發(fā)射率波譜形狀等先驗知識的依賴，將地表溫度與發(fā)射率的分離看做是消除觀測噪聲的過程，并利用發(fā)射率是連續(xù)性信號的特點，將小波變換引入地表溫度和發(fā)射率的分離過程，通過在頻域對發(fā)射率光譜進行分解與重構，成功減少了未知數(shù)個數(shù)，消除了輻射傳輸方程的病態(tài)性問題，并以循環(huán)迭代的方式逐步消除觀測噪聲，最終實現(xiàn)高光譜熱紅外地表溫度與發(fā)射率的高精度分離。

本發(fā)明不依賴于任何先驗知識，例如現(xiàn)有算法中均存在對光譜發(fā)射率曲線形狀的先驗約束，如光譜平滑或分段的線性表示等，而本發(fā)明是將發(fā)射率光譜看做是一個連續(xù)性的譜信號，通過小波變換的方式將其分解為高頻和低頻兩部分信息，并基于小波變換的特點，認為高頻部分對應的信息主要為噪聲，可僅通過低頻部分來實現(xiàn)原始發(fā)射率光譜的恢復，即通過分解與重構的方式，減少輻射傳輸方程的未知數(shù)個數(shù)，實現(xiàn)了地表溫度與發(fā)射率的分離。算法本身不依賴于任何的先驗假設或知識，簡單易行，具有更好的適用性。此外，本發(fā)明受奇異值影響小。當大氣等效溫度與地表溫度近似相等時，現(xiàn)有的光譜平滑算法和分段線性算法都會出現(xiàn)奇異值的問題，而本發(fā)明，通過譜域與頻域的轉換，以分解和重構的方式，僅利用發(fā)射率光譜的主要信息(低頻)實現(xiàn)其光譜的重構，對波段范圍內的噪聲進行了綜合消除，在一定程度上消弱了奇異值的影響，最大限度的保證了光譜形狀的合理性。

上述實施方式并非是對本發(fā)明的限制，本發(fā)明也并不僅限于上述舉例，本技術領域的技術人員在本發(fā)明的技術方案范圍內所做出的變化、改型、添加或替換，也均屬于本發(fā)明的保護范圍。