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

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

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

0 引言

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

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

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

1.1 試驗(yàn)區(qū)概況

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

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

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

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

圖1 Sentinel影像Fig.1 Sentinel images

配準(zhǔn)完成后，裁剪出研究區(qū)域的Sentinel影像.如圖1所示，(a)為無(wú)云的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)換為二值影像，再手動(dòng)添加未被探測(cè)的被云污染的區(qū)域，如圖2.

1.3 研究方法

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

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

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

(1)

(2)

(3)

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

(4)

(5)

(6)

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

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

(7)

(8)

(9)

2 結(jié)果與分析

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

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

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

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

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

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

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

3 結(jié)論與討論

3.1 結(jié)論

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

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

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

3.2 討論

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