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

背景技術(shù)：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

附圖說(shuō)明

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

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

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

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

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

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

具體實(shí)施方式

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表3云檢測(cè)后精度評(píng)價(jià)

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

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

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

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

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

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

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

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

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

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

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

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

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

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