王曉慧,譚炳香,李世明,馮林艷
(中國林業(yè)科學(xué)研究院資源信息研究所,國家林業(yè)和草原局林業(yè)遙感與信息技術(shù)重點實驗室 北京 100091)
遙感技術(shù)以其觀測范圍廣、信息量大、精度高和速度快等特點為森林資源變化檢測提供了便利手段,已經(jīng)在森林資源調(diào)查監(jiān)測中發(fā)揮著重要作用[1-6]。通過前后時相影像的變化檢測可以提取森林資源變化信息[7-11],通過時間系列影像能夠識別森林干擾[12-14]。隨著遙感影像空間分辨率的提高,影像中地物的結(jié)構(gòu)、形狀、紋理和細(xì)節(jié)等信息也更加突出,同種地物內(nèi)部光譜變異性增大,空間結(jié)構(gòu)特征更加明顯。面向?qū)ο蟮膱D像處理方法通過影像分割技術(shù)把影像分解成同質(zhì)像元組,即影像對象。影像對象充分利用像元間的空間特征關(guān)系,避免椒鹽現(xiàn)象,因此面向?qū)ο蠖嗵卣魅诤系姆椒ǔ蔀楫?dāng)前遙感影像變化檢測的重要手段[15-16]。
面向?qū)ο蟮淖兓蛄糠治龇ㄊ且环N基于圖斑的特征向量空間變化檢測方法,由各個特征值的變化量組成變化向量,以歐氏距離計算變化強度,通過選擇合適的閾值判斷地物是否變化[17-19]。該方法中的特征值通常包括影像對象的光譜特征和紋理特征,其中光譜特征又包括了各個波段的光譜均值、光譜標(biāo)準(zhǔn)差等,紋理特征包括了灰度共生矩陣的均值、同質(zhì)性、對比度、差異性和熵等[20-21]。由于特征值眾多,減少數(shù)據(jù)冗余,增強特征之間的理解,選擇有效識別地物特征的信息,是基于變化向量分析法的變化檢測中的一項重要內(nèi)容。廣西壯族自治區(qū)上思縣桉樹(Eucalyptus robustaSmith)速生豐產(chǎn)林生長快,輪伐期短,森林采伐和造林活動多,森林資源變化頻率高。本研究在上思縣設(shè)立研究區(qū),基于面向?qū)ο笞兓蛄糠治龇ㄟM行森林資源變化檢測,應(yīng)用隨機森林平均精確率減少的方法進行變化特征的選擇,通過選取的不同特征向量和常規(guī)的基于光譜均值、光譜均值和標(biāo)準(zhǔn)差的變化向量分析法進行比較,以期獲得較好的森林資源變化檢測方法和結(jié)果。
研究區(qū)位于廣西壯族自治區(qū)上思縣東北部,108°11′~108°16′ E,22°12′~22°14′ N 之間,面積為2 667.45 hm2(圖1)。屬于丘陵地形,平均海拔220 m。屬南亞熱帶季風(fēng)氣候,氣候溫暖,無霜期長,雨量適中,干濕季明顯。年平均日照時數(shù)為1 896.1 h,年平均氣溫21.7 ℃,年平均降水量1 217.3 mm。
圖1 研究區(qū)地理位置及影像(紅:band 4;綠:band 3;藍(lán):band 2)Fig.1 Geographical location and images of study area (R:band 4;G:band 3;B:band 2)
獲取了研究區(qū)兩期影像,分別為2015 年8 月24 日和2016 年10 月11 日的國產(chǎn)高分二號多光譜影像,包含藍(lán)、綠、紅、近紅外4 個波段,空間分辨率為4 m(圖1)。
應(yīng)用中國資源衛(wèi)星應(yīng)用中心高分二號衛(wèi)星絕對輻射定標(biāo)系數(shù)和多光譜相機的歸一化響應(yīng)度文件,對兩景影像分別進行輻射定標(biāo)和大氣校正,獲取大氣表觀反射率。利用30 m 分辨率的DEM 數(shù)據(jù)并結(jié)合高分二號影像自帶的RPC 文件對影像進行正射校正,消除地形引起的幾何畸變。以2016 年影像作為基準(zhǔn)影像,通過從兩景影像上選擇控制點的方法來配準(zhǔn)2015 年影像,使用最近鄰法進行重采樣,配準(zhǔn)誤差在0.5 個像元以內(nèi)。以2016 年影像作為基準(zhǔn)影像,在兩景影像上選取不變地物點,對相應(yīng)波段分別建立線性函數(shù)關(guān)系,對2015 年影像做相對輻射校正。從經(jīng)過校正的影像裁剪出研究區(qū)。
據(jù)2014 年度森林資源數(shù)據(jù)更新調(diào)查結(jié)果,上思縣林地和非林地的面積比例分別為65.5%和34.5%。主要樹種馬尾松(Pinus massonianaLamb.)、桉樹(Eucalyptus robustaSmith)和一般闊葉樹分別占林地面積的41.2%、27.3%和23.9%。森林資源的主要變化是桉樹的造林更新,以及桉樹和馬尾松的采伐,林地變化頻繁,范圍廣,面積大。桉樹輪伐期短,一般6~8 a,有的甚至縮短到4~5 a,桉樹采伐后萌芽更新[22]。由于桉樹速生豐產(chǎn)林經(jīng)濟效益高,有些地區(qū)采伐馬尾松,當(dāng)年春夏或者第二年春天通過植苗造林更新為桉樹,桉樹種植面積和比例提高很多。
本研究中無林地包括采伐跡地、未成林地和宜林荒山荒地,有林地指馬尾松林、桉樹林和一般闊葉林。針對于高分辨率遙感影像的森林資源變化類型,包括有林地轉(zhuǎn)變?yōu)闊o林地,無林地轉(zhuǎn)變?yōu)橛辛值?,以及馬尾松林轉(zhuǎn)變?yōu)殍駱淞?,桉樹幼林轉(zhuǎn)變?yōu)槌闪郑▓D2)。農(nóng)田、水體、城鎮(zhèn)等非林地不在本研究的森林資源變化檢測范圍內(nèi)。
圖2 森林類型影像特征Fig.2 Image characteristics of forest types
對前后兩時相高分二號影像進行疊加組合,形成多時相組合影像。應(yīng)用eCognition Developer 9.4軟件中的多尺度分割算法對組合影像進行分割。從像元開始進行自下而上的合并,連續(xù)地合并像元或現(xiàn)有的影像對象,把小對象組合成較大的對象。確定合適的分割參數(shù),包括分割尺度、各波段權(quán)重、光譜因子、形狀因子和緊致度因子。理想的影像分割結(jié)果是,分割得到的圖斑內(nèi)部具有相似的光譜和紋理特性,相鄰的空間位置關(guān)系,具有良好的同質(zhì)性,圖斑與鄰接圖斑之間具有良好的異質(zhì)性,保證圖斑之間的可分性強。
前后兩時相影像圖斑的各個特征值組成特征矢量。將基準(zhǔn)影像的特征矢量記為X,檢測影像的特征矢量記為Y。X和Y分別如下:
其中,n為特征值數(shù),變化矢量 ?X為:
采用歐氏距離計算變化強度|?C|:
其中,|?C|包含了從基準(zhǔn)時期到檢測時期各個特征值變化的總和,表示兩期影像的特征差異。一般來說,當(dāng)|?C|較大時,兩期影像之間的差異較大,對應(yīng)圖斑的地物類型發(fā)生變化的概率較大,反之則較小。
根據(jù)采伐地塊、造林地塊和各種森林類型在高分二號影像的表征,針對前后時相影像圖斑分別構(gòu)建初始特征空間,每期影像圖斑分別包括29 個特征,其中光譜信息包括4 個波段的均值和標(biāo)準(zhǔn)差,以及NDVI 值,紋理信息包括4 個波段的灰度共生矩陣的能量、對比度、相關(guān)性、熵和同質(zhì)性。
采用線性函數(shù)轉(zhuǎn)換方法將影像的各個特征數(shù)據(jù)歸一化到[0,1]之間。Y=(X?Xmin)/(Xmax?Xmin),其中,X為歸一化前的圖斑特征值,Y為歸一化后的圖斑特征值,Xmax和Xmin為歸一化前的圖斑特征值的最大值和最小值。
面向?qū)ο蟮淖兓蛄糠治龇ㄖ校總€圖斑用多個特征組成的多維特征向量進行描述。為了減少數(shù)據(jù)相關(guān)性和冗余度,增強算法的可解釋性,選擇出相關(guān)性小、識別地物變化好的特征參與變化檢測。
變化檢測中,主要關(guān)注2 個時相特征值的差異部分。本研究采用歐氏距離計算變化強度,變化特征是對兩個時相影像歸一化后的特征做差值。試驗區(qū)內(nèi)的300 個樣本點,包括55 個變化地類,245個未變化地類。變化地類中23 個為無林地造林后成為桉樹林,27 個為馬尾松、混交林或者桉樹采伐后成為無林地,5 個為桉樹幼林地轉(zhuǎn)變?yōu)槌闪值亍?yīng)用隨機森林平均精確率減少的方法對29 個變化特征進行篩選。主要是打亂每個變化特征順序,度量順序變動對模型精確率的影響。該方法穩(wěn)定性較好,通過直接度量每個變化特征對模型精確率的影響來確定變化特征的重要性。隨機森林中樹的數(shù)量采用默認(rèn)值100,按照7∶3 將數(shù)據(jù)分為訓(xùn)練集和測試集。
應(yīng)用ROC 曲線來確定變化強度閾值,判斷圖斑地物屬性是否發(fā)生變化。ROC 曲線又稱受試者工作特征曲線,廣泛應(yīng)用于各種診斷試驗性能的評價[23-24]。
對于二分問題,本研究中設(shè)定變化圖斑為正類,未變化圖斑為負(fù)類,會出現(xiàn)4 種情況。如果一個圖斑發(fā)生變化并且也被檢測為變化,則為真正類(TP),未變化圖斑被檢測為變化則為假正類(FP),未變化圖斑被檢測為未變化則為真負(fù)類(TN),變化圖斑被檢測為未變化則為假負(fù)類(FN)。真正類率(TPR)表示變化檢測算法所識別出的變化圖斑占所有變化圖斑的比例,計算公式為TPR=TP/(TP+FN)。假正類率(FPR)表示變化檢測算法錯分為變化的未變化圖斑占所有未變化圖斑的比例,計算公式為FPR=FP/(FP+TN)。
ROC 曲線法根據(jù)一系列不同的臨界值,以真正類率(敏感度)為縱坐標(biāo),假正類率(1-特異度)為橫坐標(biāo)繪制曲線。最靠近左上角ROC 曲線的點是錯誤最少的最佳閾值,其真正類率盡可能高和假正類率盡可能低。ROC 曲線下的面積(AUC)是重要的試驗準(zhǔn)確度指標(biāo)。AUC的取值范圍在0.5和1 之間。AUC越接近1.0,檢測方法真實性越高;AUC等于0.5 時,則真實性最低,無應(yīng)用價值。
生成ROC 曲線和確定閾值分為以下3 步:(1)采用特征選取時建立的樣本點作為測試樣本;(2)構(gòu)建變化向量后獲得的變化強度作為檢驗變量;(3)SPSS 統(tǒng)計軟件提供所有臨界值的表,每一個臨界值對應(yīng)不同的敏感度和特異度,以約登指數(shù)即(敏感度+特異度?1)最大時所對應(yīng)的點為最佳閾值。
依據(jù)選擇的變化特征以及常規(guī)的變化特征,分別構(gòu)建變化向量,根據(jù)閾值獲得變化檢測結(jié)果,即區(qū)分為變化地類或者未變化地類。試驗區(qū)內(nèi)按照隨機分布建立737 個驗證樣本點,一個圖斑只對應(yīng)一個驗證樣本點。對照兩期高分二號影像以及Google Earth 上的高空間分辨率影像判別樣本點的變化情況,其中138 個變化地類,599 個未變化地類。采用同一套驗證樣本點分別對不同特征組合得出的變化檢測結(jié)果進行精度驗證,計算總體精度、Kappa系數(shù)、誤檢率和漏檢率。
對前后兩時相高分二號影像形成的8 波段組合影像進行多尺度分割,經(jīng)過試驗測試和人工目視評判,確定分割參數(shù),獲得前后時相空間位置一致的圖斑作為變化檢測的基本單元(表1)。多尺度分割自下而上通過分割創(chuàng)建出影像對象圖層,圖層1 是像元分割后形成的圖層,圖層2 是圖層1 對象執(zhí)行分割后形成的圖層。
表1 影像分割參數(shù)Table 1 Parameters of image segmentation
鑒于試驗區(qū)采伐和造林地塊大小不一,有些比較破碎,為了能完全識別出變化地塊,分割時考慮適度的過分割。從分割圖可以看出,多時相影像組合分割結(jié)果用一個或多個圖斑表示發(fā)生變化的地物,分割線清楚地勾繪出了不同地物的邊界(圖3)。
圖3 局部分割影像和變化檢測Fig.3 Partial segmentation images and change detection
應(yīng)用隨機森林平均精確率減少的方法進行變化特征的選擇,按照29 個變化特征得分從高到低,獲得變化特征重要性排序(圖4),特征代碼對照見表2。前10 個變化特征為:第3 波段光譜均值、第1 波段光譜均值、NDVI 值、第2 波段光譜均值、第1 波段光譜標(biāo)準(zhǔn)差、第3 波段光譜標(biāo)準(zhǔn)差、第2 波段光譜標(biāo)準(zhǔn)差、第1 波段熵、第4 波段光譜均值和第1 波段能量。
圖4 變化特征重要性Fig.4 Change feature importance
表2 特征代碼Table 2 Feature code
按照變化特征排序的前1 個、前2 個一直到29 個變化特征構(gòu)建變化向量,分別應(yīng)用ROC 法確定變化閾值(圖5),獲得變化檢測結(jié)果(圖3)。同時,進行對比研究,用常規(guī)的由4 個波段光譜均值、4 個波段光譜均值和4 個波段光譜標(biāo)準(zhǔn)差構(gòu)建變化向量,以及NDVI 差值法進行變化檢測,應(yīng)用ROC 法確定變化閾值。
圖5 前12 個變化特征向量法的ROC 曲線分析Fig.5 ROC curve of the first twelve change feature vectors
圖6 不同變化特征數(shù)量的面向?qū)ο笊仲Y源變化檢測的精度Fig.6 Accuracy of object-oriented forest change detection of various change feature numbers
對不同數(shù)量變化特征參與的變化向量法獲取的森林資源變化檢測結(jié)果進行比較(圖6)。變化特征數(shù)量從1 增加到3,總體精度和Kappa 系數(shù)均較大,總體精度大于90%,Kappa 系數(shù)約為0.7。變化特征數(shù)量為4 和5 時,總體精度和Kappa 系數(shù)均增大,總體精度大于92%,Kappa 系數(shù)大于0.76。變化特征數(shù)量從6 增加到9,總體精度大于91%,Kappa 系數(shù)約為0.7。變化特征數(shù)量增加到10 以后,總體精度和Kappa 系數(shù)呈現(xiàn)下降趨勢。并不是參與變化向量分析法的變化特征越多,識別森林資源變化的效果越好。
前4 個變化特征參與的變化向量法,即高分二號多光譜影像藍(lán)、綠、紅波段光譜均值和NDVI 值識別森林資源變化精度最高,總體精度92.94%,Kappa 系數(shù)0.763 0,變化地類誤檢率15.63%,漏檢率22.86%,未變化地類誤檢率5.25%,漏檢率3.35%(圖6 和表3)。
表3 面向?qū)ο笊仲Y源變化檢測的精度驗證Table 3 Accuracy assessment of object-oriented forest change detection
特征選擇后變化向量法識別森林資源變化的精度普遍高于常規(guī)的4 個波段光譜均值、4 個波段光譜均值和4 個光譜標(biāo)準(zhǔn)差的變化向量分析法的檢測結(jié)果(表3)。4 個變化特征到9 個變化特征參與變化向量分析法的總體精度在90%以上,另外兩種常規(guī)方法的總體精度分別為83.58%和77.88%。4 個變化特征向量法識別森林資源變化檢測的總體精度高于NDVI 差值法的總體精度,變化地類漏檢率由NDVI 差值法的34.29%降低到22.86%。針對4 個變化特征向量法,對照影像核對驗證樣本,一部分檢測錯誤的圖斑是由于圖斑中的局部發(fā)生了地類變化,這部分圖斑標(biāo)識為變化地塊,但是變化檢測算法應(yīng)用的是圖斑的光譜均值,局部變化和未變化部分光譜信息平衡后,削弱了變化強度,因此這部分變化圖斑未能被正確識別出來。該問題是由于分割尺度引起的,減小分割尺度參數(shù),增強圖斑內(nèi)部的同質(zhì)性,同時保持圖斑之間的異質(zhì)性,將有助于提高由于圖斑分割尺度參數(shù)稍大引起的變化檢測精度偏低的問題。
對于森林資源發(fā)生變化的4 種類型,應(yīng)用4 個變化特征向量法識別其中3 種變化類型效果好,即桉樹或者馬尾松采伐、無林地造林后成為有林地、桉樹從幼林到成林。有些變化地類誤識別為未變化地類,對照兩期影像目視分析,一部分檢測錯誤的圖斑在前期影像上是馬尾松,馬尾松采伐后種植桉樹,桉樹進入成熟期,在后期影像上是桉樹成林,在前后期影像上均表現(xiàn)為有林地的光譜特征。在變化特征選取時未包括馬尾松林轉(zhuǎn)變?yōu)殍駱涑闪值臉颖军c,選取的變化特征未能反映該變化信息,影響到變化檢測結(jié)果。有些未變化地類誤識別為變化地類,主要是受桉樹除草的影響,前后時相影像上桉樹的光譜發(fā)生變化,變化檢測受到同物異譜的干擾。
(1)本研究應(yīng)用高分二號多光譜數(shù)據(jù),以廣西壯族自治區(qū)上思縣林地為例,基于面向?qū)ο蟮淖兓蛄糠ǎ?jīng)過變化特征選擇,采用ROC 曲線劃分變化閾值,獲得了較好的森林資源變化檢測結(jié)果。
(2)應(yīng)用隨機森林平均精確率減少的方法進行變化特征選擇,4 個變化特征向量法,即藍(lán)、綠、紅波段光譜均值和NDVI 值參與的變化向量分析法,識別森林資源變化情況精度最高,總體精度92.94%,Kappa 系數(shù)0.763 0。
(3)特征選擇后變化向量分析法對森林資源的變化檢測結(jié)果,明顯優(yōu)于常規(guī)的4 個波段光譜均值、4 個波段光譜均值和4 個光譜標(biāo)準(zhǔn)差的變化向量分析法的檢測結(jié)果,特征選擇出的前4 個和前5 個變化特征向量法識別森林資源變化檢測的總體精度高于NDVI 差值法的總體精度。
(4)本研究方法用于森林采伐和造林活動多、林地變化頻率高的地區(qū)森林資源變化識別時,具有一定的優(yōu)勢和潛力。