|
|
|
发布时间: 2016-06-25 |
遥感图像处理 |
|
|
|
赵辽英1, 范明阳1, 厉小润2, 陈辰2 |
收稿日期: 2015-07-09; 修回日期: 2016-01-08
基金项目: 浙江省自然科学基金项目(LZ14F030004);国家自然科学基金项目(61571170)
第一作者简介: 赵辽英(1970-),女,教授,2004年于浙江大学获电气工程专业博士学位,主要研究方向为模式识别、遥感图像分析。 E-mail:zhaoly@hdu.edu.cn
中图法分类号: TP751
文献标识码: A
文章编号: 1006-8961(2016)06-823-11
|
摘要
目的 高光谱图像混合像元的普遍存在使得传统的分类技术难以准确确定地物空间分布,亚像元定位技术是解决该问题的有效手段。针对连通区域存在孤立点或孤立两点等特例时,通过链码长度求周长最小无法保证最优结果及优化过程计算量大的问题,提出了一种改进的高光谱图像亚像元定位方法。 方法 以光谱解混结合二进制粒子群优化构建算法框架,根据光谱解混结果近似估计每个像元对应的亚像元组成,通过分析连通区域存在特例时基于链码长度求周长最小无法保证结果最优的原因,提出修改孤立区域的周长并考虑连通区域个数构造代价函数,最后利用二进制粒子群优化实现亚像元定位。为了减少算法的时间复杂度,根据地物空间分布特点,采用局部分析代替全局分析,提出了新的迭代优化策略。 结果 相比直接基于链码长度求周长最小的优化结果,基于改进的目标函数优化后,大部分区域边界更明显,并且没有孤立1点和孤立两点的区域,识别率可以提高2%以上,Kappa系数增加0.05以上,新的优化策略可以使算法运算时间减少近一半。 结论 实验结果表明,本文方法能有效提高亚像元定位精度,同时降低时间复杂度。因为高光谱图像中均匀混合区域不同地物的分布空间相关性不强,因此本文方法适用于非均匀混合的高光谱图像的亚像元定位。
关键词
二进制粒子群优化, 高光谱图像, 亚像元定位, 空间相关性, 光谱解混
Zhao Liaoying1, Fan Mingyang1, Li Xiaorun2, Chen Chen2 Abstract
Objective Traditional classification technologies cannot easily or accurately determine the spatial distribution of ground features for hyperspectral images because mixed pixels are widespread throughout the image. Sub-pixel mapping technology is an effective tool to solve this problem. The existing sub-pixel mapping methods that are based on linear optimization encounter two issues in their practical implementation: their inexact objective functions and their excessive computation. Method This paper proposes a new sub-pixel mapping method to solve the aforementioned problems. The algorithm framework is constructed by combining spectral unmixing with binary particle swarm optimization. The numbers of sub-pixels for each pixel are estimated according to the results of spectral unmixing. The regional perimeter is modified by analyzing the influence on the perimeter and region number as induced by some special cases, such as isolated point or regions that include only two points. The cost function is formulated by considering the regional perimeter and number of connected regions. To reduce the running time of the algorithm, global analysis is replaced with local analysis according to the feature space distribution characteristics, and a new iterative optimization strategy is proposed. Result Compared with directly minimizing the region circumference based on the image chain code, the modified object function emphasizes the boundary of most regions and does not yield any isolated points or regions that include only two points. The method also improves the recognition rate by more than 2% and the Kappa coefficient by more than 0.05. Moreover, the new iterative optimization strategy nearly halves the CPU time. Conclusion The experimental results show that the proposed algorithm can improve the mapping accuracy and that the proposed optimization strategies can accelerate the mapping. Given the weak spatial correlation in areas where the end members are uniformly mixed, the proposed algorithm is suitable for hyperspectral images without uniformly mixed areas.
Key words
binary particle swarm optimization, hyperspectral imagery, sub-pixel mapping, space correlation, hyperspectral unmixing
0 引 言
高光谱图像包含了从可见光到近红外几百个连续的狭窄波段光谱信息，克服了宽波段遥感探测的局限，可以对物质表面信息进行完整的光谱描述，被广泛应用于多种领域，成为对地观测重要的信息源之一。但由于成像原理与制造技术等因素的限制，高光谱图像的空间分辨率相对较低，混合像元普遍存在于图像中。对于土地覆盖制图、海岸线提取、变化检测和景观指数估计等应用来说，混合像元内地物的空间细节信息极其重要，如果按照传统的硬分类方法将混合像元判定为任一类都是不正确的。亚像元定位正是解决该问题的有效手段。
亚像元定位的概念是由Atkinson教授[1]于1997年提出的，经过近20年的研究，已形成了基于神经网络的方法[2-4]、基于地统计学的方法[5-6]、基于马尔可夫随机场的方法[7-8]和线性优化的方法[9-17]等几类典型的方法。其中线性优化方法近几年研究比较多，这类方法的思想是首先通过混合像元分解或模糊分类技术确定各像元中地物的比例，定义一个描述空间相关性的数学模型，构造目标函数，将亚像元定位转化为一个线性优化问题。比利时学者 Verhoeye等人[9]提出以邻域像元与中心混合像元内各亚像元之间的空间相关性最大为目标，然后用单纯形等线性优化技术求解。Mertens等人2006年提出亚像元—像元空间引力模型(SPSAM)[10]，假设同类地物之间相互吸引，有着明确的物理意义，对空间相关性理论进行了有效的实现。针对SPSAM忽略了像元内部亚像元相互之间的相关性的缺点，王群明等人[11]提出一种基于修正的亚像元—像元空间引力模型(MSPSAM)的亚像元定位方法和同时考虑像元间和像元内相关性的混合空间引力模型。钟燕飞等人先后采用免疫克隆选择算法[12]和差分进化算法[13]，对文献[9]提出的数学模型进行智能化的搜索求解。吴柯等人[14] 提出基于进化Agent实现亚像元定位方法。上述方法的目标函数计算相对比较复杂。为了更直观地反映空间相关性，同样基于同类地物之间相互吸引的原理，2013年，Villa等人[15]提出以亚像元定位后各类区域周长最小为目标，分别用模拟退火和Atkinson[16]提出的像元交换技术(PSA)进行迭代优化。2014年，Erturk等人[17]利用二进制粒子群算法对文献[15]的目标函数进行优化，该算法并行性较强。
各种基于线性优化的亚像元定位方法的区别在于目标函数和优化方法的不同，其中以各类区域周长最小为目标来刻画空间相关性的方法具有直观并且计算简单的优点，但当采用链码长度求周长时，该目标函数没有考虑孤立点和孤立两点等特例，有时不能准确度量空间相关性。另一个问题是如果每个像元每次迭代时，计算全图周长构造代价函数[17]，其计算量较大，因此如何根据地物空间分布，采用局部分析来代替全局分析以减少运算量有待进一步研究。针对上述问题，根据地物空间分布，对目标函数进行调整，并提出了新的局部优化策略，实验结果表明，本文方法即提高了算法速度，又提高了定位精度。
1 算法理论基础
1.1 高光谱混合像元分解
线性光谱混合模型[18]具有简单实用的特点，目前被广泛地应用于混合像元分解问题的研究中。假设r∈RL×1为L个波段高光谱图像中单个像元的光谱向量，M∈RL×P为图像中P种地物对应的端元光谱矩阵，s∈RP×1为该像元中所有地物对应的组分比(丰度)，ε∈RL×1为误差项，则线性光谱混合模型表示为[18]
| $ r{\rm{ = }}Ms{\rm{ + }}\varepsilon $ | (1) |
| $ {M_{lp}} \ge 0,1 \ge {s_p} \ge 0,\sum\limits_{p = 1}^P {{s_p}} = 1 $ | (2) |
式(2)表示非负性约束和全加性约束，即端元光谱和丰度分布是非负的，同时每一个像素内所含端元的丰度和为1。
1.2 二进制粒子群算法
在二进制粒子群算法(BPSO)中，每个优化问题的解称为“粒子”，每个粒子有两个特征量：速度和位置。BPSO初始化为D维目标搜索空间中的一群随机粒子，然后通过粒子速度和位置的迭代找到最优解。设m个粒子构成一个群体，其中第i个粒子的位置可表示为D维的矢量
| $ v_i^{k + 1} = {c_1} \otimes ({P_i} \oplus z_i^k) + {c_2} \otimes ({P_\varepsilon } \oplus z_i^k) $ | (3) |
| $ z_i^{k + 1} = z_i^k \oplus v_i^{k + 1} $ | (4) |
式中，
1.3 基于链码的图像周长计算
常用的链码方法有4向链码和8向链码，8向链码更能准确反映原图像边界[20]。如图 1(a)为8向链码编号。从图像边界上任意选取一个像素点作为起始点，顺时针沿着边界按照方向编号规则，给每一对边界像素间的线段编号，依次将方向编号连接即可得到图像边界的链码表示。如图 1(b)所示为图像边界的8向链码。
图像中某连通区域的周长可以根据区域边界链码计算[20]，即
| $ p{\rm{ = }}h{\rm{ + }}\sqrt {\rm{2}} s $ | (5) |
式中，h和s分别为边界上编号为偶数和奇数的编码个数，即边界上非对角连线数和对角连线数。
需要说明的是周长计算是对封闭边界进行的，图 2给出了一些特殊区域的链码周长[21]。
2 算法描述
2.1 线性优化的亚像元定位算法框架
算法框架如图 3所示，首先确定高光谱图像中地物种类数目，通过光谱解混确定每个像元包含的地物个数及其面积比例。每个纯像元对应的亚像元全部归类为相应的纯地物类别，各混合像元的亚像元定位过程为：根据丰度确定混合像元中每种地物的亚像元个数，根据空间相关性理论构造代价函数，经过线性优化方法确定各亚像元最终类别。各种线性优化方法都适用于该算法框架，本文选用二进制粒子群优化算法。
2.2 粒子表示及更新处理
设每个像元对应一个由r2个亚像元构成的亚像元窗。设第i个像元包含d种地物，其丰度按
| $ \begin{array}{l} {n_{ik}} = round({s_{ik}} \times {r^2}),k = 1,2,\cdots ,d - 1\\ {n_{id}} = {r^2} - \sum\limits_{k = 1}^{d - 1} {{n_{ik}}} \end{array} $ | (6) |
式中，round(·)表示四舍五入。
因为只要前d-1种地物在亚像元窗中的分布确定，则第d种地物的分布就随之确定了，所以第i个像元对应的粒子位置和速度分别只需用一个(d-1)×r2的二进制矩阵表示。每个粒子位置矩阵zi满足两个约束条件：1)矩阵的第k行有且仅有nk个值为1；2)矩阵的每列最多只有一个值为1。
粒子的初始化和更新过程中，zi都必须满足上述两个条件，因此粒子更新后需要进一步判断和处理，具体过程为：
1)如果zi第k行1的个数不等于nk，则随机修改一位以保持1的个数不变。
2)如果zi第k列为1的行数超过1，则找到这些行中1个数最少的行(设为第m行)，置zi(k,m)为0，并修改一个全0列的第m行值为1。
2.3 目标函数
依据地物分布的空间相关性，通常假设一个像元中所含的某种地物应该与其相邻像元中所含的同种地物空间相邻，基于这种假设，定义代价函数[15-17]
| $ C = \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^{{C_i}} {{P_{ij}}} } $ | (7) |
式中,M表示类别数，Ci表示第i类的连通区域个数，Pij为第i类第j个连通区域的周长。
图像连通区域周长计算通常有区域和背景缝隙的长度之和、链码长度和边界点数之和3种方法[20]。文献[15-17]都以式(7)为代价函数实现亚像元定位，并用链码长度求Pij。为了便于描述，将文献[15-17]的目标函数称为目标函数1。
下面分析目标函数1的问题。定义图像空间分辨率增强比为每个像元对应的亚像元区域大小。以图像空间分辨率增强比等于3×3为例分析，可能存在以下两种特殊情况：
1)可能存在只有一个亚像元的连通区域(孤立单点)。以图 4(a)所示的场景1为例，该场景中含有两种地物，设每3×3区域对应一个像元，则两种地物的丰度如图 4(b)所示，图 4(c)为满足该丰度的含有孤立单点的特殊分布。图 4(a)(c)的8连通区域数分别为3和9，每个连通区域用链码长度计算周长[20]相加得到图 4(a)(c)的周长分别为
2)可能存在只有两个亚像元的连通区域(孤立两点)。以图 4(d)所示的场景2为例，该场景中含有两种地物，设每3×3区域对应一个像元，则两种地物的丰度如图 4(e)所示，图 4(f)为满足该丰度的含有孤立单点的特殊分布。图 4(d)(f)的8连通区域数分别为3和5，周长分别为
分析目标函数1在上述两种情况下失效的原因是：通常一个孤立区域与其他区域合并后会增加总的链码长度，从而使目标函数1的值增大。因为用链码长度求孤立单点的周长为0，孤立两点区域周长为2，所以两个孤立单点总周长为0，区域合并后，总周长为2；1个孤立单点与1个孤立两点区域合并前总周长为2，区域合并后总周长为4或
根据上述分析，基于链码长度求区域周长Pij，提出如下优化的目标函数
| $ C = \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^{{C_i}} {{{\hat P}_{ij}} + kM} } $ | (8) |
| $ {{\hat P}_{ij}} = \left\{ {\begin{array}{*{20}{l}} {{P_{ij}} + \beta \;\;{P_{ij}} = 0或{P_{ij}} = 2}\\ {{P_{ij}}\;\;\;\;\;\;\;其他} \end{array}} \right. $ | (9) |
式中，系数k>0，增量β≥0。
一方面孤立单点(Pij=0)或孤立两点区域(Pij=2)的周长加β≥0，但有时这样仍无法保证最优结果对应目标函数最小，如图 4(c)，8个孤立区域，每个孤立区域的周长加1后，总周长仍然等于图 4(a)的周长，因此，在目标函数中考虑区域个数，因为小的孤立区域的存在，会使得整个图中区域个数增加。
2.4 目标优化策略
亚像元定位的目标优化是全局迭代结合局部迭代的过程。局部迭代是逐像元迭代，所有像元遍历一遍为一次全局迭代。下面主要分析局部迭代的具体策略。
为便于描述，将文献[17]的局部迭代优化策略称为策略1,即逐像元迭代，并且每个像元每次迭代时，计算全图周长构造代价函数。由于每个像元迭代优化时都要计算整个图的周长，策略1的时间复杂度比较高。为了减少时间复杂度，设计新的优化策略2。
策略2：逐像元迭代，但每个像元每次迭代时，只计算其对应的亚像元窗及周围一圈构成的局部区域周长构造代价函数。
图 5给出了策略2对应的逐像元迭代优化区域示意图，其中图 5(b)对应图 5(a)点划线区域的3×3像素分辨率增强后的亚像元定位图，图 5(a)中圆圈标注的像元为图 5(b)中正处理的像元。逐像元迭代时，因为只改变了该像元对应的亚像元分布，其对应的亚像元窗口周围一圈的类别不变(图 5(b)中虚线所示)，因此可以只计算虚线区域的周长构造代价函数。
2.5 算法时间复杂度
因为迭代过程是逐像元进行的，每次计算目标函数时只有少数亚像元位置发生了变化，所以，单个像元迭代过程中，一个粒子迭代一次的时间近似相等。设图像中混合像元个数为N，迭代次数为m，粒子数为p，以单个像元一个粒子迭代一次的时间为单位时间，整个算法迭代过程的时间复杂度近似为O(Nmp)。两种策略相比，用策略1时的时间复杂度更高。
3 实验结果
3.1 仿真数据实验
首先，假设丰度100%准确的情况下，比较不同目标函数和不同优化策略下，亚像元定位的结果。为了构造准确的丰度图，以美国印第安纳州的AVIRIS Indian pines高光谱数据真实地物分布为基础，该地物真实中共给出了17类不同的地物(包括背景)，将其中的某些类合并后得到9类，从USGS光谱库中选9种地物光谱按照真实地物分布生成仿真的高光谱图像，按3×3比例均值滤波，构造相应的混合像元丰度，得到如图 6所示的高光谱图像(第30波段)。其中图 6(a)为[1∶144,1∶144]区域按3×3缩小得到的仿真图，大小为48×48像素，图 6(b)是截取[1∶60,70∶144]区域缩小后大小为20×25像素的子图。
通过文献[22]方法确定地物种类个数，顶点成分分析法(VCA)[18]提取光谱曲线，全约束最小二乘法(FCLS)[23]求得丰度，再亚像元定位。实验环境为Intel(R) CoreTMi5-3210M CPU。
本部分共5个实验，除了第4个实验，其他实验的空间分辨率增强比均为3×3。主要以最后定位的识别率、Kappa系数及运行时间(不包括光谱解混时间和各地物在亚像元分布图中亚像元个数确定时间)为客观指标评价算法性能。为了避免粒子群算法的初值敏感性对结果的影响，实验1和实验2中的每个实验全部用了同样的初始粒子群(粒子个数为50)，初始识别率为93.83%。所有结果都是运行5次的平均值。为了便于描述，将目标函数1对应的算法命名为CBPSO，优化目标函数对应的算法称为改进的CBPSO，命名为MCBPSO。
实验1不同迭代优化策略对结果的影响。分别以目标函数1和优化的目标函数进行亚像元定位，采用2种不同的迭代优化策略，最大迭代次数为20次。表 1给出了各种目标函数和优化策略组合的亚像元定位结果。从表 1中可以看出用策略2可以明显减少计算时间，并且识别率及Kappa系数都得到了提高。图 7给出了两种算法采用策略2时的亚像元定位结果图。从图 7可以看出，CBPSO的结果有很多孤立的点，MCBPSO的结果中大部分类别边界比较明显，并且没有孤立1点和孤立两点的区域。图 8给出了目标函数值随迭代次数的变化情况。
表 1
不同目标函数和迭代优化策略组合的结果比较
Table 1
Results of different objective functions and iterative optimization strategies
| 算法 | 迭代策略 | 识别率/% | Kappa系数 | 时间/s |
| CBPSO1 | 策略1 | 94.27 | 0.848 2 | 411 |
| CBPSO2 | 策略2 | 94.91 | 0.865 0 | 216 |
| MCBPSO1 | 策略1 | 96.78 | 0.915 2 | 340 |
| MCBPSO2 | 策略2 | 97.88 | 0.944 5 | 218 |
| 注：MCBPSO算法中系数k=2，增量β=1。 | ||||
设置迭代停止条件为：迭代次数大于最大迭代次数或前后两次迭代的目标函数差值小于10-6。从图 8可以看出，所有算法都迭代4次就收敛，其实后面所有的实验都是迭代4次就收敛了，这可能是由BPSO算法决定的。CBPSO的目标函数值小于MCBPSO的目标函数值，这是因为MCBPSO将链码长为0或2的连通区域周长各加了1(β=1)，并且目标函数值加了连通区域数的2倍(k=2)。
由于策略1计算全图周长，策略2只计算局部子图的周长，所以策略1的时间复杂度相对高一些。从表 1可以看出，策略1的时间近似为策略2时间的两倍。
相比目标函数1，改进的目标函数计算每个连通区域周长后还要多一次判断，对于满足条件的区域还要多一次加法运算，所以MCBPSO算法的运行时间通常略大于CBPSO的运行时间，但由于粒子迭代过程中两个约束条件是否满足是随机的，相应的处理时间具有不确定性，所以有时MCBPSO反而比CBPSO快，表 1中MCBPSO和CBPSO时间关系与此分析是一致的。
实验2优化的目标函数中参数对结果的影响。为了验证优化的目标函数中参数β和k不同取值对结果的影响，对优化的目标函数采用策略2进行了几组参数取值的实验比较，最大迭代次数为20次。
因为用链码长度求孤立单点的周长为0，孤立两点区域周长为2，所以两个孤立单点总周长为0，区域合并后，总周长为2；1个孤立单点与1个孤立两点区域合并前总周长为2，区域合并后总周长为4或
根据上述特例分析，参数β的取值为
表 2
MCBPSO取不同参数的结果比较
Table 2
Results of MCBPSO with different parameters
| 参数取值 | 识别率/% | Kappa系数 | 时间/s |
| β=1,k=1 | 97.53 | 0.933 6 | 219 |
| β=2,k=1 | 97.45 | 0.931 9 | 221 |
| β=1,k=2 | 97.88 | 0.944 5 | 218 |
| β=2,k=2 | 97.56 | 0.936 4 | 218 |
实验3 粒子数的影响。 用优化的目标函数结合策略2，参数β=1，k=2，最大迭代次数为20次。表 3给出了粒子数分别为30、40、50、60和100的结果。
表 3
不同粒子数的结果比较
Table 3
Results of different number of the particles
| 算法 | 粒子数 | |||||
| 30 | 40 | 50 | 60 | 100 | ||
| CBPSO | 识别率/% | 94.69 | 94.79 | 94.91 | 94.88 | 94.82 |
| Kappa系数 | 0.861 9 | 0.862 2 | 0.865 0 | 0.863 9 | 0.861 2 | |
| 时间/s | 121 | 179 | 216 | 259 | 429 | |
| MCBPSO(β=1,k=2) | 识别率/% | 97.49 | 97.72 | 97.88 | 97.86 | 97.73 |
| Kappa系数 | 0.931 9 | 0.941 5 | 0.944 5 | 0.943 8 | 0.942 9 | |
| 时间/s | 135 | 177 | 218 | 259 | 430 | |
从表 3可以看出，随着粒子数的增加，运行时间越来越长，这与算法的时间复杂度分析一致。另外，识别率和Kappa系数随粒子数增加先增大后减小，对于优化的目标函数，粒子数为60时识别率和Kappa系数开始减小，目标函数1则时粒子数100时识别率和Kappa系数反而变小了。从定位精度和时间综合考虑，粒子数50比较合适。
实验4优化的目标函数的适用性。为了验证优化的目标函数对图像空间分辨率增强比大于3×3时的亚像元定位的适用性，截取[1∶60,70∶144]区域进行5×5滤波，然后按图像增强比5×5进行亚像元定位。分别用CBPSO和MCBPSO采用策略2迭代优化，粒子数为分别为50和100，最大迭代次数为20次。
因为5×5滤波后的图像中纯像元变少，类似于文献[17]，需要对丰度值做调整，将最大丰度值与其他丰度值之差的最小值大于0.9的像素作为纯像元。所有算法用相同的初始粒子群，初始识别率为90.35%。表 4给出了定位结果的评价指标值，图 9给出了粒子数为50时的定位结果图。从表 4可以看出：1)MCBPSO参数取值不同结果差别不大，并且所有结果明显优于CBPSO的结果；2)从识别率和Kappa系数指标来看，粒子数100时稍微优于粒子数50时，但前者运行时间是后者的3倍多，因此从时间和精度综合考虑，选粒子数50也是比较合适的。
表 4
空间分辨率增强比5×5的亚像元定位结果
Table 4
results of sub-pixel mapping when spatial resolution enhancement is 5×5
| 算法 | 粒子数50 | 粒子数100 | ||||||
| 识别率/% | Kappa系数 | 时间/s | 识别率/% | Kappa系数 | 时间/s | |||
| CBPSO | 90.93 | 0.775 5 | 120 | 91.02 | 0.780 5 | 375 | ||
| MCPBSO | β=1,k=1 | 92.57 | 0.819 1 | 135 | 92.86 | 0.826 5 | 498 | |
| β=1,k=2 | 92.68 | 0.822 1 | 139 | 93.17 | 0.833 0 | 382 | ||
| β=2,k=1 | 92.42 | 0.814 9 | 135 | 92.82 | 0.825 8 | 459 | ||
| β=2,k=2 | 92.29 | 0.810 6 | 136 | 92.54 | 0.817 3 | 451 | ||
类似于图 7，图 9的结果同样表明CBPSO的结果有很多孤立的点。图 9的结果说明相比CBSPO，MCBPSO对于图像增强比大于3×3的情况也能提高定位精度。但从图 9还可以看出，结果图中还有其他的孤立情况，因此，实际应用中，对于图像增强比大于3×3的情况再进行有针对性的分析和处理，也许能得到更好的结果。
实验5图像大小的影响。为了验证论文所提算法对影像的适用性，根据Indian pines整个图的真实地物分布设计仿真图，如图 6(a)所示。用优化的目标函数结合策略2，参数β=1，k=2，粒子数为50，最大迭代次数为20次。图 10为定位结果，表 5给出了评价指标。从图 10可以看出CBPSO的结果中各类之间有很多孤立的点，MCBPSO的结果中大部分类别间边界比较明显。表 5的结果也说明MCPBSO的定位结果更好，但MCPBSO的运行时间比CPBSO多了100多秒。比较表 5和表 1，可以看出图像大小对定位精度没什么影响。图像越大，混合像元数越多，经统计图 6(a)中的混合像元数近似为图 6(b)中混合像元数的5.5倍，比较表 5和表 1中的运行时间，考虑迭代过程中的两个约束条件的不确定性，可以认为运行时间与图像中混合像元数近似成正比。
3.2 真实图像实验
实验数据采用ROSIS高光谱遥感影像，采集地点是意大利北部帕维亚大学。整幅图大小为610×340像素，包含9类，波段数是103，几何分辨率为1.3 m。截取117×75大小的区域，该区域包含5类地物，其真实地物分布如图 11(a)所示，3×3均值滤波，图 11(b)为滤波后第37波段的灰度图，VD确定端元个数5，用VCA+FCLS得到丰度图。同文献[17]，对丰度值做调整，将最大丰度值与其他丰度值之差的最小值大于某个阈值的像素作为纯像元，因为金属区域比较明显，只统计金属区域的定位识别率。采用策略2迭代优化。
表 6给出了阈值为0.4和0.5时各种算法的识别率和运行时间比较。图 11(c)为随机定位结果。图 11(d)(e)分别为粒子数为阈值0.4，50个粒子时的CBPSO和MCBPSO的结果。
表 6
PaviaU子图亚像元定位结果比较
Table 6
PaviaU subgraph sub-pixel mapping results
| 阈值 | 评价指标 | 随机定位 | CBPSO | MCBPSO | |||||
| n=50 | n=100 | n=50 | n=100 | n=50 | n=100 | ||||
| 0.4 | 识别率/% | 94.28 | 94.87 | 95.81 | 95.89 | 97.79 | 98.13 | ||
| 时间/s | — | — | 329 | 832 | 427 | 887 | |||
| 0.5 | 识别率/% | 93.46 | 93.61 | 94.55 | 94.49 | 97.55 | 97.49 | ||
| 时间/s | — | — | 440 | 1 130 | 455 | 123 2 | |||
从表 6可以看出：1)所有情况下，MCBPSO的结果都明显比CBPSO的结果好； 2)阈值0.4的结果比阈值0.5的结果好，这与文献[17]的结果一致；3)从识别率看，总体上粒子数100时稍微好一点，但粒子数越大，时间代价越大。从图 11可以看出，CBPSO的结果中有较多的孤立区域，导致有些类边界不清，而MCBPSO的结果中孤立区域很少，说明MCBPSO的结果优于CBPSO的结果。
4 结 论
为了降低利用链码长度求区域周长最小为目标的亚像元定位方法的不确定性，提高定位精度并降低算法复杂度，提出了一种改进的亚像元定位新方法。考虑到连通区域为孤立点和孤立两点等特例，修改区域周长并结合连通区域个数构造代价函数，目标函数的优化采用全局迭代结合局部迭代的方法，提出了一种新的迭代优化策略。通过仿真和真实遥感图像实验，验证了本文方法相比原线性优化方法能获得视觉上更符合空间相关性分布的定位结果，Kappa系数、识别率和运行时间等指标的定量评价亦表明新方法的有效性。因而本文方法为提升基于线性优化的亚像元定位提供了一种新的可行思路。将来的研究工作为空间分辨率增强比例增大时的特例分析及处理，进一步提高亚像元定位的精度。
参考文献
-
[1] Atkinson P M. MappingSub-pixel Boundaries from Remotely Sensed Images[M]//Kemp Z. Innovations in GIS IV. London: Taylor and Francis, 1997: 166-180.
-
[2] Tatem A J, Lewis H G, Atkinson P M, et al. Super-resolution target identification from remotely sensed images using a Hopfield neural network[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001 ,39 (4) : 781 –796. [DOI:10.1109/36.917895]
-
[3] Shi W Z, Zhao Y L, Wang Q M. Sub-pixel mapping based on BP neural network with multiple shifted remote sensingimages[J]. Journal of Infrared and Millimeter Waves, 2014 ,33 (5) : 527 –532. [ 史文中, 赵元凌, 王群明. 多偏移遥感图像的BP神经网络亚像元定位[J]. 红外与毫米波学报, 2014 ,33 (5) : 527 –532.] [DOI:10.11972/j.issn.1001-9014.2014.05.012]
-
[4] Wu K, Niu R Q, Li P X, et al. Sub-pixel mapping of remote sensing images based on fuzzy ARTMAP neural network model[J]. Geomatics and Information Science of Wuhan University, 2009 ,34 (3) : 297 –300. [ 吴柯, 牛瑞卿, 李平湘, 等. 基于模糊ARTMAP神经网络模型的遥感影像亚像元定位[J]. 武汉大学学报: 信息科学版, 2009 ,34 (3) : 297 –300.] [DOI:10.13203/j.whugis2009.03.009]
-
[5] Boucher A, Kyriakidis P C. Super-resolution land cover mapping with indicator geostatistics[J]. Remote Sensing of Environment, 2006 ,104 (3) : 264 –282. [DOI:10.1016/j.rse.2006.04.020]
-
[6] Boucher A, Kyriakidis P C. Integrating fine scale information in super-resolution land-cover mapping[J]. Photogrammetric Engineering & Remote Sensing, 2007 ,73 (8) : 913 –921. [DOI:10.14358/PERS.73.8.913]
-
[7] Kasetkasem T, Arora M K, Varshney P K. Super-resolution land cover mapping using a Markov random field based approach[J]. Remote Sensing of Environment, 2005 ,96 (3-4) : 302 –314. [DOI:10.1016/j.rse.2005.02.006]
-
[8] Ardila J P, Tolpekin V A, Bijker W, et al. Markov-random-field-based super-resolution mapping for identification of urban trees in VHR images[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2011 ,66 (6) : 762 –775. [DOI:10.1016/j.isprsjprs.2011.08.002]
-
[9] Verhoeye J, De Wulf R. Land-cover mapping at sub-pixel scales using linear optimization techniques[J]. Remote Sensing of Environment, 2002 ,79 (1) : 96 –104. [DOI:10.1016/S0034-4257(01)00242-5]
-
[10] Mertens K C, de Basets B, Verbeke L P C, et al. A sub-pixel mapping algorithm based on sub-pixel/ pixel spatial attraction models[J]. International Journal of Remote Sensing, 2006 ,27 (15) : 3293 –3310. [DOI:10.1080/01431160500497127]
-
[11] Wang Q M. Research on sub-pixel mapping and its related techniques for remote sensin gimagery[D]. Harbin: Harbin Engineering University, 2012. [王群明. 遥感图像亚像元定位及相关技术研究[D]. 哈尔滨: 哈尔滨工程大学, 2012.]
-
[12] Zhong Y F, Zhang L P, Li P X, et al. A sub-pixel mapping algorithm based on artificial immune systems for remote sensing imagery[C]// IEEE International Geoscience and Remote Sensing Symposium. Cape Town: IEEE, 2009, 3: III-1007-III-1010. [DOI: 10.1109/IGARSS.2009.5417948]
-
[13] Zhong Y F, Zhang L P. Sub-pixel mapping algorithm based on adaptive differential evolution for remote sensing imagery[C]// IEEE International Geoscience and Remote Sensing Symposium. Vancouver, BC, Canada: IEEE, 2011: 1724-1727. [DOI: 10.1109/IGARSS.2011.6049568]
-
[14] Wu K, Li P X, Zhang L P. Sub-pixel toying of remote sensed images based on evolutionary agentalgorithm[J]. Journal of Remote Sensing, 2009 ,13 (1) : 60 –66. [ 吴柯, 李平湘, 张良培. 一种基于进化 Agent 的遥感影像亚像元定位方法[J]. 遥感学报, 2009 ,13 (1) : 60 –66.] [DOI:10.3321/j.issn:1007-4619.2009.01.008]
-
[15] Villa A, Chanussot J, Benediktsson J A, et al. Unsupervised methods for the classification of hyperspectral images with low spatial resolution[J]. Pattern Recognition, 2013 ,46 (6) : 1556 –1568. [DOI:10.1016/j.patcog.2012.10.030]
-
[16] Thornton M W, Atkinson P M, Holland D A. Sub-pixel mapping of rural land cover objects from fine spatial resolution satellite sensor imagery using super-resolution pixel-swapping[J]. International Journal of Remote Sensing, 2006 ,27 (3) : 473 –491. [DOI:10.1080/01431160500207088]
-
[17] Erturk A, Gullu M K, Cesmeci D, et al. Spatial resolution enhancement of hyperspectral images using unmixing and binary particle swarm optimization[J]. IEEE Geoscience and Remote Sensing Letter, 2014 ,11 (12) : 2100 –2104. [DOI:10.1109/LGRS.2014.2320135]
-
[18] Nascimento J M P, Bioucas Dias J M. Vertex component analysis: a fast algorithm to unmix hyperspectral data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005 ,43 (4) : 898 –910. [DOI:10.1109/TGRS.2005.844293]
-
[19] Afshinmanesh F, Marandi A, Rahimi-Kian A. A novel binary particle swarm optimization method using artificial immune system[C]//Proceedings of the International Conference on Computer as a Tool. Belgrade: IEEE, 2005, 1: 217-220. [DOI: 10.1109/EURCON.2005.1629899]
-
[20] Xie F Y. Digital Image Processing and Application[M]. Beijing: Publishing House of Electronics Industry, 2014 . [ 谢凤英. 数字图像处理及应用[M]. 北京: 电子工业出版社, 2014 .]
-
[21] Kulpa Z. More about areas and perimeters of quantized objects[J]. Computer Vision, Graphics, and Image Processing, 1983 ,22 (2) : 268 –276. [DOI:10.1016/0734-189X(83)90069-5]
-
[22] Chang C I, Du Q. Estimation of number of spectrally distinct signal sources inhyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004 ,42 (3) : 608 –619. [DOI:10.1109/TGRS.2003.819189]
-
[23] Heinz D C, Chang C I. Fully constrained least squares linear spectral mixture analysis method for material quantification inhyperspectral imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001 ,39 (3) : 529 –545. [DOI:10.1109/36.911111]
