改进型PSI算法及其在高空间分辨率遥感影像分割中的应用

孙小丹

(福州职业技术学院，福州 350108)

1 引 言

在高空间分辨率遥感影像(如：IKONOS、QuickBird)中，各类地物的形状、纹理、空间位置等特征均能得到清晰体现，这不仅为地物信息的判读提供了丰富的特征，同时也弥补了该类影像在光谱分辨率方面的不足。人工地物为城市的关键地类，是反映城市生态环境状况的一个重要因素。结合高空间分辨率遥感影像，实现人工地物信息的自动识别和提取仍是一大研究热点[1-3]。

空间分析技术可依据像元间的空间上下文关系，实现空间特征(如：纹理、形状)的提取或影像分割/分类。由于该技术能充分挖掘高空间分辨率遥感影像的优势，因此被广泛用于该类影像的分析与处理。目前，空间分析技术的应用主要体现在：①灰度共生矩阵(Gray Level Co-occurrence Matrix，GLCM)，灰度共生矩阵通过描述移动窗口内像元灰度值空间分布状况，从而获得一系列纹理特征值[4-6]；②小波分析技术(Wavelet Transform，WT)，该技术能提取不同空间尺度下地物的纹理、几何形状等特征，为实现多尺度影像分析提供重要依据[7-9]；③高斯马尔可夫随机场(Gaussian Markov Random Field，GMRF)从概率统计分析的角度，分析像元或对象之间的空间上下文关系，从而实现影像分割/分类[10-12]；④形态学理论可在多个空间尺度下提取地物形态结构特征[13-14]。这些技术虽然能从不同的角度描述影像中包含的各种空间特征信息，但是计算过程均比较繁琐，不利于进一步推广使用。2003年，由Shackelford和Davis提出了长度-宽度提取算法(length-width Extraction Algorithm，LWEA)，该算法依据像元之间的空间上下文关系，从长度和宽度等两个方向，描述局部区域的形状特征，因其操作步骤简单得以推广[15]。2006年，黄昕对长度-宽度提取算法做了进一步拓展，提出了像元形状指数PSI(Pixel Shape Index，PSI)算法[16]。相比LWE算法，PSI算法能从多个方向描述以某一像元为中心的邻域内的形状特征，同时计算得到的中心像元PSI指数与相位无关，具有较高的鲁棒性。由于在PSI派生波段下，人工地物的外形能得到清晰刻画，因此在影像分割时，联合该波段，能有效提高人工地物的分割质量。然而，在分割结果中，部分包含人工地物信息的图斑的边缘与人工地物的真实边缘不能完全吻合(即“边缘效应”现象)，这一现象降低了影像的分割精度。本文针对PSI算法存在的弊端，对其加以改进，并提出了一种改进型PSI算法。最后，通过影像分割对比实验，验证了新算法的可行性和优越性。

2 原理和方法

2.1 像元形状指数PSI及其算法[16]

像元形状指数PSI即是在以某一像元S为中心的邻域内，从D(D=8、12、16、20)个方向，对相邻像元与中心像元的上下文关系进行评估，其中像元间的上下文关系用像元间光谱特征同质性的高低来衡量。根据中心像元与不同方向上相邻像元上下文关系的大小，生成D条长短不一的方向线，通过该组方向线来刻画该邻域的形状和轮廓特征；接着，将D条方向线的长度通过简单的代数运算(如：求和等)得到邻域中心像元S的形状指数。例如：计算下图中建筑物内像元S的形状指数PSIS时，在以该像元为中心的圆形邻域内，沿着8个方向评估光谱特征的同质性，得到8条长短不一的方向线，通过这8条方向线来刻画像元S邻域范围内的形状特征，如图1所示。

图1 关于中心像元S的8条方向线示意图

求解像元S的PSIS算法包括方向线生成、方向线长度估算、计算PSI值等三个步骤(图2)。具体为：

(1)

其中，|·|为取绝对值操作。若PHi小于或等于阈值T1，则延长第i条方向线，反之则停止该方向线的延伸，直到方向线上像元的总数DLPNi大于阈值T2为止。

②方向线长度估算。第i方向线长度di定义为中心像元S和该方向线末端像元之间在水平、垂直等两个方向距离的最大值

di=max{|xs-xend|,|ys-yend|}

(2)

其中，xs、ys分别表示中心像元S在X轴、Y轴的坐标值；xend、yend分别表示方向线末端像元在X轴、Y轴的坐标值；max{·}为取最大值操作；

③计算PSI值。即将D条方向线的长度通过简单的代数运算，如求和运算(式(3))，得到中心像元S的PSIS值。最后，将该值映射到影像的光谱特征空间中。

(3)

2.2 PSI派生波段在影像分割中的应用

根据PSI算法，计算影像中每一个像元的像元形状指数PSI，即获得PSI派生波段数据，联合该派生波段，有助于提高影像分割中各类地物(尤其是外形规则的人工地物)的分割精度。本文联合QuickBird多光谱影像数据(包括蓝、绿、红、近红外等4个波段)，采用Definiens软件提供的区域生长分割技术，对小块的实验区做分割对比实验。分割前，先执行PSI算法，以获得实验区影像PSI派生波段数据，其中阈值T1=130、T2=50，方位角总数D=8；接着，执行两次影像分割，第一次基于原始多光谱波段数据，第二次基于原始多光谱波段和PSI派生波段数据，分割尺度均为70；最后，用蓝色线条勾勒出分割生成的图斑边缘，并叠加在由近红外、红、绿波段构成的假彩色影像上，结果如图2所示。对比两次分割结果，可以看出：相比基于原始多光谱波段数据的分割结果(图2(a))，在联合PSI派生波段数据的分割结果中(图2b)，建筑物的分割精度明显提高了，所生成的图斑边缘与建筑物边缘基本吻合，图斑中包含的建筑物信息比较完整，这将有助于提高后续建筑物信息自动识别和提取的精度。然而，在第二次分割结果中存在明显的“边缘效应”，即部分斑块的边缘与建筑物的真实边缘不能完全吻合(图2(b))，存在着偏差，即使调整阈值T1、T2、方位角总数D的数值大小，分割效果也得不到明显改进。在实验区内，以随机采样的方式，分别选取建筑物、建筑物边缘、植被等三类像元各20个，提取这三类像元的PSI波段值，并以散点图形式加以显示(图3)。从三类像元PSI值的分布状况可以看出：建筑物边缘处像元的PSI值接近植被像元的PSI值，而远离建筑物内部像元的PSI值，这导致了影像分割时建筑物边缘处的像元与植被被分在了一个图斑中，形成了“边缘效应”现象。

图2 分割结果对比图

为了分析产生“边缘效应”现象的原因，选取实验区中出现“边缘效应”的像元作为中心像元S，沿着水平方向线，提取6个相邻像元a1～a6，分别计算S与相邻像元a1～a6在Blue、Green、Red、Nir等4个波段光谱特征的差异和PH值。用折线图的方式显示计算结果，发现：PH值的变化趋势与3个可见光光谱特征差异的变化趋势较为相似，但与Nir光谱特征差异的变化趋势存在较大区别。因此PH值不能客观地反映各个波段光谱特征差异的变化趋势，而该值决定了方向线长度，这导致了方向线的长度无法客观地反映方向线上各像元在每一波段光谱特征同质性的高低，使得中心像元S的PSIs值不能客观地描述以像元S为中心邻域内的形状特征。PSI算法的弊端直接影响了上述实验区内建筑物边缘处像元PSI值的合理性，于是产生了“边缘效应”现象。

图3 关于建筑物/边缘/植被等PSI值的散点图

图4 中心像元S与水平方向线上相邻像元4个波段光谱征差异绝对值

2.3 改进型PSI算法

本文针对PSI算法存在的弊端，对其做了改进，并提出改进型PSI算法。新算法包括：多层次方向线生成、多层次方向线长度的估算、计算PSI值等三个步骤。除了第三个步骤与原算法一致外，前两个步骤均做了改进，具体如下所述。

2.3.1 多层次方向线生成

(4)

图5 在4个波段数据层上方向线扩展的示意图

2.3.2 多层次方向线长度的估算

在进行中心像元S第i方向的方向线长度估算时，先按原方法(式(2))估算各个波段数据层的方向线长度。接着，为了体现邻域内第i方向线上不同波段数据层光谱特征同质性存在的差异，将该方向上每个波段数据层方向线的长度做加权和，从而获得中心像元S第i方向的方向线长度。

(5)

为了初步验证改进型PSI波段的优越性，在分割参数一致的前提下，重新对上述小块实验区执行影像分割，结果如图6所示。分割前，先利用新算法获得改进型PSI派生波段数据，由于方向线的扩展是各波段独立进行，因此将阈值T1调整为35，T2和D保持不变。对比联合PSI派生波段的分割结果(图2(b))，可以看出：分割精度得到进一步提高，“边缘效应”现象明显减少，包含建筑物信息的图斑边缘与建筑物边缘基本吻合。

图6 联合多光谱波段和改进型PSI派生波段的分割结果

在实验区中，在出现边缘效应现象的区域，以随机采样的方式分别选取建筑物、建筑物边缘、植被等三类像元各20个，提取这三类像元的改进型PSI波段值，并以散点图形式加以显示(图7)。从三类像元的改进型PSI波段值的分布状况可以看出，在改进型PSI波段中，建筑物边缘处像元的PSI值不仅与建筑物像元的PSI值比较接近，还与植被像元的PSI值之间存在较明显的分界，数值的分布状况得以改善。因此，联合改进型PSI波段执行影像分割，“边缘效应”现象能得到明显减少，这充分说明了，在估算像元的PSI值时，考虑其邻域内方向线上各像元在不同波段光谱特征同质性存在的差异，能有效提高像元PSI值的合理性。

图7 关于建筑物/建筑物边缘/植被等改进型PSI值的散点图

3 对比实验及分析

3.1 实验数据及结果

为了避免上述实验中存在的偶然性，本文分别联合两种高空间分辨率遥感影像，选取两个较大面积的实验区，来验证改进型PSI算法的可行性和优越性。实验1区为IKONOS影像数据，大小907×338像素(图8(a)～图8(b))；实验2区为QuickBird影像数据，大小221×97像素(图9(a)～图9(b))，这两个区域均包含了城市中关键的地物信息。实验中，先分别采用PSI算法和改进型PSI算法，得到PSI、改进型PSI等两个派生波段数据；接着，在分割参数一致的前提下，分别联合这两个派生波段和原始多光谱波段数据，执行两次影像分割；最后，对比分割结果中外形规则的人工地物分割质量，对PSI、改进型PSI等两个派生波段的性能做分析和评价。

3.2 目视对比

首先，采用目视对比的方式，分析对比实验区的两种分割结果。对于实验1区来说，计算PSI派生波段时，阈值T1=110，T2=10，方位角总数D=8；计算改进型PSI派生波段时，阈值T1=20，T2和方位角总数D保持不变；影像分割时，尺度为65，各个波段的分割权值均为1，其余的分割参数均为默认值，两次影像分割的结果如图8(a)～图8(b)所示。从图可以看出：联合PSI派生波段，执行影像分割，外形规则的人工地物分割精度(图8(a))不高。例如：①方框I内的包含建筑物信息的图斑的边缘与建筑物的边缘不吻合，存在一些“边缘效应”现象；②方框II内的建筑物信息和相邻道路信息被划分到一个图斑，误分割现象较为严重。而在联合改进型PSI派生波段的分割结果中，包含人工地物的图斑边缘与建筑物、道路信息的实际边缘基本吻合，尤其在第I、II方框内，人工地物信息的分割质量得到明显提高(图8(b))。

图8 实验1区分割结果对比

对于实验2区来说，计算PSI派生波段时，阈值T1=135，T2=50，方位角总数D=8；计算改进型PSI派生波段时，阈值T1=37，T2和方位角总数D保持不变；影像分割时，分割尺度为70，各个波段的分割权值均为1，其余的分割参数均为默认值，两次影像分割的结果如图9(a)～图9(b)所示。从图可以看出：联合PSI派生波段，执行影像分割，其结果也不理想(图9(a))。主要有：①方框I内包含道路信息图斑的边缘与道路的真实边缘不吻合，存在“边缘效应”现象；②方框II、III内的道路与绿化带的信息被划分在一个图斑，存在误分割现象；③方框IV内包含建筑物信息图斑的边缘与建筑物的真实边缘不吻合，存在“边缘效应”现象。而在联合改进型PSI派生波段的分割结果中，上述方框内存在的误分割及“边缘效应”现象均得到明显地减少，分割质量得到明显改善(图9(b))。

图9 实验2区分割结果对比

3.3 分割精度验证

接着，对实验区的两种分割结果做精度验证。验证方法为：分别在两个实验区内建筑物、道路和其他地物(即植被、水体、裸土等)信息等三类信息中。以分区随机采样的方式选取250个验证像元，分别以IKONOS、QuickBird的全色波段影像为参照标准，进行人机交互式的精度验证，结果如表1、表2所示，精度验证数据表明：相比联合PSI波段的分割结果，联合改进型PSI波段的分割精度有了较大程度的提高，总精度提高到90%以上，Kappa系数能达到0.9以上。

4 结束语

本文针对PSI算法存在的弊端，对该算法的关键步骤做了相应地改进，并提出了改进型PSI算法。新算法的改进之处主要有：①在方向线生成阶段，充分考虑了不同波段光谱特征之间同质性存在的差异，每个波段数据层的方向线生成分开独立进行，从而提高了方向线生成的合理性；②提出了每一条方向线长度为各波段数据层的方向线长度的加权和，以进一步体现不同波段光谱特征之间同质性存在的差异，从而提高像元(尤其是边缘处的像元)PSI值的合理性。最后，通过实验证明：相比原算法，通过新算法获得的改进型PSI派生波段，能有效提高外形规则的人工地物的分割精度，“边缘效应”和误分割现象得到明显减少。然而，改进型PSI算法也存在一些不足：①阈值T1尚不能依据邻域内各波段光谱特征同质性的不同而动态地变化；②阈值T2尚不能根据方向线上各像元光谱特征同质性的高低而机动地做相应地调整。因此，在下一阶段的研究工作中，将从以上两点不足出发，进一步对算法加以改进。

表1 实验1区的精度验证数据

表2 实验2区的精度验证数据

参考文献：

[1] JENSEN J R.Introductory digital image processing：A remote sensing perspective[M].Beijing：China Machine Press，2007：155-161.

[2] ABRAHAM L，SASIKUMAR M.Automatic building extraction from satellite images using artificial neural networks[J].Procedia Engineering，2012，(50)：893-903.

[3] AWRANGIEB M，ZHANG C S，CLIVE S F.Automatic extraction of building roofs using LiDAR data and multispectral imagery[J].ISPRS Journal of Photogrammetry and Remote Sensing，2013，(83)：1-18.

[4] HEROLD M，LIU X H，CLARKE K C.Spatial metrics and image texture for mapping urban land use[J].Photogrammetric Engineering & Remote Sensing，2003，69(9)：991-1001.

[5] SU H，WANG Y P，XIAO J，et al.Improving MODIS sea ice detectability using gray level co-occurrence matrix texture analysis method：A case study in the bohai sea[J].ISPRS Journal of Photogrammetry and Remote Sensing，2013，(85)：13-20.

[6] 何珏，赵鹏，李浩.基于纹理的林区影像匹配窗口设置方法探讨[J].遥感信息，2013，28(4)：85-89.

[7] LU C S，CHUNG P C，CHEN C F.Unsupervised texture segmentation via wavelet transform[J].Pattern Recognition，1997，30(5)：729-742.

[8] XUE Z J，MING D，SONG W，et al.Infrared gait recognition based on wavelet transform and support vector machine[J].Pattern Recognition，2010，43(8)：2904-2910.

[9] GOH S S，GOODMAN T N T，LEE S L.Singular integrals，scale-space and wavelet transforms[J].Journal of Approximation Theory，2013，(176)：68-93.

[10] PANJWANI D K，HEALEY G.Markov random field models for unsupervised segmentation of textured color images[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1995，17(10)：939-954.

[11] 赵银娣，张良培，李平湘.广义马尔可夫随机场及其在多光谱影像分类中的应用[J].遥感学报，2006，10(1)：123-129.

[12] YE X F，ZHANG Z H，LIU P X，et al.A sonar image segmentation based on GMRF and level-set models[J].Ocean Engineering，2010，37(10)：891-901.

[13] PLAZA A，MARTINEZ P，PEREZ R，et al.A new approach to mixed pixel classification of hyper spectral imagery based on extended morphological profiles[J].Pattern Recognition，2004，37(6)：1097-1116.

[14] LOU S，JIANG X Q，SCOTT P J.Correlating motif analysis and morphological filters for surface texture analysis[J].Measurement，2013，46(2)：993-1001.

[15] SHACKELFORD A K，DAVIS C H.A hierarchical fuzzy classification approach for high-resolution multispectral data over urban areas[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41(9)：1920-1932.

[16] ZHANG L P，HUANG X，HUANG B，et al.A pixel shape index coupled with spectral information for classification of high spatial resolution remotely sensed imagery[J].IEEE Transactions on Geoscience and Remote Sensing，2006，44(10)：2950-2961.