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

背景技術(shù)：

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

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

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

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

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

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

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

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

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

1)獲取離地輻亮度R_surf

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

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

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

2)估計溫度初值T₀

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

(b)然后，利用輻射傳輸方程分別計算兩個地表發(fā)射率初值所對應(yīng)的各波長位置的地表溫度，并將其地表溫度的通道溫度最大值的均值作為初始溫度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進行小波變換，獲取相應(yīng)的高頻和低頻小波系數(shù)：

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

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

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

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

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

5)構(gòu)建代價函數(shù)

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

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

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

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

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

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

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

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

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

附圖說明

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

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

具體實施方式

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

1)獲取離地輻亮度R_surf

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

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

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

2)估計溫度初值T₀

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

(b)然后，利用輻射傳輸方程分別計算兩個地表發(fā)射率初值所對應(yīng)的各波長位置的地表溫度，并將其對應(yīng)的波長位置的地表溫度的通道溫度最大值的均值作為初始溫度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)過多次試驗表明，當(dā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進行小波變換，獲取相應(yīng)的高頻和低頻小波系數(shù)：

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

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

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

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

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

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

5)構(gòu)建代價函數(shù)

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

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

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

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

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

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

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

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

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

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

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

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