劉春亭，馮權(quán)瀧，金鼎堅(jiān)，史同廣，劉建濤，朱明水

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

0 引言

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

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

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

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

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

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

圖1 研究區(qū)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 數(shù)據(jù)源及預(yù)處理

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

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

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

(1)

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

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

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

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

2 研究方法

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

圖2 技術(shù)流程Fig.2 Flowchart of research

2.1 特征提取

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

(2)

(3)

(4)

(5)

(6)

(7)

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

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

2.2 樣本選擇

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

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

圖3 樣本點(diǎn)分布Fig.3 Distribution of sampling points

2.3 隨機(jī)森林算法理論

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

2.4 精度評估

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

3 結(jié)果與分析

3.1 分類結(jié)果及精度驗(yàn)證

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

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

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

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

3.2 不同特征組合分類結(jié)果比較

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

表4 Sentinel數(shù)據(jù)組合方案Tab.4 Combination schemes of Sentinel data

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

表5 分類精度統(tǒng)計(jì)Tab.5 Statistics of classification accuracy

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

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

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

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

3.3 與其他方法比較

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

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

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

由此可見，特征集合與分類方法是影響不透水層提取精度的關(guān)鍵因素。基于多源Sentinel數(shù)據(jù)，運(yùn)用隨機(jī)森林算法提取不透水層信息可以獲得較高的精度，說明Sentinel衛(wèi)星光譜數(shù)據(jù)與雷達(dá)數(shù)據(jù)結(jié)合在提取不透水層提取方面的應(yīng)用潛力，也驗(yàn)證了隨機(jī)森林算法在分類方面具有的優(yōu)勢。

4 結(jié)論

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

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

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

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

本研究驗(yàn)證了隨機(jī)森林算法協(xié)同多源Sentinel-1和Sentinel-2數(shù)據(jù)在黃河三角洲不透水層提取中的優(yōu)勢，為光學(xué)和雷達(dá)影像的數(shù)據(jù)融合提供參考。