孫舜蓉
（工業和資訊化部電信設備認證中心， 北京 100088）
：為了滿足遙感衛星即時雲檢測要求， 提出了一種基於光譜角的雲檢測演算法。 首先在多光譜譜段構成的歐式空間中構造參考雲向量；然後計算圖像中各像元向量與參考雲向量的光譜角， 對光譜角進行高斯濾波生成雲掩膜；最後計算得到圖像中雲量百分比。 利用該演算法對156景人工標注過的LandSat7有雲圖像進行雲檢測， 結果表明本文演算法雲檢測整體精度能到達83.6%， 平均執行速度為360 ms。 本文演算法原理簡單， 執行速度快， 檢測效果好， 具有星上即時應用潛力。
：雲污染；多光譜遙感；光譜角；雲檢測
：P407.8文獻標識碼：ADOI： 10.19358/j.issn.1674-7720.2017.06.006
引用格式：孫舜蓉. 一種基於光譜角原理的多光譜遙感圖像雲檢測演算法［J］.微型機與應用， 2017,36（6）：16-18， 21.
0引言
未來的遙感衛星將向著更高空間解析度、更高時間解析度、更高光譜解析度以及更高輻射解析度的方向發展， 這將導致圖像資料量的成倍增長， 給星上硬體資源帶來巨大壓力。 中國科學院地理空間資料雲系統［1］對2014年1月1日到2015年1月1日的Landsat 7 ETM+ 湖北省的135景圖像進行雲量統計， 結果見表1所示。 若定義雲量小於5%的圖像為乾淨圖像， 則該比例僅為20%， 即80%的圖像都或多或少地被雲圖污染。 雲量在80%以上的圖像55景， 占總圖像的40.78%， 這部分圖像屬於雲圖嚴重污染圖像， 基本上沒有應用價值， 因此可以不需要從衛星下傳到地面。
除氣象衛星外， 雲圖在高解析度遙感圖像中有害無益， 一是污染圖像， 對地物識別造成困難；二是佔用了大量寶貴的運算和存儲資源；三是佔用有限的傳輸頻寬。 這三方面使得相機有效成像時間縮短， 極大地限制了衛星的成像能力。 因此， 開發星上即時高效的雲檢測演算法勢在必行。 星上可根據雲量的檢測結果對原始圖像進行針對性的處理， 可以有效改善衛星成像效率， 然而星上雲檢測遠不如地面雲檢測成熟。
雲檢測演算法在文獻［2］、［3］中得到很好的歸納總結。
（1)雲判演算法多為地面應用， 發展比較成熟， 但星上即時雲判研究較少。
（2)單譜段雲圖檢測基於圖像閾值、圖像紋理進行雲判別。 由於雲千變萬化， 雲檢測分類器往往難以區分與雲具有相似紋理的冰雪下墊面。
（3)雲檢測演算法針對特定的感測器進行量身定制， 還沒有一種雲檢測演算法能適應所有類型的感測器。
（4)地面演算法較為複雜， 不適合直接移植到星上應用， 基於光譜和紋理的分類思想具有借鑒意義。
本文在分析雲和典型地物的光譜反射率特性的基礎上， 給出利用光譜角進行雲檢測的演算法原理和步驟， 然後對比自動雲覆蓋評估演算法， 給出雲檢測結果以及代碼執行效率，
1演算法原理及實現流程
不同地物由於物理特徵的不同， 各自具有獨特的光譜反射特性。 圖1為雲和典型地物的光譜反射率曲線。 由圖1可見， 在0.4μm ~0.9μm可見光區域， 雲具有相對高的反射率， 但是比冰雪低；在1.4μm以上的短波紅外譜段， 雲的反射率比冰雪、土壤和植被的反射率都高。 多光譜探測器的光譜響應在其譜段範圍內響應是近似線性的， 那麼， 典型地物反射率曲線等價於探測器對地物的光譜回應曲線， 也就是說地物反射率越高， 表現為圖像的量化值越大。 而且， 譜段越多， 光譜解析度越高， 描繪的光譜曲線越完整， 利用發射率資訊對地物的區分度就越高。 本文提出一種利用光譜向量表徵地物， 進而利用光譜向量與標準地物參考向量之間的相似性，
對於具有n個譜段的多光譜圖像， 圖像像元可用n維向量來表示， 例如圖2中的Pi。 假設Pr為某地物參考向量， 若Pi與Pr的指向相近， 且長度也相近， 則可以認為Pi與Pr相似。Pi的長度可以通過對Pi的求模操作即|Pi|來得到，而Pi與Pr間的光譜角可用如下數學運算式得到:
αi=arccos(Pi·Pr/|Pi|·|Pr|)(1)
式中:αi為圖像像元向量與參考像元向量之間的夾角，即光譜角，代表了光譜回應的相似性，其變化範圍是［0，π/2］。Pi·Pr為向量點乘操作，|Pi|表示Pi的長度，|Pr|表示Pr的長度。光譜角度越小，同時向量長度越接近，則被估計像元的光譜回應與參考像元的光譜回應就越相似，參考物與被估計像元為同一類型的概率就越高。如圖2所示，假設圖像有3個譜段，分別用i，j，k表示，αi就是像元Pi與參考像元Pr之間的光譜角。
基於光譜角的雲檢測流程如圖3所示。
基於以上分析，提出基於光譜角的雲檢測演算法如下：
(1)計算遙感圖像某像元Pi與雲參考像元Pr的之間的光譜角，見式（1）。
(2)為了利用向量的長度資訊，構造一個包含像元Pi向量與雲參考像元Pr向量長度的高斯濾波器，對αi進行濾波，即乘以一個高斯權重係數GsCoef，該權重係數由式（2）得到，範圍是（0,1 ］。高斯係數的作用是給那些向量長度與參考像元向量長度接近的像元賦予較大的權重，直觀解釋就是向量長度相近的像元更有可能是同一類地物。濾波結果為Cgsi，其範圍是（0，1］。
GsCoef=e-((|Pi|-|Pr|)2/(0.35·|Pr|2))(2)
Cgsi=GsCoef*（π-αi）/π(3)
(3)取適當的閾值對Cgsi進行二值化，即當minThreshold
if (Cgsi>minThreshold && Cgsi
Ci=1;//該像元判斷為雲
else
Ci=0;//該像元判斷為非雲
(4)重複步驟(1)~(3)，依次完成所有像元檢測，生成雲掩膜。統計雲像元所占比例，為後續處理做好準備。
2實驗結果
為了便於比較，本文利用開源遙感影像處理軟體Grass7.0提供的自動雲評估ACCA演算法C++代碼，與文中的光譜角雲檢測演算法C++代碼，從整體檢測精度（Cloud Overall Accuracy, COA）以及執行速度上進行比較。計算平臺為：Lenovo Windows 7 64 bit PC; Intel Core i52520 @2.5 GHz CPU; 4.00 GB RAM。
從美國地質調查局網站得到156景經過人工標注的Landsat7圖像，作為本文雲檢測精度評價的標準圖像。此156景資料分佈在9個不同的緯度區域，對不同地域的雲分佈有一定統計意義，同時它們是文獻［3］中所採用的212幅人工標注圖像的子集。在這156景8 bit量化的遙感圖像中，圖像專業分析人員利用Photoshop中的圖像套索工具一點點地分離雲和非雲像元。
綜合考慮檢測精度和運算速度，本文取Landsat7中的Band1(450~515 nm)，Band3(630~690 nm)，Band4(750~900 nm)，Band5(1 550~1 750 nm)作為試驗對象，利用這四個譜段構造光譜角進行雲檢測。同時，為了更好地模擬在軌即時處理過程，利用從原始圖像截取的1 024像元×1 024像元大小的含不同雲量的圖像進行試驗。雲參考向量Pr取(225,215,182,168)，minThreshold取0.6，maxThreshold取1，比較結果如圖4~6所示。
對於參考圖像，ACCA演算法耗時較長，平均耗時為2 299.6 ms。耗時曲線中顯示對不同的圖像，耗時相差較大，這是因為ACCA演算法針對不同的地物採用不同的演算法分支，而某些分支中包含有迴圈處理。SA執行速度很快，平均耗時僅為361.1 ms。
圖7是帶雲的原始圖像，圖8是人工標注的雲圖像，圖中，雲標注為全白，非雲地物為全黑。典型的ACCA雲檢測結果如圖9所示，SA檢測結果如圖10所示。整體檢測精度分佈如表2所示，可知ACCA的整體檢測精度達到89.05%，高於SA演算法的83.6%。對檢測誤差較大的圖像進行分析發現，高緯度地區覆蓋的冰和雪常被誤檢為雲，因為冰雪反射率與雲接近，且具有較高亮度；對薄雲的漏檢在兩種演算法中均存在，SA演算法相對較嚴重。
3結論
分析可知，SA的檢測精度略低於ACCA，但是仍然高於83%。最為關鍵的是，SA演算法執行時間僅為ACCA演算法執行時間的1/6，這在即時影像處理中尤為重要。對於高緯度地區SA演算法的過檢問題，可以針對高緯度地區的冰雪特點，適當調高閾值，提高檢測精度。在軌應用時，一些實際問題例如圖像的輻射校正、幾何校正、圖像匹配等問題需要加以考慮。
參考文獻
［1］ 侯舒維, 孫文方, 鄭小松. 遙感圖像雲檢測方法綜述［J］. 空間電子技術, 2014, 11(3).6877.
［2］ TAPAKIS R, CHARALAMBIDES A G. Equipment and methodologies for cloud detection and classification: a review［J］. Solar Energy, 2013,95(5): 392430.
則可以認為Pi與Pr相似。Pi的長度可以通過對Pi的求模操作即|Pi|來得到，而Pi與Pr間的光譜角可用如下數學運算式得到:αi=arccos(Pi·Pr/|Pi|·|Pr|)(1)
式中:αi為圖像像元向量與參考像元向量之間的夾角，即光譜角，代表了光譜回應的相似性，其變化範圍是［0，π/2］。Pi·Pr為向量點乘操作，|Pi|表示Pi的長度，|Pr|表示Pr的長度。光譜角度越小，同時向量長度越接近，則被估計像元的光譜回應與參考像元的光譜回應就越相似，參考物與被估計像元為同一類型的概率就越高。如圖2所示，假設圖像有3個譜段，分別用i，j，k表示，αi就是像元Pi與參考像元Pr之間的光譜角。
基於光譜角的雲檢測流程如圖3所示。
基於以上分析，提出基於光譜角的雲檢測演算法如下：
(1)計算遙感圖像某像元Pi與雲參考像元Pr的之間的光譜角，見式（1）。
(2)為了利用向量的長度資訊，構造一個包含像元Pi向量與雲參考像元Pr向量長度的高斯濾波器，對αi進行濾波，即乘以一個高斯權重係數GsCoef，該權重係數由式（2）得到，範圍是（0,1 ］。高斯係數的作用是給那些向量長度與參考像元向量長度接近的像元賦予較大的權重，直觀解釋就是向量長度相近的像元更有可能是同一類地物。濾波結果為Cgsi，其範圍是（0，1］。
GsCoef=e-((|Pi|-|Pr|)2/(0.35·|Pr|2))(2)
Cgsi=GsCoef*（π-αi）/π(3)
(3)取適當的閾值對Cgsi進行二值化，即當minThreshold
if (Cgsi>minThreshold && Cgsi
Ci=1;//該像元判斷為雲
else
Ci=0;//該像元判斷為非雲
(4)重複步驟(1)~(3)，依次完成所有像元檢測，生成雲掩膜。統計雲像元所占比例，為後續處理做好準備。
2實驗結果
為了便於比較，本文利用開源遙感影像處理軟體Grass7.0提供的自動雲評估ACCA演算法C++代碼，與文中的光譜角雲檢測演算法C++代碼，從整體檢測精度（Cloud Overall Accuracy, COA）以及執行速度上進行比較。計算平臺為：Lenovo Windows 7 64 bit PC; Intel Core i52520 @2.5 GHz CPU; 4.00 GB RAM。
從美國地質調查局網站得到156景經過人工標注的Landsat7圖像，作為本文雲檢測精度評價的標準圖像。此156景資料分佈在9個不同的緯度區域，對不同地域的雲分佈有一定統計意義，同時它們是文獻［3］中所採用的212幅人工標注圖像的子集。在這156景8 bit量化的遙感圖像中，圖像專業分析人員利用Photoshop中的圖像套索工具一點點地分離雲和非雲像元。
綜合考慮檢測精度和運算速度，本文取Landsat7中的Band1(450~515 nm)，Band3(630~690 nm)，Band4(750~900 nm)，Band5(1 550~1 750 nm)作為試驗對象，利用這四個譜段構造光譜角進行雲檢測。同時，為了更好地模擬在軌即時處理過程，利用從原始圖像截取的1 024像元×1 024像元大小的含不同雲量的圖像進行試驗。雲參考向量Pr取(225,215,182,168)，minThreshold取0.6，maxThreshold取1，比較結果如圖4~6所示。
對於參考圖像，ACCA演算法耗時較長，平均耗時為2 299.6 ms。耗時曲線中顯示對不同的圖像，耗時相差較大，這是因為ACCA演算法針對不同的地物採用不同的演算法分支，而某些分支中包含有迴圈處理。SA執行速度很快，平均耗時僅為361.1 ms。
圖7是帶雲的原始圖像，圖8是人工標注的雲圖像，圖中，雲標注為全白，非雲地物為全黑。典型的ACCA雲檢測結果如圖9所示，SA檢測結果如圖10所示。整體檢測精度分佈如表2所示，可知ACCA的整體檢測精度達到89.05%，高於SA演算法的83.6%。對檢測誤差較大的圖像進行分析發現，高緯度地區覆蓋的冰和雪常被誤檢為雲，因為冰雪反射率與雲接近，且具有較高亮度；對薄雲的漏檢在兩種演算法中均存在，SA演算法相對較嚴重。
3結論
分析可知，SA的檢測精度略低於ACCA，但是仍然高於83%。最為關鍵的是，SA演算法執行時間僅為ACCA演算法執行時間的1/6，這在即時影像處理中尤為重要。對於高緯度地區SA演算法的過檢問題，可以針對高緯度地區的冰雪特點，適當調高閾值，提高檢測精度。在軌應用時，一些實際問題例如圖像的輻射校正、幾何校正、圖像匹配等問題需要加以考慮。
參考文獻
［1］ 侯舒維, 孫文方, 鄭小松. 遙感圖像雲檢測方法綜述［J］. 空間電子技術, 2014, 11(3).6877.
［2］ TAPAKIS R, CHARALAMBIDES A G. Equipment and methodologies for cloud detection and classification: a review［J］. Solar Energy, 2013,95(5): 392430.