CN102789631A - 基于图像分割和引力模型的遥感图像分辨率提高处理方法 - Google Patents

Abstract

本发明涉及一种基于图像分割和引力模型的遥感图像分辨率提高处理方法，该方法包括：对图像进行预处理；对图像进行多尺度图像分割；利用正交子空间投影方法遴选出待选端元；构建端元分层的混合像素模型；利用混合像素分解方法获得端元在各混合像素中的百分比值；逐层筛选混合像素模型；得到最佳混合像素分解图；将各端元的混合像素分解图转换为亚像素图像；不断调整亚像素的位置，使得所有亚像素间的总引力达到最大；重复调整所有像素中的亚像素直到遍历图像中的所有像素至处理结束；处理完成后，即得到空间分辨率提高的亚像素级图像。与现有技术相比，本发明具有方法简单、不依赖高空间分辨率，且抗噪音性好、能节约大量时间等优点。

Description

基于图像分割和引力模型的遥感图像分辨率提高处理方法

技术领域

本发明涉及一种遥感图像分辨率提高处理方法，尤其是涉及一种基于图像分割和引力模型的遥感图像分辨率提高处理方法。

背景技术

对于光学遥感器系统而言，图像的空间分辨率和光谱分辨率是一对矛盾，在给定信噪比的条件下，较高的光谱分辨率(窄光谱波段)往往意味着以降低空间分辨率为代价，从而使图像目标识别的精确度和识别程度都受到了很大的限制。在保留光谱信息的前提下，提高空间分辨率具有重要意义。提高硬件技术是提高图像空间分辨率的最直接的方法，而在既有的硬件条件限制下，算法和软件技术成为提高图像空间分辨率的主要途径。归纳起来，常用的图像空间分辨率提高技术包括重采样、图像融合方法、混合像素模型方法等三大类。

(1)重采样方法

利用采样方法提高图像空间分辨率主要从光谱域和空间域两方面的技术。Tsai和Huang(1984)最早提出一种光谱反假频方法从光谱域提高图像的空间分辨率的。技术研究证明可以从一系列的多幅较低空间分辨率遥感图像中计算得到较高空间分辨率的图像。Kim和Bose等(1990，1993)改进了这种模型方法，利用最小二乘和总体最小二乘方法和离散傅里叶变换，从较低分辨率的采样过疏的图像中重构得到高分辨图像，这些低分辨率图像可为不同的空间分辨率，有噪音混淆信号影响。Atkinson(1988)，Schultz和Stevenson(1996)，Patti等(1997)和Eren(1997)，Hao(1998，1999)提出了一些空间域采样方法提高空间分辨率，其中Atkinson(1988)使用重采样方法将较低空间分辨率的图像与较高空间分辨率图像和空间矢量数据进行配准。通过与不经过重采样的图像模式识别进行比较，分析几种不同的重采样方法对图像模式识别的不同影响，发现了重采样与图像模式识别之间的内在关系。

(2)图像融合方法

重采样方法是较为基本的空间分辨率提高技术，但随着越来越多的高空间分辨率传感器出现，将图像与高空间分辨率图像进行融合处理，可以大大改善成像图像空间分辨率，同时保持其光谱物理特性和波段形态。采用图像融合方法提高图像空间分辨率的技术有很多，其中Albertz和Zelianeos(1990)提出一种数据集聚方法增强图像的空间分辨率。这种方法以地貌特征为辅助信息，通过傅里叶变换合并同一区域的多个单幅图像，合并后的图像具有更高的空间分辨率。Garguet-Duport等(1996)利用小波变换融合SPOT图像的全色波段和多光谱波段数据，并将小波变换方法与IHS和P+XS等融合方法进行了比较，试验证明小波变换使得融合图像中光谱特征信息失真最小。Wald等(1997)提出一系列准则评价提高图像空间分辨率的图像融合方法的精度，试验表明对不同级别的光谱分辨率和空间分辨率的图像，应该使用不同的融合方法。Ranchln和Wald(2000)设计了一种ARSIS图像多分辨率模型，这种模型结合精确的光谱信息，利用图像融合方法提高图像空间分辨率。Aiazzi等(2006)提出了一种多分辨率图像融合方法，这种方法依靠一种多尺度过采样的广义拉普拉斯金字塔信号还原滤波器，利用调制传递函数将低空间分辨率的光谱数据计算扩展到较高空间分辨率。Winter等(2007)和Capobianco等(2007)分别利用Ikonos多光谱图像和ALI全色图像跟Hyperion高光谱图像进行图像融合提高Hyperion高光谱图像的空间分辨率。Rao等(2008)提出了一种小波变换图像融合方法，同时提高图像的空间分辨率和光谱分辨率。

(3)混合像素模型方法

由于图像的空间分辨率的限制，图像中的像素，很少是由单一均匀的目标物类型组成，而一般都是几种类型目标物的混合。这种像素称为“混合像素”。混合像素现象给图像的处理和实际应用带来了极大的困难，已经成为图像处理和模式识别定量化深入发展的重要障碍。以图像的模式分类为例，对图像的像素级模式分类以传统的硬分类技术为主。而对于混合像素，将其按照一般硬分类方法归属为任何一类都是不准确的，因为它不可能完全属于这种单一类型。在实际应用中分析者常常需要更为准确的信息，包括混合像素内包含哪些类别成分、各成分所占的比例是多少、各成分在混合像素内的空间分布是怎样的，等等。随着对各应用领域对图像处理技术要求的提高，图像处理从像素级发展为亚像素级已成必然，由此而带来的许多技术难点急待改进或解决。

因此，除了以上两类技术可用来提高图像的空间分辨率，目前也有很多的科学家利用混合像素模型方法提高图像的空间分辨率。从空间特征上，混合像素模型方法基于图像中的每个像素中由可分解的亚像素以一定的空间位置排列构成，且每种目标物类型在每个混合像素中占一定的面积比例的前提条件，先通过光谱库或图像本身获取必要的先验信息，得到所有图像目标物类型的纯净像素(端元)，这一过程称为端元提取，每一端元对应于一种目标物类型；然后求解混合像素中各端元所占的面积比，这一过程称为混合像素分解；在混合像素分解后，将原始混合像素划分为单位更小的亚像素，使混合像素中端元所占亚像素数目的比例与端元在该混合像素中的面积比相等，利用目标物类型的空间相关性最终得到混合像素中各端元所属亚像素的空间分布情况，从而获得空间分辨率提高的图像，为科学家提供更可靠的视觉和数量依据，也为图像更高精度上的后续应用提供可能。

混合像素分解现在已成为提高高图像空间分辨率的一种重要手段。迄今为止，应用于混合像素分解的方法已经很多，主要有基于最小二乘混合模型、多层感知器、人工神经网络、遗传算法、最近邻分类器、独立成分分析、支持向量机和最小体积约束的非负矩阵分解等。

在混合像素分解的基础上，亚像素制图技术能够得到较原始图像更高空间分辨率的图像，从而为研究者和决策者提供更可靠的视觉和数量依据，也为高光谱图像更高精度上的后续应用提供可能。Atkinson最早引入亚像素制图概念，在软分类的基础上获得亚像素的位置信息，并最早提出了像素置换亚像素成图方法。随后，许多科学家利用各种方法进行亚像素制图，从而提高图像的空间分辨率，例如基于图像锐化的模糊软分类方法、基于图像融合的图像锐化方法、神经网络和遗传算法等智能方法、小波、基于卫星传感器点扩散函数的方法、基于半方差函数线性优化、马尔可夫随机场和协同克里格的地统计方法、元胞自动机进化模型、基于MAP正则化模型和进化Agent技术等。

目前已有的技术方法，在实际应用中，都存在以下不足：

(a)重采样方法只是简单地对图像的灰度值进行数值计算，既不考虑目标物的空间相关性，更不利用目标物的图像光谱特性，只是将图像分辨率硬性地拉大，其精度和准确性属于粗略形式，适用性很差。

(b)图像融合能在保留图像光谱特征的同时，很好地提高空间分辨率。然而，这种方法依赖于已有的更高空间分辨率的图像，对既定条件下的图像，则无法实现空间分辨率的提高。

(c)混合像素模型分析方面，已有的方法都主要从图像光谱特征的角度进行分析，通常都忽略了图像中目标物的空间相关性等信息，且已有的方法中，都是使用一个固定不变的端元集对整个图像的所有像素进行混合像素分解，这样不仅缺乏灵活性，且稳健性也很低。

(d)亚像素制图方面，已有的亚像素制图方法大多只能针对两个端元的情况进行亚像素制图，针对多个端元的亚像素制图中亚像素制图不稳定；已有的亚像素制图方法受初始化影响大；已有的亚像素制图技术一次迭代只进行一次亚像素调整，花费太多的时间来调整亚像素的位置，方法的时间开销大；已有亚像素制图方法都只考虑与亚像素直接相邻接的亚像素和像素间的空间相关性，忽略更大邻域中的亚像素和像素的空间相关性，实际上目标物的空间相关性是确定亚像素位置的重要因素。本发明能通过最大化亚像素之间的总引力有效地重构得到空间分辨率提高的图像。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种方法简单、不依赖高空间分辨率，且抗噪音性好、能节约大量时间的基于图像分割和引力模型的遥感图像分辨率提高处理方法。

本发明的目的可以通过以下技术方案来实现：一种基于图像分割和引力模型的遥感图像分辨率提高处理方法，其特征在于，该方法包括以下步骤：1)对图像进行波段选择、辐射纠正和几何纠正的预处理；2)对图像进行基于区域生长的多尺度图像分割，将图像分割成多个分割块；3)在每个分割块中，利用正交子空间投影方法获取端元，遴选出多个待选端元；4)利用待选端元进行分层，构建端元分层的混合像素模型，不同端元组合对应了不同的混合像素模型；5)利用带约束的最小二乘混合像素分解方法对所有的端元组合进行计算比较，具有最小的分解误差的端元组合为选定的最佳端元集合，并同时获得端元组合相应端元在各混合像素中的百分比值；6)对端元分层的混合像素模型进行逐层筛选；7)由混合像素分解最终得到在每个像素中各种端元的面积百分比形成的最佳混合像素分解图；8)将各端元的混合像素分解图转换为亚像素图像；9)基于引力模型的亚像素制图技术，计算得到亚像素的引力大小和亚像素的引力标识值，不断调整亚像素的位置，使得所有亚像素间的总引力达到最大；10)开始处理下一像素中的亚像素，重复步骤1)～9)处理所有像素中的亚像素直到遍历图像中的所有像素；11)所有像素处理完成后，即得到空间分辨率提高的亚像素级图像。

所述的步骤2)中的对图像进行基于区域生长的多尺度图像分割包括以下步骤：

21)首先设定一空分割块集合

将每个像素初始化值置为0，看成是待处理的像素，设定分割块序号k＝1；

22)从图像的左上角位置开始遍历未处理过的像素，从未处理过的像素中任意选定一个，把它看做一个初始区域s_k，同时选定的第一个像素的值设为k；

23)依次检测区域s_k中每一像素v_x，y的8邻域内的未处理的像素v_x′，y′与该区域s_k之间的相关性，从而判断像素v_x′，y′是否能够被添加到区域s_k中，此处像素v_x′，y′和v_x，y都为列向量表示的像素值，像素v_x′，y′与区域s_k之间的相关性可由三个指标来衡量：

(i)像素v_x′，y′与像素v_x，y之间的光谱角，光谱角

可由以下公式计算得出；

(式-1)

式中，||v||表示向量v的二阶范数；

(ii)最大光谱角

并由以下公式计算得出：

(式-2)

式中，

由式(式-1)计算得到，表示区域sk中的每个像素；

(iii)均方根误差(RMS Error)

${\mathrm{\ϵ}}_{{s_k}^{*}}=\frac{1}{N}\sqrt{{\mathrm{\Σ}}_{t=1}^N{(v_{x^{*},y^{*}}-{\mathrm{\Σ}}_{t=1}^N{v}_{x^{*},y^{*}}/N)}^2}$ (式-3)

式中，s_k ^*＝v_x′，y′∪s_k，即假定像素v_x′，y′添加到当前区域sk中而生长成新的区域sk*；

N是新区域sk*的像素数目；

24)算法中为三个指标

和都设定一个阈值，当三个指标的值都小于预设阈值时，像素v_x′，y′将被添加到区域s_k中，得到新的区域s_k，像素v_x′，y′也将被标识为已处理；

25)循环步骤23)～25)检测区域中每个像素，直到没有新的邻域像素被添加到区域s_k中，区域s_k不再生长；

26)添加区域s_k到分割块集合S中，使得S＝{s_i|i＝1，2，...，k}，并且设k＝k+1，然后算法返回到步骤22)，再次遍历图像中未处理过的像素，从未处理过的像素中任意选定一个，把它看做一个新的初始区域s_k，直到图像中所有像素都被处理完，算法即结束，由此得到基于区域生长的图像分割图像，每个区域表示了一个分割块。

所述的步骤3)中的利用正交子空间投影方法获取端元具体包括以下步骤：

31)选定一任意像素向量t₀，设θ为误差阈值；设选定的端元集合

且迭代次数i＝0；

32)对图像中所有的像素v运用正交投影方法计算像素集的正交投影矩阵，正交投影是一种空间投影算子，可由下式得到：

P_E＝I-EE^#＝I-E(E^TE)^-1E^T                        (式-4)

式中E^#＝(E^TE)^-1E^T为端元集E的伪逆矩阵，然后再计算基于凸多面锥的子空间投影值：

O(E，v)＝(P_Ev)^T(P_Ev)                        (式-5)

式中E为矩阵形式的端元集，v表示图像中任一像素，P_E为(式-4)式给定的正交子空间投影算子；

33)搜寻具有最大凸多面锥的子空间投影值的像素，记作t₁：

$t_1=\arg \left\{\underset{v}{\max }\right[{\left(P_{t_0}v\right)}^T\left(P_{t_0}v\right)\left]\right\},$ (式-6)

式中，

表示使O(t₀，v)达到最大值的v的值，并且设置E₁＝[t₁]且迭代次数i＝1；

34)如果

则算法跳至第38)步进行；否则，设迭代次数i＝i+1，并往下进行；

35)搜寻具有最大凸多面锥的子空间投影值的像素，记作t_i：

$t_i=\arg \left\{\underset{v}{\max }\right[{\left(P_{\left[t_0{E}_{i-1}\right]}v\right)}^T\left(P_{\left[t_0{E}_{i-1}\right]}v\right)\left]\right\},$ (式-7)

式中E_i-1＝[t₁ t₂...t_i-1]为第(i-1)次迭代过程中生成的端元集；[t₀ E_i-1]＝[t₀ t₁...t_i-1]；

36)设E_i＝[t₁ t₂...t_i]，并计算并比较η_i和误差阈值θ的值；

37)如果η_i＞θ，算法返回到第35)步；否则就往下进行，在此，如果端元集E中的端元数目i达到了预设的每个分割块组的最大端元数目，算法也往下进行；

38)到此步骤，正交子空间投影方法进行完毕，至此，端元集E包含有i个端元向量，其中不包括初始向量t₀。

所述的步骤9)具体包括以下步骤：

41)根据(式-9)计算得到亚像素的引力大小p和亚像素的引力标识值p′，不断调整亚像素的位置，使得下式成立：

$\max P=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}$ (式-8)

式中，N为图像的像素总数，M为每个像素中的亚像素数目(M＝n²)，当所有亚像素间的总引力达到最大时，得到的图像即最终空间分辨率提高的亚像素图像；

42)在亚像素基础上，按照同质相吸、异质相斥的规则，计算该像素的邻域范围内亚像素之间的引力，亚像素之间的引力根据物体间引力模型公式计算：

$p_{12}=\mathrm{\ρ}\·\frac{m_1{m}_2}{r_{12}^2}$ (式-9)

式中p₁₂为亚像素间的引力大小；m₁、m₂为亚像素的质量；r₁₂为亚像素之间的距离，引力大小与亚像素质量成正比，与亚像素之间的距离的平方成反比关系，ρ为引力系数，相同类型亚像素间ρ＝1，不同类型亚像素间ρ≤0，考虑到图像中目标物的分布服从近似的高斯正态分布，所以不同类型目标物的亚像素间ρ＝-1/e，即：

(式-10)

式中，e为指数函数的底，计算引力大小时，按照相同类型亚像素间ρ＝1，不同类型亚像素间ρ＝-1/e的规则；计算引力标识值时，按照相同类型亚像素间ρ＝1，不同类型亚像素间ρ＝0的规则；

43)搜索每种目标物类型对应的亚像素中引力标识值最大和最小值所对应的亚像素位置x_imax和x_imin：

x_imax＝{x_i|p′_imax＝max(p′_i)}                        (式-11)

x_imin＝{x_i|p′_imin＝min(p′_i)}                        (式-12)

并搜索所有亚像素中引力标识值最大和最小值所对应的亚像素位置x_max和x_min：

x_max＝{x_imax|p′_max＝max(p′_imax)}                    (式-13)

x_min＝{x_imin|p′_min＝min(p′_imin)}                    (式-14)

44)由此，如果p(x_min)＜p(x_max)，则置换x_min和x_max位置上的亚像素，并返回到42)子步骤；否则不进行置换，进行步骤45)；

45)比较第二最大和第二最小引力标识值亚像素上的引力大小值以进行置换，依次类推，当置换(或不置换)循环次数小于端元类别数C时，返回到步骤42)进行下一次引力计算；如果C次比较过程中，都没有亚像素置换，则表示已完成亚像素调整。

与现有技术相比，本发明具有以下优点：

1、重采样方法只是简单地对图像的灰度值进行数值计算，既不考虑目标物的空间相关性，更不利用目标物的图像光谱特性，只是将图像分辨率硬性地拉大，其精度和准确性属于粗略形式，适用性很差。而本发明充分利用目标物像素的空间相关性，在保留图像光谱特征的同时，能很好的提高图像的空间分辨率。

2、图像融合依赖于已有的更高空间分辨率的图像，对既定条件下的图像，则无法实现空间分辨率的提高。而本发明不受条件的限制，无需有更高空间分辨率图像，只需要本图像包含的多个图像波段即可实现保留图像光谱特性的同时，很好地提高图像的空间分辨率。

3、已有的方法都主要从图像光谱特征的角度进行混合像素模型分析，通常都忽略了图像中目标物的空间相关性等信息，且已有的方法中，都是使用一个固定不变的端元集对整个图像的所有像素进行混合像素分解，这样不仅缺乏灵活性，且稳健性也很低。本发明使用一种基于图像分割的混合像素模型方法。不同于以往的方法，基于图像分割的混合像素模型方法先对图像进行基于区域生长的多尺度图像分割，然后在每个分割块内获取统计端元，同时利用了像素的光谱特征和空间特征信息，不仅提取到更加精确的端元，而且该方法较其它方法具有更好的抗噪音性。

4、已有的亚像素制图方法大多只能针对两个端元的情况进行亚像素制图，本发明针对多个端元的亚像素制图，具有较好的制图精度和稳定性；已有的亚像素制图方法受初始化影响大，本发明可采用顺序初始化和随机初始化，且制图精度和效果不受初始化的影响；已有的亚像素制图技术一次迭代过程中只进行一次亚像素调整，花费太多的时间来调整亚像素的位置，本发明一次迭代过程可以将一个像素中所有符合条件的亚像素进行调整，节省了大量的时间消耗；已有亚像素制图方法都只考虑与亚像素直接相邻接的亚像素和像素间的空间相关性，本发明考虑更加大邻域范围的亚像素和像素的空间相关性，更加符合目标物空间相关性的实际情况，从而制图结果也更加精确和可靠。

附图说明

图1为本发明的流程图；

图2为本发明中亚像素制图过程示意图；

图3为本发明中亚像素制图初始化过程示意图，其中(31)为顺序初始化放置，(32)为随机初始化放置；

图4为本发明中端元分层的混合像素模型原理框架图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例

如图1～图4所示，该发明的具体实施过程包括以下步骤：

(1)对图像进行波段选择、辐射纠正和几何纠正等预处理；

(2)对图像进行基于区域生长的多尺度图像分割；

(3)在每个分割块中，利用正交子空间投影方法获取端元，遴选出多个待选端元；

(4)利用待选端元进行分层，构建端元分层的混合像素模型，不同端元组合对应了不同的混合像素模型；

(5)利用带约束的最小二乘混合像素分解方法对所有的端元组合进行计算比较，具有最小的分解误差的端元组合为选定的最佳端元集合，并同时获得端元组合相应端元在各混合像素中的百分比值；

(6)对端元分层的混合像素模型进行逐层筛选；

(7)由混合像素分解最终得到在每个像素中各种端元的面积百分比(0％～100％)形成的最佳混合像素分解图，如图2中的(a)和图3中的(d)所示；

(8)将各端元的混合像素分解图转换为亚像素图像。此处，比例尺度为n，即每个像素被分为n×n＝n²个亚像素，每个亚像素对应等于1/n2个像素的面积。初始时随机放置或者顺序放置亚像素，如图2中的(b)和图3中的(e)所示；

(9)基于引力模型的亚像素制图技术：根据式(2)计算得到亚像素的引力大小p和亚像素的引力标识值p′，不断调整亚像素的位置，使得下式成立：

$\max P=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}---\left(1\right)$

式中，N为图像的像素总数，M为每个像素中的亚像素数目(M＝n²)。当所有亚像素间的总引力达到最大时，得到的图像即最终空间分辨率提高的亚像素图像，如图2中的(f)所示。

(9-1)在亚像素基础上，按照同质相吸、异质相斥的规则，计算该像素的邻域域范围内亚像素之间的引力，亚像素之间的引力根据物体间引力模型公式计算：

$p_{12}=\mathrm{\ρ}\·\frac{m_1{m}_2}{r_{12}^2}---\left(2\right)$

式中p₁₂为亚像素间的引力大小；m₁、m₂为亚像素的质量；r₁₂为亚像素之间的距离，引力大小与亚像素质量成正比，与亚像素之间的距离的平方成反比关系。ρ为引力系数，相同类型亚像素间ρ＝1，不同类型亚像素间ρ≤0，考虑到图像中目标物的分布服从近似的高斯正态分布，所以不同类型目标物的亚像素间ρ＝-1/e，即：

式中，e为指数函数的底。计算引力大小时，按照相同类型亚像素间ρ＝1，不同类型亚像素间ρ＝-1/e的规则；计算引力标识值时，按照相同类型亚像素间ρ＝1，不同类型亚像素间ρ＝0的规则。

(9-2)搜索每种目标物类型对应的亚像素中引力标识值最大和最小值所对应的亚像素位置x_imax和x_imin：

x_imax＝{x_i|p′_imax＝max(p′_i)}                        (4)

x_imin＝{x_i|p′_imin＝min(p′_i)}                        (5)

并搜索所有亚像素中引力标识值最大和最小值所对应的亚像素位置x_max和x_min：

x_max＝{x_imax|p′_max＝max(p′_imax)}                    (6)

x_min＝{x_imin|p′_min＝min(p′_imin)}                    (7)

(9-3)由此，如果p(x_min)＜p(x_max)，则置换x_min和x_max位置上的亚像素，并返回到(9-1)子步骤；否则不进行置换，往(9-4)子步骤；

(9-4)比较第二最大和第二最小引力标识值亚像素上的引力大小值以进行置换，依次类推，当置换(或不置换)循环次数小于端元类别数C时，返回到(9-1)子步骤进行下一次引力计算；如果C次比较过程中，都没有亚像素置换，则至(10)步骤；

(10)返回至(9)步骤，从(9-1)子步骤开始调整下一像素中的亚像素，重复调整直到遍历图像中的所有像素至算法结束。

(11)算法成功执行完成后，即得到空间分辨率提高的亚像素级图像。

在此，对(2)～(6)步骤中的本发明中的图像分割、端元提取和混合像素分解还做以下详细说明：

图像分割：

本发明运用一种基于区域生长的图像分割方法将图像分割成许多分割块。每一个分割块可由多种特征来综合描述，包括图像像素值和空间特征，如纹理、形状等，其基本原理可以概括如下：

(a)首先设定一空分割块集合

将每个像素初始化值置为0，看成是待处理的像素，设定分割块序号k＝1；

(b)从图像的左上角位置开始遍历未处理过的像素，从未处理过的像素中任意选定一个，把它看做一个初始区域sk，同时选定的第一个像素的值设为k；

(c)依次检测区域s_k中每一像素v_x，y的8邻域内的未处理的像素v_x′，y′与该区域s_k之间的相关性，从而判断像素v_x′，y′是否能够被添加到区域s_k中。此处像素v_x′，y′和v_x，y都为列向量表示的像素值。像素v_x′，y′与区域sk之间的相关性可由三个指标来衡量：

(i)像素v_x′，y′与像素v_x，y之间的光谱角，光谱角

可由以下公式计算得出。

式中，||v||表示向量v的二阶范数；(ii)最大光谱角

并由以下公式计算得出：

式中，由式(8)计算得到，

表示区域sk中的每个像素；(iii)均方根误差(RMS Error)

${\mathrm{\ϵ}}_{{s_k}^{*}}=\frac{1}{N}\sqrt{{\mathrm{\Σ}}_{t=1}^N{(v_{x^{*},y^{*}}-{\mathrm{\Σ}}_{t=1}^N{v}_{x^{*},y^{*}}/N)}^2}---\left(10\right)$

式中，x_k ^*＝v_x′，y′∪s_k，即假定像素v_x′，y′添加到当前区域sk中而生长成新的区域sk*；

N是新区域s_k ^*的像素数目。

(d)算法中为三个指标

和

都设定一个阈值，当三个指标的值都小于预设阈值时，像素v_x′，y′将被添加到区域s_k中，得到新的区域s_k，像素v_x′，y′也将被标识为已处理。

(e)循环(c)～(e)步骤检测区域中每个像素，直到没有新的邻域像素被添加到区域s_k中，区域s_k不再生长；

(f)添加区域s_k到分割块集合S中，使得S＝{x_i|i＝1，2，...，k}，并且设k＝k+1，然后算法返回到步骤(b)，再次遍历图像中未处理过的像素，从未处理过的像素中任意选定一个，把它看做一个新的初始区域s_k。直到图像中所有像素都被处理完，算法即结束，由此得到基于区域生长的图像分割图像，每个区域表示了一个分割块。

端元提取：

本发明使用的正交子空间投影方法遴选端元属于非监督迭代方法，不需要端元的先验信息。方法可以按以下过程进行：

(a)选定一任意像素向量t₀，可以为均值向量，也可以为具有最大模的向量；设θ为误差阈值；设选定的端元集合且迭代次数i＝0。

(b)对图像中所有的像素v运用正交投影方法计算像素集的正交投影矩阵。正交投影是一种空间投影算子，可由下式得到：

P_E＝I-EE^#＝I-E(E^TE)^-1E^T                    (11)

式中E^#＝(E^TE)^-1E^T为端元集E的伪逆矩阵。然后再计算基于凸多面锥的子空间投影值：

O(E，v)＝(P_Ev)^T(P_Ev)                                (12)

式中E为矩阵形式的端元集，v表示图像中任一像素，P_E为(11)式给定的正交子空间投影算子。

(c)搜寻具有最大凸多面锥的子空间投影值的像素，记作t₁：

$t_1=\arg \left\{\underset{v}{\max }\right[{\left(P_{t_0}v\right)}^T\left(P_{t_0}v\right)\left]\right\},---\left(13\right)$

式中，表示使O(t₀，v)达到最大值的v的值。并且设置E₁＝[t₁]且迭代次数i＝1。

(d)如果

则算法跳至第(h)步进行；否则，设迭代次数i＝i+1，并往下进行。

(e)搜寻具有最大凸多面锥的子空间投影值的像素，记作ti：

$t_i=\arg \left\{\underset{v}{\max }\right[{\left(P_{\left[t_0{E}_{i-1}\right]}v\right)}^T\left(P_{\left[t_0{E}_{i-1}\right]}v\right)\left]\right\},---\left(14\right)$

式中E_i-1＝[t₁ t₂...t_i-1]为第(i-1)次迭代过程中生成的端元集；[t₀ E_i-1]＝[t₀ t₁...t_i-1]。

(f)设E_i＝[t₁ t₂...t_i]，并计算并比较η_i和误差阈值θ的值。

(g)如果η_i＞θ，算法返回到第(e)步；否则就往下进行。在此，如果端元集E中的端元数目i的达到了预设的每个分割块组的最大端元数目，算法也往下进行。

(h)到此步骤，正交子空间投影方法进行完毕。至此，端元集E包含有i个端元向量，其中不包括初始向量t₀。

混合像素分解：

本发明使用一种基于端元分层的混合像素分解模型，在该模型中，端元是可变的，模型在图像分割端元提取的基础上，对不同分割块组中的每个像素采用不同的端元集进行混合分解。该模型强调每一个端元均由一组向量而不是一个固定的唯一的向量来表示，当进行混合像素分解的时候，对每个像素，从代表各自端元的向量组中各选取一个最合适的向量，使得均方根误差最小，其中端元可以从图像中或者从该区域的光谱库中选取。其模型方法的基本思想为：

(a)方法使用带约束的最小二乘分解方法对像素进行混合像素分解，采用混合分解的RMS误差来评价待选的端元模型，并为混合分解的RMS误差设定适当的阈值。

(b)使用基于图像分割的端元提取方法获取出原始端元集，在此基础上，建立端元分层的混合像素模型，图4所示为端元分层的混合像素模型原理框架图。图中，

表示在第m层中的端元，Nm为第m层中的端元数目，

表示在第m层中的m-端元模型，nm为第m层中的端元模型数目，且有N₁＞N₂＞...＞N_m和n₁＞n₂＞...＞n_m。

依据图4，按土地覆盖类型的端元划分，第1层模型为原始端元集。基本有效的混合像素分解从第2层开始建立，则从第2层到第7层的端元模型可以定义为：

2-EM模型：水体-非水体2-端元模型；

3-EM模型：水体-植被-不透水层3-端元模型；

4-EM模型：水体-植被-不透水层-土壤4-端元模型；

5-EM模型：水体-植被-不透水层-土壤-阴影5-端元模型；

6-EM模型：水体-阴影-高含水植被-低含水植被-不透水层-土壤6-端元模型；

7-EM模型：水体-阴影-高含水植被-低含水植被-高反射率不透水层-低反射率不透水层-土壤7-端元模型。

本发明的基于端元分层的混合像素分解模型分析方法，首先考虑水体-非水体2-EM模型。利用所有遴选得到的原始端元集，构成端元组合，组成2-EM模型，进行混合像素分解，计算混合分解的RMS误差，选取符合误差阈值条件的2-EM模型组成2-EM模型集S2。

(c)选择任意一个2-EM模型

然后选取与模型

中端元不同的端元作为第3个端元补充到模型

中，由此重复，生成一系列的水体-植被-不透水层3-EM模型。对所有生成的3-EM模型，利用混合像素分解，计算混合分解的RMS误差，选取符合误差阈值条件的3-EM模型组成3-EM模型集S3。

(d)如此重复选择和计算，依次组成水体-植被-不透水层-土壤4-EM模型、水体-植被-不透水层-土壤-阴影5-EM模型等。

(e)对各类模型，利用利用带约束的最小二乘混合像素分解方法，得到各种模型下具有最小分解误差的端元组合为该模型情况下最佳的端元集。比较各种模型情况下的分解误差，具有全局最小分解误差的模型和相应的端元集为最终求解的最佳端元模型。

Claims (4)

1.一种基于图像分割和引力模型的遥感图像分辨率提高处理方法，其特征在于，该方法包括以下步骤：

1)对图像进行波段选择、辐射纠正和几何纠正的预处理；

2)对图像进行基于区域生长的多尺度图像分割，将图像分割成多个分割块；

3)在每个分割块中，利用正交子空间投影方法获取端元，遴选出多个待选端元；

4)利用待选端元进行分层，构建端元分层的混合像素模型，不同端元组合对应了不同的混合像素模型；

5)利用带约束的最小二乘混合像素分解方法对所有的端元组合进行计算比较，具有最小的分解误差的端元组合为选定的最佳端元集合，并同时获得端元组合相应端元在各混合像素中的百分比值；

6)对端元分层的混合像素模型进行逐层筛选；

7)由混合像素分解最终得到在每个像素中各种端元的面积百分比形成的最佳混合像素分解图；

8)将各端元的混合像素分解图转换为亚像素图像；

9)基于引力模型的亚像素制图技术，计算得到亚像素的引力大小和亚像素的引力标识值，不断调整亚像素的位置，使得所有亚像素间的总引力达到最大；

10)开始处理下一像素中的亚像素，重复步骤1)～9)处理所有像素中的亚像素直到遍历图像中的所有像素；

11)所有像素处理完成后，即得到空间分辨率提高的亚像素级图像。

2.根据权利要求1所述的一种基于图像分割和引力模型的遥感图像分辨率提高处理方法，其特征在于，所述的步骤2)中的对图像进行基于区域生长的多尺度图像分割包括以下步骤：

21)首先设定一空分割块集合将每个像素初始化值置为0，看成是待处理的像素，设定分割块序号k＝1；

22)从图像的左上角位置开始遍历未处理过的像素，从未处理过的像素中任意选定一个，把它看做一个初始区域s_k，同时选定的第一个像素的值设为k；

23)依次检测区域s_k中每一像素v_x，y的8邻域内的未处理的像素v_x′，y′与该区域s_k之间的相关性，从而判断像素v_x′，y′是否能够被添加到区域s_k中，此处像素v_x′，y′和v_x，y都为列向量表示的像素值，像素v_x′，y′与区域s_k之间的相关性可由三个指标来衡量：

(i)像素v_x′，y′与像素v_x，y之间的光谱角，光谱角

可由以下公式计算得出；

(式-1)

式中，||v||表示向量v的二阶范数；

(ii)最大光谱角

并由以下公式计算得出：

(式-2)

式中，

由式(式-1)计算得到，

表示区域sk中的每个像素；

(iii)均方根误差(RMS Error)

${\mathrm{\ϵ}}_{{s_k}^{*}}=\frac{1}{N}\sqrt{{\mathrm{\Σ}}_{t=1}^N{(v_{x^{*},y^{*}}-{\mathrm{\Σ}}_{t=1}^N{v}_{x^{*},y^{*}}/N)}^2}$ (式-3)

式中，s_k ^*＝v_x′，y′∪s_k，即假定像素v_x′，y′添加到当前区域sk中而生长成新的区域sk*；

N是新区域sk*的像素数目；

24)算法中为三个指标

和

都设定一个阈值，当三个指标的值都小于预设阈值时，像素v_x′，y′将被添加到区域s_k中，得到新的区域s_k，像素v_x′，y′也将被标识为已处理；

25)循环步骤23)～25)检测区域中每个像素，直到没有新的邻域像素被添加到区域s_k中，区域s_k不再生长；

26)添加区域s_k到分割块集合S中，使得S＝{s_i|i＝1，2，...，k}，并且设k＝k+1，然后算法返回到步骤22)，再次遍历图像中未处理过的像素，从未处理过的像素中任意选定一个，把它看做一个新的初始区域s_k，直到图像中所有像素都被处理完，算法即结束，由此得到基于区域生长的图像分割图像，每个区域表示了一个分割块。

3.根据权利要求1所述的一种基于图像分割和引力模型的遥感图像分辨率提高处理方法，其特征在于，所述的步骤3)中的利用正交子空间投影方法获取端元具体包括以下步骤：

31)选定一任意像素向量t₀，设θ为误差阈值；设选定的端元集合且迭代次数i＝0；

32)对图像中所有的像素v运用正交投影方法计算像素集的正交投影矩阵，正交投影是一种空间投影算子，可由下式得到：

P_E＝I-EE^#＝I-E(E^TE)^-1E^T                    (式-4)

式中E^#＝(E^TE)^-1E^T为端元集E的伪逆矩阵，然后再计算基于凸多面锥的子空间投影值：

O(E，v)＝(P_Ev)^T(P_Ev)                            (式-5)

式中E为矩阵形式的端元集，v表示图像中任一像素，P_E为(式-4)式给定的正交子空间投影算子；

33)搜寻具有最大凸多面锥的子空间投影值的像素，记作t₁：

$t_1=\arg \left\{\underset{v}{\max }\right[{\left(P_{t_0}v\right)}^T\left(P_{t_0}v\right)\left]\right\},$ (式-6)

式中，

表示使O(t₀，v)达到最大值的v的值，并且设置E₁＝[t₁]且迭代次数i＝1；

34)如果

则算法跳至第38)步进行；否则，设迭代次数i＝i+1，并往下进行；

35)搜寻具有最大凸多面锥的子空间投影值的像素，记作t_i：

$t_i=\arg \left\{\underset{v}{\max }\right[{\left(P_{\left[t_0{E}_{i-1}\right]}v\right)}^T\left(P_{\left[t_0{E}_{i-1}\right]}v\right)\left]\right\},$ (式-7)

式中E_i-1＝[t₁ t₂...t_i-1]为第(i-1)次迭代过程中生成的端元集；[t₀ E_i-1]＝[t₀ t₁...t_i-1]；

36)设E_i＝[t₁ t₂...t_i]，并计算

并比较η_i和误差阈值θ的值；

37)如果η_i＞θ，算法返回到第35)步；否则就往下进行，在此，如果端元集E中的端元数目i达到了预设的每个分割块组的最大端元数目，算法也往下进行；

38)到此步骤，正交子空间投影方法进行完毕，至此，端元集E包含有i个端元向量，其中不包括初始向量t₀。

4.根据权利要求1所述的一种基于图像分割和引力模型的遥感图像分辨率提高处理方法，其特征在于，所述的步骤9)具体包括以下步骤：

41)根据(式-9)计算得到亚像素的引力大小p和亚像素的引力标识值p′，不断调整亚像素的位置，使得下式成立：

$\max P=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}$ (式-8)

式中，N为图像的像素总数，M为每个像素中的亚像素数目(M＝n²)，当所有亚像素间的总引力达到最大时，得到的图像即最终空间分辨率提高的亚像素图像；

42)在亚像素基础上，按照同质相吸、异质相斥的规则，计算该像素的邻域范围内亚像素之间的引力，亚像素之间的引力根据物体间引力模型公式计算：

$p_{12}=\mathrm{\ρ}\·\frac{m_1{m}_2}{r_{12}^2}$ (式-9)

式中p₁₂为亚像素间的引力大小；m₁、m₂为亚像素的质量；r₁₂为亚像素之间的距离，引力大小与亚像素质量成正比，与亚像素之间的距离的平方成反比关系，ρ为引力系数，相同类型亚像素间ρ＝1，不同类型亚像素间ρ≤0，考虑到图像中目标物的分布服从近似的高斯正态分布，所以不同类型目标物的亚像素间ρ＝-1/e，即：

(式-10)

式中，e为指数函数的底，计算引力大小时，按照相同类型亚像素间ρ＝1，不同类型亚像素间ρ＝-1/e的规则；计算引力标识值时，按照相同类型亚像素间ρ＝1，不同类型亚像素间ρ＝0的规则；

43)搜索每种目标物类型对应的亚像素中引力标识值最大和最小值所对应的亚像素位置x_imax和x_imin：

x_imax＝{x_i|p′_imax＝max(p′_i)}                            (式-11)

x_imin＝{x_i|p′_imin＝min(p′_i)}                            (式-12)

并搜索所有亚像素中引力标识值最大和最小值所对应的亚像素位置x_max和x_min：

x_max＝{x_imax|p′_max＝max(p′_imax)}                    (式-13)

x_min＝{x_imin|p′_min＝min(p′_imin)}                    (式-14)

44)由此，如果p(x_min)＜p(x_max)，则置换x_min和x_max位置上的亚像素，并返回到42)子步骤；否则不进行置换，进行步骤45)；

45)比较第二最大和第二最小引力标识值亚像素上的引力大小值以进行置换，依次类推，当置换(或不置换)循环次数小于端元类别数C时，返回到步骤42)进行下一次引力计算；如果C次比较过程中，都没有亚像素置换，则表示已完成亚像素调整。

Images