基于Sentinel-2A和Landsat8的城市不透水面的提取

趙怡，許劍輝，鐘凱文，王云鵬，胡泓達，吳萍昊

(1.廣東省科學院廣州地理研究所/廣東省遙感與地理信息應用重點實驗室/廣東省地理時空大數據工程實驗室，廣州 510070；2.南方海洋科學與工程廣東省實驗室(廣州)，廣州 511458；3.中國科學院廣州地球化學研究所，廣州 510640；4.中國科學院大學，北京 100049)

0 引言

城市不透水面是快速城市化的典型產物[1-3]。在城市土地利用中，植被和土壤等自然地表逐漸減少，取而代之的是不透水面。不透水面絕大部分是由人工材料制成，水體無法滲透，例如道路(水泥、瀝青等材料)和建筑(水泥、巖石等材料)[4]。許多研究表明，不透水面是城市重要的土地覆蓋類型，也是城市生態環境相關問題的重要參考指數，例如城市內澇和暴雨天氣等[5]。因此，不透水面的遙感提取已經成為城市遙感的研究熱點[6-13]。

在城市不透水面遙感提取中，線性光譜混合分解(linear spectral mixture analysis,LSMA)是常用的方法[14-15]，該方法基于亞像元尺度和植被-不透水面-土壤(vegetation-impervious-surface,V-I-S)模型[16]，假設遙感影像中混合像元的反射率是該像元中所有地物端元反射率的線性組合[17]。首先從混合像元內各類端元的純凈像元獲取端元的光譜信息，然后通過光譜解混計算，可以得到混合像元中高反射率和低反射率不透水面的面積比例，得到不透水面的面積比例，即不透水面蓋度。由于環境影響，以及在特定波段上不透水面的光譜特征與其他地物相似，導致不透水面的提取精度較低[18]。在LSMA解混過程中，端元選擇是至關重要的一步[19]，原因在于LSMA是基于某一端元光譜特性，判斷混像元內該類端元的面積比例，因此，所選取端元的光譜信息是進行解混混合像元的基礎。

端元的純凈像元光譜信息來源有兩個渠道：一是從標準地物光譜庫中直接選取目標地物的光譜曲線，二是從實際目標影像中選擇純凈像元，得到不同地物的光譜曲線。由于地物的光譜特征曲線會受到環境因素的影響，且不同地物之間也會產生相互作用，因此，標準地物光譜庫中的光譜曲線并不能較好地表示目標地物，故而大部分研究選擇第二種方式獲取端元的光譜曲線[20]。傳統的LSMA在端元中采用最小噪聲分離變換(minimum noise fraction rotation，MNF)[21]將遙感影像的信息集中到前幾個波段，達到去除噪聲的目的；然后通過計算純凈像元指數(purity pixel index，PPI)[22]篩選出純凈度較高的像元；最后通過N維可視化選擇純凈像元。此外，N-FINDAR[23]算法也常用于純凈端元的選擇。

為了描述城市進程發展，長時間序列的遙感影像必不可少。在目前的遙感發展中，能滿足長時間序列要求的影像大都是中等空間分辨率的。然而在中等空間分辨率的遙感影像(例如30 m空間分辨率的Landsat8 OLI影像)中通常存在大量的混合像元，純凈像元選擇存在不確定性，其精度直接影響不透水面的提取精度。而在較高空間分辨率的遙感影像中(例如10 m空間分辨率Sentinel-2A MSI影像)，可以選擇更多的純凈像元。在同一時期，目標地物在特定波段上的反射率是一定的，這就意味著，高空間分辨率中端元的光譜特性與中等空間空間分辨率中相似。為了獲取更為可靠的純凈像元，本研究提出結合Sentinel-2A MSI影像改進從Landsat8 OLI 選擇端元的光譜曲線。

此外，由于地物的光譜曲線受外界影響，且端元選擇存在誤差，因此，經過LSMA解混的不透水面蓋度仍然存在錯分現象，即在土壤分量中可能含有少量的不透水面，而在不透水面分量中也存在土壤或植被蓋度。為了解決以上問題，本研究提出解混結果優化方案，即在像元尺度上，基于閾值分割法，利用歸一化植被指數(normalized differential vegetation index，NDVI)和干旱裸土指數(dry bare-soil index，DBSI)識別不透水面、植被和土壤，從而分離出錯分的像元，提高不透水面蓋度提取精度。

1 研究區概況與數據源

1.1 研究區概況

研究區影像如圖1所示。廣州市位于E112°57′～114°30′，N22°26′～23°56′之間，是廣東省的省會，屬于國際中心城市，是國際商貿中心和綜合交通樞紐，是粵港澳大灣區、泛珠江三角洲經濟區的中心城市以及一帶一路的樞紐城市。本研究選擇廣州市6個主城區，包括越秀區，海珠區，天河區，白云區，黃埔區，荔灣區(圖1)。研究區域內土地利用類型豐富，包括大面積的植被區域(森林、農田)，裸土區域，不同密度的居民用地和商業用地。

圖1 研究區Landsat8 B4(R)，B3(G)，B2(B)合成影像Fig.1 Study area image combined with Landsat8 B4(R)，B3(G)，B2(B)

1.2 數據源

Sentinel-2A 衛星發射于2015年6月23日，是“全球環境與安全監測”計劃的第二顆衛星。該衛星的重訪周期為十天，幅寬為290 km,攜帶的多光譜成像儀(multispectral instrument，MSI)有13個波段，從可見光和近紅外到短波紅外，具有不同的空間分辨率[24-25]。本研究選取2017年11月1日的Sentinel-2A MSI遙感影像的B2，B3，B4，B8，B11，B12作為研究目標，采用Sen2cor處理器進行大氣校正。此外，通過哨兵系列衛星產品綜合應用平臺(Sentinel Application Platform，SNAP)，將B11，B12的遙感影像進行重采樣，使其空間分辨率為10 m。

Landsat 8衛星發射于2013年2月11日，是美國陸地衛星計劃的第8顆衛星。該衛星的重訪周期為16 d，攜帶的陸地成像儀(operational land imager,OLI)有9個波段，其中包括7個常用波段(B1—B7，空間分辨率為30 m)，一個全色波段(B8，空間分辨率為15 m)和一個云檢測波段(B9，空間分辨率為30 m)。本研究選取2017年10月23日的Landsat8 OLI遙感影像的B1—B7作為研究目標，通過輻射定標和大氣校正，將影像的DN值轉換為地表反射率。地表發射率產品是在美國地質調查局(United States Geological Survey，USGS)地球資源觀測(Earth Resources Observation and Science，EROS)和科學中心科學處理機構(Center Science Processing Architecture，ESPA)(https://espa.cr.usgs.gov/)[26]的定制頁面進行定制得到。

本研究采用的Sentinel-2A MSI 和Landsat8 OLI數據參數表1所示。

表1 Landsat8 OLI 和Sentinel-2A MSI 數據介紹Tab.1 Landsat8 OLI and Sentinel-2A MSI data

此外，采用2017年WorldView-2遙感影像(空間分辨率為1.8 m)，通過對樣本區域內的不透水面進行矢量化，計算不透水面面積，作為精度驗證數據。

2 研究方法

本研究所提出的方法分為3個步驟：首先，從Sentinel-2A和Landsat8 OLI 遙感影像中選取純凈像元，獲取線性光譜解混所需的端元光譜曲線；然后利用傳統線性光譜解混方法對兩種數據進行光譜解混；最后，利用NDVI和DBSI兩種光譜指數，對解混結果進行后處理，得到精度較高的城市不透水面比例。本文整體流程如圖2所示。

圖2 流程圖Fig.2 Flow chart

2.1 端元選擇優化方案

端元選擇是LSMA過程中的重要步驟。理論上，作為光譜解混的基礎，所選取的端元應是遙感影像上的純凈像元(只包括植被，裸土，高反射率不透水面與低反射率不透水面4類地物中的一類)，然后利用所選純凈像元的光譜信息，對混合像元進行光譜解混，計算得到混合像元內每一種地物的面積比例。因此，實驗中所選端元的純凈度直接影響解混的精度。Sentinel-2A MSI影像的分辨率較高，像元大小為10 m×10 m，與Landsat8 OLI影像(像元尺寸為30 m×30 m)相比，其像元的純凈度較高，純凈像元的選取也比較容易。此外，在城市環境中，土壤分布比較破碎，很少有大面積的裸土呈現，因此，像元面積越小，裸土純凈像元就越多，也就是說在Sentinel-2A影像上裸土的純凈像元更多。

為了提高端元光譜的選去精度，本研究提出端元選取優化方案：結合Sentinel-2A MSI影像優化從Landsat8選取的端元光譜信息。具體過程如下：

1)采用傳統方法對Landsat8 OLI影像進行端元選取，將所選純凈像元的光譜信息以文件的形式輸出。

2)利用選擇感興趣區域(region of interest，ROI)的方式對Sentinel-2A影像進行端元選取，選擇植被、裸土、高反射率不透水面與低反射率不透水面4類地物的純凈像元，并統計其光譜信息，生成端元光譜曲線，并以文件的形式輸出。

3)利用從Sentinel-2A MSI 影像中選取的的純凈像元在每個波段(B2,B3,B4,B8,B11,B12)上的反射率替換Landsat8 OLI影像中選取的純凈像元在相應波段(B2—B7)的反射率，生成一個新的端元光譜信息文件，作為線性光譜混合分解的基礎(從表1可以看出，Sentinel-2A MSI影像的6個波段與Landsat8 OLI的前6個波段的波譜范圍比較接近，可以視為相同的波段)。

2.2 線性光譜混合分解(LSMA)

諸多研究表明，LSMA具有較強的理論基礎。該方法假設遙感影像中混合像元在某一波段上的反射率是組成混合像元的各類端元在該波段上的反射率的線性組合，像元反射率越接近某類端元的反射率，就說明這種端元在該像元中的面積比例越大[27]。LSMA公式為[28]：

(1)

(2)

式中：i=1,2…,M,M為波段數；n為端元數目；Ri為混合像元在i波段的反射率；fk為端元k的面積比例；ERi為i波段的計算殘差；RMS為計算殘差，以此評價解混結果的優劣[29]。

由于混合像元中多類端元的光譜之間會發生相互作用，此外，環境等因素也會影響端元的光譜信息。因此在解混過程中，會存在不透水面被錯誤估計的問題，例如在不透水面覆蓋度低于20%的區域，采用LSMA解混混合像元時，不透水面比例會被嚴重高估[3]。

由于不透水面的低反射率部分與水體部分反射率都比較低，十分接近，為了減少解混誤差，在數據預處理階段，利用改進的歸一化水體指數(modified normalized difference water index,MNDWI)[30]將水體去除，計算公式為：

(3)

2.3 解混結果優化方案

在LSMA結果中，解混誤差普遍存在，具體表現為：①在土壤蓋度中含有不透水面，即不透水面被錯分為土壤；②在不透水面蓋度中(高反射率與低反射率蓋度之和)含有植被和土壤，即少量植被和土壤被錯分為不透水面。為了降低以上誤差，提高城市不透水面的提取精度，本研究結合像元尺度，利用DBSI和NDVI，通過二者的閾值，提取解混過程中被錯分的像元，將其歸還到相應的蓋度中去。其中，DBSI和NDVI的計算公式分別為：

(4)

(5)

式中：BNIR為影像的近紅外波段(Landsat8 B5)；BR為影像的紅色波段(Landsat8 B4)。

優化過程分為3步，具體如下：

1)從土壤蓋度分量中提取錯分的不透水面蓋度。利用DBSI閾值篩選土壤蓋度中的不透水面像元，借助DBSI影像的直方圖，通過實驗確定0.1作為該步DBSI的閾值，即在土壤蓋度中，DBSI值小于0.1的像元被認為是不透水面，將這部分像元的解混值認為是不透水面蓋度，對應圖2中的不透水面蓋度(1)；

2)將LSMA解混結果中的高反照率和低反照率不透水面蓋度相加，得到LSMA解混遙感影像得到的不透水面蓋度，然后將上一操作中的不透水面蓋度(1)與之相加，得到不透水面蓋度(2)；

3)去除上一操作中的不透水面蓋度(2)中含有的少量植被和土壤蓋度，得到最終的不透水面蓋度。這部分誤差來源于LSMA本身的解混誤差。理論上，NVDI值越接近1，說明目標像元越接近植被，負值代表水體，0代表土壤或者沙漠化土地[32]。為了保證識別的像元為純凈的植被和土壤，借助DBSI和NDVI影像的直方圖，確定DBSI和NDVI閾值分別為0.2和0.3，即DBSI<0.2且NDVI>0.3的像元被認為是植被，DBSI>0.2且NDVI<0.3的像元被認為是土壤。

在以往研究中，DBSI>0.26的像元被認為是裸土，NDVI值接近于1被認為是植被。在步驟1)中，從土壤蓋度中提取出不透水面，為了盡量保留LSMA結果的土壤蓋度中土壤的真實面積比例，DBSI值應稍小于理論值，因此DBSI值選擇為0.1。而在步驟3)中，從不透水面中去除植被和土壤，為了盡量保留真實不透水面的面積比例，DBSI和NDVI閾值應稍大于理論值，因此NDVI閾值被確定為0.3；然而在實際實驗過程中，經過反復試驗得知，DBSI閾值為0.2能較好地區分不透水面與其他地物。

2.4 精度驗證

為了驗證本文所提出的方法對提取城市不透水面的有效性，隨機選取了171個樣本區域，樣本大小為480 m×480 m，并且采用WorldView-2 影像對樣本區域進行矢量化，得到不透水面的樣本真實面積比例。通過系統誤差(system error，SE)、平均絕對誤差(mean absolute error，MAE)、均方根誤差(root mean square error，RMSE)和決定系數R2來衡量實驗結果的精度。前3個指標的計算公式分別為：

(6)

(7)

(8)

“喝西鳳、吃泡饃、吼秦腔”是被概括成了“關中人的形象”，當然，這也是游人必不可少的原生態體驗。性烈醇勁的西鳳酒，勁道味重大碗泡饃，街市上都隨處可見。烈酒好面，在北方并不罕見，唯有這秦腔卻是西安人和關中漢所獨有特征?！鞍税倮锴卮▔m土飛揚，三千萬楞娃齊吼秦腔?！甭犌厍?、吼秦腔至今仍是西安人日常生活中的一大樂事，和兵馬俑一樣，秦腔也可被視為陜西的象征。

3 結果與分析

3.1 端元選取

圖3展示了不同方法選擇的端元光譜曲線。在端元選取過程中，由于影像自身的原因，以及環境因素，采用傳統端元選擇方法選取的土壤和高反射率不透水面這兩類端元容易混淆，且精度不高，如圖3(a)所示。結合Sentinel-2A影像選擇較純凈的像元，優化傳統方法選擇的端元光譜信息，如圖3(b)所示，降低了各類地物端元光譜的混淆問題，且與實際地物光譜信息更接近，精度更高。

(a)傳統端元選擇方法 (b)端元選取優化方案

3.2 結果比較分析

通過采集到的端元光譜曲線，利用LSMA對Landsat8 OLI遙感影像進行像元解混，得到如圖4所示的不透水面蓋度分布圖。通過對比實驗分析，結果顯示在城市密度較高的區域，圖4(b)中不透水面蓋度比圖4(a)明顯增加。圖4(c)圖采用了解混結果優化方案，結果顯示在森林等透水區域，不透水面蓋度大幅降低，接近于0，與實際不透水面分布十分接近。

(a)傳統LSMA (b)結合端元優化方案的LSMA (c)結合端元優化方案和解混結果優化方案的LSMA

通過精度評價，如圖5和表2所示，結果表明，通過端元選擇優化方案，將高空間分辨率較高的Sentinel-2A MSI影像上選取的純凈像元的光譜信息代替相應波段上從Landsat8 OLI影像中選取的端元的光譜信息，經過LSMA解混后，不透水面蓋度的精度得到提高，從線性擬合結果來看，圖5(b)中的觀測值(利用不同LSMA提取的不透水面蓋度)與真值比圖5(a)更接近，擬合線周圍的樣本點更為密集。如表2，結合優化方案的LSMA所提取的不透水面蓋度的R2比傳統LSMA高16%。

(a)傳統LSMA (b)結合端元優化方案的LSMA (c)結合端元優化方案和解混 結果優化方案的LSMA

表2 不同方法提取的城市不透水面蓋度精度分析Tab.2 accuracy assessment of urban IS fraction extracted by different methods

綜上所述，本研究所提出的端元優化方案可以使LSMA更加準確地識別不透水面，降低了不透水面蓋度被低估的問題(在不透水面覆蓋度較高的區域)，提高城市化程度較高的區域中不透水面蓋度的提取精度。

但是圖5(a)和(b)同時存在一個問題，即在不透水面實際蓋度較低的區域(橫軸接近于0的位置)，觀測值明顯比真值大。也就是說，在不透水面覆蓋度較低的區域，傳統LSMA和結合端元光譜優化方案的LSMA會高估不透水面蓋度。圖5(c)顯示了結合端元優化方案和解混結果優化方案的LSMA提取的水面蓋度與真值的線性擬合結果，該圖表明，在不透水面蓋度較低的區域，觀測值明顯降低。如表2所示，結合端元優化方案和解混結果優化方案的LSMA所提取的不透水面蓋度的R2比傳統LSMA高20%，比僅結合端元優化方案的LSMA高3%，且SE為-0.066，說明整體系統誤差較小，該方法所得到的不透水面蓋度與實際情況接近。以上結果說明通過解混結果優化方案，大大改善了不透水面被高估的問題(在不透水面蓋度較低的區域)。

圖6展示了通過不同方式獲取的不透水面蓋度細節圖，通過對比圖6(b)和(c)，圖6(f)和(g)來看，在城市化較快的區域，結合端元優化方案可以提高不透水面蓋度，從而提高提取精度。從圖6(d)和(h)來看，結合解混結果優化方案，可以明顯去除森林和土壤等透水區域內不透水面蓋度，同時保持原有的不透水面蓋度，有效解決不透水面被高估的問題。

(a)Sentinel-2AMSI影像 (b)傳統LSMA提取的不透水面蓋度 (c)結合端元選擇優化方案的LSMA提取的不透水面蓋度 (d)結合端元選擇優化方案和解混結果優化方案的LSMA提取的不透水面蓋度

4 結論

1)本研究針對Landsat8 OLI遙感影像的端元選擇問題提出了端元選擇優化方案；針對LSMA后處理問題，結合解混結果優化方案，以此提高不透水面的提取精度。

2)與Landsat系列影像相比，Sentinel-2A MSI影像的空間分辨率較高，可以快速選擇相對比較純凈的像元作為LSMA的端元，然后將選取像元的光譜信息代替對應在波段上從Landsat8 OLI影像上選擇的端元的光譜信息，優化了光譜解混所必需的端元光譜信息，從而達到了提高LSMA對混合像元解混精度的目的。

3)對于LSMA后的各地物蓋度分量，結合解混結果優化方案，利用DBSI和NDVI兩個指數作為識別植被和土壤的重要參數，可以比較準確地從土壤蓋度分量中提取被錯分的不透水面蓋度分量，且剔除了不透水面蓋度分量中被錯分的植被和土壤蓋度分量，整體提高了不透水面蓋度的提取精度。

以上兩種優化方案解決了傳統LSMA在不透水面蓋度較高的區域(例如繁華的商業區，聚集的居民地等區域)被低估的問題，以及在不透水面蓋度較低的區域(例如大面積的森林區域等)被高估的問題。然而，本研究在DBSI和NDVI閾值選擇過程中，大部分操作依賴統計實驗，未來可以從兩個指數的閾值入手，為不透水面蓋度提取提供可靠的理論方法，從而達到提高不透水面蓋度提取精度的目的。