胡鐵瀧,蔣良群,王 杰

(西華師范大學(xué) 國土資源學(xué)院,四川 南充 637009)

植被是連接土壤、大氣和水分的紐帶,具有減少雨滴擊濺、減緩地表徑流、增加土壤保固等功能[1]。植被覆蓋度是指植被的葉、莖、枝等在地面的垂直投影面積占統(tǒng)計(jì)區(qū)總面積的百分比[2]，通常用植被覆蓋度來表示地表植被的覆蓋情況。目前,獲取植被覆蓋度的方式有地面監(jiān)測和遙感監(jiān)測,其中地面人工監(jiān)測耗時(shí)耗力,因此通過遙感數(shù)據(jù)反演植被覆蓋信息的方法被廣泛應(yīng)用[3]。南方丘陵地區(qū)常年多云霧,光學(xué)影像受到云霧的污染,很難準(zhǔn)確反演植被覆蓋度,因此利用光學(xué)影像進(jìn)行植被覆蓋度的反演時(shí),通常要考慮云霧的影響。合成孔徑雷達(dá)能適應(yīng)各種惡劣的天氣全天候成像,同時(shí)由于波長較長能穿透云霧,可較清晰地呈現(xiàn)地表的景觀格局,被廣泛地應(yīng)用于生產(chǎn)實(shí)踐之中[4]。雖然合成孔徑雷達(dá)具有較大的優(yōu)勢,但相對于目前在軌的多光譜資源衛(wèi)星如Landsat 8[5]、Sentinel 2A[6,7]等,極化組合的波段較少、斑點(diǎn)噪聲多、人工目視解譯難,涉及較復(fù)雜的輻射傳輸模型[8],大大限制了合成孔徑雷達(dá)的廣泛利用。光學(xué)與SAR影像各自具有優(yōu)勢,如何有效地集合它們的優(yōu)勢是目前研究的熱點(diǎn)[9,10]。通常情況下,這類算法需要多時(shí)相的合成孔徑雷達(dá)數(shù)據(jù),但會涉及到變化檢測的問題,增加植被覆蓋度提取的難度和不確定性。單時(shí)相的SAR與光學(xué)數(shù)據(jù)的融合研究非常具有科研價(jià)值,它簡化了融合算法的復(fù)雜程度,可兼顧時(shí)效性。植物的光譜曲線隨著時(shí)間不斷變化,同時(shí)由于地形、光照、水分等因素的影響,植被光譜的變化較大,增加了植被覆蓋度提取的難度。如何有效地減弱地物光譜的變化,近幾年大量的光譜變化混合分析模型被提出[11-17]，大致分為兩類,即確定性模型與統(tǒng)計(jì)性模型。雖然光譜變化混合分析模型能解決端元間的變化,但是對相似度較高的端元內(nèi)的光譜變化大都很難解決[18,19]。

針對上述研究中的問題,本文基于一景Sentinel 1A和Sentinel 2A影像,避免多時(shí)相數(shù)據(jù)的復(fù)雜性與時(shí)效性,利用全局光譜最小距離法去除Sentinel 2A影像的云,同時(shí)恢復(fù)Sentinel 2A云覆蓋下的地物光譜曲線,然后采用光譜歸一化技術(shù)減弱地物內(nèi)的光譜變化,在此基礎(chǔ)上采用擴(kuò)展線性混合模型反演植被覆蓋度。

1 研究區(qū)概況

本文選擇夏季植被覆蓋茂盛、地形復(fù)雜的四川丘陵區(qū)作為研究區(qū)。南充市位于四川省東北部,地處嘉陵江中游,屬于亞熱帶濕潤季風(fēng)氣候區(qū),四季分明、雨熱同期。研究區(qū)位于四川省南充市,屬丘陵地貌,介于30°74′—31°9′N、106°08′—106°35′E。研究區(qū)夏季植被覆蓋茂盛、類別復(fù)雜,光學(xué)遙感影像常年多云覆蓋,適合研究光學(xué)遙感影像的去云和植被覆蓋度的反演。

2 數(shù)據(jù)源及預(yù)處理

數(shù)據(jù)源采用Sentinel 1A和Sentinel 2A影像,其中Sentinel 1A影像為干涉寬模式(IW)下的地距產(chǎn)品(GRDH),而Sentinel 2A影像為L1C格式,需要進(jìn)行大氣校正。從歐洲航天局(ESA，https://scihub.copernicus.eu/dhus/)獲取南充市附近2018年7月15日的Sentinel 1A、2018年7月21日的Sentinel 2A遙感影像。由于Sentinel-1A和Sentinel-2A兩顆衛(wèi)星過境時(shí)間不一致,夏季沒有時(shí)間點(diǎn)完全一致且能利用的影像,加之在7月份相差6天的時(shí)間對影像的影響較小,最終選定了這兩景影像。本文的數(shù)據(jù)處理流程見圖1。圖1的右部分為Sentinel 2A影像去云處理(結(jié)合Sentinel 1A影像),左部分為Sentinel 2A影像的光譜歸一化與光譜混合分析,從而獲取植被覆蓋度,并進(jìn)行精度評價(jià)。

圖1 數(shù)據(jù)處理流程圖

Sentinel 2A影像的預(yù)處理主要通過ENVI 5.3和ESA提供的SNAP 6.0(Sentinel Application Platform)軟件完成。首先，運(yùn)用SNAP軟件的Sen2cor模塊對Sentinel 2A L1C數(shù)據(jù)進(jìn)行大氣校正處理,輸出10m、20m、60m分辨率的11個(gè)波段(Band 1、Band 2、Band 3、Band 4、Band 5、Band 6、Band 7、Band 8、Band 8a、Band 11和Band 12)的L2A地表反射率數(shù)據(jù)，利用空間分辨率增強(qiáng)工具(Sentinel-2 super-resolution)將L2A所有波段增強(qiáng)到10m,使影像的波段保持相同的分辨率;然后將增強(qiáng)的影像重采樣至20m,以保持與Sentinel 1A預(yù)處理后的空間分辨率一致;最后利用ENVI5.3將11個(gè)波段數(shù)據(jù)合成多光譜影像。

運(yùn)用SNAP軟件對Sentinel 1A影像進(jìn)行預(yù)處理。由于地形導(dǎo)致輻射的偏差,本次試驗(yàn)對Sentinel 1A數(shù)據(jù)采用地形輻射校正。預(yù)處理步驟:①應(yīng)用軌道文件(Apply Orbit File)自動下載精確的軌道信息文件。②對影像進(jìn)行輻射校正(Calibrate)。輸出beta0波段,因后續(xù)地形輻射校正算法需要輸入beta0波段;對影像進(jìn)行斑點(diǎn)濾波(Frost Filter),消除影像上的噪聲。③進(jìn)行多視(Multilooking)處理。將影像重采樣到20m分辨率,在保持?jǐn)?shù)據(jù)有效的情況下,同時(shí)提高數(shù)據(jù)處理的效率,為Sentinel 2A去云處理做準(zhǔn)備。④進(jìn)行輻射地形校正(Radiometric Terrain Flattening)。自動下載30m的數(shù)字高程模型。⑤進(jìn)行距離多普勒地理編碼(Range-Doppler Terrain Correction),獲取準(zhǔn)確的地理位置。為提高后續(xù)的匹配運(yùn)算精度,分別生產(chǎn)VH與VV極化的輻射強(qiáng)度(db)波段。

為減弱兩類影像的配準(zhǔn)誤差,完成Sentinel 2A和Sentinel 1A的預(yù)處理后,利用SNAP軟件的地理位置配準(zhǔn)工具,以Sentinel 1A作為主影像,Sentinel 2A為輔影像,將兩類影像進(jìn)行精確地配準(zhǔn),為后續(xù)的去云算法做好數(shù)據(jù)準(zhǔn)備。在Sentinel 2A影像上裁剪出研究區(qū),選擇12、8、2這三個(gè)少云的波段合成影像，見圖2；其他波段大都呈現(xiàn)大量的云,見圖3顯示的波段1。Sentinel 2A影像經(jīng)過大氣校正后,產(chǎn)生云覆蓋的概率圖,將其轉(zhuǎn)換為二值影像,手動添加未被探測到的云污染區(qū)域,見圖4。同時(shí)，裁剪出相應(yīng)區(qū)域的Sentinel 1A影像,見圖5(VH,VV_db,VV合成)。從圖5可見,雖然Sentinel 1A影像經(jīng)過光譜濾波與多視處理后紋理非常清晰,但是影像斑噪仍然非常明顯。

圖2 Sentinel 2A影像彩色圖

圖3 Sentinel 2A波段1

圖4 Sentinel 2A影像云掩模圖

圖5 Sentinel 1A彩色合成圖

3 研究方法

3.1 Sentinel 2A影像去云處理

由于光學(xué)影像Sentinel 2A上存在大量的云,而微波影像Sentinel 1A無云,因此利用Sentinel 1A影像的光譜特征,將Sentinel 2A的云霧去除。具體算法流程為:首先，通過Sen2cor大氣輻射校正模型獲取Sentinel 2A影像的地表反射率和云掩模產(chǎn)品,在此基礎(chǔ)上,用云掩模產(chǎn)品(二值影像)將相應(yīng)的Sentinel 1A影像標(biāo)記成云與非云兩類像元；然后在整個(gè)Sentinel 1A影像空間采用歐式距離計(jì)算離每個(gè)云像元距離最短的非云像元的位置；最后根據(jù)在Sentinel 1A獲取的非云像元位置,在Sentinel 2A上找到相同位置的非云像元,再用找到的Sentinel 2A非云像元位置的光譜曲線代替距離最短的云像元的光譜。以上處理可有效結(jié)合兩類數(shù)據(jù)優(yōu)勢,達(dá)到Sentinel 2A影像去云的目的,為后續(xù)的端元變化解混準(zhǔn)備數(shù)據(jù)源。去云處理主要通過IDL代碼來實(shí)現(xiàn),由于篇幅的限制,文中沒有羅列去云的相關(guān)代碼。

去云的核心部分見圖6。首先在Sentinel 2A找到云像元的位置,在Sentinel 1A找到相同位置的像元,然后在整個(gè)Sentinel 1A影像空間采用歐氏距離計(jì)算出距離像元最小距離的非云像元的位置,再在Sentinel 2A上找到與位置相同的非云像元,最后用Sentinel 2A非云像元的光譜曲線替代有云像元的光譜。對Sentinel 2A有云的像元,運(yùn)用程序依次循環(huán)以上的過程,最終達(dá)到去云的效果。

圖6 去云流程圖

3.2 光譜歸一化

常見的光譜解混模型是首先固定端元組,然后應(yīng)用解混算法獲取各個(gè)地類的光譜曲線。如果端元組內(nèi)光譜的相關(guān)性過高、變化較大,那么選擇合理的端元難度增大,從而解混精度會降低[18,19],因此采用光譜歸一化技術(shù)減弱端元組內(nèi)的變化,具有十分重要的意義。光譜歸一化就是各個(gè)端元光譜除以它們各自的范數(shù),在不改變端元之間的相關(guān)性和它們在高維特征空間的相對位置的情況下,壓縮數(shù)據(jù)空間,減弱端元變化的絕對量,減弱端元組內(nèi)的差異。

3.3 光譜混合分析

由于空間分辨率的限制,常常導(dǎo)致植被與其他地物以混合像元形式呈現(xiàn)在影像上。通過前期的光譜歸一化技術(shù),雖然能降低植被內(nèi)的光譜變化,但是不能從根本上解決植被與其他地物的混合狀態(tài),因此混合像元分解技術(shù)被采用。在此次試驗(yàn)中,擴(kuò)展線性混合模型[15,20]被用于反演植被的覆蓋度。假設(shè)包含N個(gè)像元的L波段遙感圖像包含P個(gè)端元,用x∈RL×N、A∈RP×N、S∈RL×N、E∈RL×N分別表示像元矩陣、豐度矩陣、端元矩陣和誤差項(xiàng)。那么對單個(gè)像元,線性混合模型可表示為:

(1)

式中,xk(k=1…N)為一個(gè)像元;ek(k=1…N)為誤差向量。通過對豐度進(jìn)行限制,使其具有較強(qiáng)的物理意義,再增加兩個(gè)限制條件:

(2)

本文將式(1)與式(2)稱為全約束線性最小二乘(FCLS)算法[21]。當(dāng)考慮每個(gè)像元的端元變化時(shí),式(1)可改寫為:

(3)

當(dāng)考慮到光照因子為主要影響因素時(shí),有fpk(s0p)=Ψpksop。式中,sop為參考端元,式(3)可改寫為:

(4)

式中,φk=diag(φk),φk∈RP×P為值的對角矩陣,相當(dāng)于部分約束最小二乘法(SCLSU)[17]。式(4)只是將光照因子作為影響端元變化的主要因素,當(dāng)同時(shí)考慮像元間的空間鄰域信息和光譜變異問題時(shí),在保留ELMM框架下可寫為:

(5)

=λA(PHh(A)P2,1+PHv(A)P2,1)+LRP×N+(A)+μT(ATlP-lN)

(6)

(7)

4 結(jié)果及分析

利用Sentinel 1A和Sentinel 2A結(jié)合去云,得到的結(jié)果見圖7。從圖7可見,在上段河流中出現(xiàn)了大量的土壤像元,再結(jié)合圖5,發(fā)現(xiàn)Sentinel 1A預(yù)處理后的河流中也有很多強(qiáng)反射率的土壤像元,說明去云算法有效利用了Sentinel 1A影像的特征。雖然去云算法可能會導(dǎo)致恢復(fù)的圖像出現(xiàn)部分錯(cuò)誤,但是未影響植被覆蓋度的提取。從圖7可見,去云后的Sentinel 2A影像植被出現(xiàn)斑點(diǎn)狀且不連續(xù),主要是由Sentinel 1A影像的斑點(diǎn)狀引起。為了減弱Sentinel 1A影像的斑點(diǎn)對植被覆蓋度提取的影響,再對去云后的Sentinel 2A影像進(jìn)行中值濾波平滑,以減小Sentinel 2A影像噪聲的影響(圖8)。從圖8可見,中值濾波的影像更趨于平滑,壓制了部分噪聲。

圖7 去云后的Sentinel 2A影像

圖8 中值濾波后的Sentinel 2A影像

在進(jìn)行去云與中值濾波處理后,對Sentinel 2A影像進(jìn)行光譜歸一化,以減弱地物內(nèi)的光譜差異,結(jié)果見圖9。對比圖8和圖9可見,建筑物、土壤、水體、植被呈現(xiàn)大面積的一致性,說明光譜歸一化技術(shù)能減弱地類內(nèi)的光譜變化,但其中建筑物與土壤的光譜還存在一定的光譜變化。完成影像光譜歸一化后,進(jìn)行端元提取與光譜混合分析。經(jīng)過對研究區(qū)的實(shí)地調(diào)研,發(fā)現(xiàn)多個(gè)水塘有水生植物,可能會影響植被覆蓋度反演的精度,為此結(jié)合ENVI5.3軟件與以往的土地利用分類圖提取研究區(qū)的水體,掩膜掉Sentinel 2A影像上的水體,只評價(jià)陸地上的植被覆蓋度。端元提取算法采用了非降維的體積迭代[22],共提取了33個(gè)端元,并獲取了端元的位置,最終確定了5類端元的光譜曲線,見圖10。由于篇幅的限制,本文未對大量的植被光譜歸一化前后進(jìn)行對比,可參考相關(guān)文獻(xiàn)[18,19]。

圖9 光譜歸一化后的Sentinel 2A影像

圖10 歸一化前后所提端元的光譜曲線

采用ELMM模型定量植被的覆蓋度,其中ELMM算法需要數(shù)據(jù)的緊湊度參數(shù)、豐度系數(shù)參數(shù)和尺度比例參數(shù)。為了減弱隨機(jī)輸入?yún)?shù)的誤差,我們設(shè)置了20組參數(shù)的組合,取RMSE最小值對應(yīng)的參數(shù)作為最終參數(shù),將其解混結(jié)果作為ELMM算法的最終結(jié)果。試驗(yàn)發(fā)現(xiàn),其組合(λs=0.1,λA=0.001,λψ=0.0)最優(yōu),反演的結(jié)果見圖11(右)。為了驗(yàn)證植被覆蓋度解混的精度,以近幾年的土地利用分類圖和時(shí)間點(diǎn)相近的無云Sentinel 2A影像作為參考,使用Sentinel 2A數(shù)據(jù)(10m分辨率影像,SNAP軟件Sentinel 2A空間分辨率增強(qiáng)運(yùn)算)進(jìn)行神經(jīng)網(wǎng)絡(luò)分類,其中云與陰影覆蓋的區(qū)域采用手工矢量化改正,將分類的結(jié)果作為參考影像對解混的結(jié)果進(jìn)行精度評價(jià)(同時(shí)去掉建筑物與裸地等非植被類型),結(jié)果見圖11(左)。評價(jià)指標(biāo)采用均方根誤差(RMSE)、相關(guān)系數(shù)(R)。精度評價(jià)后其RMSE約為0.14,R的值約為0.95,由此可知去云后影像植被覆蓋度的反演結(jié)果精度較高,接近地表植被覆蓋的真實(shí)情況。

對框架各部分的程序運(yùn)行效率進(jìn)行分析,在電腦處理器為Inter(R) Xeon(R) CPU E5-2690 v4@2.6GHZ,運(yùn)行內(nèi)存128G和64位Windows10的操作系統(tǒng),以及MATLAB版本為2016a 64bit(主要ELMM模型的運(yùn)行)、ENVI 5.3 64bit(主要進(jìn)行去云模塊的實(shí)現(xiàn))下,對781×1305大小的Sentinel 2A影像進(jìn)行處理,其中ELMM模型用時(shí)大約470秒,而去云模塊用時(shí)大約為7.5h,說明逐像元循環(huán)去云模塊的效率極低。今后應(yīng)考慮多臺電腦進(jìn)行集群與分塊計(jì)算,以提高處理效率。

圖11 植被覆蓋度的參考影像與反演結(jié)果

5 結(jié)論與展望

本文針對南方光學(xué)影像多云的特征,提出了結(jié)合SAR去除光學(xué)影像上的云;針對丘陵地區(qū)復(fù)雜的地形、水分、光照等因素對植被的影響,對影像和端元進(jìn)行歸一化,以減弱因端元變化導(dǎo)致的解混誤差,并采用擴(kuò)展線性模型進(jìn)行植被覆蓋度反演和精度評價(jià)。通過以上分析,得出如下結(jié)論:基于去云的Sentinel 2A影像反演植被覆蓋度的效果較好,接近地表植被覆蓋的真實(shí)情況。由于Sentinel 1A影像的特點(diǎn),部分地物反演的面積不準(zhǔn)確,特別是較破碎化的裸地,但對植被的覆蓋度影響不大。本次試驗(yàn)僅采用了一景Sentinel 1A影像,以后可考慮多景Sentinel 1A數(shù)據(jù)進(jìn)行變化檢測,同時(shí)進(jìn)行多波段合并與配準(zhǔn),增加極化波段,提高Sentinel 1A的光譜匹配能力。為了提高程序的運(yùn)行效率,以后將數(shù)百臺計(jì)算機(jī)進(jìn)行集群與分塊計(jì)算。