基于樸素貝葉斯方法的FY-4A/AGRI云檢測模型

鄢俊潔， 郭雪星， 瞿建華， 韓 旻

(北京華云星地通科技有限公司，北京 100081)

0 引言

云在地球大氣系統的輻射收支平衡中的作用是由云的光學性質決定的[1]。云層的遮擋會導致衛星影像出現噪音，從而使遙感反演的參數出現誤差，干擾后續的圖像分類或目標識別。因此，對云的準確檢測是遙感產品精準反演的前提。

傳統的云檢測方法主要利用云的高反射率和云頂的低溫等特點[2]，結合物理閾值法[3]或同態濾波去云法[4]等進行檢測。物理閾值法包括國際衛星云氣候計劃(international satellite cloud climatology project，ISCCP)云檢測方法[5-8]、高級甚高分辨率輻射計云、陸、海處理方案(AVHRR processing scheme over cloud, land and ocean，APOLLO)檢測算法[9-11]、美國國家海洋和大氣管理局云高級甚高分辨率輻射計(the NOAA cloud advanced very high reso-lution radiometer，CLAVR)云檢測算法[12-13]和二氧化碳(CO2)薄片法[14-15]。ISCCP法基于可見光和紅外窗區波段的數據，把每一個像元的觀測輻射值與晴空輻射值比較，若兩者的差大于晴空輻射值本身的變化時，判定該像元點是云點; APOLLO法利用高級甚高分辨率輻射計(advanced very high resolution radiometer, AVHRR)5個全分辨率探測通道資料，包括5個閾值檢測因子; CLAVR法關注的是極區，利用光譜及空間變化特征檢測云的存在; CO2薄片法主要針對全球高云覆蓋，包括薄卷云，利用對CO2敏感的紅外光譜輻射，把穿透性云從不透明的云和晴空中區分出來。閾值法簡單易行，但是對閾值的敏感性高，而閾值的選擇具有一定的主觀性，同時閾值法以像元為主要處理單元，沒有考慮遙感影像的結構信息[16]。同態濾波去云法將頻率過濾與灰度變化結合起來，分離云與背景地物，最終從影像中去除云的影響。這種方法由于涉及到濾波器以及截止頻率的選擇，在濾波過程中有時會導致一些有用信息的損失。隨著模式識別理論的發展應用，正演模式被引入到云檢測研究，如大氣輻射傳輸模式[17]，考慮了大氣溫濕和云微物理狀況，能夠正演模擬云頂亮溫[18]。傳統的云檢測算法需要大量的先驗知識，經過大量的人工判讀和調整，才能實現全球范圍的云判識。針對這一問題，機器學習方法被逐漸引入到云檢測研究中[19-20]。機器學習方法利用訓練樣本進行特征訓練，將先驗知識和人工判識的經驗提供給計算機，由計算機來綜合相關制約因素，完成回歸學習過程，同時使算法的適應性增強[21]。利用機器學習進行云檢測時，檢測效果通常取決于選取的特征描述能力的強弱，因此需要尋找一種好的描述特征[20]。

以往氣象衛星光學載荷的云檢測主要還是采用與美國地球觀測系統/中分辨率成像光譜儀(earth observing satellites/moderate resolution imaging spectroradiometer, EOS/MODIS)載荷類似的方法，通過多光譜在云上的反射和輻射特性進行云的識別，需要大量人為干預。據此，本文基于風云四號A星(FY-4A)多通道掃描成像輻射計(advanced geosynchronous radiation imager，AGRI)數據，提出了一種基于樸素貝葉斯方法的全自動云檢測模型，使用樸素貝葉斯算法作為核心結構，基于光學載荷基本云檢測原理選擇合適的紅外通道作為特性分類器參數，同時針對不同的地表類型和不同月份分別分類訓練構建，最終得到基于樸素貝葉斯算法的全自動云檢測模型，實現光學載荷影像云像元的高效、精準識別。

1 云檢測方法流程

基于樸素貝葉斯的云檢測方法從類別上分屬于統計模型，因而整體上分為模型訓練和產品推斷2個階段(圖1)。

圖1 樸素貝葉斯云檢測方法流程Fig.1 Methodology of naive bayes cloud detection algorithm

模型訓練階段的目標是生成云檢測特征概率的查找表，具體步驟如下： ① 先執行預處理，主要內容包括初始化準備工作，按照文件時間確定模型的月度標志，讀取L1b數據，按照下墊面歸集云和晴空像元的位置； ② 基于第一步歸集到各類下墊面的云和非云的像元，計算用于與云檢測相關的特征，在本研究中使用了6類特征，分別為T11，Tstd，Emiss4，Tmax-T，TD_11-85和GeoColor，3.3節對各類特征進行了詳細介紹。

產品推斷階段的目標是基于時次觀測的L1b產品文件推斷其云檢測產品。與模型訓練階段相同，推斷階段同樣執行了預處理操作，與之不同的是推斷階段依賴訓練階段的云檢測特征查找表來判識當前觀測的云像元，對每個像元的6種云屬特征，查找與特征值對應的機率，按照條件概率公式推斷云屬概率值，并按照下墊面不同云屬閾值得到最終的云產品。

2 創建數據集

FY-4氣象衛星是新一代靜止軌道定量遙感氣象衛星，其載荷AGRI通道數由風云二號氣象衛星的5個增加為14個(表1)，覆蓋了可見光、近紅外、短波紅外、中波紅外和長波紅外等波段。星上輻射定標精度為0.5 K、靈敏度為0.2 K、可見光空間分辨率為0.5 km，能實現云、氣溶膠、水汽、陸地表面特性、海洋水色等大氣、陸地、海洋參量的高精度定量反演[22]。

表1 FY-4A AGRI載荷通道 [23]Tab.1 Channel setting of FY-4A AGRI

3 樸素貝葉斯云檢測算法

3.1 算法原理

樸素貝葉斯是基于貝葉斯定理與特征條件獨立假設的分類方法。對于給定的訓練數據集，首先基于特征條件獨立假設學習輸入輸出的聯合概率分布。然后基于此模型對給定的輸入利用貝葉斯定理求出后驗概率最大的輸出。

具體地，設云特征為n維向量集合X={x1,x2,x3,...,xn}，輸出為云類別標記Y={0,1}。則云檢測問題描述為：

(1)

基于特征條件獨立假設，上式轉化為：

(2)

式中：X(i)為向量集合X中的第i維；xi為X中第i維的取值。

本文利用極大似然估計計算學習聯合概率分布P(X,Y)。假定訓練樣本集合容量為J，云屬先驗概率的極大似然估計為：

(3)

特征條件概率分布的極大似然估計為：

(4)

式中ajk為第j個樣本的第k個可能取值。

由此，對給定像元的云屬特征x，通過學習到的模型計算其后驗概率分布，根據后驗概率最大化準則判別其云屬類型。另外應云檢測產品的要求，本文對檢測結果按照概率閾值劃分了可能云和可能晴空2種中間類別。

3.2 地表覆蓋分類

由于云在不同下墊面類型上的特征差異顯著，本文將地球表面分為7類型(深海、淺水、陸地、積雪、北極、南極、沙漠)。對不同地表類型進行分類的目的是為了得到對晴空條件認識存在的系統性偏差，這種偏差因地表類型的不同而有很大差異。分類的輸入數據來自FY-4氣象衛星導航文件，圖2顯示了2019年1月和6月FY-4A全圓盤中下墊面類型的分布。

圖2 全球地表覆蓋分類

從圖2中可以看出，地表類型變化中最大的驅動因素是積雪的變化，本文中的積雪覆蓋使用近實時積雪覆蓋產品(near-real-time ice and snow extent, NISE)，積雪覆蓋產品包括由微波成像儀(SSM/I)生成的北半球和南半球海冰、積雪覆蓋范圍。積雪覆蓋產品以南、北半球2個25 km方位角、等面積投影(EASE grid)方式提供。首先，將EASE grid投影轉換為全球等經緯度投影，然后根據經緯度信息與FY-4A全圓盤4 km網格數據進行匹配，得到全圓盤范圍的積雪覆蓋信息。圖3顯示了2019年具體各月份積雪和海冰的空間覆蓋變化，從圖中可以看出，積雪和海冰的空間覆蓋度隨季節變化差異顯著。

圖3-1 2019年全球各月積雪覆蓋

圖3-2 2019年全球各月積雪覆蓋

3.3 云檢測分類器

提取良好的云屬特征是構建云檢測分類模型的關鍵步驟。本文在調研現有云檢測算法的基礎上，針對FY-4A衛星AGRI載荷的觀測特性選取了6種云相關特征：

1)T11： 11 μm亮溫是經由AGRI載荷觀測到的11 μm譜段輻亮度計算得到的亮度溫度xbt_11。因中、高層云在11 μm譜段的亮溫相較于地面低很多，因此可以直接使用該特征檢測云[24]。

2)Tstd： 11 μm亮溫局地標準差是為以目標像元為中心，局部窗口范圍內的所有像元亮溫的標準差。使用該特征能夠較為準確地區分中、高層云的邊緣，薄云和晴空背景。

3)Tmax-T： 11 μm亮溫局地最大值差被定義為以目標像元為中心局部窗口范圍內亮溫最大值與目標像元亮溫的差。與Tstd類似，中、高層云邊緣處的該特征值顯著高于晴空。

4)TD_11-85： 11 μm與8.5 μm亮溫差定義為目標像元11 μm亮溫與8.5 μm亮溫之差。由于冰和水的吸收峰值位于在窗區通道不同的波長[25]，有云像元處該特征值一般為正值，而晴空像元處該特征值則普遍小于0。

5)Emiss4： 偽4 μm發射率定義為：

(5)

式中： e為自然對數； c2為常數；xbt_11為11 μm亮溫；xbt_4為4 μm亮溫。圖4表現了FY-4A衛星AGRI載荷2020年1月1日晝夜時刻云與晴空發射率隨太陽天頂角的關系。分析圖4可以發現，云的偽4 μm發射率隨晝夜變化十分明顯。當太陽天頂角小于85°時(圖4(a))，云的偽4 μm發射率多數大于1.8，晴空則集中于1附近，可能云、可能晴空存在于兩者的過渡區間。而當太陽天頂角大于85°時(圖4(b))，云的偽4 μm發射率的下界降低到1以下，與晴空分布存在較大的重疊。因此，本研究選擇計算偽4 um發射率2種特征Emiss4_Day和Emiss4_Night，分別在太陽天頂角小于85°和大于85°時參與后驗概率計算。

圖4 偽4 μm發射率隨太陽天頂角的變化

6)GeoColor： GeoColor是一種由AGRI載荷多通道觀測值合成的3通道彩色圖像，具有良好的視覺表現力[26]。圖5是一幅彩色合成示例圖，其制作方法為： 在太陽天頂角小于90°時，采用0.46 μm、0.64 μm和0.86 μm反射率經過拉伸、變換得到白晝圖像； 當太陽天頂角大于90°時，采用3.9 μm和10.4 μm亮溫經由拉伸、組合、變換得到夜間圖像；按照太陽天頂角計算晨昏線位置處的漸變掩模，合并晝夜圖像得到最終的彩色合成圖。云和晴空在GeoColor合成圖中具有較為明顯的視覺差別，具體地，云表現為白色，晴空則表現為下墊面的色彩。由于云在GeoColor中的色彩是RGB通道共同決定的結果，且晝夜較為均勻，本研究決定在最大估計時學習通道值特征的聯合概率分布，并將其引入到樸素貝葉斯的框架中作為一種新的獨立特征。>

圖5 2020年1月1日10時 UTC FY-4A/AGRI GeoColor彩色合成圖Fig.5 Color composite image of FY-4A/ AGRIGeoColor at 10: 00 UTC January 1, 2020

3.4 模型構建

研究使用2019年12個月的FY-4A衛星AGRI載荷L1產品，定位(Geolocation， GEO)產品和云檢測(cloud mask，CLM)產品。其中，L1記錄譜段觀測值，GEO提供時間和地理相關信息，CLM作為云檢測的真實標記。每個月選取20%的時次構成訓練集，80%的時次構成測試集合，圖6是帶有集合時次數量的數據劃分統計情況。

圖6 數據劃分統計Fig.6 Data partition statistical chart

圖7為特征分布與云屬概率曲線修正圖。為了降低載荷觀測特征隨時間變化對云屬判斷產生的影響，更加準確地建立其云屬關系，本文以月度作為模型求取的最小時間尺度。研究求取云檢測單一特征月度模型的步驟描述如下：

1)在收集數據集中隨機選取20%的時次作為樣本集，按云屬和非云屬統計每個時次各下墊面上特征的頻數直方圖，如圖7(a)。

2)根據步驟1)得到頻數直方圖計算特征的云屬和非云屬條件概率分布，隨后按公式(4)求取特征的貝葉斯模型參數。

3)在統計樣本特征的分布時，區間邊緣的樣本總數較小，晴空與云的比例易受噪聲的干擾，云屬概率值有時會表現異常。因此需要對步驟2)得到的模型參數進行手動的修正，圖7(b)展示了修正后11 μm亮溫特征對應的云屬概率曲線。

3.5 云檢測類別劃分

在根據求解得到的模型檢測產品時，首先需要判斷像元的下墊面類型，然后計算算法定義的6類特征值，最后照式(2)融合多個特征的條件概率分布系數，由此可以得到該像元的云屬概率值。

現有主流的云檢測算法提供的是4類檢測結果，晴空、可能晴空、可能云和云，如EOS/MODIS，FY3D/MERSI和NOAA/AVHRR云檢測產品。為了生成類別產品，本研究選定0.5作為云與非云的概率閾值。為了得到較為準確的云和晴空類別，設定0.1的閾值劃分晴空和可能晴空，0.9的閾值劃分為可能云和云。

4 結果驗證分析

4.1 FY-4A/AGRI數據云檢測結果

以FY-4A業務CLM產品作為參照，本文采用分別針對云和晴空的召回率(probability of detection,POD)、誤判率(false alarm ratio,FAR)并綜合Kuipers評分(Kuiper’s skill score,KSS)來從多個角度反映算法的精度。假設參照為云的像元且實際被判別為云的像元數量為a，參照為云的像元但實際被判別為晴空的像元數量為b，參照為晴空的像元但實際被判別為云的像元數量為c，參照為晴空的像元且實際被判別為晴空的像元數量為d，算法對云判識的POD反映的是本文算法判別出的云占實際的云的比例，云的POD公式為：

(6)

云的FAR反映的是本文算法判別出的錯誤的云占判出的云的比例，公式為：

(7)

晴空的POD反映的是本文算法判別出的晴空占實際的云的比例，公式為：

(8)

晴空的FAR反映的是本文算法判別出的錯誤的晴空占實際的晴空的比例，公式為：

(9)

KSS是一種對錯誤分類敏感的補充性綜合指標，常用于算法傳統指標相近時更嚴格的評價，計算公式為：

(10)

2019年12個月測試集驗證統計結果如表2。從表中可以看出，本文方法云檢測的POD最高為6月份的98.9%，最低為9月份的90.4%。總體上，在2019年全年驗證集中，云的平均POD為97%，平均FAR為6.3%，晴空的平均POD為89.0%，平均FAR為2.9%，方法的平均KSS為87.4%。

表2 2019年12個月樸素貝葉斯以業務CLM為真值的交叉比對結果

圖8展示了測試集1月1日10時的云檢測結果。整體而言，貝葉斯云檢測算法與業務算法的結果大體一致。對比同時次的彩色合成圖(圖5)，細節上，貝葉斯算法在白天部分判斷出的晴空多于業務算法，可能晴空少于業務算法，與合成圖觀察到的結果更一致。在耀斑區，貝葉斯算法出現了部分中間狀態，業務算法表現更好。在澳大利亞中部，業務算法出現了錯判，而貝葉斯算法表現更好。在夜晚部分，業務算法在中國南海海域判云較多，而貝葉斯算法可能云較多，更符合實際云圖分布。

圖8 2020年1月1日10時UTC云檢測結果對比

4.2 模型分析

利用統計手段計算模型連續的參數值是本文算法區別于經典云檢測方法的最大不同之處。根據不同的特征量和下墊面信息從模型中獲得對應的有云概率值，并將不同特征量的概率值歸一化到0～1之間獲得云屬概率并根據分類閾值獲得云檢測結果。下面針對T11(11 μm亮溫)特征參數模型進行分析。

圖9為1月和6月7類下墊面類型所構建模型的T11特征概率曲線。其中，橫坐標AGRI載荷在11 μm譜段的觀測亮溫，縱坐標為經由式(2)計算得出的后驗云屬概率。11 μm通道為長波紅外通道，對于此通道，大氣組成成分對輻射的影響可以忽略，水汽只有微弱的吸收和再發射作用，衛星接收到的輻射主要是云和下墊面的發射輻射。根據普朗克定律，通常情況下高云比低云溫度低，低云比下墊面溫度低。但是由圖9可知，不同下墊面和季節，云的T11特征概率差別很大，在使用單一閾值進行區別時就需要特別注意。當下墊面為南極、北極和積雪時，云和下墊面差別較小，溫度很低的高厚云才容易被識別； 在不同季節，當下墊面為深海和淺海時，云與晴空的曲線變化十分陡峭，說明云的區分明顯，T11在此下墊面比較適用； 由于沙漠和積雪的下墊面主要分布在北半球，在不同季節時，沙漠和積雪的云屬概率差異很明顯，在1月時出現明顯的逆溫現象，這是因為此時沙漠和積雪的溫度有可能比云低。另外，綜合不同季節的模型和檢驗結果，在南極下墊面的有效情況樣本點過少，此下墊面的模型正確性待確認。

圖9 1月和6月各下墊面T11特征概率曲線

觀察圖10左側上部的耀斑區，可以發現本文算法與FY-4A CLM產品表現接近一致，都檢測到了成片的云區。但本文算法在亮度反應較為強烈的海洋上會產生一定程度地誤判，其原因是GeoColor合成圖中耀斑處。觀察左側中部的薄云區，可以發現本文算法在薄云邊緣處的判識精度仍有待提升。觀察右側中下部的海陸邊界處，可以發現FY-4A CLM產品沿著邊界線出現了較為明顯的誤判，而本文算法表現良好。同時，由于FY-4A CLM產品使用了12 μm亮溫，所以當譜段觀測出現缺失時會形成明顯的條帶。本文算法使用的譜段較少，且基于概率的條件獨立假設在信息缺失條件下仍能夠得出后驗云屬概率，因而可以有效避免條紋現象。另外，本文對12個月份分別構建模型，因而需要評估不同月度模型之間的差異對云檢測精度的影響。

圖10 2019年9月26日5時 UTC云檢測局地結果對比

圖11直觀展示了不同月度模型對云檢測精度的影響。在圖像的上部，6月份模型檢測結果出現了大量的錯判晴空。這主要是不同月度模型的下墊面劃分不同，出現晴空錯判的下墊面是積雪，而6月份積雪下墊面上云的特征與1月份有較大不同。2個模型在其他區域上的表現接近一致。

圖11 不同月度模型檢測差異

為了研究模型的月度差異，研究在1月的測試數據上定量對比基于1月和6月標記數據得到的2組模型，得到如表3的百分比混淆矩陣。分析表格可以發現，1月模型的云判識結果要高出6月模型3.84百分點，晴空甚至高出7百分點； 而明顯錯判的比例也明顯較低，例如把云判為晴空的比例低6.91百分點，把晴空判為云的比例低4.02百分點。說明使用同期月度模型推斷的云屬類別會更加準確。

表3 不同月度模型云檢測混淆矩陣統計Tab.3 Model cloud detection confusion matrix statistics for different months (%)

5 結論

本文針對FY-4A/AGRI云檢測問題，提出了一種基于樸素貝葉斯算法的全自動云檢測模型，實現了FY-4A/AGRI影像云像元的高效、精準識別。與已有的云檢測方法相比，該模型引入了樸素貝葉斯理論，將云檢測的閾值判斷轉化為概率評價，即不局限于是否有云，而是以云出現的概率做評價，評價結果更加合理。同時，新算法中涉及的數據只有紅外通道，較好地解決了傳統算法同時使用紅外通道和可見光數據造成的晨昏區云檢測分割線的問題。此外，該模型還引入了新的合成圖特征，學習聯合概率，取得了較好的云檢測效果。

總體而言，基于樸素貝葉斯算法的全自動云檢測模型展現出良好的識別光學載荷影像云像元的潛力，具有一定的理論意義和應用價值。不足之處在于使用的訓練集真值來自FY-4業務產品，模型精度依賴于產品數據的精度，因此未來可以引入測云衛星CALIPSO的數據作為真值進行模型學習訓練。另外，現有方法沒有考慮到特征隨太陽高度角和耀斑角的變化，未來可以加入M估計的方法，以改善分類評價效果。我們還計劃將該模型應用于風云三號E星光學載荷的云檢測，為解決氣象衛星晨昏軌道云檢測提供新的方法。