熊元康,張鴻輝*,梁宇哲,羅偉玲,洪 良
(1.廣東國地規劃科技股份有限公司,廣東 廣州 510650;2.廣東省土地調查規劃院,廣東 廣州 510075)
土地利用/覆被變化(LUCC)記錄了人類在地球表面的空間格局活動,是導致生物多樣性減少、氣候變化、生態環境演變、生物化學循環乃至全球變化的主要因素[1-2]。快速、準確、全方位地獲取一個地區的土地覆蓋變化信息,可為該地區的社會經濟發展、生態環境建設、國土空間規劃等提供重要的基礎數據。近年來,遙感技術因其具有監測范圍廣、實時性強以及使用成本低等優點已被廣泛應用于土地覆蓋變化監測中[3]。
目前,根據變化監測對象的不同,基于遙感影像的土地覆蓋變化監測方法可分為3大類:①基于像素級的土地覆蓋變化監測方法,直接利用像元的光譜特征值,通過多時相的遙感數據進行影像間的差值、比值等代數運算,進而獲取差異影像,再利用經驗或通過樣本集獲得的閾值進行土地覆蓋變化監測,如于冰[4]等通過對GlobeLand30數據進行重編碼等操作,提出了一種基于像元轉換的土地覆蓋變化監測方法;②基 于特征級的土地覆蓋變化監測方法,通過不同地類所表現出來的光譜特征、紋理特征等進行土地覆蓋變化監測,與基于像素級的土地覆蓋變化監測方法相比,該方法不易受遙感影像時相變化的影響,差異信息更突出、監測精度更高,是目前土地覆蓋變化監測的主流方法之一,如Jose D A[5]等結合Landsat數據以及L波段SAR數據計算得到各類指數(NDVI、EVI等),再利用隨機森林分類器成功提取了緬甸南部Tanintharyi地區的土地覆蓋動態變化信息;③基于對象級的土地覆蓋變化監測方法,主要針對高分辨率的遙感影像,利用多尺度分割生成不同的基元對象,再結合基元對象的特征進行土地覆蓋變化監測,與前 兩類方法相比,該方法能獲取更加豐富的特征信息,便于提高分類精度,如吳田軍[6]等基于SPOT 4的多時相數據,利用對象級的土地覆蓋變化監測方法成功提取了廣東省東莞市東北區域2005-2008年的土地覆蓋變化信息。
目前,利用中高分辨率遙感影像在我國南方地區進行大范圍土地覆蓋變化動態監測的研究較少(四季多雨,無云數據更加難以獲取),且多集中于利用統計數據與土地利用數據相結合的方法進行土地覆蓋變化動態監測[7]。為了解決大范圍土地覆蓋變化監測所面臨的低分辨率遙感影像混合像元嚴重、中高分辨率遙感影像監測范圍小和重返周期長等問題,本文以廣東省為研究區,基于Google Earth Engine (GEE)云平臺,結合中等分辨率的Sentinel-1/2和 Landsat 7/8等多源遙感數據,通過構建多元時間序列影像的方式進行土地覆蓋變化動態監測。GEE云平臺是一個能在大尺度范圍下進行遙感數據處理和分析的云平臺,提供了一個完備的集成環境。目前,國內外多個研究組已基于該平臺開展了各種對地觀測研究,如水稻遙感制圖[8]、耕地遙感提取[9]、農作物種植結構 提取[10]等。
廣東省地處我國大陸最南部,東鄰福建,北接江西、湖南,南臨南海,西連廣西;屬于熱帶或亞熱帶季風氣候區。廣東省以廣州市為中心,東起潮州市、西至雷州市,橫跨21個地級市(圖1);陸地面積為17.97萬km2;地勢表現為北高南低;地貌類型復雜,北部、東北部和西部均有較高山脈,中部和南部沿海地區多為低丘、臺地或平原[11]。隨著廣東省各地區城市化和工業化的快速推進,耕地和林地不斷被蠶食,土地可墾率降低,土地后備資源嚴重不足,導致人地矛盾日益突出[12]。因此,快速、準確、全方位地獲取廣東省的土地覆蓋變化信息,可為其社會經濟發展、生態環境建設、國土空間規劃等提供重要的基礎數據。

圖1 研究區位置與驗證樣本點的空間分布 (審圖號:GS(2016)2929)
1.2.1 遙感數據
由于廣東省四季多雨,無云數據難以獲取,且中高分辨率遙感影像監測范圍小、重返周期長,采用單一數據源的中高分辨率影像不能滿足大范圍土地覆蓋變化監測的需求,因此本文結合中高分辨率的 Sentinel-1/2和Landsat7/8等多源遙感數據進行研究區內的土地覆蓋變化監測。遙感影像數量如表1所示。

表1 本文采用的遙感影像數量
1.2.2 地面參考數據
本文結合野外實地考察、航拍影像(0.2 m)、Google高清影像、Planet高清影像(5 m)以及土地利用數據,將研究區內的土地覆蓋類型分為6大類,如表2所示,并根據上述數據集獲得訓練數據集。

表2 土地覆蓋類型分類表
本文利用Sentinel-2、Landsat 8、Landsat 7自身的質量評價波段(QA)進行相應的去云處理。同時,由于不同傳感器在設計上的差異,導致不同傳感器具有不同的光譜響應函數(SRF),如Sentinel-2和Landsat 7數據在紅光波段的均方根誤差(RMSE)超過了8%[13],因此為了減小不同光譜響應函數帶來的誤差,本文以Landsat 8數據為基礎,對Sentinel-2和Landsat 7數據相應波段進行相關的線性轉換[14-15],并將所有數據的分辨率重采樣至30 m。
土地覆蓋變化監測與分析主要的依據是遙感影像時間序列的變化特征信息,因此在構建時間序列時,時序特征的選擇顯得尤為重要。通過對研究區內土地覆蓋變化類型進行分析發現,研究區內主要的土地覆蓋變化過程為:耕地或林地轉化為建設用地、耕地轉化為魚塘或魚塘轉化為耕地、耕地轉化為園林地等。鑒于此,本文選擇的時序特征為:
1)歸一化植被指數(NDVI)。NDVI常被作為特征參數來評估地表植被的生長狀況[16],能反映紅光波段(植物吸收強烈)與近紅外波段(植被反射強烈)之間的關系,因此能良好地區分植被與建筑用地等。其計算公式為[17]:

式中,ρNIR為近紅外波段的反射率;ρRed為紅光波段的反射率。
2)歸一化建筑指數(NDBI)。NDBI常被用于城鎮建設用地提取,能反映城鎮建設用地在近紅外和中紅外波段具有高反射率值,而其他地物具有低反射率值的趨勢。其計算公式為[18]:

式中,ρMIR為中紅外波段的反射率;ρNIR為近紅外波段的反射率。
3)改進的歸一化差異水體指數(MNDWI)。MNDWI是徐涵秋[19]為了改進NDWI[20]提取水體信息的局限性(影像中有建筑陰影的水體提取效果不佳)而提出的,能明顯增強水體與建筑物的反差,減少背景誤差,因此能良好地區分水體與建筑等。其計算公式為:

式中,ρGreen為綠光波段的反射率;ρMIR為中紅外波段的反射率。
4)SAR影像時序特征。針對研究區內長時間序列無云數據難以獲取等問題,結合合成孔徑雷達(SAR)不受云雨霧等自然條件影響以及全天候等特性,利用不同地類的不同散射機理[21]來區分水體、植被以及建筑等土地覆蓋類型。本文將Sentinel-1的VV極化數據作為研究區土地覆蓋變化監測的雷達數據。為了將VV極化數據與NDVI、NDBI以及MNDWI數據更好地結合起來,本文利用式(4)對其進行歸一化處理。

式中,maxρVV為研究區內SAR長時序數據上的最大VV極化散射系數;minρVV為研究區內SAR長時序數據上的最小VV極化散射系數;ρVV為研究區內長時序上某一時期的VV極化散射系數。
根據選取的時序特征,通過Google的高清影像,本文選取研究區內典型的土地覆蓋變化類型進行時序光譜響應分析,結果如圖2所示,可以看出,在耕地轉化為建設用地的過程中,NDVI和VV后向散射系數在不斷下降,而MNDWI和NDBI在地類表現為耕地時,隨時間的變化出現了較大的浮動,在地類轉化為建設用地時,隨時間的變化表現較為平穩;在林地轉化為建設用地的過程中,NDVI和MNDWI隨時間的變化呈下降趨勢,而NDBI和VV后向散射系數隨時間變化呈逐步上升趨勢;在林地轉化為裸地的過程中,NDVI、MNDWI和VV后向散射系數隨時間變化呈逐步下降趨勢,而NDBI呈逐步上升趨勢;在水體轉化為裸地的過程中,NDVI、NDBI和VV后向散射系數隨時間變化呈上升趨勢,而MNDWI呈下降 趨勢。

圖 2 各種土地覆蓋類型的時間序列
本文分別利用隨機森林分類(RF)、支持向量機(SVM)以及回歸分類決策樹(CART)進行土地覆蓋類型提取,其中RF和SVM已被廣泛應用于土地覆蓋類型分類中,并取得了良好的分類結果,而CART易于表達且能很好地解釋某些特定的規則[22]。
RF是一種集成學習方法,采用Bootstrap抽樣技術從原始數據集中抽取訓練集,并通過訓練集構建CART決策樹。影響RF分類器性能和效率的主要參數為決策樹數量、候選特征子集以及葉節點最小樣本數。本文將這3個參數設置為:RF分類器的決策樹數量=100;候選特征子集=4;葉節點最小樣本數=1。SVM是一種非參數機器學習算法,核心是找到一個最優的超平面作為高維空間中的決策函數,進而將輸入向量分成為不同的類別。影響SVM性能和效率的主要參數包括核函數的選擇和cost參數。本文選擇線性函數作為SVM的核函數,并將cost參數設置為100。CART是一個在二進制遞歸分區過程中成長的決策樹,通過子集內變量的最大方差和最小方差將訓練數據集分成不同的類別。樹的最大深度參數決定了CART模型的復雜性,大的深度可能具有更高的精度,但也會增加過度擬合的風險。綜合以往的研究和實際需求,本文將樹的深度參數設置 為10。
對遙感影像分類結果進行精度評價是一個十分重要的工作[10]。為了保證精度驗證的準確性,本文選擇總體精度(OA)和F-Score作為土地覆蓋變化監測結果的精度評價指標。這兩個指標均來源于混淆矩陣,其中OA用以評價整體算法的有效性,而F-Score用以評價每一類的分類精度。
為了對土地覆蓋類型變化監測結果進行精度評價,本文通過分層隨機采樣在研究區內生成了1 000個驗證樣本點,如圖1所示,并通過人工目視解譯判讀其是否為土地覆蓋類型發生變化的區域,進而獲得587個非變化樣本點和413個變化樣本點。由精度評價結果可知(表3),RF方法變化監測的OA為83.5%,F-Score為0.82;SVM方法變化監測的OA為73.5%,F-Score為0.58;CART方法變化監測的OA為70.9%,F-Score為0.51。因此,RF表現最優,SVM次之,而CART結果最不理想。

表3 分類精度統計
根據多元分類特征,在GEE云平臺上分別融合生成2017年、2018年研究區內的多元分類特征影像(49個特征波段、12個NDVI時間序列波段、12個MNDWI時間序列波段、12個NDBI時間序列波段、12個VV后向散射系數時間序列波段,1個Slope波段),并同時利用收集到的訓練數據集進行RF、SVM以及CART訓練,以建立研究區內不同土地覆蓋類型的監測模型。在此基礎上,利用訓練好的RF、SVM以及CART分類模型進行2017年、2018年研究區內的土地覆蓋類型監測,結果如圖3所示。

圖3 廣東省2017-2018年土地覆蓋類型監測結果
根據2017年、2018年廣東省土地覆蓋類型分類結果,本文采用基于空間疊置分析和統計原理分析的分類后變化檢測方法進行以a為單位的土地覆蓋變化監測,結果如圖4、5所示。在此基礎上,本文對研究區內的土地覆蓋變化區域的空間分布狀況進行統計發現,研究區內土地覆蓋類型發生變化的區域主要集中在珠三角地區(廣州市等9個地區),其次為粵東地區(潮州市等4個地區),最后為粵西地區(湛江市等3個地區)和粵北山區(清遠市等5個地區)。由于珠三角地區人口的集聚和工業的發展,其土地覆蓋變化類型主要為耕地或林地轉化為建設用地或工業用地;由于粵東部分地區人口眾多,其土地覆蓋變化類型主要為耕地或林地轉化為建設用地;由于粵西地區和粵北山區是重要的農業生產空間,其土地覆蓋變化類型主要為耕地轉化為魚塘或魚塘轉化為耕地,以及耕地轉化為園林地。

圖4 廣東省2017-2018年土地覆蓋變化監測結果

圖5 珠江口2017-2018年土地覆蓋變化監測結果
為了發掘影響土地覆蓋變化的潛在因素,結合2004-2018年廣東省的各類統計數據,本文對廣東省的土地覆蓋變化進行分析發現,影響廣東省土地覆蓋變化空間分布的主要原因為人口增長、區域社會經濟發展以及政策變化,如圖6所示。在人口增長方面(圖6a),珠三角地區人口年均增長99.2萬人,粵東地區人口年均增長4.6萬人,粵西地區人口年均增長 7.8萬人,粵北山區人口年均增長8.8萬人;在區域社會經濟發展方面(圖6b),珠三角地區的年均生產總值增長為5 031億元,粵東地區的年均生產總值增長為386億元,粵西地區的年均生產總值增長為 438億元,粵北山區的年均生產總值增長為281億元; 在政策變化方面,主要體現為政府通過制定發布相關的政策來干預、調整土地的使用(如加強基礎設施建設等),由研究區內各地區的固定資產投資(基本建設投資、更新改造投資、房地產開發投資)情況可知(圖6c),珠三角地區的年均固定資產投資增長為 1 833億元,粵東地區的年均固定資產投資增長為 397億元,粵西地區的年均固定資產投資增長為290億元, 粵北山區的年均固定資產投資增長為242億元,進而導致研究區內建設用地與耕地面積比值的不斷增加,即建設用地面積不斷增加,而耕地面積不斷減 少(圖6d)。

圖6 影響廣東省土地覆蓋變化的潛在因素
本文以廣東省為研究區,通過分析研究區內的主要土地覆蓋變化類型,選取NDVI、MNDWI、NDBI、VV后向散射系數的時序影像以及DSM影像作為研究區土地覆蓋類型監測的特征波段影像,同時結合訓練數據集、RF、SVM以及CART對研究區內 2017年、2018年的土地覆蓋類型進行了監測,并在此基礎上進行了研究區內的土地覆蓋變化監測。通過對土地覆蓋變化監測結果進行分析,本文得出的主要結 論為:
1)構建的49個分類特征波段以及RF分類器能較好地適應于研究區內的土地覆蓋類型監測。較SVM和CART而言,RF具有最高的土地覆蓋變化監測精度,其土地覆蓋變化監測OA為83.5%,F-Score為 0.82。
2)對研究區內土地覆蓋變化監測結果進行分析發現,研究區內土地覆蓋發生變化的區域主要集中在珠三角地區,其土地覆蓋變化類型主要為耕地或林地轉化為建設用地或工業用地;其次為粵東地區,其土地覆蓋變化類型主要為耕地或林地轉化為建設用地;最后為粵西地區和粵北山區,其土地覆蓋變化類型主要為耕地轉化為魚塘或魚塘轉化為耕地,以及耕地轉化為園林地。
3)結合統計年鑒等相關數據,對研究區內的土地覆蓋變化進行分析發現,影響研究區內土地覆蓋變化空間分布的主要原因為人口增長、區域社會經濟發展和政策變化。