基于OTSU 在MODIS 數據中的云檢測應用?

解本巨孫 巖于龍振

（1.青島科技大學信息科學技術學院 青島 266100）（2.青島科技大學經濟與管理學院 青島 266100）

1 引言

知識的積累，推動人類對新領域的探索，不斷擴展數據獲取的渠道，人們通過發射衛星來得到更多的地球信息，以便更好地了解自然規律。MODIS是搭載在Terra 和Aqua 兩顆衛星上的中分辨率成像光譜儀［1］，通過觀察地球表面獲取相關信息，傳輸到地面進而對數據反演利用［2］。然而地球表面超過60%的部分都被云層覆蓋［3］，因此通過衛星取得的信息，就不得不面對云層帶來的干擾［4］，所以云檢測是充分利用衛星遙感信息的前提［5］，為地質觀察、氣象預報以及自然災害預警［6~9］等都提供極大的幫助。

我國地處中低緯度地區，受季風洋流影響較大，尤其在沿海區域，雨量充沛，上空更是常年籠罩著大量的云，對遙感信息的利用和轉化產生了極大阻力。目前，云檢測的方法主要分為兩大類［10］，一是比較常見的物理方法［11~12］，也稱之為閾值法，如何全軍等［13］提出了基于MODIS 數據的閾值的云檢測識別，該方法簡單易操作，但在不同季節不同地點，檢測效果有很大出入，并且在薄云和低云效果很差。又如馬浩等［14］在2018 年提出基于某一地區的多波段固定閾值云檢測，在對冬季的黃渤海地區進行云檢測能達到一個較好的效果，但是適用的時相性較差。另一類為模式識別方法，常見的有SVM統計學方法，人工神經網絡等，如瞿建華等［15］在2019 年提出基于深度學習的神經網絡云檢測方法，該方法準確率高，普適性強，能夠達到較好的效果，但該方法模型訓練周期長，使用的數據量多，并且對于數據集的選擇比較困難。本文在多波段閾值法研究的基礎上，融合改進后的OTSU 模型，實現閾值的自動獲取，同時也提高了檢測精度。

2 資料處理

MODIS 共有36 個中等分辨率的光譜波段［16］，每一到兩天對地球表面觀測一次，獲取陸地、海洋、汽溶膠、以及云層［17］等眾多信息。MODIS完全免費公開，可以從美國國家地質調查局地球探索者數據庫網站上獲取，根據自己的實際需求，選擇相應的產品類型。為了方便研究，本文采用MODIS 1B 級別產品，選擇在35.0°~42.0°緯度和117.0°~124.6°經度之間的區域數據，地理位置處于山東半島以及黃渤海區域，該地區四季分明，且下墊面種類豐富，具有典型的代表性，其中用到的波段信息如表1 所示。

表1 選用波段的相關信息

前四個波段，只需要計算其反射值，而第29 波段則需要計算每個像元的亮溫。本文選取2019 年11 月7 日、2020 年1 月5 日、2020 年2 月2 日以及2020 年4 月6 日 的 數 據 為 例，分 別 以20201107、20200105、20200202 和20200406 作為圖像在文中的編號，經過預處理后得到的影像如圖1所示。

圖1 預處理后的影像

3 云檢測方法

云檢測主要原理是根據云和下墊面如水、土地、雪、植被等在各個的波段中所產生的輻射率以及輻射亮度的不同，來進行區別檢測。云通常在可見光的范圍反射率較其他下墊面物體高很多，而在紅熱波段亮溫卻很低，充分利用云這一特性，就是云檢測的關鍵。以往的多波段閾值云檢測方法，通過人眼觀察，確定相應段數據的閾值，再綜合分析每一幅影像檢測效果得到統一數值，雖然方法實現簡單，但是在不同時間不同高度，云的輻射率和亮溫是有很大差異的，閾值的固定會使云檢測魯棒性差，具有極大缺陷。

云檢測的本質同樣是從背景信息中提取目標，不同的是在遙感數據中背景物體很多，并且各波段數據分布情況不同，因此需要進行一些轉換和處理。通過融合OTSU 思想，將遙感里的數據里的云看成前景，把非云看成背景，選擇合適的波段，滿足算法模型的需求，在多波段固定閾值法的基礎上改進，從而達到閾值的自動確定，減少主觀性誤差，適應不同季節不同類型的數據信息，以此實現影像分割。

3.1 多波段閾值方法

云檢測算法對閾值判斷分析如下：

1）在可見光波段波段1 中，植被反射率最低，然后是水、陸地、冰雪，它們的反射率在0.2上下，且數值之間相接近，分布狀況相似，可視為同類型數據，其與云的反射率有明顯差距，并且隨著云層的高度厚度，差距將更明顯，因此將該波段反射率的閾值設為T1，當B1大于此閾值時判定為云。

2）在MOODIS 數據的26 波段中，由于陸地，植被、水等下墊面的輻射將會被水汽吸收，衛星能接收的輻射數據只有冰雪和云，而高云和冰雪的輻射率值差距很大，因此可以用來檢測高云，將該波段反射率的閾值設為T2，當B26大于此閾值時判定為云。

3）云和雪在可見光范圍內的反射性比較相似，而在第6 波段處，云對太陽輻射的吸收能力較弱，反射率較高，而雪卻恰好相反。計算波段1 和波段6 的歸一化指數V1=（B1-B6）/（B1+B6），雪、云以及水這一指數是正數，而植被土壤是負數，利用V1>0 這個特性將它們分離出來，同時歸一化也可以消除一定的大氣輻射以及儀器影響。在波段29中，其他物體相比冰雪的亮溫較低，利用V2>C29將冰雪分離出來。在波段1 中，由于云的反射率明顯高于其他下墊面很多，且水的輻射率不會超過0.2，所以可以設定B1>0.2將水去除，經過以上操作后樣本中只剩云像元即可判斷為云。由于歸一化的這個設定可能會將海岸線、湖岸判定為云，因此再利用在波段8 中對云的出現很敏感的特性，設置B8>T3來進行對海岸線的誤判進行修正。

3.2 動態閾值實現

OTSU 算法也被稱為最大類間方差法，被譽為是圖像分割中確定閾值的最佳算法［19~20］。其根據圖片里像素點的灰度值不同，將圖片分成兩個部分，分別為前景（目標）和背景，當這兩部分的類間方差越大時，則可說明這兩部分的差別越大，以此進行圖像分割。

對于一個圖像I（a，b），其大小記為W*L，其灰度范圍為0～M-1，記i 為某點的灰度值，把分離前景和背景的閾值記為T，將i小于閾值T的所有像素占圖像的比例記為p0，平均灰度記為u0，大于閾值T的像素占圖像全部像素點的比例為p1，記其平均灰度為u1，圖像的總平均灰度記為u，類間方差記為V則有：

對于上述公式，當V 取得最大值時，能夠使得分割出的目標和背景都與圖像中心距離最遠，對應的T便是得到的閾值。但由于MODIS數據量大，云的邊界不明顯，傳統OTSU 算法對于這種情況分割后的結果不理想。除了均值外，還有灰度平均方差也能夠反映灰度分布，邊界附近的灰度值往往波動較大，而平均方差值更能反映其離散程度，因此這里采用平均方差來代替傳統OTSU 中的均值［21~22］，改進后的閾值T應滿足：

其中g2、g02、g12如下，Pi代表圖像中灰度值i 出現的概率。

由于圖像的灰度值是在0 到255 區間內的整數，在計算圖像的類間方差時只需直接帶入公式計算就可以，對于遙感數據里的反射率和亮溫并不在這個整數值區間內，所以需要進行標準化處理，從而映射到相應區間，適應改進后的OTSU 模型，對應公式如下，其中min是樣本中的最小值，max則為樣本中最大值。

4 實現過程

云檢測實現過程借助VS2016 開發工具，通過C#語言搭建云檢測平臺，使用IDL語言對遙感數據的進行處理，基于組件技術模型技術調用COM_IDL_connect 進行讀取IDL 編譯后的SAV 文件，引入ArcGIS 的插件MapControl 將處理后的數據轉換為為圖片顯示在頁面中，無需再通過GIS 查看地圖，具體實現步驟如下。

1）首先將遙感數據文件通過C#程序讀取，進行預處理，分別為輻射校正、幾何校正和區域裁剪。

2）數據經過預處理后，將生成TIF 文件，使用MapControl 插件載入該文件即可看到原始影像圖。使用IDL 讀取該文件，其數據主要是各個波段中的像元的輻射率或亮溫，通過HISTOGRAM 函數對這些數據進行密度的分析計算，均等分成256 個小區間，將標準化后的數據注入。

3）觀察統計后的數據，在分成的256 個小區間中，會發現在數據的首端或者末尾存在極少量像元“獨占”多個小區間情況，由于整個數據量很大，幾萬甚至幾十萬個，而區間大小是由數據的最值決定，這些點將會擴大區間的范圍，從而影響結果的精度，因此可以將這些點視為噪點進行過濾，從而減少誤差。經過多次實驗分析，當前端或者尾部小區間里落入的像元個數連續少于或者等于3，記錄這些像元，將其去除，會使得去云效果更好。采用設計的DIVIDI 函數進行噪點去除工作，具體內容為，第一步對數據首端遍歷，直到遇到的小區間內元素的個數值大于3，則終止，并且統計噪點個數為N1；第二步對數據的尾部進行倒遍歷，統計噪點個數為N2；第三步去掉數據前N1個，尾端的后N2個，再次將新數據通過HISTOGRAM 函數統計，重復上述步驟，經過多次迭代，直至N1和N2值均為0，此時為云檢測使用數據。

4）通過上述操作就可以得到不同波段里每一個像元對應反射率或者亮溫的分布情況，將數據輸入到融合OTSU 算法模型的OTSU_Threshold 函數中，動態計算得出每個波段對應的閾值，結果如表2。

表2 各個圖像閾值

5）將閾值輸入到云檢測模型中，利用其自帶的掩膜工具分離云和其他下墊面物體數據，將生成的掩膜圖像中云像元的值設為0，其他地方像元值不變，從而達到去云效果。

5 結果分析

5.1 檢測結果

為了更好地對比和檢驗效果，將固定閾值［14］與改進后OTSU 模型產生的動態閾值帶入云檢測模型處理，影像如圖2所示。

圖2 云檢測結果對比

共有四組不同日期的分圖，其中每個分圖包含的3 個小圖從從左到右分別為經過預處理后的影像、固定閾值和OTSU 自動閾值的云識別結果。觀察圖像可以看出，固定閾值法在不同日期檢測效果差異較大，其中在20200406 影像中漏云情況較大。固定閾值法中閾值確定是由人目視觀察多幅影像的去云效果得出，無法覆蓋全部類型云，導致其的桎梏性，對于邊界區域的薄云、碎云漏判現象嚴重。

5.2 精度驗證

本文從定量的角度進行評估檢測結果，分別對兩種方法的總體精度、精確率和召回率進行評價，其定義如表3所示。

表3 混淆矩陣

其中，正樣本代表云像素點，負樣本則代表非云點。總體精度計算方法如式（8），用來表示檢測結果中判斷正確云像元的比例，精確率計算如式（9），表示預測正樣本占所有預測樣本的比例，召回率計算如式（10），表示正樣本中被預測出的比例。

對于檢測后的圖像進行計算，得出的結果如表4所示。

通過表4 可以看出，經過改良后OTSU 動態閾值法較傳統的固定閾值法總體精度提升較大，對于漏判現象明顯改善，結合之前影像效果圖，可以看出同日期的低云薄云檢測效果明顯加強，尤其在云邊界地區。其根據不同的影像中所包含云類型狀況的不同計算出對應的閾值，更好地將云與非云區分出，在環渤海地區具有良好的適應性。

6 結語

該云檢測方法對每一個波段閾值的確定都能根據不同的遙感數據自動計算出其最佳閾值，不需人為觀察調整，能夠快捷準確地得出閾值，且經過改進后的OTSU 算法，能夠更好地將薄云區分開來，極大地提高了在黃渤海地區適用性，同時也滿足了遙感數據反演需求，為后續的海冰識別，天氣預警，環境監測等提供重要的基礎。同時本文也存在一些不足，當遙感信息數據極大時，可能計算時間過長，并且模型只使用到五個波段的遙感信息，沒有更大限度的利用遙感信息，且當出現極端的情況，即云像元占比極少的時，誤差比較大，這些在日后的科研工作中，都將繼續深入分析，不斷完善。