隨機森林協同Sentinel-1/2的東營市不透水層信息提取

劉春亭，馮權瀧，金鼎堅，史同廣，劉建濤，朱明水

(1.山東建筑大學測繪地理信息學院，濟南 250101；2.中國農業大學土地科學與技術學院，北京 100083；3.中國自然資源航空物探遙感中心，北京 100083)

0 引言

不透水層是指諸如瀝青或水泥道路、屋頂、停車場等阻止水下滲的天然或人工地表覆蓋[1]。當前城市化正處在快速發展的階段，隨著城市化進程的不斷加快，城市空間擴張和人口的快速增長，城市熱島、環境惡化、洪澇災害、交通擁堵等一系列問題愈發嚴重，給生態環境及城市發展帶來巨大挑戰。城市化進程的顯著表現之一就是不透水層的擴張。在城市化過程中，土壤、水域、植被等自然景觀被大范圍的不透水層替代，阻礙了地表的蒸散作用，破壞了地表熱平衡，城市熱島效應使得城市增溫顯著。不透水層的廣泛存在還導致水循環的失衡，造成地下水補給減少、水質下降、水污染、暴雨流量和洪水內澇災害頻發[2-3]，嚴重影響人類生產生活。研究表明，不透水層的面積和分布是城市化對環境造成影響的重要因子[4]，因此及時準確地獲取不透水層信息對于城市化及生態環境保護具有重要意義[5]。

遙感技術由于大范圍同步觀測、重復周期短、動態、迅速的優勢，成為從局部到全球尺度估算不透水層的主要技術手段。近年來不同數據源、不同尺度下提取不透水層的研究已經取得了很多成果，利用遙感數據提取不透水層的方法主要有混合光譜分析[6-8]、指數模型[9-11]、基于決策樹[12-13]、支持向量機[14-17]、神經網絡[18-19]等方法的影像分類，以及多源遙感數據結合的方法[20-22]等。但由于城市景觀的復雜性，準確及時地估算不透水層空間分布信息仍然是一項具有很強挑戰性的任務[23-25]，如建筑物及樹木的陰影容易與反射率較低的不透水層相混淆，裸土、鹽田等也易錯分為不透水層。

通過文獻綜述發現目前關于黃河三角洲地區的不透水層的信息提取多依靠光學遙感影像。隨著遙感技術的發展，多平臺、多傳感器、高空間分辨率、高光譜分辨率的遙感數據不斷涌現，融合多源數據被認為是改善不透水層提取精度的一種有效方法。光學遙感容易受到云、雨天氣的影響，而雷達數據具有全天時、全天候的優勢，能夠獲得穩定的周期性數據，可以彌補光學影像特殊天氣情況下無法及時監測的缺陷[26-28]。協同使用光學和雷達數據進行地物識別已經成為近年來的研究熱點之一，部分研究表明使用合成孔徑雷達遙感數據補充光學遙感數據，可在一定程度上提高不透水層的提取精度[29-32]。

近幾十年來，石油開采[33]、旅游開發[34]、人工建筑[35-36]等人類活動導致黃河三角洲地區生態系統發生明顯改變，生態環境逐漸惡化，引起國內外學者的廣泛關注。但是通過文獻調研發現結合主被動遙感對黃河三角洲地區進行不透水層提取，尤其是高空間分辨率不透水層提取的研究較少，因此本研究以東營市為研究區域，探索機器學習方法協同Sentinel-2光學數據和Sentinel-1雷達數據提取不透水層信息的有效性。具體來說，本文協同使用地表反射率、紋理、后向散射系數等特征，利用隨機森林算法對研究區不透水層進行提取，探究和定量分析光學遙感與主被動遙感結合在信息提取中的精度差異，為黃河三角洲區域遙感監測、開發與保護提供科學依據。

1 研究區概況與數據源

1.1 研究區概況

黃河三角洲有著我國暖溫帶地區最年輕、最廣闊、生物多樣性最豐富的濕地生態系統，在維護區域生態安全和可持續發展方面有重要意義[37]。根據《黃河三角洲高效生態經濟區發展規劃》中劃定的黃河三角洲范圍，本文以山東省東營市為研究區域(圖1)。東營市是黃河三角洲地區的中心城市，地理位置橫跨N36°55′～38°10′，E118°07′～119°15′，屬于典型的溫帶大陸性季風氣候。其東部和北部瀕臨渤海，境內河流較多，水資源豐富，黃河在東營市墾利縣流入渤海，泥沙淤積現象嚴重，后備土地資源量大，開發空間充足；土地類型多，可利用性較強。受河流、海洋的影響，鹽堿地分布的范圍較大，土壤發育不穩定，在一定程度上影響了當地的發展。東營市有豐富的石油資源、旅游資源，其獨特的環境條件使得該地區土地利用/覆蓋變化劇烈。

圖1 研究區Sentinel-2 B4(R),B3(G),B2(B)波段合成影像Fig.1 Image of Sentinel-2 B4(R),B3(G),B2(B)bands in the study area

1.2 數據源及預處理

本文所使用的影像來自歐州航天局“哥白尼計劃”地球觀測衛星系列的Sentinel系列衛星，數據通過歐州航天局數據共享網站(https://scihub.copernicus.eu/dhus/#/home)獲得。為了提高不透水層與其他地物之間的可分性并結合研究區影像質量狀況選取了2019年9月份的遙感影像開展不透水層的提取研究。

Sentinel-1號衛星星座擁用A和B 2顆衛星，可以達到每6 d對同一地點進行重復觀測。本研究中選用的2景Sentinel-1數據為雷達干涉寬刈幅(interferometric wide swath，IW)成像模式的Level-1級別地距多視影像(ground range detected，GRD)，極化模式為VV和VH，影像的獲取時間為2019年9月26日，空間分辨率10 m。

Sentinel-1數據預處理主要在SNAP Desktop軟件中完成，處理步驟主要包括軌道校正、熱噪聲去除、輻射定標、斑點濾波、地理編碼、影像拼接等。處理的具體流程為：①使用精密軌道文件進行細化，以獲取準確地衛星位置和速度信息，提高影像配準和地理編碼的精度；②根據元數據集提供的噪聲查找表消除雷達系統內在的熱噪聲；③雷達數據存在輻射偏差，通過輻射定標后使雷達影像的像素值真正表示后向散射系數σ0，輻射定標公式為：

(1)

式中：DNi與Ai分別為像元i的灰度值與定標參數；④斑點噪聲隨機分布于影像中，使用Refined Lee濾波器進行斑點濾波，在降噪的同時能夠更好地保留輻射和紋理信息[38-40]；⑤地理編碼是將投影轉換為地理坐標投影，并且糾正雷達成像過程中由于透視收縮、疊掩等引起的幾何形變，采用30 m的SRTM DEM改正地物在雷達影像中發生的畸變，輸出10 m空間分辨率的重采樣數據。

Sentinel-2號衛星星座同樣也是由A和B 2顆衛星組成，重訪周期為5 d。本文選用的2景Sentinel-2數據是經過大氣校正的地表反射率產品L2A，獲取時間為2019年9月29日。Sentinel-2 MSI數據包含13個波段，具體參數如表1所示。在開展實驗時，本文僅使用了10 m和20 m空間分辨率的波段。

表1 Sentinel-2波段信息Tab.1 Band information of Sentinel-2

為保證不同空間分辨率數據在空間上的匹配性，本文將20 m空間分辨率的數據重采樣至10 m，且所有數據都統一轉換為WGS-84 UTM 50N投影。

2 研究方法

不透水層信息提取的技術流程如圖2所示，主要包括4個步驟：①影像下載及預處理，包括Sentinel-1影像獲取后進行軌道校正、輻射定標、斑點濾波等處理得到后向散射系數，Sentinel-2影像重采樣以及影像配準、鑲嵌、裁剪等預處理工作；②特征選擇，基于灰度共生矩陣計算紋理特征，協同地表反射率和后向散射系數構建多維特征空間；③確定分類體系及樣本選擇；④由選好的訓練樣本訓練隨機森林，執行分類；⑤根據驗證樣本得到混淆矩陣進行精度驗證，并與不同特征組合、不同方法進行精度對比。

圖2 技術流程Fig.2 Flowchart of research

2.1 特征提取

已有的研究表明紋理特征有助于改善地表信息提取的精度[41-44]。本文使用基于灰度共生矩陣(gray level co-occurrence matrix，GLCM)的紋理特征，選取了均值(MEA)、方差(VAR)、同質性(HOM)、相異性(DIS)、熵(ENT)、角二階矩(ASM)6種彼此相關性弱的紋理特征，計算公式分別為：

(2)

(3)

(4)

(5)

(6)

(7)

式中：N為灰度級數目；P為N×N維歸一化灰度共生矩陣；P(i,j)為矩陣中第i行、第j列的標準化灰度值。計算處理窗口設置為7像素×7像素[45]。另外考慮到光譜波段之間的相關性，為減少計算量，只使用紅光波段計算紋理特征。

將紋理特征、后向散射系數與反射率特征相融合，共同構造不透水層提取的輸入特征，即包含10個反射率特征波段、6個紋理特征波段、2個雷達數據的后向散射系數波段(VV和VH)。

2.2 樣本選擇

樣本的選擇直接影響識別模型的穩定性和分類精度。使用ENVI軟件中ROI提取工具建立感興趣區實現樣本的選擇，結合研究區的特點，將其土地利用/覆蓋類型劃分為：亮不透水層、暗不透水層、有作物耕地、空閑耕地、大棚用地、林地、水域、灘涂、鹽田、未利用地共10類，其分類標準及數量如表2所示，訓練樣本和驗證樣本空間分布如圖3所示。

表2 分類體系及樣本數量Tab.2 Classification scheme and number of samples (個)

圖3 樣本點分布Fig.3 Distribution of sampling points

2.3 隨機森林算法理論

隨機森林算法是Breiman于2001年提出的[46]，目前已經廣泛應用于遙感應用研究[47-50]。隨機森林算法對數據的適應能力強，數據集無須規范化，隨機森林分類受噪聲和異常值的影響較小，且不會過擬合[46]，能夠直接處理高維數據，對于缺省值也可以獲得較好的結果。該算法的基本思想是首先利用Bootstrap方法從原始訓練集中有放回的隨機抽取n個樣本并構建n個決策樹構成森林，n個樣本的個數約為總樣本個數的2/3。假設訓練樣本數據中有m個特征，每次分裂時根據最優特征進行分裂，每棵樹都充分生長，直到該節點的所有訓練樣本都屬于同一類，讓每棵決策樹在不做任何修剪的前提下最大限度無限生長，最后將生成的多棵分類樹組成隨機森林，根據多棵樹分類器投票決定最終分類結果。其中每次隨機抽取的n個樣本小于樣本量總數用于替換訓練樣本，其余用來進行交叉驗證，以此來提高隨機森林方法的泛化能力。在使用隨機森林分類時，通常需要設置生長樹的數量和節點分裂時輸入的特征變量數量。經過一系列實驗發現，當N=500時，袋外誤差(out of bag，OOB)趨于穩定，因此選擇500棵樹用于隨機森林分類，特征變量數則由參與分類的特征數的平方根計算得到[51]。

2.4 精度評估

分類完成后，采用混淆矩陣對分類結果進行精度驗證，由混淆矩陣派生的評價指標主要有總體分類精度、制圖精度(producer’s accuracy，PA)、用戶精度(user’s accuracy，UA)和Kappa系數。根據目視解譯獲得驗證樣本，在ENVI軟件中利用Confusion Matrix Using Ground Truth ROIs 工具進行精度驗證。

3 結果與分析

3.1 分類結果及精度驗證

采用本文方法分類結果的混淆矩陣如表3所示，分類結果總體精度達到93.37%，Kappa系數達到0.925 8。從表3中可見，空閑耕地、灘涂分別有11個和9個樣本點錯分為亮不透水層，這些土地利用類型表面存在的鹽堿顆粒與不透水層光譜性質相似，造成了錯分，導致不透水層的高估現象。研究區內的未利用地主要是油田附近被壓實的裸土以及正在建設的區域，因此未利用地易與亮不透水層或暗不透水層混淆，其中有36個未利用地樣本點錯分為亮不透水層，13個暗不透水層錯分為未利用地。有作物耕地與林地也發生了混淆現象，其中68個林地樣本被錯分為有作物耕地，同時21個有作物耕地錯分為林地，這類錯分可能與樣本的選擇有關，同時也受到林地密度、樹木高度及季節等因素的影響。

表3 分類結果混淆矩陣Tab.3 Confusion matrices for classification

將分類結果的暗不透水層和亮不透水層合并為不透水層，其他地類合并為非不透水層，得到東營市不透水層分布情況(圖4)。結合多光譜影像及Google Earth影像發現建城區、農村居民點、道路等不透水層與現實分布狀況吻合，取得了較好的提取效果?？傮w來看，研究區內不透水層多為暗不透水層，主要分布在東營區城區、墾利區城區，河口區、廣饒縣、利津縣的建成區和分散式分布的農村居民地及道路。

圖4 不透水層提取結果Fig.4 Extracted result of impervious surface

3.2 不同特征組合分類結果比較

為探索Sentinel衛星的光譜和雷達數據對不透水層提取精度的影響，本文設計了3種組合方案進行對比(表4)。

表4 Sentinel數據組合方案Tab.4 Combination schemes of Sentinel data

通過對上述3種特征組合方案進行隨機森林分類(表5,圖5)，發現僅使用后向散射系數(方案B)分類時精度較低，總體精度為40.76%，Kappa系數為0.335 8。使用地表反射率特征及紋理特征(方案A)進行隨機森林分類時已經可以得到較為理想的分類精度(總體精度為92.04%、Kappa系數為0.911 0)，而且暗不透水層的PA和UA均大于90%，亮不透水層的PA和UA也達到87%以上。在方案A基礎上引入后向散射系數輔助特征后，總體精度和Kappa系數分別達到93.37%和0.925 8，暗不透水層的PA和UA分別提高了1.80和2.72百分點，亮不透水層的UA提高0.47百分點，說明協同Sentinel-1和Sentinel-2數據可以有效改善不透水層的估算精度。

表5 分類精度統計Tab.5 Statistics of classification accuracy

(a)方案A (b)方案B (c)方案C

由于研究區域較大且使用10 m空間分辨率，圖5中不同特征組合方案分類結果差異不太明顯，尤其是方案A和C。因此結合混淆矩陣發現，后向散射系數對各類地物的提取精度都不高，但能提供光譜數據的互補信息，使得方案C分類精度有所提高。結合后向散射系數后，除大棚用地和鹽田的分類精度在正常波動范圍內略有降低外，其他地類提取精度均有不同程度的提高。數據組合A方案中，暗不透水層、未利用地和空閑耕地的混淆程度較高，有作物耕地與林地也存在一定程度的混淆。在方案C中這種情況被明顯改善，后向散射系數提高了不透水層與空閑耕地、未利用地等裸露土地的區分度，地物錯分減少，提取精度獲得提高。

同時可以看到，后向散射系數(方案B)的分類結果除了精度較低以外，還有明顯的椒鹽現象，圖像上存在大量的離散斑塊和孤立像元?；诠庾V和紋理(方案A)的分類結果椒鹽現象有所緩解，孤立像元數據有所減少。光學、紋理和雷達結合(方案C)的分類結果中分類精度提高，地物邊界分類清晰，椒鹽現象被較好抑制。

總的來說，利用Sentinel-2地表反射率、紋理特征和Sentinel-1雷達影像后向散射系數能夠提高不同地物類型的可分性，尤其在減少不透水層與裸露土地的混淆程度上有顯著作用，并且對于基于像素分類常出現椒鹽現象有明顯抑制作用，因此這種特征組合在獲取高空間分辨率不透水層信息方面有較大的優勢。

3.3 與其他方法比較

為了進一步評估隨機森林方法在提取不透水層的優勢，本文選取了支持向量機、決策樹等遙感圖像分類領域較常用的機器學習方法進行了對比實驗。為了保證不同方法所得結果的可比性，本文使用了相同的數據(方案C)和樣本，3種方法的比較結果如表6所示。

表6 隨機森林、支持向量機、決策樹分類精度比較Tab.6 Comparison of classification accuracy of RF,SVM and CART

從表6中可以看出，隨機森林分類的精度最高，支持向量機分類精度略低(總體精度為93.19%，Kappa系數為0.923 8)，決策樹分類精度最低(總體精度為87.79%，Kappa系數為0.863 5)。隨機森林與支持向量機分類方法精度都較好，體現了Sentinel數據對于不同方法的良好適用性。但是與隨機森林分類算法相比，支持向量機分類方法需要設置多個參數來執行分類，計算量大，耗費時間長，這也正是隨機森林得到廣泛應用的原因之一。

由此可見，特征集合與分類方法是影響不透水層提取精度的關鍵因素?；诙嘣碨entinel數據，運用隨機森林算法提取不透水層信息可以獲得較高的精度，說明Sentinel衛星光譜數據與雷達數據結合在提取不透水層提取方面的應用潛力，也驗證了隨機森林算法在分類方面具有的優勢。

4 結論

本文基于多源Sentinel-1和Sentinel-2數據，利用光譜特征、紋理特征和后向散射系數，設計了3組不同的特征組合方案，運用隨機森林算法在黃河三角洲東營市開展了不透水層的遙感提取研究。通過對不同特征組合與分類算法進行對比，得到如下結論：

1)隨機森林協同反射率光譜特征、紋理特征以及反向散射系數在東營市不透水層提取中能夠取得最高的精度；單獨使用Sentinel-2反射率光譜數據所提取的不透水層精度低于上述的特征組合。

2)在光譜特征和紋理特征的基礎上引入后向散射系數能增強地物之間的可分性，特別是對提高暗不透水層與裸露土地的提取精度有積極作用，說明了Sentinel-1雷達數據與Sentinel-2光譜數據協同在提取不透水層具有較好的應用前景。

3)通過與支持向量機和決策樹分類的對比實現，隨機森林算法的提取精度和效率較好。

本研究驗證了隨機森林算法協同多源Sentinel-1和Sentinel-2數據在黃河三角洲不透水層提取中的優勢，為光學和雷達影像的數據融合提供參考。