一種融合光學與SAR的植被指數(shù)提取算法

林用智，胡鐵瀧，李映潭，賀海明，王 杰

(西華師范大學地理科學學院，四川南充 637009)

0 引言

植被指數(shù)是表征地表植被覆蓋的簡單有效的度量參數(shù)，被廣泛用于植被類型監(jiān)測[1]、葉面積指數(shù)監(jiān)測[2-3]、覆蓋度監(jiān)測[4]和生物量估算[5]等研究.常見植被指數(shù)有差值植被指數(shù)(DVI)[6]、歸一化植被指數(shù)(NDVI)[7]、增強型植被指數(shù)(EVI)[8-9]和土壤調(diào)節(jié)植被指數(shù)(SAVI)[10]等.植被指數(shù)主要是通過光學影像的近紅外波段和紅波段反射率而獲取[11]，Sentinel-2影像被廣泛應用于植被指數(shù)的提取與作物信息監(jiān)測.Sentinel-2影像為多光譜影像，有13個波段，空間分辨率分為10、20、60 m，時間分辨率為5 d，Sentinel-2影像在紅邊范圍有3個波段(Band 5、Band 6、Band 7)，3個紅邊波段能有效地監(jiān)測植被信息[12]，然而，Sentinel-2影像易受云霧污染，云層的存在導致傳感器無法獲取地物反射信息，使植被指數(shù)獲取較為困難，因此獲取植被指數(shù)時，要考慮去除云霧的影響.傳統(tǒng)的去云方法通常分為兩大類：多時相的光學影像去云和SAR與光學融合去云.國內(nèi)外學者在去除光學影像中的云霧、獲取云層下的植被指數(shù)方面做了大量研究，Cerra等人首次將基于稀疏解混的去噪應用于Sentinel-2影像的去云，該算法的速度較快、效率較高[13]；合成孔徑雷達具有全天候、全天時監(jiān)測的特點，其波長較長，受云霧影響較小，且微波能穿透植被，獲取植被的種植結(jié)構(gòu)信息，使其在植被監(jiān)測中有一定優(yōu)勢.光學和SAR數(shù)據(jù)各自具有優(yōu)勢，有效地結(jié)合它們的優(yōu)勢能提高數(shù)據(jù)可用性和協(xié)作性，因此有學者利用SAR與光學融合去云[14-15]，以機器學習去云為主，機器學習去云包括回歸方法和深度學習去云算法，Meraner等人設計了一種融合SAR和光學數(shù)據(jù)的深度殘差神經(jīng)網(wǎng)絡構(gòu)架去除Sentinel-2影像中的云霧，并最大限度地保留了原始信息[16]；Filgueiras等人利用8種回歸算法，結(jié)合Sentinel-1和Sentinel-2影像，驗證了4種建模方法，得出8種回歸算法中，隨機森林方法獲取的精度最好[17].

以上研究中，基于多時相的光學影像去云算法需要多景影像數(shù)據(jù)，影像之間存在的時間間隔，加上降雨、光照等因素的影響，植被的光譜曲線會發(fā)生變化，導致多時相去云算法的誤差；機器學習算法雖然能最大程度減小光學與SAR的成像時間間隔，但是其結(jié)果通常會丟失部分光譜特征，或者獲取的光譜曲線發(fā)生巨大的變化，從而導致植被指數(shù)獲取呈現(xiàn)巨大的誤差.如何有效地去除光學影像的云像元并保留云下像元的光譜特征是計算有云影像植被指數(shù)的關(guān)鍵.針對上述問題，本文利用Sentinel-1影像構(gòu)建出6種匹配模型，分別對6種匹配模型采用全局光譜最小距離法去除Sentinel-2影像上的云與陰影像元，獲取6景去云后的Sentinel-2影像，恢復云下地物的光譜曲線，然后獲取植被指數(shù)，最后驗證去云與陰影區(qū)域的植被指數(shù)提取的精度.

1 數(shù)據(jù)與方法

1.1 試驗區(qū)概況

研究區(qū)位于新疆維吾爾自治區(qū)烏蘇市(44°27′～44°39′N，84°29′～84°45′E)，屬溫帶大陸性氣候，四季分明，地形平坦，主要農(nóng)作物為棉花、小麥和玉米等，且種植面積廣，植被覆蓋面積較大.因而本研究選擇該區(qū)域作為研究區(qū).

1.2 數(shù)據(jù)獲取及處理

本文的研究數(shù)據(jù)來自ESA數(shù)據(jù)共享網(wǎng)站(https://scihub.copernicus.eu/dhus/#/home)，包括SAR數(shù)據(jù)和光學數(shù)據(jù).SAR數(shù)據(jù)為2020年9月11日的Sentinel-1影像GRD數(shù)據(jù)；光學影像為2020年9月3日無云的Sentinel-2B和2020年9月8日有云的Sentinel-2A影像.

本文選擇的Sentinel-2影像為L1C數(shù)據(jù).使用SNAP 8.0和ENVI 5.3對數(shù)據(jù)進行預處理，步驟如下：(1)用Sen2cor模塊對影像進行大氣校正；(2)再用Sentinel-2 super-resolution模塊進行空間分辨率增強；(3)在ENVI中合成多光譜影像.由于Sentinel-1A、Sentinel-2A兩顆衛(wèi)星過境時間不一致，且研究區(qū)的夏季沒有時間點完全一致且能使用的影像，而有云的Sentinel-2影像時間為2020年9月8日，SAR影像和光學影像相差3天的時間對影像影響較小，所以本文采用2020年9月11日的Sentinel-1影像.用SNAP軟件對SAR影像進行預處理：(1)應用軌道文件；(2)輻射校正，輸出sigma0波段；(3)斑點濾波；(4)多普勒地理編碼，輸出本地入射角(LIA，Local Incidence Angle)；(5)同時生產(chǎn)雙極化的dB波段；(6)最后將Sentinel-1A、Sentinel-2A影像進行配準.

圖1 Sentinel影像Fig.1 Sentinel images

配準完成后，裁剪出研究區(qū)域的Sentinel影像.如圖1所示，(a)為無云的Sentinel-2B影像(RGB：8A，4，3)，(b)為有云的Sentinel-2A影像(RGB：8A，4，3)，(c)為Sentinel-1影像(VV，VH，VV/VH合成).Sentinel-2A影像經(jīng)大氣校正后，產(chǎn)生云覆蓋的概率圖，將其轉(zhuǎn)換為二值影像，再手動添加未被探測的被云污染的區(qū)域，如圖2.

1.3 研究方法

1.3.1 Sentinel-2A影像去云處理 胡鐵瀧[18]等人通過Sentinel-2A影像去云方法有效地結(jié)合了SAR數(shù)據(jù)和光學數(shù)據(jù)的優(yōu)勢，達到Sentinel-2A影像去云與陰影的目的，如圖3所示，其基本流程如下：通過預處理后的Sentinel-2A影像獲取云掩膜圖，利用云掩膜圖將研究區(qū)的Sentinel-1A影像分為云和非云兩類像元，然后采用歐式距離獲取Sentinel-1A影像中離每個云像元(X′)距離最近的非云像元(Y′)的位置，最后根據(jù)獲取的非云像元位置，在Sentinel-2A上找到相同位置的非云像元(Y)，再用非云像元的光譜曲線代替距離最近的云像元(X)的光譜.通過以上的算法，尋找Sentinel-1A影像上的最相近的位置來替代Sentinel-2A影像的云光譜曲線，能夠有效地恢復地表物質(zhì)的光譜曲線，相比機器學習算法，其更具有物理意義，因為其影像上大量存在相似光譜曲線.

圖2 云與陰影的感興趣區(qū)域Fig.2 Clouds and shadows region of interest圖3 去云與陰影流程圖Fig.3 Flow figure of removal cloud and shadow

以上算法的核心是實現(xiàn)Sentinel-1A影像上像元的正確匹配，Sentinel-1A影像只有兩個極化波段，分別為VV與VH，由于類內(nèi)光譜曲線呈現(xiàn)較大的變化[19-21]，這樣會增大匹配的誤差，從而導致后期的植被指數(shù)的計算精度較低.為了增加Sentinel-1A影像像元的匹配精度，需要考慮增加波段信息，為此，本文通過不同的特征運算，增加了數(shù)個特征，同時對這些特征進行組合，產(chǎn)生六種不同的匹配模型，驗證不同波段模型組合的匹配精度.其中，增加的特征運算包括：[1]VV與VH極化波段的dB值計算，如式(1)所示；[2]Sentinel-1A影像的比值植被指數(shù)(RVI)[22-23]，如式(2)所示；[3]Sentinel-1A影像歸一化比值指數(shù)(NRPB)[17]，如式(3)所示；[4]Sentinel-1A影像本地入射角(LIA)，這個波段是對Sentinel-1A影像地理編碼后生成的，將其導入即可[17].

(1)

(2)

(3)

1.3.2 植被指數(shù) 獲取去云與陰影影像之后，計算云與陰影覆蓋區(qū)域的植被指數(shù)，包括歸一化植被指數(shù)(NDVI)、土壤調(diào)節(jié)植被指數(shù)(SAVI)和增強型植被指數(shù)(EVI).NDVI是根據(jù)植被在近紅外波段和紅波段的梯度差得出的公式，NDVI的優(yōu)勢在于能消除大部分地形和大氣產(chǎn)生的影響，但易受土壤等因素的干擾，因此，Huete[10]等人在此基礎上考慮了土壤的影響，加入了土壤調(diào)節(jié)參數(shù)，提出了SAVI.同時，LIU[9]等人為了解決NDVI存在的問題，提出了EVI，該指數(shù)引入的土壤調(diào)節(jié)參數(shù)和大氣修正參數(shù)減少了土壤和大氣的影響.上述植被指數(shù)計算公式為

(4)

(5)

(6)

其中：ρNIR和ρred分別為近紅外波段反射率和紅波段反射率；ρblue為藍波段反射率；C1為大氣修正紅光校正參數(shù)，取6.0；C2為大氣修正藍光校正參數(shù)，取7.5；L為土壤調(diào)節(jié)參數(shù)，取1[22].

1.3.3 精度評價 本文采用相關(guān)系數(shù)(R)、平均絕對誤差(MAE)和均方根誤差(RMSE)3個指標作為模型精度的驗證[24-26].

(7)

(8)

(9)

2 結(jié)果與分析

圖4 無云的Sentinel-2B影像計算的3種植被指數(shù)Fig.4 Cloud-free Sentinel-2B image calculated three kinds of vegetation indices

為了驗證6個匹配模型所計算的植被指數(shù)的精度，使用2020年9月3日無云的Sentinel-2B影像計算出的EVI、NDVI和SAVI作為真值，其結(jié)果如圖4.結(jié)合去云算法和6個匹配模型的光譜特征去除Sentinel-2影像上的云與陰影像元，獲取6景去云與陰影后的Sentinel-2影像，如圖5所示，(a)～(f)分別與匹配模型的01～06編號相對應.采用相關(guān)系數(shù)、均方根誤差和平均絕對誤差對植被指數(shù)的整體精度和6個匹配模型的精度進行評價.首先評價去云算法，對比圖1(a)和圖5(a)～(f)可以看出，去云算法能較好地修復受云與陰影污染的像元，恢復地物的光譜曲線，清晰地呈現(xiàn)出地塊的基本形狀.觀察發(fā)現(xiàn)，圖5綠色方框中的水體區(qū)域恢復出了部分土壤的像元，這與實際情況不符合，分析圖1的Sentinel-1A影像之后，發(fā)現(xiàn)其有部分土壤的信號，主要原因是水體較淺且渾濁(最近有一次強的降雨)，導致后向散射系數(shù)增大所致，后期可結(jié)合前期的Sentinel-2A影像，探測出相關(guān)的水體，然后掩膜掉，以提高恢復影像的質(zhì)量.

圖5 去云與陰影后的Sentinel-2A影像Fig.5 The Sentinel-2A image of removal cloud and shadow圖6 6種匹配模型計算的EVIFig.6 EVI calculated by six kinds of matching models

對比圖4(a)、圖4(c)、圖6、圖8和表1分析發(fā)現(xiàn)，在3種植被指數(shù)的整體精度中，云與陰影覆蓋區(qū)域所計算的EVI和SAVI與真值較為接近，且圖像色差較小.與真值相比，SAVI整體的均方根誤差和平均絕對誤差分別小于0.092 71和0.071 40，EVI整體的均方根誤差和平均絕對誤差分別小于0.124 31和0.093 24，SAVI和EVI的相關(guān)系數(shù)均大于0.8，說明SAVI的整體精度高于EVI.結(jié)合圖4(b)和圖7可以看出，云與陰影覆蓋區(qū)域計算的NDVI與真值的圖像相比，顏色略有差異，同時分析表1可以得出，去云與陰影后計算的NDVI的RMSE和MAE分別大于0.14和0.11，研究氣象資料發(fā)現(xiàn)，9月3日至9月8日有降雨，對植被的影響較小，但降雨對土壤影響較大，因此，沒有引入土壤調(diào)節(jié)參數(shù)的NDVI的精度要略低于SAVI和EVI.

圖7 6種匹配模型計算的NDVIFig.7 NDVI calculated by six kinds of matching models圖8 6種匹配模型計算的SAVIFig.8 SAVI calculated by six kinds of matching models

通過表1，綜合對比六個匹配模型的評價指標，可以得出04_VVdB_VHdB_RVI匹配模型在三種植被指數(shù)中的精度都要高于其余五個匹配模型，在SAVI中的精度最高，其RMSE和MAE分別為：0.09236和0.07117.在六個匹配模型中，01_VV_VH_VVdB匹配模型的精度較低，但與其余五個匹配模型相比，精度相差較小，且六個匹配模型在三個植被指數(shù)中的相關(guān)系數(shù)都大于0.8，說明六個匹配模型均適用于植被指數(shù)的計算，且能獲取較好的結(jié)果.

表1 6種匹配模型計算的植被指數(shù)精度對比Tab.1 Comparison of precision of vegetation index calculated by 6 matching models

3 結(jié)論與討論

3.1 結(jié)論

為了有效地恢復受云與陰影污染的光學影像，獲取云與陰影覆蓋區(qū)域的植被指數(shù)，本文利用Sentinel-1影像構(gòu)建6種匹配模型，結(jié)合去云算法，獲取6景去云與陰影后的Sentinel-2影像，再分別利用6景去云與陰影的影像計算3種植被指數(shù)，并進行精度評價.通過上述分析，得出結(jié)論如下：

(1)去云算法能有效地去除云與陰影對光學影像產(chǎn)生的影響，恢復地表物質(zhì)的光譜曲線，為計算受云與陰影污染的光學影像的植被指數(shù)提供算法支持.

(2)六個匹配模型均適用于植被指數(shù)的計算，且能獲取較好的結(jié)果，其中04_VVdB_VHdB_RVI模型在三種植被指數(shù)中的精度均為最高.

3.2 討論

六個匹配模型中，01、02模型以VV和VH極化為主，結(jié)合dB值構(gòu)建模型，而03～06模型以VVdB和VHdB為主，結(jié)合雷達指數(shù)、本地入射角等特征構(gòu)建模型，用三個特征構(gòu)建模型時，其精度已趨于穩(wěn)定，模型包含四個特征時，精度無明顯提升.隨著特征的增加，算法運行的時間也會增加.該算法適用于SAR與光學影像行列數(shù)與分辨率相同的遙感數(shù)據(jù)，基于SAR和光學影像融合的去云算法計算植被指數(shù)能獲取較好的效果，但仍有改進空間，以后可以考慮對影像進行歸一化處理[19-21]，降低類內(nèi)的光譜變化.去云的過程中，本文采用全局搜索進行匹配，雖然精度較高，但是效率極低，后續(xù)考慮地表精細化分類(相鄰時間無云影像數(shù)據(jù))，以減小搜索的范圍.