基於小波變換的高光譜熱紅外地表溫度與發射率分離方法與流程
2023-11-02 07:03:52 2
本發明涉及一種分離方法，尤其涉及基於小波變換的高光譜熱紅外地表溫度與發射率分離方法。
背景技術：
地表發射率和地表溫度作為定量遙感領域的兩個關鍵參數，已被廣泛的應用於水文，生態，環境和生物地球化學等領域。近年來，越來越多的學者開始關注地表溫度和發射率在各領域中的關鍵作用，相關的研究也在不斷的開展。而這其中，關於如何通過遙感反演的方式精確的獲取高光譜熱紅外地表溫度與發射率一直是人們關注的研究熱點之一。
目前，利用遙感方式進行熱紅外地表溫度與發射率分離的最大的難點在於輻射傳輸方程病態性問題的求解。即便是在已獲取精準大氣信息的前提下，對於給定的N個通道的星上輻亮度，其可建立的方程組仍有N+1個未知數(N個通道發射率及1個地表溫度)，不可解。針對這一問題，現有的高光譜熱紅外地表溫度與發射率的分離算法從消除反演病態的出發點上來看，大致可分為有代表性的兩類：一類是引入約束增加方程組個數的反演算法，如光譜平滑算法；另一類是減少未知數個數的反演算法，如分段線性算法。
光譜平滑算法是在地表發射率具有平滑光譜的前提下，引入了平滑度函數的概念，並以此建立約束方程，消除輻射傳輸方程的病態性，實現地表溫度與發射率的分離。光譜平滑算法在一定程度上依賴於其平滑度函數，當實際發射率光譜無法滿足平滑條件時，尤其是對於一些特定的人造材料，其算法會導致一定誤差。此外，當大氣等效溫度與地表溫度接近時，採用光譜平滑算法往往會存在奇異值的問題。
分段線性地表溫度與發射率分離算法則是通過合理的設置波長區間，將區間內的地表發射率作為波長的函數，並用以線性的方式進行表示。通過這樣的方式，實現了減少輻射傳輸方程中未知數個數的目的，使得方程可解。通常以分段線性的方式對地表發射率進行表示所帶來的誤差可以忽略，其精度與光譜平滑算法相當。但是，受限於分段區間設定的不唯一性，對於地表發射率波動明顯的光譜曲線，過大的分段區間必然會導致光譜信息的丟失，而過小的分段區間則會損失線性子函數對發射率曲線的代表性，同時也會增加算法運算的複雜度。此外，分段線性算法將發射率光譜曲線進行了類似於離散化的處理(波長子區間)，這也導致了其無法保證各子區間的連續性，與實際發射率光譜連續的客觀事實相違背。
技術實現要素：
為了解決上述技術問題中的不足之處，本發明提供了一種基於小波變換的高光譜熱紅外地表溫度與發射率分離方法。
根據維恩位移定律可以得知，地表的輻射能量主要集中在熱紅外波段，並且此波段內的太陽輻射能量很小，可忽略不計。因此，晴空條件下，在滿足局地熱平衡時的熱紅外波段的輻射傳輸方程可以表示為：
Rsensor＝(ε·B(Ts)+(1-ε)·Rdown)·τ+Rup (1)
其中，Rsensor為衛星接收的星上輻亮度；ε為地表發射率；Ts為地表溫度；B(Ts)為溫度Ts的普朗克函數；Rdown大氣下行輻亮度；Rup為大氣上行輻亮度；τ為大氣透過率。
從式(1)可知，在已完成精準大氣校正的條件下，當估算的地表溫度Ts與真實地表溫度存在差異時，其反演得到的地表發射率光譜曲線將會含有大氣特徵，尤其是在大氣吸收線的位置，會出現明顯的峰谷現象。從獲取真實地表發射率的角度來看，此時在地表發射率光譜中的大氣特徵可以看做是一種大氣噪聲的表現。另一方面，由於衛星數據中往往會存在儀器的噪聲影響，而這種影響也會直接體現在地表發射率的光譜曲線之上。也就是說，獲取真實地表發射率可以看做是從含有噪聲的觀測信號中分離出地表溫度和發射率信息的過程。基於此，本專利提出了一種基於小波變換的高光譜熱紅外地表溫度與發射率分離算法。
為了解決上述技術問題，本發明採用的技術方案是：一種基於小波變換的高光譜熱紅外地表溫度與發射率分離方法，主要由以下步驟組成：
1)獲取離地輻亮度Rsurf
通過輻射傳輸模型，在準確的配套大氣溫溼度剖面廓線數據的輔助下，估算大氣上/下行輻亮度(Rup/Rdown)和大氣透過率(τ)；隨後，利用輻射傳輸方程計算離地輻亮度Rsurf，即
Rsurf＝(Rsensor–Rup)/τ (2)
其中，Rsurf為離地輻亮度；Rsensor為星上觀測輻亮度；Rup為大氣上行輻亮度；τ為大氣透過率。
2)估計溫度初值T0
(a)首先，分別設定兩個地表發射率的初值，ε1和ε2；
(b)然後，利用輻射傳輸方程分別計算兩個地表發射率初值所對應的各波長位置的地表溫度，並將其地表溫度的通道溫度最大值的均值作為初始溫度T0，即：
T0＝(max(B-1((Rsurf-(1-ε1)·Rdown)/ε1))+max(B-1((Rsurf-(1-ε2)·Rdown)/ε2)))/2 (3)
其中，T0為溫度初值；Rsurf為離地輻亮度；Rdown為大氣下行輻亮度；ε1和ε2分別為地表發射率初值；B-1表示普朗克反函數。
3)將T0代入下式，計算得到地表發射率εs
εs＝(Rsurf-Rdown)/(B(T0)-Rdown) (4)
其中，Rsurf為離地輻亮度；Rdown為大氣下行輻亮度；B表示普朗克函數；T0為溫度初值。
4)對εs進行小波變換，獲取重建後的地表發射率估值ε′s
(a)首先，對εs進行小波變換，獲取相應的高頻和低頻小波係數：
[CA,CD]＝WT(εs,db) (5)
其中，CA,CD分別表示低頻小波係數與高頻小波係數；εs為地表發射率；db表示所採用的小波基；WT表示小波變換。
(b)然後，進行逆小波變換，將高頻小波係數所對應的部分全部視為噪聲予以排除，僅利用低頻小波係數重建信號，獲取重建後的地表發射率ε′s：
ε′s＝WT-1(CA,db) (6)
其中，WT-1表示逆小波變換；CA表示低頻小波係數；db表示小波基。
5)構建代價函數
在已知地表溫度T0與重構後的發射率ε′s後，我們可利用輻射傳輸方程估算其對應的星上輻亮度值R′sensor:
R′sensor＝(ε′s·B(T0)+(1-ε′s)·Rdown)·τ+Rup (7)
其中，R′sensor為估算的星上輻亮度；ε′s為重構的地表發射率；B(T0)為溫度T0下的普朗克函數；Rdown，Rup分別為大氣下行及上行輻亮度；τ為大氣透過率。此時，若地表溫度與發射率接近真值，則其估算的星上輻亮度應與測量值近似相等。因此，可建立如下的代價函數：
C＝∑((R′sensor-Rsensor)/mean(Rsensor))2 (8)
其中，C為代價函數；Rsensor為星上觀測輻亮度；R′sensor為估算的星上輻亮度；當代價函數C達到最小時，其對應的地表溫度T0與發射率ε′s便最為接近真實值。
6)將ε′s與T0帶入式(8)，計算此時的代價函數C。
7)利用牛頓迭代法，並結合退火算法對地表溫度進行擾動，迭代計算地表溫度，並將其作為下一次迭代的溫度初值T0，重複步驟(3)至(7)，直至滿足代價函數達到最小值，結束迭代，獲取最終的地表溫度和發射率。
其中，步驟2)中ε1和ε2的取值分別為0.9和1.0。
本發明克服了對發射率光譜形狀假設的依賴，通過小波變換將地表發射率從譜域信號轉為頻域信號，並通過獲取的地表發射率的低頻信息重構原始發射率光譜，這一過程合理減少了反演中未知數個數，消除了遙感反演的病態性，實現了高光譜熱紅外地表溫度和發射率的高精度遙感分離和準確獲取。
附圖說明
下面結合附圖和具體實施方式對本發明作進一步的詳細說明。
圖1是本發明方法的整體流程圖。
具體實施方式
如圖1所示，本發明主要由以下步驟組成：
1)獲取離地輻亮度Rsurf
通過輻射傳輸模型，在準確的配套大氣溫溼度剖面廓線數據的輔助下，估算大氣上/下行輻亮度(Rup/Rdown)和大氣透過率(τ)；隨後，利用輻射傳輸方程計算離地輻亮度Rsurf，即
Rsurf＝(Rsensor–Rup)/τ (2)
其中，Rsurf為離地輻亮度；Rsensor為星上觀測輻亮度；Rup為大氣上行輻亮度；τ為大氣透過率。
2)估計溫度初值T0
(a)首先，分別設定兩個地表發射率初值，如ε1＝0.9，ε2＝1.0；
(b)然後，利用輻射傳輸方程分別計算兩個地表發射率初值所對應的各波長位置的地表溫度，並將其對應的波長位置的地表溫度的通道溫度最大值的均值作為初始溫度T0，即：
T0＝(max(B-1((Rsurf-(1-ε1)·Rdown)/ε1))+max(B-1((Rsurf-(1-ε2)·Rdown)/ε2)))/2 (3)
其中，T0為溫度初值；Rsurf為離地輻亮度；Rdown為大氣下行輻亮度；ε1＝0.9,ε2＝1.0為地表發射率初值(對於典型地類，地表發射率在熱紅外譜段的取值範圍通常介於0.9到1之間，經過多次試驗表明，當其初值分別取值為0.9和1.0時，算法可具有較好的執行效率)；B-1表示普朗克反函數。
3)將T0代入下式，計算得到地表發射率εs
εs＝(Rsurf-Rdown)/(B(T0)-Rdown) (4)
其中，Rsurf為離地輻亮度；Rdown為大氣下行輻亮度；B表示普朗克函數；T0為溫度初值。
4)對εs進行小波變換，獲取重建後的地表發射率估值ε′s
(a)首先，對εs進行小波變換，獲取相應的高頻和低頻小波係數：
[CA,CD]＝WT(εs,db) (5)
其中，CA，CD分別表示低頻小波係數與高頻小波係數；εs為地表發射率；db表示所採用的小波基；WT表示小波變換。
由於小波變換的特點決定其不具有單一性，理論上存在無限多種的小波基，這使得本發明提出的方法適用性較強。因此，可根據不同的需求來選擇不同的小波基，如更關注算法的整體執行效率，則可選擇Haar小波基，而若更關注算法的反演精度，則可選用Daubechies或Symlets小波基等。
(b)然後，進行逆小波變換，將高頻小波係數所對應的部分全部視為噪聲予以排除，僅利用低頻小波係數重建信號，獲取重建後的地表發射率ε′s：
ε′s＝WT-1(CA,db) (6)
其中，WT-1表示逆小波變換；CA表示低頻小波係數；db表示小波基。
5)構建代價函數
在已知地表溫度T0與重構後的發射率ε′s後，我們可利用輻射傳輸方程估算其對應的星上輻亮度值R′sensor:
R′sensor＝(ε′s·B(T0)+(1-ε′s)·Rdown)·τ+Rup (7)
其中，R′sensor為估算的星上輻亮度；ε′s為重構的地表發射率；B(T0)為溫度T0下的普朗克函數；Rdown，Rup分別為大氣下行及上行輻亮度；τ為大氣透過率。此時，若地表溫度與發射率接近真值，則其估算的星上輻亮度應與測量值近似相等。因此，可建立如下的代價函數：
C＝∑((R′sensor-Rsensor)/mean(Rsensor))2 (8)
其中，C為代價函數；Rsensor為星上觀測輻亮度；R′sensor為估算的星上輻亮度。基於非線性最優化理論，當代價函數C達到最小時，即為非線性方程的最優解，因此可以認為此時對應的地表溫度T0與發射率ε′s便最為接近真實值。
6)將ε′s與T0帶入式(8)，計算此時的代價函數C。
7)利用牛頓迭代法，並結合退火算法對地表溫度進行擾動，迭代計算地表溫度，並將其作為下一次迭代的溫度初值T0，重複步驟(3)至(7)，直至滿足代價函數達到最小值(即前後兩次迭代過程中的代價函數差值小於一定閾值，而在綜合考慮了算法精度和執行效率的基礎上，認為10-6是一個較為合理的數值)，結束迭代，獲取最終的地表溫度和發射率。
本發明充分考慮了地表發射率光譜曲線具有連續性的特點，無需發射率波譜形狀等先驗知識，將小波變化引入地表溫度與發射率的同步分離過程，充分利用了小波變化在處理高頻細節(噪聲)信息上的優勢，將地表發射率光譜轉換至頻域進行去噪處理，即僅利用低頻信息實現發射率光譜的重構，在不損失反演精度的同時，最大程度的消除了各類噪聲並恢復發射率光譜信息，更加簡單易行，具有更廣泛的適用性。
研究表明，在小波域，一個空間上(或時間域)具有一定連續性的有效信號，其小波係數模值往往較大，而對於具有隨機性的噪聲來說，其模值則很小。也就是說，當對含有噪聲的信號進行小波分解時，含噪聲部分主要包含在高頻小波係數中，而信號信息則主要包含在低頻小波係數中。另一方面，從地表發射率的光譜曲線上我們可以清晰的發現，地表發射率是一個隨著波長連續變化的曲線，其正好符合小波變換在消除噪聲中的應用特點，具有可行性，也為本發明的提出奠定了理論基礎。
小波去噪方法包括三個基本的步驟：對含噪聲信號進行小波變換；對變換得到的小波係數進行某種處理，以去除其中包含的噪聲；對處理後的小波係數進行小波逆變換，得到去噪後的信號。因此，本發明主要在此基礎上，克服了現有算法中對發射率波譜形狀等先驗知識的依賴，將地表溫度與發射率的分離看做是消除觀測噪聲的過程，並利用發射率是連續性信號的特點，將小波變換引入地表溫度和發射率的分離過程，通過在頻域對發射率光譜進行分解與重構，成功減少了未知數個數，消除了輻射傳輸方程的病態性問題，並以循環迭代的方式逐步消除觀測噪聲，最終實現高光譜熱紅外地表溫度與發射率的高精度分離。
本發明不依賴於任何先驗知識，例如現有算法中均存在對光譜發射率曲線形狀的先驗約束，如光譜平滑或分段的線性表示等，而本發明是將發射率光譜看做是一個連續性的譜信號，通過小波變換的方式將其分解為高頻和低頻兩部分信息，並基於小波變換的特點，認為高頻部分對應的信息主要為噪聲，可僅通過低頻部分來實現原始發射率光譜的恢復，即通過分解與重構的方式，減少輻射傳輸方程的未知數個數，實現了地表溫度與發射率的分離。算法本身不依賴於任何的先驗假設或知識，簡單易行，具有更好的適用性。此外，本發明受奇異值影響小。當大氣等效溫度與地表溫度近似相等時，現有的光譜平滑算法和分段線性算法都會出現奇異值的問題，而本發明，通過譜域與頻域的轉換，以分解和重構的方式，僅利用發射率光譜的主要信息(低頻)實現其光譜的重構，對波段範圍內的噪聲進行了綜合消除，在一定程度上消弱了奇異值的影響，最大限度的保證了光譜形狀的合理性。
上述實施方式並非是對本發明的限制，本發明也並不僅限於上述舉例，本技術領域的技術人員在本發明的技術方案範圍內所做出的變化、改型、添加或替換，也均屬於本發明的保護範圍。