基于多源遙感數(shù)據(jù)的石漠化信息提取研究

陳 希，蘭安軍，熊康寧，廖建軍

(貴州師范大學 喀斯特研究院/國家喀斯特石漠化防治工程技術(shù)研究中心，貴州 貴陽 550001)

0 引言

近年來，石漠化問題已經(jīng)成為制約喀斯特地區(qū)發(fā)展的最關鍵問題．石漠化是指人類所實施的違反自然規(guī)律的社會生產(chǎn)經(jīng)濟活動所導致的喀斯特地區(qū)人地矛盾尖銳、植被破壞、水土流失、巖石逐漸裸露、土地生產(chǎn)力衰退喪失[1]，地表在視覺上呈現(xiàn)類似于荒漠化景觀的演變過程[2]．其中，貴州省的喀斯特石漠化程度在全國范圍內(nèi)最為嚴重．目前，貴州省的石漠化土地面積達到4.2×104km2，占貴州省總面積的比例為22.9%，占貴州省喀斯特面積的37.1%[3]．因此，為了治理石漠化以及解決該現(xiàn)象所衍生出的諸多生態(tài)經(jīng)濟等難題，各界需要進一步加深對貴州省喀斯特石漠化的研究力度．

石漠化治理需要大量的生態(tài)環(huán)境數(shù)據(jù)、自然背景數(shù)據(jù)、社會經(jīng)濟數(shù)據(jù)等作為支撐，依托精準的遙感影像數(shù)據(jù)進行石漠化地區(qū)土地利用類型圖、土壤侵蝕圖等信息解譯和石漠化等級判定是石漠化治理過程中必不可少的工作環(huán)節(jié)．但是貴州省喀斯特石漠化地區(qū)多云多雨的氣候特征導致了在該地區(qū)難以獲得有用的高質(zhì)量的光學遙感數(shù)據(jù)．與光學遙感相比較，合成孔徑雷達具備不受大霧或者云層影響的優(yōu)點，穿透性和抗干擾性都較強，對紋理和地形的解譯效果較好[4]．現(xiàn)階段，國內(nèi)對于石漠化信息提取的研究多為嘗試利用光學遙感數(shù)據(jù)去開展，但很少涉及雷達遙感數(shù)據(jù)的利用，如2010年李麗等[5]利用中巴衛(wèi)星光學遙感數(shù)據(jù)提出了一種基于植被覆蓋度的石漠化信息提取方法，對云南隆林各族自治縣完成了石漠化信息提??；2016年朱大運等[6]利用GF-1與Landsat-OLI光學遙感數(shù)據(jù)，對不同影像的植被指數(shù)分類能力進行了定量評價．

同時，國內(nèi)對于利用雷達遙感數(shù)據(jù)去提取石漠化信息多基于單一時相的數(shù)據(jù)，解譯效果并不理想．利用多時相雷達遙感數(shù)據(jù)的后向散射系數(shù)去構(gòu)建時間序列的研究多集中于農(nóng)作物信息提取、地表檢測等，嘗試其進行石漠化信息提取的研究較少，如2009年錢峻屏等[7]利用Radarsat雷達數(shù)據(jù)時間序列，通過密度異常檢測的方法對南方多云雨地區(qū)的地面形變進行了檢測；2015年鐘禮山等[8]利用SAR時間序列對江蘇省徐州市進行了耕地信息的提取研究.所以，目前利用雷達遙感數(shù)據(jù)的后向散射系數(shù)構(gòu)建時間序列去提取石漠化信息和利用光學遙感數(shù)據(jù)與雷達遙感數(shù)據(jù)的融合進行石漠化強度定級的研究都尚屬于探索階段．基于上述，本研究嘗試利用合成孔徑雷達數(shù)據(jù)，通過構(gòu)建后向散射系數(shù)值時間序列的方法進行研究區(qū)石漠化信息的提取，然后結(jié)合Landsat8遙感數(shù)據(jù)以及石漠化等級分類標準進行石漠化強度定級.

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

1.1 研究區(qū)概況

本研究選取了“畢節(jié)撒拉溪喀斯特高原山地輕—中度石漠化綜合治理示范區(qū)”和“關嶺—貞豐花江喀斯特高原峽谷中—強度石漠化綜合治理示范區(qū)”兩個區(qū)域作為研究區(qū)，研究區(qū)分別位于貴州省的西北部和西南部地區(qū)，都具備了喀斯特地貌發(fā)育典型且喀斯特面積分布廣泛的特點[9]，且各自代表了喀斯特高原山地輕-中度石漠化和喀斯特高原峽谷中-強度石漠化這兩種典型的地貌，兩個示范區(qū)的位置分布見圖1.

圖1 研究區(qū)位置圖Fig. 1 Location of study areas

畢節(jié)撒拉溪喀斯特高原山地輕—中度石漠化綜合治理示范區(qū)(簡稱“撒拉溪示范區(qū)”)位于畢節(jié)市西部的撒拉溪鎮(zhèn)和野角鄉(xiāng)境內(nèi)，經(jīng)緯跨度為105°1′00″E至105°09′30″E、27°39′11″N至27°18′00″N．研究區(qū)行政區(qū)劃涉及2個鄉(xiāng)鎮(zhèn)，9個行政村，52個村民組，包括撒拉溪鎮(zhèn)的朝營村、鐘山村、沖鋒村、永豐村、龍鳳村、沙樂村、楊柳、撒拉溪和野角鄉(xiāng)的茅坪村[10]．研究區(qū)屬北亞熱帶濕潤季風氣候，植被類型以旱地、荒草地及稀疏灌木林地為主，土壤類型以沙壤和黃壤為主．區(qū)內(nèi)土地利用類型以耕地和林地為主，水土流失整體狀況以微度和中度為主．多年平均降雨量為984.40 mm，全年平均氣溫約為12 ℃．海拔高度為1500至2300 m，整體地勢東高西緩，地表呈現(xiàn)斜坡丘陵自然延伸的形態(tài)，海拔差距較大，是典型的喀斯特高原山地輕度石漠化地區(qū)．

關嶺—貞豐花江喀斯特高原峽谷中-強度石漠化綜合治理示范區(qū)(簡稱“花江示范區(qū)”)位于貴州安順市關嶺縣與黔西南自治州貞豐縣交界處的北盤江峽谷花江段，經(jīng)緯跨度為105°36′30″至105°46′30″E、25°39′13″至25°41′00″N．研究區(qū)行政區(qū)劃涉及3個鄉(xiāng)鎮(zhèn)，8個行政村，50個村民組，包括北盤江鎮(zhèn)的水淹壩、查爾巖、云洞灣和板貴鄉(xiāng)的木工、壩山、三家寨、孔落箐以及花江鎮(zhèn)的五里村[11]．研究區(qū)屬亞熱帶季風濕潤氣候，植物類型多以次生灌木林為主，土壤類型以黃壤、黃色石灰土為主．區(qū)內(nèi)土地利用類型以林地和耕地為主，水土流失整體狀況以微度和輕度為主．多年平均降雨量為1100 mm，全年平均氣溫約為18.4 ℃，海拔高度為450至1450 m，整體地勢陡峭，是典型的喀斯特高原峽谷強度石漠化地區(qū)．

1.2 雷達數(shù)據(jù)選擇

本研究的SAR數(shù)據(jù)選擇Sentinel-1所載的C波段SAR．Sentinel-1是GMES的重要組成部分，它由2顆衛(wèi)星構(gòu)成，重訪周期為6 d，空間分辨率較高．它是搭載了C波段(5.405 GHz)的SAR傳感器，它有4種成像模式：條帶模式、干涉寬幅模式、超寬幅模式以及波模式．本研究中撒拉溪示范區(qū)的SAR數(shù)據(jù)獲取時間分別為2016年1月12日、2月3日、3月15日、4月26日、5月23日、6月7日、7月13日、8月19日、9月30日、10月14日、11月13日和12月22日；花江示范區(qū)的SAR數(shù)據(jù)獲取時間分別為2016年1月9日、2月14日、3月9日、4月14日、5月11日、6月15日、7月18日、8月21日、9月13日、10月17日、11月8日和12月11日.本研究的SAR數(shù)據(jù)均采用VV極化方式．

1.3 光學數(shù)據(jù)選擇

本研究的光學數(shù)據(jù)選擇Landsat8衛(wèi)星所成的影像．Landsat8由美國國家航空航天局研發(fā)，它攜帶有兩個主要荷載，分別是陸地成像儀(operational land imager)和熱紅外傳感器(thermal infrared sensor)．Landsat8的陸地成像儀一共有9個波段，除全色的波段8空間分辨率是15 m外，其他波段的空間分辨率均為30 m．本研究中撒拉溪示范區(qū)和花江示范區(qū)分別選取2016年6月14日和10月4日的Landsat8陸地成像儀所成的影像作為光學遙感數(shù)據(jù)源．

2 提取方法

2.1 數(shù)據(jù)預處理

由于24景SAR數(shù)據(jù)的圖像屬性以及時相特征不一致，所以不能夠直接把他們的圖像信息有規(guī)律地聯(lián)系起來．為了讓各圖像的原始信號幅度值能夠準確量化地表達出來，使得不同時相的SAR數(shù)據(jù)具備可比性，必須要對影像進行輻射定標，將原始振幅信息轉(zhuǎn)換為后向散射系數(shù)[12]．其公式為：

A·σo+Pn,

(1)

地形輻射校正的目的是通過消除遙感影像的幾何變形，使多時相的SAR遙感數(shù)據(jù)能夠互相匹配，從而提高數(shù)據(jù)的可對比性和可用性以及分析的準確性．本研究中的研究區(qū)生態(tài)環(huán)境復雜，地形破碎，加上SAR采用傾斜觀測的成像方式，導致影像中地形畸變較大．本研究采用ASTER GDEM(先進星載熱發(fā)射和反射輻射儀全球數(shù)字高程模型)進行地形輻射校正，其空間分辨率為30 m．本研究采用薛東劍等2017年提出的校正模型[13]，其公式為：

(2)

式中：β0為成像面散射單元對應的后向散射系數(shù)，α1為距離向坡度角，α2為方位向坡度角，δ為參考入射角，θ為當?shù)厝肷浣牵?/p>

本研究采用精致Lee算子方法對SAR進行斑點噪聲的濾波，這是一種基于邊緣檢測的算法．其原理是若檢測范圍內(nèi)像元值較為一致，則可以認為該地區(qū)為勻質(zhì)區(qū)域，此時利用低通濾波方法．反之，則采用高通濾波方法，這樣可以最大程度上保留影像的紋理細節(jié)信息．

2.2 SAR后向散射系數(shù)時間序列的構(gòu)建

SAR遙感數(shù)據(jù)后向散射系數(shù)時間序列構(gòu)建是基于影像中同一個像元在不同時相數(shù)據(jù)中像元值的排列，為了能使各時相SAR數(shù)據(jù)的配準精度能夠控制在0.5個像元誤差之內(nèi)[14]，需要在之前數(shù)據(jù)預處理的基礎上進行幾何精校配的步驟．幾何精校配在ArcMap中進行，該操作以人工配準的方法，從兩個示范區(qū)的12時相影像中選取質(zhì)量較好的影像作為主影像，其余11時相影像作為輔影像，依次完成基于主影像的配準．設第m個點的配準誤差為R，R的計算公式為[20]：

(3)

式中：x和y為第m個控制點的在主影像上的經(jīng)度和緯度坐標，xi和yi是完成幾何精校配之后第m個控制點在第i個時相在輔影像上的經(jīng)度和緯度坐標．將配準后的24景影像按照示范區(qū)分區(qū)和時間順序合成兩幅影像，使得不同時相影像中的像元值一一對應．合成的影像中，X和Y軸分別代表經(jīng)緯度，Z軸作為時間軸，那么每個像元則構(gòu)成了后向散射系數(shù)的時間序列．SAR時間序列表達式為：

T={[(x1,y1),V1],[(x1,y1),V1],

…,[(xm×n,ym×n),Vm×n]},

(4)

Vi={((v1,t1),(v1,t1),…,(vn,tn)},

(5)

式中：(xi，yi)為第i個像元的行列坐標位置，m和n為SAR影像的行數(shù)和列數(shù)，Vi是基于第i個像元構(gòu)建的時間序列，Vj為第i個像元處在時刻為tj的SAR影像后向散射系數(shù)值．

2.3 動態(tài)時間調(diào)整算法提取石漠化信息

DTW(dynamic time warping，動態(tài)時間調(diào)整)算法，能夠?qū)r間序列中的時間軸實現(xiàn)動態(tài)彎曲，對像元的后向散射系數(shù)在各個時間節(jié)點的變化具有很高的識別和匹配精度，它在解決尺度偏移的問題上具備一定的優(yōu)勢，而且在篩選和過濾畸變數(shù)值方面有較好的效果[15]．DTW算法提取石漠化信息的步驟為：首先，利用典型石漠化地物像元的后向散射系數(shù)時間序列作為石漠化地物標準時間序列；接著計算研究區(qū)內(nèi)每個待分類像元的時間序列曲線與石漠化地物標準時間序列曲線之間的DTW距離，其算法公式為[16]：

D(1,1)=α11,

(6)

D(i,j)=αij+min{D(i-1,j-1),D(i,j-1),

D(i,j-1),D(i-1,j)}.

(7)

式中：i=2，3，…，m；j=2，3，…，n；D(i，j)為待分類像元的時間序列曲線與標準時間序列曲線距離的最小累加值，即DTW值，又稱為最佳彎曲路徑．然后，以目視解譯的方法選取研究區(qū)石漠化地區(qū)邊緣位置的若干像元，我們稱之為混合像元[17]．接著計算混合像元在不同時間節(jié)點的后向散射系數(shù)值的平均值，這樣就構(gòu)建了一條新的基于12時相的混合像元時間序列曲線．最后計算混合像元時間序列和石漠化地物標準時間序列的DTW值，以該值作為研究區(qū)待分類地物類型判斷的閾值[18]．若待分類像元的時間序列曲線與石漠化地物標準時間序列曲線的DTW值小于該閾值，則可以認為該像元屬于石漠化地物類型；若待分類像元的時間序列曲線與石漠化地物標準時間序列曲線的DTW值大于該閾值，則可以認為該像不屬于石漠化地物類型．

2.4 SAR與遙感數(shù)據(jù)融合的石漠化強度定級

將SAR與光學遙感數(shù)據(jù)在ENVI中進行融合操作，然后以3S(遙感技術(shù)、地理信息系統(tǒng)、全球定位系統(tǒng))技術(shù)作為支撐，結(jié)合坡度、基巖裸露率、土壤信息、植被覆蓋等作為主要判斷因子，輔以平均降水量、土層厚度、水土流失狀況等因子創(chuàng)建石漠化強度等級分類標準(表1)，最后利用ArcMap的空間分析功能對研究區(qū)的石漠化強度進行基于圖斑識別的定級．

表1 石漠化強度等級分類標準Tab. 1 Classification standard of rock desertification intensity grade

3 結(jié)果與討論

3.1 基于像元的后向散射系數(shù)分析

對于已構(gòu)建的SAR時間序列，在影像中的石漠化地區(qū)選取訓練樣本，訓練樣本的選擇應該遵循均勻地落在研究區(qū)內(nèi)，且盡量能滿足對整個研究區(qū)的代表性的選擇原則．本研究在各示范區(qū)都選擇了50個訓練樣本，均為6像元×6像元的正方形區(qū)域．下圖為撒拉溪示范區(qū)部分訓練樣本在12時相(T1-T12)的后向散射系數(shù)值時間序列曲線圖(見圖2，篇幅有限不能全部列出)．

圖2 部分訓練樣本后向散射系數(shù)值時間序列曲線圖Fig. 2 Parts of training samples’backscatteringcoefficients time series trace

從數(shù)據(jù)上分析，撒拉溪示范區(qū)訓練樣本12個時相中的后向散射系數(shù)值大概集中在-24 db到-18 db之間，花江示范區(qū)訓練樣本12個時相中的數(shù)值集中在-27 db到-20 db之間，整體上數(shù)值大小范圍為1月到6月數(shù)值較低，7月到9月數(shù)值較高，10月到12月數(shù)值再次呈現(xiàn)降低的趨勢，其中，7月和8月數(shù)值最大，1月和12月達到最小值．究其原因，SAR后向散射系數(shù)值大小的主要影響因子主要有土壤水分、地表粗糙度以及植被覆蓋等，而訓練樣本所在區(qū)域基巖裸露情況嚴重，植被覆蓋較少，導致石漠化區(qū)域后向散射系數(shù)值一整年都比較平穩(wěn)，變化趨勢較小．花江示范區(qū)訓練樣本后向散射系數(shù)值較之撒拉溪示范區(qū)偏高一點的原因是該地區(qū)石漠化強度以中度和強度為主．而數(shù)值變化的原因在于數(shù)值高的月份太陽輻射較強，植被覆蓋有所增強，土壤濕度也隨之上升，后向散射系數(shù)值偏高；數(shù)值低的月份太陽輻射較弱，溫度降低，植被稀疏，土壤水分凍結(jié)，導致土壤介電常數(shù)下降，后向散射系數(shù)值偏低．

3.2 SAR提取石漠化信息結(jié)果

在24時相的影像中目視選擇質(zhì)量最好的2016年7月13日的撒拉溪示范區(qū)雷達數(shù)據(jù)和2016年9月5日的花江示范區(qū)示范區(qū)雷達數(shù)據(jù)與選擇的光學遙感數(shù)據(jù)進行影像融合，分別以51和42為DTW算法的閾值對撒拉溪示范區(qū)以及花江示范區(qū)進行石漠化信息提取，并對提取的結(jié)果進行二值化處理，結(jié)果見圖3和圖4.石漠化信息統(tǒng)計見表2.

表2 基于SAR的撒拉溪示范區(qū)和花江示范區(qū)石漠化信息提取結(jié)果Tab. 2 Information extraction of rocky desertification basedon SAR data in Salaxi and Huajiang demonstration areas

圖3 基于SAR的撒拉溪示范區(qū)石漠化信息提取結(jié)果Fig. 3 Information extraction of rocky desertificationbased on SAR data in Salaxi demonstration area

圖4 基于SAR的花江示范區(qū)石漠化信息提取結(jié)果Fig. 4 Information extraction of rocky desertification ofHuajiang demonstration area based on SAR data

3.3 數(shù)據(jù)融合的石漠化強度定級結(jié)果

基于SAR和光學遙感數(shù)據(jù)融合的研究區(qū)石漠化強度定級結(jié)果見圖5和圖6.

圖5 基于SAR和光學遙感數(shù)據(jù)融合的撒拉溪示范區(qū)石漠化強度定級結(jié)果Fig. 5 Result of defining the levels of rock desertification intensitybased on SAR and optical data fusion in Salaxi demonstration areas

圖6 基于SAR和光學遙感數(shù)據(jù)融合撒拉溪示范區(qū)石漠化強度定級結(jié)果Fig. 6 Result of defining the levels of rock desertification intensity ofHuajiang demonstration area based on SAR and optical data fusion

基于SAR和光學遙感數(shù)據(jù)融合的研究區(qū)各等級石漠化面積提取結(jié)果見表3.

表3 基于SAR和光學遙感數(shù)據(jù)融合的撒拉溪示范區(qū)各等級石漠化面積提取結(jié)果Tab. 3 Area extraction of different levels of rock desertification intensitybased on SAR and optical data fusion in Salaxi demonstration areas

3.4 DTW算法及提取結(jié)果精度驗證

為了驗證DTW算法的精度，本研究利用皮爾森相關系數(shù)算法和歐式距離算法對像元同樣進行閾值分割處理，然后提取石漠化和非石漠化面積.為了驗證提取結(jié)果的精度，本研究利用光學遙感數(shù)據(jù)對研究區(qū)石漠化強度進行等級確定，數(shù)據(jù)和定級方法與前文一致，分別選擇2016年6月14日、10月4日的Landsat8數(shù)據(jù)，定級方法同樣采用基于圖斑識別的方法，定級結(jié)果見表4.

表4 基于光學遙感數(shù)據(jù)的撒拉溪示范區(qū)和花江示范區(qū)各等級石漠化面積提取結(jié)果Tab. 4 Area extraction of different levels of rock desertificationintensity based on optical data in Salaxi and Huajiang demonstration areas

對于DTW算法及提取結(jié)果精度的驗證，本研究采用野外實地驗證的方法．首先，在兩個示范區(qū)分別選取200個精度驗證訓練樣本．訓練樣本的選擇遵循均勻分布且能盡量能代表整個研究區(qū)的原則．選取方式可以在地圖上以目視的方式進行．然后利用GPS記錄下整個野外實地驗證的行程路徑，用數(shù)碼相機記錄訓練樣本的照片并完成照片編號和該點石漠化等級的確定．最后，利用HOLUX ezTour軟件，把行程路徑和拍攝的照片文件轉(zhuǎn)化為kml文件，再在谷歌地球中打開kml文件，作為DTW算法及提取結(jié)果精度驗證的依據(jù)[19]．本研究選取石漠化界定正確率(該訓練樣本是否屬于石漠化地物)和石漠化強度等級正確率兩個指標作為精度驗證的結(jié)果．結(jié)果見表5和表6.

從精度驗證結(jié)果中分析，DTW算法在兩個示范區(qū)的石漠化界定精度均要高于皮爾森相關系數(shù)算法和歐式距離算法；撒拉溪示范區(qū)基于SAR后向散射系數(shù)時間序列的提取方法的石漠化界定正確率為91%，比光學遙感數(shù)據(jù)提取方法的精度高7.69%；花江示范區(qū)基于SAR后向散射系數(shù)時間序列的提取方法的石漠化界定正確率為91%，比光學遙感數(shù)據(jù)提取方法的精度高9.41%；撒拉溪示范區(qū)SAR與光學遙感數(shù)據(jù)融合提取方法的石漠化界定正確率為88.5%，比光學遙感數(shù)據(jù)提取方法的精度高4.73%，石漠化強度等級界定正確率為83.5%，比光學遙感數(shù)據(jù)提取方法的精度高5.99%；花江示范區(qū)SAR與光學遙感數(shù)據(jù)融合提取方法的石漠化界定正確率為92.5%，比光學遙感數(shù)據(jù)提取方法的精度高8.82%，石漠化強度等級界定正確率為88%，比光學遙感數(shù)據(jù)提取方法的精度高5.39%.

表5 撒拉溪示范區(qū)和花江示范區(qū)不同算法的精度驗證Tab. 5 Accuracy verification of different algorithms inSalaxi and Huajiang demonstration areas

表6 撒拉溪示范區(qū)和花江示范區(qū)石漠化信息提取結(jié)果精度驗證Tab. 6 Accuracy verification of information extraction of rockydesertification in Salaxi and Huajiang demonstration areas

4 結(jié)論

本研究利用研究區(qū)多時相SAR數(shù)據(jù)集合，利用構(gòu)建后向散射系數(shù)時間序列曲線的方式進行石漠化信息的提取，該方法使用DTW算法確定閾值從而對石漠化地物像元和非石漠化像元進行分割，利用SAR和光學遙感數(shù)據(jù)，結(jié)合石漠化強度等級分類標準表進行石漠化強度定級．實驗結(jié)果表明：1)石漠化地物的后向散射系數(shù)值一年內(nèi)趨于穩(wěn)定，總體變化趨勢較小，這也反證了利用SAR的后向散射系數(shù)時間序列提取石漠化信息的可行性與精準性；2)對比基于光學遙感數(shù)據(jù)的提取方法，基于SAR后向散射系數(shù)時間序列的提取方法在石漠化地物和非石漠化地物分割上精度更高；3)對比基于光學遙感數(shù)據(jù)的提取方法，基于多時相SAR與光學遙感數(shù)據(jù)融合的提取方法在石漠化強度定級上精度更高．

研究方法在一定程度上克服了使用單一時相SAR數(shù)據(jù)的缺陷，同時集雷達遙感和光學遙感之所長，兼顧雷達遙感成像中蘊含的豐富的紋理地形信息以及光學遙感成像中蘊含的豐富的光譜信息，使得石漠化信息提取精準度有所提高．