本发明涉及一种结合相位与光谱的高分遥感影像分割方法，属于遥感影像处理技术领域。
背景技术：
随着遥感影像空间分辨率的不断提高，同类地物内部的光谱差异也更加显著，使传统像素级分析方法已越来越难以满足实际应用的需求。对象级图像分析object-basedimageanalysis(obia)技术则其中geographicalobject基本分析单元，为高分辨率遥感影像的信息提取提供了新的思路。开展obia首先需要有效的图像分割，图像分割也是后续处理中实现准确的特征提取及目标识别的基本前提与重要保证。
与中、低分辨率遥感影像相比，高分辨遥感影像具有更加丰富的光谱、纹理、几何等空间细节特征。特别是在城市场景中，杂乱的图像背景、多样的地物种类以及结构复杂的人造地物等使高分辨率遥感影像呈现出结构化、层次性特征，因此需要从多层次、多角度加以分析和理解。这些因素均导致了采用传统的图像分割方法难以取得理想的效果。
面对这些困难与挑战，学者们已经从不同角度提出了一些有效的应对策略。例如，zhang等通过引入多尺度分析策工具对传统分水岭分割方法加以改进，有助于深入剖析对象的空间细节特征，从而更好地定位对象的边缘[1]；hu等首先通过可信度指标进行边缘检测并结合商业软件ecognition来确定初始对象，再通过构建区域邻接图及分析边缘完整性来进一步精炼分割结果，该方法能够有效分割不同面积、不同尺寸的对象(改写)[2]；raffaele等提出了一种基于标记点自适应控制的多尺度分水岭分割方法，通过在最精细尺度中对边缘信息、光谱以及形态学进行融合，能够有效识别不同尺度下代表性地物的轮廓特征[3]。尽管如此，这些方法均过于依赖图像的梯度信息，对于影像中对比度较低的部分响应较弱，容易造成轮廓信息的丢失。另一方面，基于局部同质性信息来构建边缘检测模型不仅能够减少对梯度信息的依赖，而且对弱边缘非常敏感，但同时分割过程受到的噪声和伪边缘的干扰也更加显著。
参考文献：
[1]zhangb,hebb.multi-scalesegmentationofhigh-resolutionremotesensingimagebasedonimprovedwatershedtransformation[j].journalofgeo-informationscience,2014,16(1):142-150.
[2]huy,chenj,pand,etal.edge-guidedimageobjectdetectioninmultiscalesegmentationforhigh-resolutionremotelysensedimagery[j].ieeetransactionsongeoscienceandremotesensing,2016,54(8):4702-4711.
[3]gaetanor,masig,poggig,etal.marker-controlledwatershed-basedsegmentationofmultiresolutionremotesensingimages[j].ieeetransactionsongeoscienceandremotesensing,2015,53(6):2987-3004.
技术实现要素：
发明目的：针对现有技术中存在的问题与不足，为实现高精度、高可靠的高分辨率遥感影像分割，本发明利用传统彩色纹理分割算法jseg中j-value指标所具有的强大局部同质性信息检测能力与多尺度特性，与相位一致性信息相结合来进行优势互补，构建了一种参数自适应的边缘响应模型，进而提出了一种联合相位一致性及局部同质性的高分辨率遥感影像分割方法。通过不同传感器类型的高分辨遥感影像的分割实验表明，所提出的方法能够获得低依赖于梯度信息的边缘响应，准确定位对象弱边缘并减少噪声及伪边缘的影响。
技术方案：一种结合相位与光谱的高分遥感影像分割方法，主要包括四个步骤：(1)根据图像边缘点之间的相位信息具有较高的一致性，实现了基于相位一致性信息的边缘检测；(2)计算局部均匀性指数j值，提取局部均匀性信息；(3)将相位一致性与局部均匀性相结合，构建边缘响应模型，同时遍历所有可能的参数组合，寻找目标函数的最小值，得到自适应的最优参数组合和响应模型；(4)在优化模型的基础上，对高分辨率遥感影像图像进行多尺度分割，并通过结构相似性(ssim)对相邻区域进行融合。
边缘检测
首先基于相位一致性提取图像的边缘信息，具体过程如下：
设f(x)为变量为x的一维离散信号，其傅里叶展开式可表示为：
其中，an，ωn，
由于在边缘特征点处各次谐波分量的值相等，因此公式(2)可以转化为:
在此基础上，采用loggabor小波的多尺度特征对相位一致性模型进行改进；定义变量pc来表示任意一个像素z的相位一致性大小，令(a,b)为像素z的坐标值，则改进后的二维相位一致性模型为：
因此，可采用δφn表示尺度n下所有方向相位差异的和，可进一步采用如下公式表示：
e(a,b)和o(a,b)分别为loggaborfilter的奇对称和偶对称小波，n为滤波器尺度，eno(a,b)、ono(a,b)分别为loggaborfilter的奇对称和偶对称小波，
δφno(a,b)为瞬时相位差异，t0为估计噪声，en(a,b)2为尺度n下奇对称小波的平方、on(a,b)为尺度n下偶对称小波，φn(a,b)为尺度n的相位值，
j-value指标的提取
引入梯度信息建立参数自适应边缘响应模型，参数自适应边缘响应模型的建立过程如下：
step1：影像量化；在luv空间中采用fcm动态聚类对原始多光谱影像进行颜色量化，获得单一波段的包含q个灰度级的量化影像；
step2：基于局部同质性的边缘检测；在量化影像中，设任意一个像素为z，依据公式(1)～(4)计算每个像素z的相位一致性值pcz(a,b)；遍历所有像素，获得初步边缘检测结果；
step3：计算多尺度j-value集合；在量化影像中，构建以像素z为中心，尺寸为d×d像素的窗口p，其中d为尺度参数，去除窗口p中的角点；
令pq为窗口中属于灰度级q的所有像素的集合，采用像素z的坐标z(a,b)作为替换该像素的灰度值，令mq为pq内所有像素的平均值，则窗口z中属于同一灰度级像素的方差的和可表示为：
pq中所有像素的均值为m，则总体方差可表示为：
则尺度参数为d时像素z对应的j-value为：
jd(a,b)＝(st-sw)/sw(8)
边缘响应模型构建与优化
基于所提取的相位一致性及局部同质性信息，对于任意一个像素z，构建边缘响应模型为：
ier(a,b)＝pc(a,b)×j(a,b)(9)
pc(a,b)、j(a,b)分别表示z的相位一致性指标值和局部相似性值；该模型主要涉及三个参数，即计算局部同质性模型时loggabor滤波器的尺度参数lgscale和方向参数lgori，以及计算j-value时的尺度参数d。为提高算法的自动化程度，提出基于尺度间互信息最小化的目标函数优化策略，来自适应的确定最优模型参数的组合；实现过程为：
首先确定用于后续图像分割的总尺度数目k，参数lgscale，lgori以及d的取值范围。在此基础上，记两个连续尺度间的互信息为ami，构建目标函数：
amitot＝-[ami(1，2)+ami(2，3)+...+ami(k-1，k)](10)
ami(k-1，k)表示相邻尺度间的互信息；通过迭代操作在给定范围内遍历所有可能的参数组合，当amitot取得最小值时，即认为此时获得的多尺度边缘图像序列seqopt＝{scale1,scale2...scalek}能够最为充分的反映影像中具有代表性的边缘细节特征，则此参数组合所构建边缘响应模型为最优模型；同时，确定了最优的lgscale,lgori和d。
多尺度区域分割及合并
(1)多尺度分割
多尺度分割在所提取的seqopt中进行；令tk＝μk+0.2σk以确定提取区域的阈值，其中，μk和σk为尺度scalek下的灰度均值与标准差，k＝1,2...k；在最粗糙尺度scalek中，检测所有灰度值小于阈值tk的像素，进而采用4连通方法获得种子区域，其余像素(联通区域像素，作为区域增长的起点)按照该像素的j-value及该像素与距离最近种子区域的距离从小到大的顺序逐个并入相邻的种子区域，从而获得scalek下的对象集合rk＝{rk1,rk2...rkn}；在下一个较精细尺度scalek-1中，分别以尺度scalek中提取的每个对象为基本分割单元，计算tk-1从而构建种子区域并进行区域增长。以此类推，直到最精细尺度scale1计算完毕，可获得初步分割结果r1；
(2)基于ssim的区域合并
在r1中，计算任意两个相邻对象ra和rb的结构相似性s(a,b)：
其中，a、b分别ra和rb对应的特征向量，μa，μb，σa，σb，
在seqopt的所有尺度及原始影像的所有波段中分别计算s(a,b)，并取均值为
step1：依据r1生成区域邻接图rag；
step2：对任意一个对象ra，根据rag确定所有与ra相邻的对象；
step3：计算ra与相邻的一个对象rb，计算所有尺度以及原始影像所有波段下ra和rb的欧式距离均值
step4：若
step5：更新rag，重复step2至step5，遍历所有对象及其相邻对象，直到没有能够合并的区域为止，从而获得最终分割结果。
附图说明
图1是本发明方法流程图；
图2是当d＝9时的窗口示例；
图3是0.6m航空遥感dom影像；
图4是2mworldview-2影像；
图5是试验一不同策略的边缘检测结果：(a)相位一致性模型，(b)j-value检测，(c)边缘响应响度模型；
图6是试验二不同策略的边缘检测结果：(a)相位一致性模型，(b)j-value检测，(c)边缘响应响度模型；
图7是试验一不同方法的分割结果：(a)本发明方法，(b)方法2，(c)方法3，(d)方法4，(e)方法5；
图8是试验二不同方法的分割结果：(a)本发明方法，(b)方法2，(c)方法3，(d)方法4，(e)方法5；
图9是j-value定量评价结果：(a)实验一；(b)实验二；
图10是lmiinor定量评价结果：(a)实验一；(b)实验二。
具体实施方式
下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
如图1所示，结合相位与光谱的高分遥感影像分割方法，主要包括四个步骤：(1)根据图像边缘点之间的相位信息具有较高的一致性，实现了基于相位一致性信息的边缘检测；(2)计算局部均匀性指数j值，提取局部均匀性信息；(3)将相位一致性与局部均匀性相结合，构建边缘响应模型，同时遍历所有可能的参数组合，寻找目标函数的最小值，得到自适应的最优参数组合和响应模型；(4)在优化模型的基础上，对高分辨率遥感影像图像进行多尺度分割，并通过结构相似性(ssim)对相邻区域进行融合。
1、基于局部同质性模型的边缘检测
morrone的研究表明，相位信息在图像边缘特征点处具有高度一致性，并提出了相位一致性边缘检测模型。为减少分割过程对梯度信息的依赖，提高对弱边缘的检测能力，本发明首先基于相位一致性提取图像的边缘信息。具体过程如下：
设f(x)为一维离散信号，其傅里叶展开式可表示为：
其中，an，ωn，
由于在边缘特征点处各次谐波分量的值相等，因此公式(2)可以转化为:
在此基础上，采用loggabor小波的多尺度特征对相位一致性模型进行改进。定义变量pc来表示任意一个像素z的相位一致性大小，令(a,b)为像素z的坐标值，则改进后的二维相位一致性模型为：
因此，可采用δφn表示尺度n下所有方向相位差异的和，可进一步采用如下公式表示：
e(a,b)和o(a,b)分别为loggaborfilter的奇对称和偶对称小波，n为滤波器尺度，eno(a,b)、ono(a,b)分别为loggaborfilter的奇对称和偶对称小波，
δφno(a,b)为瞬时相位差异，t0为估计噪声，en(a,b)2为尺度n下奇对称小波的平方、on(a,b)为尺度n下偶对称小波，φn(a,b)为尺度n的相位值，
2、j-value指标的提取
尽管依据局部同质性信息有助于提高对弱边缘的检测能力，同时也必然使图像分割对噪声及伪边缘等高频干扰因素更加敏感。因此，引入梯度信息来进一步提高边缘检测的可靠性，建立参数自适应边缘响应模型，参数自适应边缘响应模型的建立过程如下：
step1：影像量化。在luv空间中采用fcm动态聚类对原始多光谱(或者pan-sharpenedimage)影像进行颜色量化，获得单一波段的包含q个灰度级的量化影像。
step2：基于局部同质性的边缘检测。在量化影像中，设任意一个像素为z∈z，z表示窗口，依据公式(1)～(4)计算每个像素z的相位一致性值pcz(a,b)。遍历所有像素，获得初步边缘检测结果。
step3：计算多尺度j-value集合。在量化影像中，构建以像素z为中心，尺寸为d×d像素的窗口p。其中d为尺度参数(sp)，为了保证各个方向高频信息的一致性，窗口中的角点被去除。图2为当d＝9时的窗口示例。
令pq为窗口中属于灰度级q的所有像素的集合，采用像素z的坐标z(a,b)作为替换该像素的灰度值，令mq为pq内所有像素的平均值，则窗口z中属于同一灰度级像素的方差的和可表示为：
pq中所有像素的均值为m，则总体方差可表示为：
则尺度参数为d时像素z对应的j-value为：
jd(a,b)＝(st-sw)/sw(8)
由此可以看出，j-value是基于局部区域的灰度统计特征而来的，实质上反映了原始影像所包含的梯度信息，因此对孤立点等噪声及伪边缘具有较高的鲁棒性，可以有效弥补相位一致性信息存在的局限。
3、边缘响应模型构建与优化
基于所提取的相位一致性及局部同质性信息，对于任意一个像素z，构建边缘响应模型为：
ier(a,b)＝pc(a,b)×j(a,b)(9)
pc(a,b)、j(a,b)分别表示z的相位一致性指标值和局部相似性值；该模型主要涉及三个参数，即计算局部同质性模型时loggabor滤波器的尺度参数lgscale和方向参数lgori，以及计算j-value时的尺度参数d。为提高算法的自动化程度，提出基于尺度间互信息最小化的目标函数优化策略，来自适应的确定最优模型参数的组合。实现过程为：
首先确定用于后续图像分割的总尺度数目k，参数lgscale，lgori以及d的取值范围。在此基础上，记两个连续尺度间的互信息为ami，构建目标函数：
amitot＝-[ami(1，2)+ami(2，3)+...+ami(k-1，k)](10)
ami(k-1，k)表示相邻尺度间的互信息；通过迭代操作在给定范围内遍历所有可能的参数组合，当amitot取得最小值时，即认为此时获得的多尺度边缘图像序列seqopt＝{scale1,scale2...scalek}能够最为充分的反映影像中具有代表性的边缘细节特征，则此参数组合所构建边缘响应模型为最优模型。同时，确定了最优的lgscale,lgori和d。
4、多尺度区域分割及合并
(1)多尺度分割
多尺度分割在所提取的seqopt中进行。令tk＝μk+0.2σk以确定提取区域的阈值。其中，μk和σk为尺度scalek下的灰度均值与标准差，k＝1,2...k。在最粗糙尺度scalek中，检测所有灰度值小于阈值tk的像素，进而采用4连通方法获得种子区域，其余像素按照j-value值及与种子区域的距离按照从小到大的顺序逐个并入相邻的种子区域，从而获得scalek下的对象集合rk＝{rk1,rk2...rkn}。在下一个较精细尺度scalek-1中，分别以尺度scalek中提取的每个对象为基本分割单元，计算tk-1从而构建种子区域并进行区域增长。以此类推，直到最精细尺度scale1计算完毕，可获得初步分割结果r1。
(2)基于ssim的区域合并
由于初步分割结果中难以避免的会存在过分割现象，因此依据相邻区域在各个尺度及各个波段中的总体相似性来进行区域合并。选择ssim(结构相似性)对分割结果中的相邻对象进行描述。
在r1中，计算任意两个相邻对象ra和rb的结构相似性s(a,b)：
其中，a、b分别ra和rb对应的特征向量，μa，μb，σa，σb，
在seqopt的所有尺度及原始影像的所有波段中分别计算s(a,b)，并取均值为
step1：依据r1生成区域邻接图rag；
step2：对任意一个对象ra，根据rag确定所有与ra相邻的对象；
step3：计算ra与相邻的一个对象rb，计算所有尺度以及原始影像所有波段下ra和rb的欧式距离均值
step4：若
step5：更新rag，重复step2至step5，遍历所有对象及其相邻对象，直到没有能够合并的区域为止，从而获得最终分割结果。
实验与分析
1、实验数据与比较方法
实验采用两幅城市场景下不同传感器类型的多光谱高分辨率遥感影像。实验一采用航空遥感dom(digitalorthophotomap)影像，数据包括red、green、blue三波段，采集时间为2009年3月，所在地区为中国江苏南京，空间分辨率为0.6m，尺寸为1024×1024像素，如图3所示。实验二采用worldview-2卫星遥感影像，数据包括red、green、blue以及nearinfrared四个波段，采集时间为2011年8月，所在地区为中国重庆，空间分辨率为2m，尺寸为512×512像素，如图4所示。
实验中，除了相位一致性模型(称为“方法2”)和局部同质性模型(称为“方法3”)外，还选择了其他两种先进的高分辨率遥感影像图像分割方法进行对比。一种方法为前文介绍过的zhang等提出的多尺度改进分水岭分割算法(下文简称“方法4”)，另一种是商业软件ecognition(下文简称“方法5”)。其中，德国definiensimaging公司开发的ecognition是目前公认最好的面向高分辨率遥感影像分类的商业软件，其综合考虑了对象的光谱特征形状特征，实现了基于分型网络演化(fractalnetevolutionapproach，fnea)方法的多分辨率分割。ecognition的分割过程主要有三个参数需要进行人工设定，包括：scale；shape；compactness。本文在实验中通过改变以上三个参数，采用试错法选择分割效果最佳的参数组合。
2、参数设定
根据本发明提出的基于尺度间互信息最小化的目标函数优化策略，通过迭代运算确定自适应的最优参数组合。该参数组合可以看作是通过公式(9)计算多尺度序列的最优尺度参数集。
在两组实验中，本文方法分割尺度均设定为k＝5，并根据公式(10)对相关参数进行优化。实验一中所提取的最优参数组合为：lgscale＝4，lgori＝8，d＝{5,7,8,11,15}。实验二中，所提取的最优参数组合为：lgscale＝4，lgori＝8，d＝{5,6,7,9,12}。对于方法2和方法3，分别从不同模型的最优参数组合中选择相关参数，并在两个实验中采用相同的多尺度分割和区域合并策略。对于其他两种方法，采用试错法来保证最佳的尺度参数和最佳的分割效果。基于此，方法4在两组实验中设定滤波系数c＝(15,1/7,9)，标记阈值为10，区域合并的参数为(70，0.5，0.5)。对于方法5，实验一中设定参数scale为110，shape为0.5，compactness为0.5，实验二中设定参数scale为85，shape为0.5，compactness为0.5。
3、分割结果与分析
将所提出的边缘响应模型与相位一致性模型以及基于j-value的边缘检测结果进行比较，如图5、图6所示。
通过观察图5(a)和图6(a)可以发现，仅依赖相位一致性信息尽管能够提取图像的边缘细节信息，但真实边缘却不突出，反而在一些非边缘区域出现了伪边缘现象。图5(b)和图6(b)为d＝5时j-value提取结果。可以看出，仅依赖梯度信息能够准确反映强边缘，但却损失了大量弱边缘信息。而在图5(c)和图6(c)中，通过结合相位一致性与局部同质性信息后的边缘检测结果，真实边缘尤其是弱边缘得到了突出，而伪边缘则得到了有效抑制。采用psnr(peaksignaltonoiseratio)对边缘检测结果进行进一步的定量分析，结果如表1所示。
表1不同方法边缘检测精度的峰值信噪比评价
由于psnr越高时边缘检测结果越接近真实边界，因此本发明所构建的边缘响应模型优于其他两种策略，与目视分析结果一致。在此基础上，对比5种方法的分割结果，如图7、图8所示。
如图所示，为便于分析，我们对影像中的典型地物或关键位置采用字母a～l进行了标注。在相同的多尺度分割和区域合并策略下，三种边缘响应模型的分割结果存在显著差异。本发明提出的方法以及方法2比方法3对弱边缘更为敏感，例如在位置b，只有前两种方法提取到了建筑物与其相邻道路之间的弱边缘。方法2分割结果的主要缺陷是会产生伪边缘。例如，在e位置，方法2在湖泊大面积区域有过分割现象，同样的现象也出现在荒地区域f位置和建筑物屋顶i位置。与其他两种方法相比，对于植被和荒地区域，这三种算法具有几乎相同的分割性能，如在位置f。由于方法4具有突出的过度分割现象，因此其分割结果中存在大量的碎片区域，如位置d、e、h。此外，由于只使用局部纹理，与方法3一样，方法4在弱边缘检测能力上显著低于其他两种方法，例如在位置a和g。方法5则主要存在大量类似细长形状区域的虚假边缘，例如c、i、j、k、l等。而对于较大尺寸的均质区域，如位置e的湖泊，仅有本发明提出的算法既能够保持整体水体轮廓的完整，又能够准确提取湖面上的小岛等弱小目标。因此，通过对五种分割方法的比较，可以初步得出结论：本发明提出的方法能够准确定位建筑物、城市景观道路等典型人工地理对象的边缘，提取的轮廓最完整，明显优于其他四种方法。
3、定量分析与评价
(1)评价指标
为了更客观地评价性能算法，采用同质性指数j值和异质性指数lmi(localmoran'si)对两组实验进行了定量分析。选择这两个指标的原因是，j值和lmi可以分别评价内部对象的均匀性和相邻对象之间的异质性。其定义和计算如下：
当分割精度越高时，则分割结果中对象内部的平均均质程度越高，即所有对象所对应的j-value的均值
其中，u为影像中的像素总数，v为分割结果中的对象总数，wv为对象rv内部的像素总数，b为原始影像的波段总数，jvb表示波段b下对象rv的j-value。
对于lmi，定义如下：
li是分割结果中第i个对象的平均灰度值，
v是分割结果中的对象总数，e2可通过以下方式计算：
wij为相邻节点权重矩阵，表示每个ri和rj.节点之间的邻域空间关系。
假设ri和rj没有相邻关系，wij＝1，否则wij＝0。为了便于lmii的分析，假定lmii取值范围(0，1)，并使用如下归一化公式：
lmiinor＝(lmii-lmimin)/(lmimax-lmimin)(16)
在此基础上，归一化lmi的均值可采用如下公式计算：
其中，更小的
(2)评价结果
基于j值和lmi的两个实验的精度评估结果如表2所示。
表2.定量精度评价
基于表2可以看出，本发明所提出方法分割精度明显优于其它方法，与目视分析结果一致。基于该评价准则，进一步对不同方法分割结果中的j-value和lmiinor分布情况进行统计。将每个对象rv对应的j-value和lmiinor按照大小在[0,1]区间内量化为20个单位。以纵坐标为分割结果中对应不同j-value/lmiinor值的对象在总对象中所占的比例，横坐标为j-value/lmiinor值，构建分布曲线，如图9和图10所示。图中，不同方法的j-value和lmiinor分布曲线采用不同样式的线型进行表示，三条竖线分别代表每种方法的精度评价指标j。根据j-value分布情况可以看出，在[0.05,0.15]这一对象内部的均质程度很高的区间，本发明方法获得的比值更高，也从一定程度上说明了本发明方法定位对象的边缘更为准确。同样，根据lmiinor分布情况，相邻目标之间的高度异质性对应的区间[0.05,0.2]也存在相似的情况。另外，实验二中5种算法的分割精度均较实验一有所降低，且更加集中于均值程度较低的区间，我们认为主要是由于实验二的场景中人造地物形状更加复杂多样且分布更加密集所导致的。
实验表明，本发明方法不仅具有较强的弱边缘检测能力，而且对伪边缘等干扰具有较强的鲁棒性。与传统依赖梯度信息的分割方法相比，本发明方法能更准确地定位对象边缘，保持对象轮廓更加完整，能够较好地抑制过分割现象；与econgnition相比，本发明所提出的方法不仅分割精度高，而且分割结果中不存在细长形状区域的虚假边缘。在此基础上，通过对三种不同边缘模型的对比实验，得出将一致性与局部均匀性相结合是可行和有效的。最后，由于主要参数均可以自适应确定，因此本发明方法还具有较高的自动化程度。