1. 引言
随着对地观测技术的不断发展，高光谱遥感影像已经成为土地制图、目标识别等城市信息提取研究的重要数据源之一。高光谱数据虽然具有极高的光谱分辨率，但是光谱特征存在“同物异谱，异物同谱”的特点 [1] 。因此，提取有效的纹理、结构等空间特征对于高光谱影像解译也尤为重要。黄昕等 [2] 所提出的像元形状指数(Pixel Shape Index, PSI)是一种有效的纹理提取算法，它利用光谱同质性进行方向线扩展，进而计算得到每一个像元的PSI特征指数，补充光谱信息的不足，提高光学影像数据解译精度。但是，PSI算法主要针对单波段影像进行特征提取，对于具有数百个波段的高光谱数据，PSI算法难以发挥作用。
然而，对于城市区域复杂的空间格局和高程变化，光学影像单一数据源难以满足当前应用需求 [3] 。近年来，激光雷达(Light Detecting and Ranging, LiDAR)作为一种新型对地观测技术，得到了广泛研究和应用。不同于高光谱等光学影像，LiDAR数据通过大量具有空间三维信息的激光回波脚点，记录地面目标的高度、反射率等信息，可直接获取地物高程信息。因此，在城市遥感信息提取领域，点云辅助光学影像进行多源遥感数据分类相比单一数据源具有极大的优势。
高度特征是LiDAR点云数据中最显著的信息 [4] 。利用点云数据生成的数字表面模型(Digital Surface Model, DSM)或归一化数字表面模型(normalized Digital Surface Model, nDSM)等栅格数据，可以区分光谱和空间特征难以判别的地物，例如房屋和道路、树木和草地等，弥补光谱维和空间维的信息缺失。
LiDAR点云数据与高光谱遥感影像融合用于城市土地覆盖分类已经得到了大量研究，并可以获得比单一数据源更高的分类精度 [5] 。董彦芳等 [6] 利用机载LiDAR点云数据与高光谱数据融合后分类，取得了比单一数据源更好的分类精度。谷延峰等 [7] 提出一种多核学习框架，可以充分利用两种数据源获得更好的分类结果。这些研究分别提取高光谱遥感影像的空谱特征和LiDAR点云数据的高程特征，在分类过程中将这些特征进行融合。然而，这种方法忽略了多源数据之间的内在联系。因为遥感影像和LiDAR点云数据均为地面区域的信息载体。对于同一地物，多源数据从不同的角度进行分别表达。因此，如果将多源数据直接进行融合，可以减小地物内部的异质性，增大地物之间的差异性，便于后续特征提取操作。
本文主要改进了纹理特征和高程特征的提取方法，首先针对高光谱数据的特点，利用光谱角距离优化了像元形状指数计算方法。其次，基于这种改进后的算法，我们将点云数据与高光谱数据进行波段叠加，并利用改进后的PSI算法提取得到更加稳健的空间特征。
2. 特征提取方法
在点云辅助光学影像进行城市遥感解译的过程中，光谱特征、空间特征和高度特征是三类最常用的特征信息。其中，光谱维特征一般选择所输入的高光谱影像所有波段，或经过降维处理后所得到的光谱特征；间维特征包括纹理、结构等二维影像特征，主要是利用相关算法从光学影像中提取得到；高程维的特征则是从LiDAR点云数据中提取得到的地物高程信息。本节首先对像元形状指数进行了改进，然后利用改进后的算法提取纹理特征。
2.1. 像元形状指数
像元形状指数PSI是以影像中某一像元P为中心像元，当邻近像元Pn与中心像元Pc满足一定的光谱同质性时，对其邻域进行方向线扩展，然后获取该中心像元所有方向线的统计特征（均值、最大值、最小值等）作为该像元的形状指数。PSI具体计算过程可以表示如下：
Step.1方向线扩展：对于影像中某一像素
，以其为中央像元进行方向线扩展，方向线扩展过程如图1所示，图中黄色像元C表示中央像元，蓝色像元S对应中央像元的各条方向线。邻近像元
与中央像元之间的同质性测度可以表示为：
(1)
式中，
表示第k条方向线上邻域像元与中心像元之间的同质性测度，
和
分别表示中央像素和第k条方向线上邻近像元在第i个波段上的影像灰度值。i表示影像第i个波段，N表示影像波段总数。当
小于同质性阈值条件T1时，该条方向线扩展到下一像素。当
大于T1或该条方向线长度大于空间阈值T2时，结束该条方向线扩展。
Step.2方向线长度计算：当第i条方向线扩展结束时，其长度为：
(2)
式中，
表示第i条方向线长度，
和
分别表示该条方向线的两个端点像元。
Step.3统计PSI值：当中心像元的所有方向线计算完成时，利用全部方向线长度即可计算该像元的PSI值。其中，最为常用的计算方式是计算各条方向线长度的均值作为该像元处的PSI值。
(3)
式中，
表示影像
点的PSI值，D表示方向线总数，
表示第i条方向线长度。
相比灰度共生矩阵，PSI采用方向线扩展的形式计算局部纹理信息，准确探测对象形状结构信息的同时，避免了GLCM所带来的窗口效应。但是，PSI同质性测度的计算过程主要考虑了高空间分辨率遥感影像的空间信息。而对于高光谱影像，公式(1)所采用的计算方式忽略了其光谱维信息，难以获取准确的形状指数。因此，目前常用的做法是利用主成分变换(Principal Component Analysis, PCA)等降维方法对高光谱影像进行降维，并采用其第一主分量作为基影像计算PSI特征。然而，降维影像虽然保留了高光谱数据绝大部分的空间信息，却大量损失了光谱信息。因此，本文提出一种适应高维遥感数据像元形状特征提取的方法，采用光谱角距离(Spectral Angle Distance, SAD)作为高光谱影像同质性测度。
Figure 1. The direction lines of PSI algorithm
图1. PSI算法的方向线扩展
对于高光谱影像，其中央像元
，第i条方向线上邻近像元
，其中n表示高光谱影像波段总数。中央像元与其邻近像元的光谱角距离可以表示为：
(4)
表示反三角函数，SAD取值位于[−1,1]区间。中央像元和邻近像元同质性越高，SAD越接近于1。
光谱角距离可以充分利用高光谱图像的光谱维信息，从光谱曲线的相似程度来判定两个像元的同质性程度。并且，SAD测度对像元向量的相似程度十分敏感，而与向量长度无关。相比原始的PSI算法，利用SAD替换公式(1)中的街区距离，在高维遥感数据特征提取中可以获得更加稳健的PSI特征。
2.2. 点云数据特征
LiDAR点云数据高程特征提取的相关研究中，nDSM是最为常用的高度特征 [8] 。相比数字表DSM，nDSM特征消除了地形的影响，因此其灰度值表示对应地物的真实高度。nDSM的计算过程可以表示如下：
(5)
式中，为影像中第i个像素，n为影像像素总数。DSM和DTM分别为数字表面模型和数字地面模型，前者利用点云数据可以直接采样得到，而后者是经过点云滤波去除非地面点后采样得到的。
强度信息是点云数据中另一个重要特征，其对地物具有一定的判别能力。但由于LiDAR强度信息容易受到多种条件干扰，因此强度特征影像通常具有较大的噪声。在强度信息提取的基础上，Jinha Jung等人 [9] 利用伪波模型进行高光谱影像和点云数据融合分类。伪波模型是利用点云数据的高度和强度信息，构建多波段点云特征。
伪波特征的构建过程可以表述为：
Step.1：利用DTM数据对点云激光脚点的高程转换为地面高度，去除地形对伪波模型的影响；
Step.2：利用点云的空间三维坐标将激光脚点划分为dx × dy × dz的体元，其中dx和dy表示DTM的空间分辨率，dz表示高度划分的距离，每一个高度层对应伪波模型中的一个波段；
Step.3：对每一个体元中的所有激光脚点，计算点云强度信息的均值作为对应波段和像素的像素值。
最终所获得的点云伪波模型是多波段特征影像，不同波段对应不同高度层点云回波所采样生成的特征影像。因此，伪波模型中每一波段所对应的高程划分，对特征具有较大影响。高程划分太大，不同高度的地物难以区分，高程划分太小，同一高度的地物插值到不同波段，导致特征中包含大量零散点。图2是本实验所采用的一种三波段伪波模型。其中，地面层地物指高程小于1 m的地物，例如操场、植被等；中层地物指高程位于1 m至10 m之间的地物，主要包括居民区、树木等。高层地物指高程大于10 m的地物，主要包括商业区等高层建筑。
2.3. 点云–高光谱融合特征提取
多源遥感影像的特征融合方法通常包括两种类型：多特征融合和决策融合 [10] 。前者提取多源数据的光谱、空间、高程等信息，归一化后输入单一分类器进行数据融合；后者主要利用不同特征和分类器进行组合。相比决策融合的方法，多特征融合减少分类先验知识，是一种更为简单、实用的方法。
由于LiDAR特征影像与高光谱影像在地物分布上具有互补性，因此本文的融合策略是选取LiDAR特征与光谱数据进行叠加后，采用SAD测度改进的PSI算法提取空间特征，该算法技术路线示意图如图3所示。叠加了LiDAR特征后的光谱数据，其同类地物之间的光谱同质性得到增强，不同地物之间的光谱异质性得到增强。本实验中，采用波段叠加的方式进行多特征融合，对于光谱、空间和高程特征影像，将其作为特征影像的不同波段进行组合，即可得到融合特征。
3. 实验及分析
实验主要包括两个部分，第一部分用来比较两种改进的多波段PSI计算方法与原始PSI计算方法，在提取高光谱数据特征时的效果和分类精度；第二部分用来测试利用多波段PSI算法进行先融合后提取的效果，并与其他融合方法进行比较。
Figure 2. Pseudo-waveform false color image of Huston University (R: median height, G: ground, B: large height)
图2. Huston University点云数据伪波特征假彩色影像(红：中层地物，绿：地面层，蓝：高层地物)
Figure 3. The flow chart of the multi-sources PSI features based on SAD-PSI
图3. 基于SAD-PSI的多源PSI特征提取流程图
3.1. 实验数据
本实验选用的数据是美国休斯顿大学及其邻近区域的高光谱影像和LiDAR数据，主要覆盖商业区、居民地等城市地物，研究数据由“2013 IEEE GRSS Data Fusion Contest”组委会提供 [11] 。其中，高光谱影像包含144个光谱波段，波长范围为380 mm~1080 mm，空间分辨率为2.5 m。机载LiDAR数据包括原始LiDAR点云数据和DSM影像，其中DSM影像与高光谱数据具有相同空间分辨率。
由于实验数据只提供了DSM特征影像，因此首先利用DSM影像生成nDSM特征。在nDSM生成过程中，使用数字形态学开重建，选用直径为80的结构算子对DSM影像进行处理。开重建可以消除影像中小于结构算子大小的目标，因此经过开重建操作后DSM影像中大部分地物目标被剔除，仅保留了地形高程，即DTM影像。然后，对DSM和DTM影像取差值，即可生成nDSM特征影像。
由于原始LiDAR点云数据所包含的是离散三维空间坐标，因此本实验需要对LiDAR点云数据进行处理，提取点云强度特征和伪波模型。强度特征采用自然插值算法，将激光脚点采样为高光谱数据相同空间分辨率的栅格影像。伪波模型生成过程中，空间分辨率dx，dy设为2.5 m，高程距离dz的设置为：1）第一层：高程小于2 m的回波点，代表高程维中地面地物；第二层：高程位于2 m~15 m的回波点，代表高程维中具有中等高度的地物；3) 第三层：高程大于10 m的其他点，代表高程维中具有较大高度的地物。需要注意的是，如果每一层高程范围过小，则该层栅格化后仅存在少量离散点信息。
本实验采用支持向量机(Support Vector Machine, SVM)进行监督分类，影像中地物被分为15类，包含植被、建筑、道路等，并对植被、建筑类型做了具体划分。实验所采用的数据及样本分布如图4所示。
3.2. PSI测度实验
在PSI计算过程中，对于多波段影像，通常的做法是利用主成分变换等方法进行降维处理，并取第一主成分作为基影像，进行PSI纹理特征提取。然而，该方法在计算同质性测度 时，没有充分利用高光谱数据的光谱信息。本实验选取基于SAD测度的多波段PSI计算方法分别提取空间特征，并与基于第一主成分的原始PSI计算方法相比较。
对于多波段计算方法，首先基于高光谱全波段数据计算纹理特征。经过多组对比实验，在SAD测度实验中，同质性测度T1设为0.1，空间阈值T2均设为10，方向线总数N设为10。所提取的特征如图所示。提取空间特征后，将光谱特征与空间特征叠加输入SVM分类器，使用5折叠交叉验证的方法自动寻优。在试验中，光谱特征分别选取了高光谱全波段信息和PCA前4波段与空间特征进行波段叠加。最终分类精度和分类结果图分别如表1和图5所示，其中HSI表示原始高光谱影像，PCA表示主成分变换的前4波段光谱数据。图5中由上至下分别表示原始高光谱数据、PCA降维数据、PCA + PSI数据和PCA + SAD-PSI数据的分类结果。
从表1中可以看出，PSI空间特征对于分类结果具有明显的改善，而SAD-PSI实验方案则取得了最佳实验精度。对于植被、土壤、水体等自然地物，高光谱数据的光谱特征具有较好的判别性，因此仅利用光谱信息可以达到95%以上的分类精度。但是，对于城市区域的地物目标，例如居民地、商业用地、道路、停车场等地物，其材质主要为混凝土或沥青等，具有类似的光谱特征。因此，在仅利用光谱信息进行分类时，其结果精度小于80%，难以进行准确划分。
相比光谱特征分类结果，在引入空间特征，尤其是SAD-PSI特征后，上述城市目标类别的分类精度均有一定程度的提高。其中商业用地的分类精度从40%提高至65%，停车场从70%提高至85%，因此形状结构特征对于城市人工地物具有较高的判别性。实验结果中，SAD-PSI相比PSI总体精度提高2%，各类别精度也有一定程度提升，利用全部光谱波段比仅利用第一主分量所提取的PSI特征更加稳健。
Figure 4. The false color composition of hyperspectral data, nDSM, distribution of test samples and training samples of Huston University Data, successively
图4. 由上至下分别是Huston University高光谱数据、nDSM数据、测试样本集分布和训练样本集分布
Table 1. The information of ground references and the classification accuracies of hyperspectral dataset of Huston University
表1. Huston University数据样本分布及高光谱数据分类结果精度(%)
Figure 5. The heperspectral image classification result maps of Huston University
图5. Huston University高光谱数据分类结果图
从分类结果图中也可以看出，原始光谱数据和PCA降维后数据的分类结果中不仅存在大量椒盐错分，而且由于缺少必要的空间信息，导致不同地物之间的错误也非常严重。而利用降维数据提取PSI形状结构特征后与PCA特征进行叠加，其右侧云层覆盖区域分类结果较差。这是因为采用第一主分类进行形状结构特征提取的情况下，所提取的空间特征仅与第一主分类有关。而本组数据中经过主分类变换后云层集中在第一分类中，因此严重影响了空间特征的准确性。采用SAD-PSI计算形状结构特征，其分类效果相比单一波段特征提取有显著改善。
从实验精度和效果来看，采用SAD-PSI更适用于高维数据形状结构特征提取。然而，分类效果图也表明SAD-PSI仍然无法解决单一数据源所带来的信息缺失问题。因此，后续实验将采用SAD-PSI提取融合数据的形状结构特征以改善实验结果。
3.3. 多源数据融合实验
本实验影像解译的一个难点在于影像右方存在大面积云层覆盖区域。使用光学影像单一数据源时，由于光学影像受到云层干扰，该区域难以获取正确的分类结果。LiDAR属于主动遥感技术，不受天气及云层影响，但缺失光谱信息。因此多源数据融合可以弥补单一数据源信息不足，达到更好的解译结果。
实验首先分别提取了点云数据的nDSM特征、强度特征(Intensity)和3波段伪波特征模型(PW)，以及高光谱数据的前4维主分量数据。对于PSI特征提取的对比实验，我们采取常用的方法，利用光谱数据第一主分量作为基影像，提取PSI特征作为空间信息。而在SAD-PSI特征提取的实验中，我们采取“先融合，后提取”的方法，首先将原始高光谱数据与点云特征进行融合，然后利用融合数据提取SAD-PSI特征作为空间信息。最后，将空间信息与光谱信息(前4维主分量数据)、点云特征进行波段叠加后输入SVM分类器，最终分类精度和分类结果图分别如表2和图6所示。图中由上至下分别表示原始高光谱数据、PCA与nDSM融合，PCA与强度数据分类，PCA、nDSM和强度数据融合分类以及PCA和伪波模型的分类结果。
从点云数据分类的结果可以看出，使用nDSM数据和强度数据进行融合分类的情况下，nDSM数据精度明显高于强度数据。此外，由于伪波模型将不同高度的地物划分到不同波段，对于商业用地、道路、铁路具有更好的判别性，因此使用伪波模型提取SAD-PSI特征可以获得更好的分类结果。
从表2的结果可以看出，在有点云特征辅助的情况下，数据分类精度得到大幅提高。尤其在城
市目标，例如建筑、道路、停车场等地物类别上，分类精度提高了10%~40%左右。以商业用地为例，商业用地与居民地具有极其相似的光谱特征和空间特征，而在高度特征上具有较大差异。因此，本实验利用SAD-PSI提取融合数据空间特征的方法，取得了最好的精度。
从实验结果图来看，点云数据和高光谱影像融合的分类结果明显优于高光谱单一数据源分类效果。尤其是在高光谱数据云层覆盖区域，点云数据可以补充高光谱影像的信息缺失，在一定程度上提取出云层覆盖区域的地物信息。
实验采用“先融合，后提取”的方法提取得到融合数据的SAD-PSI特征，从两组PSI和SAD-PSI的对比实验可以看出，SAD-PSI特征在云层覆盖区域可以得到更加完整的分类效果，避免了单一波段中
Table 2. The classification result (%) of the Huston university multi-source dataset
表2. Huston University多源数据融合的分类结果精度(%)
Figure 6. The classification result maps by fusing hyperspectral and LiDAR dataset of Huston University
图6. Huston University点云数据与高光谱影像融合分类结果图
异常数据的干扰。而PSI特征因为仅使用单波段数据进行空间特征提取，收到云层干扰较为严重。此外，SAD-PSI特征可以显著抑制分类结果中的椒盐噪声，改善分类精度。如图6中黑色圆圈所标识区域所示，引入SAD-PSI特征的结果，其地物目标完整性大大增强。
4. 结论
本文针对点云辅助高光谱影像融合分类所存在的一些问题，提出了一种同时利用高光谱影像和点云数据进行形状结构特征提取的方法。首先，本文利用光谱角距离SAD改进像元形状指数PSI的同质性测度，使PSI算法适用于高光谱影像特征提取。改进后的SAD-PSI算法利用数据所有波段进行同质性计算，可以融合光谱信息和空间信息计算得到更加稳健的PSI特征。其次，本文将高光谱影像和点云特征影像进行融合，然后利用SAD-PSI算法提取融合数据的形状结构特征。为了突出融合特征的类间差异性，采用伪波特征替换nDSM、强度信息作为点云信息与高光谱影像进行波段叠加。最后，将高光谱影像、点云数据和SAD-PSI特征联合，归一化后输入SVM分类器进行分类。实验结果表明，本文所采用的这种“先融合、后提取”的特征提取方法，相比传统方法取得了更好的效果。
致谢
感谢美国休斯顿大学的高光谱图像分析组和美国国家科学基金会资助的机载激光测绘中心提供本研究中使用的遥感数据。