本發(fā)明屬于遙感及影像處理技術(shù)領(lǐng)域，具體涉及一種批量遙感影像預(yù)處理的方法。

背景技術(shù)：

近年來，遙感影像在成像過程中受到天空中的云的影響，影響地物的判讀與解譯。對于遙感影像信息提取而言，含云影像將大大降低遙感影像的利用率，特別是給遙感影像預(yù)處理工作帶來極大的不便，加大了信息提取處理的難度，降低了遙感影像預(yù)處理的效率。

從上世紀(jì)九十年代以來，國內(nèi)不少遙感領(lǐng)域的科研工作者先后進(jìn)行了這方面的研究。采用閾值法、小波紋理分析法、面向?qū)ο蠓?、神?jīng)網(wǎng)絡(luò)法開展了一系列關(guān)于云檢測的工作。

如2010年采用閾值法，根據(jù)云和地物的反射率隨波長的變化而變化實(shí)現(xiàn)FORMOSAT-2,VENDS,LANDSAT and SENTINEL-2影像云地分離，但是不能保證閾值對每景影像上的云都能達(dá)到很好的識別效果。2007年采用無抽樣小波變換對全色遙感影像進(jìn)行薄云去除，針對單景遙感影像效果較好，對于批量遙感數(shù)據(jù)而言，效率不高，適應(yīng)性較弱。2011年采用小波SCM方法提取紋理特征識別云，取得了一定的成果，在云雪混淆的時候適應(yīng)性較好，但受到傳感器類型不同，同時由于影像上各種地物繁雜，空間分布并不集中，特別在影像地物類別較多的情況下，此法識別出的云的效果略差。2012年采用面向?qū)ο蟮腇mask算法檢測云，算法精度較高，但云和冰混淆的部分難以區(qū)分，特別是在批量遙感影像處理中，效率不高。采用支持向量機(jī)實(shí)現(xiàn)對遙感影像厚云的檢測與去除，并采用統(tǒng)計學(xué)補(bǔ)償?shù)姆椒ㄟM(jìn)行修復(fù)，精度相對較高，但是在薄云的檢測方面相對較弱。2015年采用BP神經(jīng)網(wǎng)絡(luò)法進(jìn)行云檢測，因其較強(qiáng)的學(xué)習(xí)能力，對樣本的理解與識別而致使云判效果較好，但是運(yùn)算量過大，不適合批量遙感影像處理。上述工作都不同程度地取得了一些成就，但工作方法相對單一，效果時好時壞。對云邊界識別精度不準(zhǔn)不僅導(dǎo)致了大量的無效影像被檢索，造成巨大浪費(fèi)。因此對批量遙感影像預(yù)處理云的檢測工作已成為當(dāng)務(wù)之急。

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

為解決上述問題，本發(fā)明公開了一種批量遙感影像預(yù)處理方法，包括以下步驟：

步驟1，批量獲取遙感影像并建立批量遙感影像數(shù)據(jù)庫，通過統(tǒng)計遙感影像的光譜特征和紋理特征獲取含云影像；

步驟2，檢測含云影像上的云域范圍及其含云量，去除含云量超過含云閾值的含云影像，得到剩余含云影像；

步驟3，去除步驟2中剩余含云影像上的云，得到無云影像；

步驟4，利用得到的無云影像，更新批量遙感影像數(shù)據(jù)庫。

進(jìn)一步地，步驟1中所述的通過統(tǒng)計遙感影像的光譜特征獲取含云影像是指：

步驟11，在批量遙感影像數(shù)據(jù)庫中任選一遙感影像作為當(dāng)前遙感影像，獲取該當(dāng)前遙感影像在紅波段的灰度頻率直方圖；

步驟12，若所述的灰度頻率直方圖中出現(xiàn)兩個峰值頻率，隨著灰度值的增大依次為第一峰值頻率和第二峰值頻率，且兩個峰值之間存在一個頻率相對平緩且低于兩個峰值的緩沖區(qū)，同時第二峰值頻率明顯高于第一峰值頻率，則將該第二峰值頻率對應(yīng)的灰度值標(biāo)記為灰度異常值且將該當(dāng)前遙感影像標(biāo)記為含云影像；否則，將該當(dāng)前遙感影像標(biāo)記為待定影像；

步驟13，重復(fù)步驟11-12，直至批量遙感影像數(shù)據(jù)庫中所有的遙感影像都被標(biāo)記為止。

進(jìn)一步地，步驟1中所述的通過統(tǒng)計遙感影像的紋理特征獲取含云影像是指：

步驟14，從待定影像中任取一個作為當(dāng)前待定影像，對當(dāng)前待定影像進(jìn)行直方圖均衡化處理，分別統(tǒng)計當(dāng)前待定影像在直方圖均衡化處理前后的均值和二階矩；

步驟15，若滿足且-0.03≤ASM_前-ASM_后≤0.03，則將該當(dāng)前待定影像標(biāo)記為含云影像，否則標(biāo)記為無云影像；

其中，Mean_后為待定影像均衡化處理后的均值，Mean_前為待定影像均衡化處理前的均值，ASM_前為待定影像均衡化處理前的二階矩，ASM_后為待定影像均衡化處理后的二階矩；

步驟16，重復(fù)步驟14、步驟15，直至所有的待定影像都被標(biāo)記為止。

進(jìn)一步地，步驟2中所述的檢測含云影像上的云域范圍及該云域范圍的含云量，去除含云量超過含云閾值的含云影像是指：

步驟21，通過建立決策樹云檢測規(guī)則依次檢測步驟2中獲取的所有含云影像中的云域范圍，所述的決策樹檢測規(guī)則如下：

其中，a為云在紅波段的灰度異常值，x為灰度值；

步驟22，分別統(tǒng)計含云影像和云域范圍的像素個數(shù)，通過下式計算云域范圍的含云量：

其中，C為含云量，N_c為云域范圍的像素個數(shù)，N_I為含云影像的像素個數(shù)；

步驟23，設(shè)定含云閾值D,若含云影像中的含云量C大于D時，去除該含云影像。

進(jìn)一步地，步驟21中所述的云在紅波段的灰度異常值a的計算方法為：

步驟211，對含云影像的灰度直方圖中的頻率做一階向前差分；

步驟212，設(shè)定閾值L，若前后頻率差分值＞L時，表明該頻率段不在緩沖區(qū)；否則，該頻率段為疑似緩沖區(qū)；

取最長的一段疑似緩沖區(qū)為頻率緩沖區(qū)；

步驟213，通過冒泡法得到含云影像的灰度直方圖上的可能異常值，計算頻率緩沖區(qū)的頻率均值w，若該可能異常值大于2w，則該可能異常值即為云在紅波段的灰度異常值a。

進(jìn)一步地，步驟3中所述的去除步驟2中剩余含云影像上的云，得到無云影像是指：

步驟31，依次選取含云影像作為含云影像A，從批量遙感影像數(shù)據(jù)庫中選取與該含云影像A的經(jīng)緯度相同但不同時相的所有無云影像,從所述無云影像中選擇質(zhì)量好且分辨率與該含云影像A相同的無云影像B；

步驟32，分別裁剪出所述含云影像A與無云影像B的云域區(qū)域C1、C2，用C2區(qū)域替換掉C1區(qū)域，鑲嵌得到中間影像A1；

步驟33，對中間影像A1采用均值濾波處理得到新的遙感影像A2，即無云影像；

進(jìn)一步地，步驟33中所述的對中間影像A1采用均值濾波處理得到新的遙感影像A2是指：

步驟331，選取中間影像A1中拼接線處的像素作為中心像素，并找到該中心像素8鄰域的8個像素，求出這8個像素和中心像素的均值作為新的中心像素；

步驟332，以該新的中心像素為中心，分別向內(nèi)外依次各緩沖5個像素；

步驟333，對中間影像A1中拼接線處所有像素都進(jìn)行步驟331、步驟332后，則得到新的遙感圖像A2。

與現(xiàn)有技術(shù)相比，本發(fā)明具有以下技術(shù)效果：

1.本發(fā)明快速高效，精度較高，提供了一種針對大批量遙感影像預(yù)處理的方法，極大程度上避免了數(shù)據(jù)過度浪費(fèi)，具有較強(qiáng)的實(shí)用價值和推廣意義；

2.本發(fā)明在影像預(yù)處理的過程中，確定了影像上的含云量，自動圈定了影像上的云域范圍并去除，對影像信息提取獲取數(shù)據(jù)給出了準(zhǔn)確的依據(jù)；

3.本發(fā)明采用了工作效率較高的決策樹方法，不僅效率高，而且人為干預(yù)少，減輕了人工云量判讀的工作量，縮短了云檢測和去除的工作時間，可用于批量選取遙感影像及遙感影像信息提取研究。

附圖說明

圖1為本發(fā)明總體流程圖；

圖2為本實(shí)施例中云在紅波段的灰度頻率直方圖；

圖3為原始初級云檢測規(guī)則；

圖4為資源三號衛(wèi)星影像常見的含云影像圖，其中(a)為含點(diǎn)云的影像圖，(b)為含卷云的影像圖，(c)為含厚云的影像圖，(d)為含薄云的影像圖；

圖5為本實(shí)施例經(jīng)過云檢測后的影像圖，其中(a)為經(jīng)過云檢測后含點(diǎn)云的影像圖，(b)為經(jīng)過云檢測后含卷云的影像圖，(c)為經(jīng)過云檢測后含厚云的影像圖，(d)為經(jīng)過云檢測后含薄云的影像圖；

圖6為本實(shí)施例中均值濾波處理的示意圖。

具體實(shí)施方式

下面結(jié)合附圖和實(shí)施例對本發(fā)明作進(jìn)一步的說明。

如圖1所示，本發(fā)明公開了一種批量遙感影像預(yù)處理方法，具體包括以下步驟：

步驟1，批量獲取遙感影像并建立批量遙感影像數(shù)據(jù)庫，通過統(tǒng)計遙感影像的光譜特征和紋理特征獲取含云影像；

根據(jù)遙感影像的軌道號及傳感器類型，分別建立遙感影像數(shù)據(jù)庫；

不同的傳感器往往導(dǎo)致遙感影像的分辨率不同，因此根據(jù)影像的軌道號及傳感器載荷分別建立遙感影像數(shù)據(jù)庫。如表1對資源三號衛(wèi)星影像按照不同的傳感器及軌道號建立關(guān)系型數(shù)據(jù)庫的方式存儲于影像數(shù)據(jù)庫中。

表1資源三號衛(wèi)星影像數(shù)據(jù)庫

不同的遙感影像上云的光譜特征及紋理特征表現(xiàn)相似。根據(jù)云在紅波段反射最明顯挑選出部分典型含云影像，并統(tǒng)計含云影像的光譜特征及紋理特征。

其中，通過統(tǒng)計遙感影像的光譜特征獲取含云影像是指：

步驟11，在批量遙感影像數(shù)據(jù)庫中任選一遙感影像作為當(dāng)前遙感影像，獲取該當(dāng)前遙感影像在紅波段的灰度頻率直方圖；

步驟12，若所述的灰度頻率直方圖中出現(xiàn)兩個峰值頻率，隨著灰度值的增大依次為第一峰值頻率和第二峰值頻率，且兩個峰值之間存在一個頻率相對平緩且低于兩個峰值的緩沖區(qū)，同時第二峰值頻率明顯高于第一峰值頻率，則將該第二峰值頻率對應(yīng)的灰度值標(biāo)記為灰度異常值且該當(dāng)前遙感影像標(biāo)記為含云影像；否則，將該當(dāng)前遙感影像標(biāo)記為待定影像；

步驟13，重復(fù)步驟11-12，直至批量遙感影像數(shù)據(jù)庫中所有的遙感影像都被標(biāo)記為止。

云的光譜特征可以通過灰度直方圖表征，云在紅波段反射是最明顯的，本實(shí)施例中，如圖2所示，無云影像的灰度直方圖呈現(xiàn)近似的標(biāo)準(zhǔn)正態(tài)分布。除薄云外，點(diǎn)云、厚云、卷云的灰度直方圖均會顯示出兩個明顯的峰值，隨著灰度值的增大依次為第一峰值頻率和第二峰值頻率，且兩個峰值之間存在一個頻率相對平緩且低于兩個峰值的緩沖區(qū)，同時第二峰值頻率明顯高于第一峰值頻率。本實(shí)施例中，如圖2所示，以2014年12月7日獲取，軌道號016170，Path/Row為894/127的資源三號衛(wèi)星影像紅波段的灰度直方圖，在通過對影像上的已有的云進(jìn)行灰度值統(tǒng)計，圖2中的圓圈圈定的峰值橫坐標(biāo)就是云在遙感影像上的光譜異常值a，平滑部分則為兩個峰值之間存在的緩沖區(qū)。

對于全色單波段影像，則可以直接獲取如圖2所示的灰度直方圖異常值，就可以標(biāo)記出含云影像。

其中，通過統(tǒng)計遙感影像的紋理特征獲取含云影像是指：

步驟14，從待定影像中任取一個作為當(dāng)前待定影像，對當(dāng)前待定影像進(jìn)行直方圖均衡化處理，分別統(tǒng)計待定影像在直方圖均衡化處理前后的均值和二階矩；

步驟15，若滿足且-0.03≤ASM_前-ASM_后≤0.03，則將該當(dāng)前待定影像標(biāo)記為含云影像，否則標(biāo)記為無云影像；

其中，Mean_后為待定影像均衡化處理后的均值，Mean_前為待定影像均衡化處理前的均值，ASM_前為待定影像均衡化處理前的二階矩，ASM_后為待定影像均衡化處理后的二階矩；

步驟16，重復(fù)步驟14、步驟15，直至所有的待定影像都被標(biāo)記為止。

利用灰度共生矩陣統(tǒng)計紋理特征是一種常見的方式，本發(fā)明為了突出遙感影像上的薄云，對待定影像進(jìn)行直方圖均衡化處理，突出影像上的云域亮度，減弱非云區(qū)域亮度?；叶裙采仃嚽昂笞兓?，分別統(tǒng)計了均值(Mean)和二階矩(ASM)，如表2所示：

表2含云影像灰度共生矩陣變化

從Mean角度來說，沙地由于較粗的紋理特征，均值逐步遞減，而云逐漸增大，利用均衡化前后的均值做比值就可以實(shí)現(xiàn)區(qū)分。如式(1)：

其中，“前”、“后”表示影像均衡化前后的均值。

同時，雪在均值和二階矩(ASM)方面前后差異較小，但從均衡化后，雪的二階矩低于云，采用如式(2)所示的差分處理可實(shí)現(xiàn)云雪易混淆的地物之間的區(qū)分。

-0.03≤ASM_前-ASM_后≤0.03 式(2)

同時滿足式(1)和式(2)的劃定為含云影像，否則為無云。

步驟2，確定含云影像上的云域范圍及該云域范圍的含云量，并劃定含云影像質(zhì)量等級；

其中，檢測含云影像上的云域范圍及該云域范圍的含云量，去除含云量超過含云閾值的含云影像是指：

步驟21，通過建立決策樹云檢測規(guī)則依次檢測步驟2中獲取的所有含云影像中的云域范圍，所述的決策樹檢測規(guī)則如下：

其中，a為云在紅波段的灰度異常值，x為灰度值；

其中，云在紅波段的灰度異常值a的計算方法為：

步驟211，對含云影像的灰度直方圖中的頻率做一階向前差分；

步驟212，設(shè)定閾值L，若前后頻率差分值＞L時，表明該頻率段不在緩沖區(qū)；否則，該頻率段為疑似緩沖區(qū)；

取最長的一段疑似緩沖區(qū)為頻率緩沖區(qū)；

步驟213，通過冒泡法得到含云影像的灰度直方圖上的可能異常值，計算頻率緩沖區(qū)的頻率均值w，若該可能異常值大于2w，則該可能異常值即為云在紅波段的灰度異常值a。

本實(shí)施例中最長的緩沖區(qū)即為圖2中的平滑部分(約為350-930)，求出這段的頻率均值為w。通過冒泡法得到曲線上的峰值異常值可能為103和938，只有938的頻率與w的差值大于2w，且938后的灰度與均值w沒有交點(diǎn)，因此938即為異常值。

步驟22，分別統(tǒng)計含云影像和云域范圍的像素個數(shù)，通過下式計算云域范圍的含云量：

其中，C為含云量，N_c為云域范圍的像素個數(shù)，N_I為含云影像的像素個數(shù)；

步驟23，設(shè)定含云閾值D，本實(shí)施例中D取50％,若含云影像中的含云量C大于D時，去除該含云影像。

如圖4、圖5所示分別是原始含云影像和檢測后的云域范圍，其中，圖5中白色的為云斑，黑色的為非云；

本實(shí)施例采用混淆矩陣可以快速評判檢測的精度，如表3所示。

表3云檢測后精度評價

從表3可知，Kappa系數(shù)均在0.93以上，說明具有較好的一致性，再次驗(yàn)證了本方案的精度。

主要根據(jù)含云量C的大小劃定質(zhì)量等級，主要包括優(yōu)秀(C≤0％)、良好(0％<C≤30％)、及格(30％<C≤50％)、廢片(C>50％)四等。對資源三號衛(wèi)星影像檢測了53軌共計4233景影像，檢測出含云影像579景影像，影像絕大部分處于良好及及格級別，說明總體成像效果較好，大部分可以直接進(jìn)行信息提取。對于含云影像，則需經(jīng)過去云才可再次使用。更新影像數(shù)據(jù)庫，得到表4。

表4更新后的影像數(shù)據(jù)庫

步驟3，去除含云影像上的云，并得到無云影像；

步驟31，依次選取劃分為優(yōu)秀、良好、及格等級的含云影像作為含云影像A，從批量遙感影像數(shù)據(jù)庫中選取與該含云影像A的經(jīng)緯度相同但不同時相的所有無云影像,從所述無云影像中選擇質(zhì)量好且分辨率與該含云影像A相同的無云影像B；

步驟32，分別裁剪出所述含云影像A與無云影像B的云域區(qū)域C1、C2，用C2區(qū)域替換掉C1區(qū)域，鑲嵌得到中間影像A1；

由于采用的影像都是具備地理坐標(biāo)的，經(jīng)過決策樹云檢測生成的矢量范圍C同樣具備地理坐標(biāo)，因此裁剪出的C1和C2具備相同的地理位置，不必進(jìn)行二次配準(zhǔn)；

步驟33，對中間影像A1采用均值濾波處理得到新的遙感影像A2，即無云影像；

通過步驟32得到的中間影像A1不能直接使用，會出現(xiàn)鑲嵌后明顯的拼接線，因此需要對A1拼接線處的灰度值進(jìn)行處理。本發(fā)明采用3×3窗口搜索拼接縫中心線處的像素周圍的9個柵格，如(x,y)處灰度值為g(x,y),取這9個柵格的像素灰度值的均值g′(x,y)為其灰度值。同時以中心線像素為準(zhǔn)，分別往內(nèi)外各緩沖5個像素，依次采用相同的方式進(jìn)行采樣，取均值作為新的像素灰度值，得到新的遙感影像A2，更新遙感影像數(shù)據(jù)庫。其中，中間邊界是拼接線，中間邊界內(nèi)外的兩個邊界是以中間邊界向內(nèi)向外5個像素進(jìn)行緩沖得到的緩沖區(qū)，表格中的322即是均值濾波后的像素灰度值,A1是替補(bǔ)后的影像，A2是均值濾波處理的圖像。

步驟4，重建批量遙感影像數(shù)據(jù)庫。

通過步驟3去云處理之后得到影像質(zhì)量相對較好的影像，按照表4所示對影像數(shù)據(jù)庫進(jìn)行重新編目，對所有影像重新劃分質(zhì)量等級,同時更新遙感影像數(shù)據(jù)庫。

本實(shí)施例中采用含有淺色礦物的遙感影像進(jìn)行去云處理，遙感影像上的淺色礦物，顏色淺，密度小。含二氧化硅和三氧化二鋁較高，如長石、白云母、石英等。影像上的云易與淺色礦物造成混淆，提取難度加大，而單景處理根本不適用于大批量遙感影像的生產(chǎn)，因此本發(fā)明先對影像上云進(jìn)行批量處理，大幅度提高淺色礦物的提取精度。

從遙感影像數(shù)據(jù)庫中提取淺色礦物，以硅化蝕變信息為例(主要是二氧化硅等氧化物)，首先利用光譜儀對野外采集到的礦物測定反射光譜曲線。采用波段比值運(yùn)算(如其中Band₄,band₃，band₂分別是近紅外波段、紅波段和綠波段)，求出各變量相關(guān)矩陣系數(shù)，得到如式(5)所示的線性回歸方程，即可提取出遙感影像上的淺色礦物。

其中，x，y分別是光譜特征向量和淺色礦物中二氧化硅的變量。