基于機載CCD和ALS偽波形數據的山區地表分類研究

黃 侃 于 強 黃華國

(北京林業大學省部共建森林培育與保護教育部重點實驗室，北京 100083)

0 引言

地表覆被是地球表面各種地物類型及其屬性特征的綜合體[1]，其空間分布變化反映了人類社會的經濟活動發展過程以及生態格局的變化趨勢。同時，地表覆被的變化會引起地球表面物質和能量的循環過程的改變，繼而對區域及全球的環境與氣候產生重要影響[2]。因此，地表覆被變化研究是當今國際科學研究的前沿領域，也是我國全球變化研究所關注的重點領域。精準的地表覆被數據對于研究社會經濟發展的健康狀況、指導土地科學利用，以及研究區域乃至全球的生態環境和氣候變化具有重要意義[3-4]。遙感技術的出現極大方便了地表覆蓋分類及其變化監測，但精度已經基本達到了上限。而機載激光雷達(ALS)作為一種新興的高分辨率、高點密度的遙感空間數據獲取技術，在土地分類中具有重要應用。目前其應用主要集中于城市建筑提取分類[5-6]、森林分類[7-8]以及地質地貌分類與識別[9-10]等方面。

然而，ALS點云數據缺乏直接表達物體表面的語義信息(例如紋理和結構)，且空間三維點云數據本身具有不連續性、不規則性以及數據的密度不均勻等，直接利用ALS點云數據進行分類，操作較為復雜，且效果并不理想[5,11]。由于地物的復雜性、城市高大建筑的陰影問題以及光譜異質性問題的存在，單獨利用光譜數據進行土地利用分類同樣存在缺陷[12]，而結合ALS與光譜數據進行分類則具有更大的潛力。

基于此，已有許多研究嘗試通過ALS數據與光譜影像數據融合進行模式識別和地物提取與分類。其中，大多數研究均是通過ALS數據衍生出的歸一化地表模型nDSM或DEM和強度參數結合光譜數據衍生出的NDVI值或其他光譜特征參數對地表進行分類和特征提取[13-20]，所采用的分類方法主要包括支持向量機SVM、決策樹、面向對象方法及對應的改進方法，對建筑、水體、綠地等地物的分類都表現出較高的精度。袁楓[21]利用統計方法，將ALS數據高程、激光強度信息和紅綠藍波段進行分布統計，提取各種土地利用類別的特征值，利用框架理論對城鎮地區土地利用類型進行了分類。

目前，大多數類似研究僅針對城市區域，區域內地面相對規整，分類精度較高。山地區域由于地形起伏的影響，光學圖像容易產生陰影、同物異譜現象，而ALS數據則存在DEM的提取精度不足、光譜紋理信息缺乏等困難，光學圖像和ALS數據地表分類均面臨挑戰。針對山區復雜地形下的地表分類問題，本文從數據融合入手，提出一種基于ALS點云數據構建大光斑偽波形，并結合點云強度信息和CCD影像的RGB 3波段信息的地表分類方法，并在安徽黃山地區進行驗證。

1 研究區概況及研究數據

研究區位于安徽省黃山市南部，地理位置為北緯29°39′58.76″～29°41′33.45″，東經118°13′48.50″～118°23′3.09″。具體位置如圖1a所示。研究數據包括ALS點云數據(三維坐標和強度)、配套的高分辨率CCD影像以及GlobeLand30全球分類產品。主要地物類型包括農田、水體、村莊及森林。

1.1 點云數據和CCD

數據采集系統為Leica ALS60，飛行時間為2011年4月。ALS數據點云密度約為0.39點/m2，CCD數據空間分辨率約為1 m，共有15景影像。圖1b所示點云數據經過ENVI LiDAR處理生成的DSM圖像，圖1c為經過ENVI軟件幾何校正、精校正和配準拼接后的CCD影像。

1.2 全球30 m地表覆蓋(GlobeLand30)分類數據

GlobeLand30分類產品數據(http:∥www.globallandcover.com)是中國國家基礎地理信息中心研制的全球30 m地表覆蓋產品，其在環境變化分析、土地資源管理等方面已經得到廣泛應用[22-27]。該產品基于30 m多光譜影像，包括TM5、ETM+多光譜影像和中國環境減災衛星(HJ-1)多光譜影像等。數據采用WGS 84坐標系，UTM投影，6°分帶，參考橢球為WGS 84橢球。本文從官方網站上下載2011年圖幅號N50-25的 GlobeLand30數據，包含了研究區全部范圍。

2 研究方法

2.1 偽波形提取

偽波形是激光雷達點云的一種轉換方式，可以更好地利用波形特征進行分類。BLAIR等[28]在對熱帶森林冠層垂直結構研究中發現，利用小光斑激光雷達高程信息模擬的垂直結構和全波形裝置記錄的全波形數據基本相似。FARID等[29]在以小光斑激光雷達模擬的波形數據估計森林葉面積指數(LAI)的研究中，通過對比發現小光斑激光雷達模擬的波形和地基激光掃描儀的回撥波形相關程度能達到0.75。尤號田等[30-31]利用0.5 m為層間距分層統計小光斑激光雷達樣方內不同高程的點云數據，模擬偽波形并提取了平均樹高、LAI、郁閉度和生物量等森林結構參數，效果優于其他方法。

為了能夠更好地將偽波形與CCD影像的像元進行融合，假設偽波形的光斑近似于邊長為L的正方形。利用CSF濾波算法[32]對研究區ALS點云進行濾波并獲取研究區5 m分辨率DEM。在此基礎上，用光斑內所有點云高度減去DEM得到每個點的歸一化高程h。對所有樣地內歸一化后的點云高程h經過統計發現，絕大多數點云高程h均位于0～20 m范圍內，約占99.8%，因此為了統一偽波形的統計范圍與橫坐標刻度值，偽波形高程范圍設為0～20 m，分為10個區間段，每個區間段間隔2 m。將h從低到高以2 m為區間統計點云數量占光斑內總點云數量百分比，以此構成模擬的波形，即邊長為L光斑對應的偽波形fw，其表達式為

(1)

式中nx——光斑內高程h在[2(x-1),2x]范圍中點的數量

N——光斑內所有點的數量

通過設置不同大小的光斑(L=10～80 m)進行波形模擬便能得到一系列不同光斑大小的偽波形。偽波形的模擬公式與計算過程均基于Matlab程序實現。

本文所用的高程歸一化方法比較簡略，僅能大致消減地形起伏對實際地物高程的影響，為避免引起因分類地形所產生較大誤差，后續可以考慮加入更多相關因子或者算法以更高的分辨率進行歸一化。另外，由于水體能夠吸收激光脈沖，返回的能量幾乎為0，點云數量極少，基于此特性能夠很容易將水體從其他地類中區分，因此分析試驗中未加入水體類型。

2.2 CCD影像處理

對15幅有重疊的CCD影像分別進行幾何校正、精校正配準后，將多幅圖像無縫拼接成一幅覆蓋整個研究區的圖像，并將研究區域從中分割，為后續研究提供數據基礎，圖像處理步驟均在ENVI軟件中實現。拼接中只在重合區進行色彩融合調整。

2.3 特征曲線提取

利用同樣的方式對同一個光斑范圍內所有點的強度信息，以及光斑垂直投影面覆蓋范圍對應CCD影像數據相同區域范圍內的各個像元RGB 3波段值，分別按從最小值到最大值分為10個等間隔區間統計每個區間點(或像元)的數量占總點(或像元)數的百分比，構成與偽波形相類似的統計點云(或像元)相對數量分布曲線。其中，強度信息fI為未經過標定的相對度量，無量綱，經過統計分析發現，大多數點云強度都位于0～10之間，因此強度信息區間范圍設為0～10。強度信息fI可表示為

(2)

式中mx——光斑中強度I在[x-1,x]范圍內點的數量

CCD影像像元的RGB波段范圍均集中于50～150之間，因此所有像元的區間統計范圍設為50～150。RGB波段信息fR/G/B可表示為

(3)

式中n′x——光斑范圍內CCD影像中像元R/G/B波段值域范圍[10(x-1)+50,10x+50]內像元的數量

N′——光斑范圍內所有像元的數量

將所有曲線相結合連接成一條曲線f，作為直徑L的光斑所屬特征曲線(類似高光譜曲線)，以此作為分類基礎。曲線f可表示為

(4)

2.4 判別分析

2.4.1光斑

為了探究偽光斑大小與分類精度的關系，找出最適合于分類的偽光斑長度L，預先試驗1將光斑邊長L以10 m為步長，從10 m到80 m設置8個區梯度。在村莊、農田、森林3種地類中分別選擇樣點進行分類試驗，每一個樣點代表一個偽光斑。樣點類別的選擇和區分均通過在高分辨率CCD影像上進行目視判讀，記錄樣點坐標與類別信息。部分樣點位置如圖2(圖中藍色圓形代表村莊；紅色三角形代表農田；綠色正方形代表森林)所示。

圖2 部分樣點分布位置

由于作為分類的特征曲線具有較多的變量，在分類之前需要對分類曲線參數利用主成分分析(PCA)進行降維處理，以滿足包含80%以上信息量為標準提取前n個主成分因子作為分類的依據。利用判別分析的方法對樣點進行分類判別。將所有樣點隨機分為訓練樣本和驗證樣本，其中訓練樣本占總樣本數量30%，驗證樣本占70%。對不同偽光斑大小的分類特征曲線判別正確率進行統計分析，得光斑大小與分類精度關系曲線。試驗流程如圖3所示。所有統計分析和計算均通過Matlab程序實現。

圖3 分類試驗流程圖

2.4.2坡向

不同坡向下由于陰影的影響太陽輻射差異較大，這會影響到通過光譜圖像進行分類的精度，較容易將不同坡向的森林分為不同的類型。進行坡向分析試驗，利用CCD影像結合偽波形進行分類試驗。共選取150個樣點，農田、村莊和森林各50個樣點，森林中陰坡、陽坡各選取15個樣點，另外局部地勢相對平緩的森林中也挑選了20個樣點。設置訓練樣本為總樣本量的30%，即45個樣點，其中每種地類中各取15個樣點。為了能更好地驗證坡向對分類的影響，森林的訓練樣本均選自地勢平緩的樣點內。其余的70%樣點作為驗證樣本進行判別分析，可以獲得判別精度。

2.4.3林分密度

同樣，林分密度不一的情況下光譜的分異性也較大，稀疏林分和密集林分在光譜圖像上差別也較大，在一定程度上影響光譜圖像的分類。進行林分密度分析試驗。3種地類分別選取50個樣點，其中位于森林的樣點疏森林、密森林各布有15個，疏密相間的混合區域也設置了20個樣點。驗證樣本和訓練樣本的選擇與坡向分析試驗類似，森林的訓練樣本選自混合區域樣點。

2.5 基于人工神經網絡的土地利用分類

以前文所得的最優光斑大小，提取研究區范圍內區域的特征曲線，并且將特征曲線上50個變量看作50個波段，生成分辨率為L、波段數為50的研究區高光譜圖像。將研究區中地物分為森林、農田、村莊、水體4種類型，分別從每種類型中選出60個分布相對均勻的像元作為訓練樣本，采用人工神經網絡(Artificial neural network，ANN)方法對地物進行分類。人工神經網絡ANN是由大量的處理單元(神經元)廣泛互相連接而成的網絡，是一種具有高度非線性的分類分析方法，其特點和優越性主要體現在3方面：自學習功能、聯想存貯功能和高速尋找優化解的能力，特別適用于復雜非線性關系的遙感分類、回歸、評價和預測[33]。

根據訓練樣本，分別利用偽波形曲線與CCD影像的RGB 3波段采用ANN對整個研究區進行分類，在此基礎上結合50個波段30 m分辨率的分類圖像地表進行研究區地表分類，用目視判讀的方式在研究區內為每種地類選取特征明顯、具有代表性的區域作為真實值，對分類結果進行驗證。為了進一步與僅采用光學數據進行分類的常規方法進行比較，選擇GlobeLand30分類結果產品對偽波形結合CCD影像的分類結果進行對比驗證。地表分類及精度驗證均利用ENVI軟件實現。

3 結果與分析

3.1 偽波形提取效果

以L為30 m為例，圖4顯示了90個樣點的偽波形累積能量強度平均值。從圖中可以看出，不同地類偽波形曲線有較明顯差異，農田地勢較為平緩，在較低高程內的偽波形累積能量強度較大，因此大多數的偽波形累積能量集中于0～2 m區間內。村莊由于有較多地面的回波，其0～2 m區間內的偽波形累積能量強度平均值比森林高，而16～20 m的高程區間內有回波的主要是森林，其他兩類均趨于0。這歸根于森林植被冠層的覆蓋，使得其地面低層回波點的數量少于建筑物零散分布的村莊和植被相對低矮或幾乎為平地的農田，而高層的回波則高于另外兩類。

圖4 3種地類偽波形平均累積能量強度

3.2 特征曲線

圖5 3種地類特征曲線

3種地類的特征曲線如圖5所示，圖中曲線均為3種地類90個樣地的均值曲線。曲線范圍1～10代表點云平均強度信息fI，11～40分別代表R、G、B3個波段的頻率統計數據fR/G/B，41～50代表偽波形信息fw。可以看出，幾種地類平均值曲線在RGB 3個波段范圍內(11～40)均呈現單峰分布，峰值居中，地類之間峰值差異不太顯著，只有村莊與其他兩類在R波段范圍內峰值略有差異。而強度平均值曲線3種地類峰值點差異比較明顯，森林強度峰值明顯偏低，村莊其次，農田強度平均值曲線峰值偏高。各段特征曲線互相之間有一定的互補性。

3.3 光斑對分類精度的影響

以選定的90個樣地作為試驗樣本，分類精度與光斑邊長L關系如圖6所示。L在10～20 m范圍內，精度相對較低，這是由于光斑過小，點云密度偏低，統計特征不明顯。L達到30～40 m時精度基本達到最高，而L超過50 m時精度有所下降，原因可能是光斑過大時一個光斑樣地中包含了多種地物，形成了混合像元。進一步增大樣本量到130和150時，趨勢基本相同，說明規律較為穩定。因此本文中選取L為30 m進行進一步的分類試驗。

圖6 判別精度隨光斑邊長的變化曲線

3.4 坡向對分類精度的影響

單獨利用CCD影像RGB 3波段信息進行判別分析，總體正確率為60.95%，Kappa系數0.414，混淆矩陣如表1所示，其中森林陽坡有較多樣點被錯分為農田，村莊與農田之間混分現象也比較多；加入偽波形信息后，誤判率顯著降低(表2)。105個驗證樣本中只有4個判別錯誤，判別總體正確率96.19%，Kappa系數0.943，其中森林中陰坡陽坡及平地樣本均判別正確，對不同坡向的森林地類基本無明顯誤分為其他地類的情況。在山區，坡向會造成同為森林，但是由于光照陰影的影響,植被生長的差異和反射的光譜信息不同的現象，進而降低分類精度。加入偽波形的特征曲線分類對地形的坡向差異適應性較好。

表1 基于CCD影像的坡向試驗混淆矩陣

表2 加入偽波形信息后的坡向試驗混淆矩陣

3.5 林分密度對分類精度的影響

單獨利用CCD影像RGB 3波段信息進行判別分析，總體正確率為67.62%，Kappa系數0.512，混淆矩陣如表3所示，其中疏林有較多樣點被錯分為農田，村莊與農田之間同樣有混分現象，同為森林，由于密度不同，可能會造成較大的光譜差異；加入偽波形信息后總體精度為91.43%，Kappa系數0.871，分類結果混淆矩陣見表4。可以看出農田沒有出現誤判，農田與疏林之間判別精度有較大提高，但是森林仍有兩個樣點被誤判為村莊，村莊有7個樣點誤判為森林。出現誤判的森林均為疏林，由此可以看出森林過于稀疏依然會對分類產生一定的影響，但是試驗中疏林的誤判率能夠維持在較低的水平。而村莊的誤判率明顯更高，原因可能是村莊中部分區域有一定的綠化種植或原生植被覆蓋，加上村莊中建筑高度與喬木植被高度相當，使其樣地中光斑內同時具備類似于森林中的光譜反射條件和高程條件，很容易被誤分為森林。

表3 基于CCD影像的林分密度試驗混淆矩陣

表4 加入偽波形信息后的林分密度試驗混淆矩陣

3.6 圖像分類精度評價

1 m分辨率CCD影像、驗證樣地位置范圍及各分類結果見圖7。圖7a中驗證樣地綠色斑塊為森林、紅色斑塊為村莊、藍色斑塊為水體、粉色斑塊為農田。基于偽波形曲線和CCD影像的分類結果總體精度分別為91.93%和68.08%，Kappa系數分別為0.866 8和0.518 6，混淆矩陣見表5和表6。

圖7 CCD影像及驗證樣地與各分類結果對比

表5 基于偽波形分類結果的混淆矩陣

表6 基于CCD影像RGB三波段信息分類結果的混淆矩陣

從偽波形分類結果的混淆矩陣可以發現，基于偽波形分類結果中森林分類精度最高，幾乎沒有錯分誤差，只有少量的漏分誤差，主要與村莊有一定的混分。村莊區域則有較大的錯分誤差，其中多數被錯分為森林，這證明了利用偽波形高程信息進行分類不利于區分部分高度與森林相類似的村莊區域。漏分誤差最嚴重的地類為水體，約有51%水體被分成了農田和村莊，這是由于水體與農田和村莊部分區域均比較平坦，在偽波形中更容易混淆難以被分辨。這些都是高程信息不足以分辨的部分。

而基于CCD影像RGB 3波段信息分類結果混淆矩陣中村莊的錯分誤差相對較低，相反，村莊的漏分誤差極高，被混分為其他地類的較多；水體也同樣與偽波形呈現相反的精度，水體的漏分誤差較低，而錯分誤差極高，錯分來源主要是森林與村莊，這也是由于3.4節中提到的部分陽坡面的森林和村莊中建筑表面光譜反射率較高，與水體相類似所引起的混分。

偽波形與CCD影像單獨分類結果的驗證表明，兩類數據源各有優勢且互相之間對于劣勢具有較高的互補性，有必要將其結合用于提升分類精度。

圖7d為結合CCD與ALS提取的50個波段30 m分辨率分類圖像地表分類結果，圖7e為GlobeLand30 分類產品在研究區范圍內的分類結果。基于驗證樣地驗證結果表明，CCD與ALS數據融合得到的50個波段分類圖像的分類結果總體精度為95.22%，Kappa系數0.919 2，混淆矩陣及各種地類的生產者精度和用戶精度見表7。對應的研究區GlobeLand30產品總體分類精度為79.56%，Kappa系數0.661 8，混淆矩陣見表8，這與其產品描述中提到的總體分類精度83.51%，Kappa系數0.78相近，檢驗結果可靠。

表7 基于50個波段分類圖像分類結果混淆矩陣

從結果中可以看出，集合CCD與ALS數據的50個波段分類圖像分類結果總體精度和Kappa系數(95.22%，0.919 2)均高于偽波形單獨分類(91.93%，0.866 8)，并且相較于CCD單獨分類(68.08%，0.518 6)有較大幅度提升，各種地類的用戶精度與生產者精度也都與兩者的結果基本持平或有顯著提高。同時，數據融合分類的結果總體精度也優于GlobeLand30產品的總體分類精度。其中農田的分類用戶精度(91.71%)相對于GlobeLand30產品(58.11%)有大幅度提升，森林(97.59%與73.89%)和村莊(90.25%與72.31%)的生產者精度也有較大幅度提升。

表8 基于GlobeLand30分類結果的混淆矩陣

4 結論

(1)ALS偽波形對于地面高程差別較大的地表差異顯著，利用其與光譜信息相結合能夠實現信息量的互補，有助于提高分類及地物判別精度。

(2)不同光斑大小的偽波形對分類有一定的影響。光斑邊長過小，可能存在信息量不足，分類精度較低；而光斑邊長過大，則光斑內可能包含多種地物，分類精度相應也會有所降低。基于每平米點密度約0.39的ALS數據所構造的偽波形及特征曲線，當光斑邊長在30 m以內時，隨著光斑邊長的增加，分類精度有明顯的提升，光斑邊長達到30～40 m時，精度達到峰值，此后隨著光斑邊長的增加，分類精度不再上升，且在一定范圍內(50～60 m)有所降低，之后在峰值附近趨于穩定。

(3)相對于單一光譜數據，利用偽波形進行地物判別能夠更好地適應復雜地形及地表覆被的影響，對于不同坡向、不同疏密度的森林，加入偽波形時判別準確率較高。但是，對于綠色植被較多或者與其他地類有重疊的村莊地區，也出現了一些將村莊誤判為其他地類的情況。因此，偽波形結合光譜數據對于村莊地區中與其他地類混合度較大的分類效果有所欠缺。

(4)通過偽波形結合強度信息和光譜RGB 3波段信息進行土地利用分類，分類總體精度達到95.22%，Kappa系數0.919 2，較同一地區的純光譜分類GlobeLand30產品(總體精度79.56%，Kappa系數0.661 8)精度有明顯的提高。