基于机载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)精度有明显的提高。