基于无监督与傅里叶核函数的遥感影像像元纯度识别方法

摘要

本发明提供一种基于无监督与傅里叶核函数的遥感影像像元纯度识别方法，包括：分割原始遥感影像，得到P个对象；对P个对象进行聚类操作，得到n个团块；对于任意一个团块ti，设团块ti包括d个对象，均执行以下操作：计算团块ti对应的傅里叶核函数包凸；计算指定对象V

说明书

技术领域

本发明属于遥感影像分析与处理技术领域，具体涉及一种基于无监督与傅 里叶核函数的遥感影像像元纯度识别方法。

背景技术

遥感影像是传感器探测元件对目标地物反射或发射电磁辐射能量的记录和 反映，高空间分辨率的遥感影像可以对地物进行更加精细的描述，提供丰富的 细节信息。

像元为遥感影像的基本单位，主要包括两类：混合像元和纯像元。其中， 混合像元是指包含了多种地物类型的像元，产生原因主要为：地物分布的复杂 性、电磁辐射传输过程中各种环境的影响以及探测元件本身的物理特性等多种 因素。而纯像元是指由单一地物组成的像元。

获得各个地区对应的像元相对纯度，具有重要现实意义，例如：对于纯度 较高的像元，说明地物类型相对单一，适于定量遥感分析；而对于纯度较低的 像元，适于做自动分类的训练数据，相关企业或科研人员可以依据像元的纯度 信息制定实地考察路线与方案，使得考察方案更具效率与理论依据，以较少的 成本和时间获得更具价值的高质量样本，提高遥感影像处理的能力和质量。

但是，现有技术中，对遥感影像像元纯度的识别，主要为定性方法，具有 精度低的不足，且不方便根据遥感影像像元纯度进行后续应用。

发明内容

针对现有技术存在的缺陷，本发明提供一种基于无监督与傅里叶核函数的 遥感影像像元纯度识别方法，为一种定量方法，能够较准确的识别出遥感影像 像元纯度，有利于根据遥感影像像元纯度进行后续应用。

本发明采用的技术方案如下：

本发明提供一种基于无监督与傅里叶核函数的遥感影像像元纯度识别方法， 包括以下步骤：

S1，读取原始遥感影像；其中，所述原始遥感影像由M行N列共M*N个 像元组成，像元集为S＝{S₁，S₂，...，S_M×N}；

S2，分割所述原始遥感影像，得到P个相对独立且大小大于等于1个像元 的对象，表示为：V＝{V₁，V₂，...，V_P}；

S3，对P个对象进行聚类操作，得到n个团块，n个团块用集合T表示， T＝{t1，t2，...，tn}；

S4，对于任意一个团块ti，其中，i∈(1、2...n)，设团块ti包括d个对象； 均执行以下操作：

S41，计算团块ti对应的傅里叶核函数包凸；其中，所述包凸由d个对象中 的b个对象组成；b小于d；

S42，当需要识别d个对象中任意一个指定对象V_f的纯度时，计算指定对象 V_f分别到b个对象中每一个对象的距离，得到距离指定对象V_f距离最小的对象 T，指定对象V_f到对象T的距离为M1；

计算指定对象V_f距离包凸中心的距离M2；

计算对象T距离包凸中心的距离M3；

计算映射距离值为M2／(M3+M1)，该映射距离值即为指定对象V_f的整体纯 度值；

S43，所述指定对象V_f的整体纯度值即为该指定对象V_f所包含的各个像元 的纯度值。

优选的，S2具体为：使用马尔科夫随机场图像模块，基于最大后验概率准 则，对所述原始遥感影像进行图像分割。

优选的，S3具体为：

S31，分别计算P个对象的波段属性向量值B＝{b₁，b₂，...，b_u}；其中，u为波段 属性维度，b₁，b₂，...，b_u分别为对象在各个维度的波段值；

具体为：设原始遥感影像中每一个像元具有u维度波段属性，则对于任意 一个对象Vi，i∈(1、2...P)，其同样具有u维度波段属性，并且，对象Vi中每一 维度波段属性为其所包含的全部像元同一维度波段属性的均值；

S32，通过K-MEANS算法对P个对象进行聚类操作。

优选的，S41具体为：

S411，构造空的初始包凸；

输入团块ti所包括的d个对象，设d个对象分别表示为V₁，V₂，...，V_d；

S412，令z＝1；读取对象V_z，将对象V_z加入到初始包凸之中；

S413，计算对象V_z到初始包凸的距离；

S414，判断该距离是否小于阈值m，如果小于，则表明对象V_z在初始包凸 之中，不重构初始包凸，转到S415；否则，表明对象V_z不在初始包凸之中，重 构初始包凸，然后将对象V_z加入到重构后的包凸之中，用重构后的包凸替代初 始包凸，然后转到S415；

S415，令z＝z+1，循环S412-S415，直到z＝d+1时，终止循环，输出最终得 到的包凸。

优选的，S413，计算对象Vz到初始包凸的距离具体为：

S4131，计算团块ti的标准差σ和均值μ；

其中，标准差σ为团块ti包括的d个对象的波段属性向量值B的标准差； 均值μ为团块ti包括的d个对象的波段属性向量值B的均值；

S4132，依公式一计算正则化因子q：

$q=1-\frac{{\sigma }^2}{{\mu }^2};$    公式一

S4133，引入傅里叶核函数K，设团块ti所包括的d个对象V₁，V₂，...，V_d

的波段属性向量值分别为B₁，B₂，...，B_d；

依公式二计算团块ti内的任意一个对象V_z距离包凸TConvS的距离：

$\mathrm{dis}\tan \mathrm{ce}(V_Z,\mathrm{convS})={\min }_{\theta }(K(B_z,B_z)-2\underset{i=1}{\overset{d}{\Sigma }}{\theta }_{i}K(B_z,B_i)+\underset{i=1}{\overset{d}{\Sigma }}\underset{j=1}{\overset{d}{\Sigma }}{\theta }_i{\theta }_{j}K(B_i,B_j)$    公式二

$s.t.\underset{i=2}{\overset{d}{\Sigma }}{\theta }_i=1,{\theta }_i\ge 0,i=\mathrm{1,2},...,k$

其中，$K(B_z,B_z)=\frac{1-q^2}{2(1-2\mathrm{qCos}\left(\right|B_z-B_z\left|\right)+q^2)};$

$K(B_z,B_i)=\frac{1-q^2}{2(1-2\mathrm{qCos}\left(\right|B_z-B_i\left|\right)+q^2)};$

$K(B_i,B_j)=\frac{1-q^2}{2(1-2\mathrm{qCos}\left(\right|B_i-B_j\left|\right)+q^2)}.$

优选的，对于b个对象中的任意一个对象V_b，设其波段属性向量值为B_b； 设指定对象V_f的波段属性向量值为B_f；

则：S42中，通过以下公式计算指定对象V_f到对象V_b的距离：

优选的，通过以下公式计算指定对象V_f距离包凸中心的距离M2：

包凸由b个对象组成，设分别表示为V₁，V₂，...，V_b，该b个对象的波段属性向 量值分别为B₁，B₂，...，B_b；则依下面公式计算M2：

$\mathrm{discenterCenter}(V_f,\mathrm{TconvS})=K(B_f,B_f)+\frac{1}{b}\underset{i=1}{\overset{b}{\Sigma }}\underset{j=1}{\overset{b}{\Sigma }}K(B_i,B_j)-\frac{2}{b}\underset{i=1}{\overset{b}{\Sigma }}K(B_f,B_i);$

用V_T表示对象T，其波段属性向量值用B_T表示，则依下面公式计算M3：

$\mathrm{discenterCenter}(V_T,\mathrm{TconvS})=K(B_T,B_T)+\frac{1}{b}\underset{i=1}{\overset{b}{\Sigma }}\underset{j=1}{\overset{b}{\Sigma }}K(B_i,B_j)-\frac{2}{b}\underset{i=1}{\overset{b}{\Sigma }}K(B_T,B_i);$

通过以下公式计算M1：

优选的，S4之后，还包括：

S5，团块ti共包括d个对象，设d个对象分别表示为V₁，V₂，...，V_d；分别计算 每一个对象的映射距离值，分别表示为L₁，L₂，...，L_d；

S6，计算L₁，L₂，...，L_d这d个映射距离值的均值μ和标准差σ；

S7，对于任意一个映射距离值L_i，其中，i∈(1、2...d)，判断其是否满足 以下关系式：μ-2×σ＜＝L_i＜＝μ+2×σ；如果不满足，则如果L_i＜μ-2×σ，令L_i＝0；如 果L_i＞μ+2×σ，令L_i＝1；如果满足，则对L_i进行转换计算，令L_i＝L_i-(μ-2×σ)／4×σ， 所得到的L_i为最终映射得到的映射距离值。

本发明的有益效果如下：

本发明提供的基于无监督与傅里叶核函数的遥感影像像元纯度识别方法， 为一种定量方法，能够较准确的识别出遥感影像像元纯度，并已通过大量实验 证明，所计算得到的纯度值与实际像元纯度情况相符，有利于根据遥感影像像 元纯度进行后续应用。

附图说明

图1为本发明提供的基于无监督与傅里叶核函数的遥感影像像元纯度识别 方法的整体流程示意图；

图2为计算团块中任意一个对象纯度的方法流程示意图；

图3为计算团块对应的傅里叶核函数包凸的流程示意图；

图4为计算映射距离值的一个具体示意图；

图5为后续计算超高维度的遥感影像的映射距离值的流程图；

图6为实际拍摄得到的一个遥感影像图；

图7为对图6遥感影像图进行处理得到的第一个团块图；

图8为对图6遥感影像图进行处理得到的第二个团块图；

图9为对图6遥感影像图进行处理得到的第三个团块图；

图10为对图6遥感影像图进行处理得到的第四个团块图；

图11为对图6遥感影像图进行处理得到的第五个团块图。

具体实施方式

以下结合附图对本发明进行详细说明：

如图1所示，本发明提供一种基于无监督与傅里叶核函数的遥感影像像元纯 度识别方法，包括以下步骤：

S1，读取原始遥感影像；其中，所述原始遥感影像由M行N列共M*N个 像元组成，像元集为S＝{S₁，S₂，...，S_M×N}；

S2，分割所述原始遥感影像，得到P个相对独立且大小大于等于1个像元 的对象，表示为：V＝{V₁，V₂，...，V_P}；

本步骤中，可以使用马尔科夫随机场图像模块，基于最大后验概率准则， 对所述原始遥感影像进行图像分割。

需要强调的是，本发明在后续步骤中，均以对象为研究对象进行研究，一 个对象中包括一个以L的像元，而不直接以像元为对象进行研究，具有简少运 算量、提高运算速度的优点。当然，根据实际精度需求，由于一个对象中也可 以只包括一个像元，因此，后续步骤中，如果以像元为研究对象，也属于本发 明的保护范围。

S3，对P个对象进行聚类操作，得到n个团块，n个团块用集合T表示， T＝{t1，t2，...，tn}；

本步骤具体包括以下步骤：

S31，分别计算P个对象的波段属性向量值B＝{b₁，b₂，...，b_u}；其中，u为波段 属性维度，b₁，b₂，...，b_u分别为对象在各个维度的波段值；

具体为：设原始遥感影像中每一个像元具有u维度波段属性，则对于任意 一个对象Vi，i∈(1、2...P)，其同样具有u维度波段属性，并且，对象Vi中每一 维度波段属性为其所包含的全部像元同一维度波段属性的均值；

S32，通过K-MEANS算法对P个对象进行聚类操作。

S4，对于任意一个团块ti，其中，i∈(1、2...n)，设团块ti包括d个对象； 如图2所示，为计算团块中任意一个对象纯度的方法流程示意图，包括以下操 作：

S41，计算团块ti对应的傅里叶核函数包凸；其中，所述包凸由d个对象中 的b个对象组成；b小于d；

如图3所示，为计算团块对应的傅里叶核函数包凸的流程示意图，包括以 下步骤：

S411，构造空的初始包凸；

输入团块ti所包括的d个对象，设d个对象分别表示为V₁，V₂，...，V_d；

S412，令z＝1；读取对象V_z，将对象V_z加入到初始包凸之中；

S413，计算对象V_z到初始包凸的距离；

可以通过以下方法计算：

S4131，计算团块ti的标准差σ和均值μ；

其中，标准差σ为团块ti包括的d个对象的波段属性向量值B的标准差； 均值μ为团块ti包括的d个对象的波段属性向量值B的均值；

S4132，依公式一计算正则化因子q：

$q=1-\frac{{\sigma }^2}{{\mu }^2};$    公式一

S4133，引入傅里叶核函数K，设团块ti所包括的d个对象V₁，V₂，...，V_d

的波段属性向量值分别为B₁，B₂，...，B_d；

依公式二计算团块ti内的任意一个对象Vz距离包凸TConvS的距离：

$\mathrm{dis}\tan \mathrm{ce}(V_Z,\mathrm{convS})={\min }_{\theta }(K(B_z,B_z)-2\underset{i=1}{\overset{d}{\Sigma }}{\theta }_{i}K(B_z,B_i)+\underset{i=1}{\overset{d}{\Sigma }}\underset{j=1}{\overset{d}{\Sigma }}{\theta }_i{\theta }_{j}K(B_i,B_j)$    公式二

$s.t.\underset{i=2}{\overset{d}{\Sigma }}{\theta }_i=1,{\theta }_i\ge 0,i=\mathrm{1,2},...,k$

其中，$K(B_z,B_z)=\frac{1-q^2}{2(1-2\mathrm{qCos}\left(\right|B_z-B_z\left|\right)+q^2)};$

$K(B_z,B_i)=\frac{1-q^2}{2(1-2\mathrm{qCos}\left(\right|B_z-B_i\left|\right)+q^2)};$$K(B_i,B_j)=\frac{1-q^2}{2(1-2\mathrm{qCos}\left(\right|B_i-B_j\left|\right)+q^2)}.$

S414，判断该距离是否小于阈值m，如果小于，则表明对象V_z在初始包凸 之中，不重构初始包凸，转到S415；否则，表明对象V_z不在初始包凸之中，重 构初始包凸，然后将对象V_z加入到重构后的包凸之中，用重构后的包凸替代初 始包凸，然后转到S415；

S415，令z＝z+1，循环S412-S415，直到z＝d+1时，终止循环，输出最终得 到的包凸。

S42，当需要识别d个对象中任意一个指定对象V_f的纯度时，计算指定对象 V_f分别到b个对象中每一个对象的距离，得到距离指定对象V_f距离最小的对象 T，指定对象V_f到对象T的距离为M1；

计算指定对象V_f距离包凸中心的距离M2；

计算对象V_T距离包凸中心的距离M3；

计算映射距离值为M2／(M3+M1)，该映射距离值即为指定对象V_f的整体纯 度值；

参考图4，为计算映射距离值的一个具体示意图，在图4中，包凸由8个对 象组成，分别为V1、V2、V3、V4、V5、V6、V7和V8，W为假想的包凸中心， V8为计算得到的距离指定对象V_f距离最小的对象。本实用新型中，计算距离过 程中不需要知道包凸中心的具体位置，仅需要知道距离包凸中心的距离即可。

在本步骤中，对于b个对象中的任意一个对象V_b，设其波段属性向量值为 B_b；设指定对象V_f的波段属性向量值为B_f；

则：通过以下公式计算指定对象V_f到对象V_b的距离：

通过以下公式计算指定对象V_f距离包凸中心的距离M2：

包凸由b个对象组成，设分别表示为V₁，V₂，...，V_b，该b个对象的波段属性向 量值分别为B₁，B₂，...，B_b；则依下面公式计算M2：

$\mathrm{discenterCenter}(V_f,\mathrm{TconvS})=K(B_f,B_f)+\frac{1}{b}\underset{i=1}{\overset{b}{\Sigma }}\underset{j=1}{\overset{b}{\Sigma }}K(B_i,B_j)-\frac{2}{b}\underset{i=1}{\overset{b}{\Sigma }}K(B_f,B_i);$

用V_T表示对象T，其波段属性向量值用B_T表示，则依下面公式计算M3：

$\mathrm{discenterCenter}(V_T,\mathrm{TconvS})=K(B_T,B_T)+\frac{1}{b}\underset{i=1}{\overset{b}{\Sigma }}\underset{j=1}{\overset{b}{\Sigma }}K(B_i,B_j)-\frac{2}{b}\underset{i=1}{\overset{b}{\Sigma }}K(B_T,B_i);$

通过以下公式计算M1：

S43，所述指定对象V_f的整体纯度值即为该指定对象V_f所包含的各个像元 的纯度值。

通过上述流程即可获得每一个像元的纯度值，对于超高维度的遥感影像数 据，计算距离包凸中心距离时，可能会出现“空心”的现象，即数据并不是分 布在[0，1]之间而是一个非常接近于1的一个较小区间，要获得一个团块中相对纯 度并映射到[0，1]之间，需要继续进行下述的操作，即：如图5所示，为后续计算 超高维度的遥感影像的映射距离值的流程图，包括：

S5，团块ti共包括d个对象，设d个对象分别表示为V₁，V₂，...，V_d；分别计算 每一个对象的映射距离值，分别表示为L₁，L₂，...，L_d；

S6，计算L₁，L₂，...，L_d这d个映射距离值的均值μ和标准差σ；

S7，对于任意一个映射距离值L_i，其中，i∈(1、2...d)，判断其是否满足 以下关系式：μ-2×σ＜＝L_i＜＝μ+2×σ；如果不满足，则如果L_i＜μ-2×σ，令L_i＝0；如果 L_i＞μ+2×σ，令L_i＝1；如果满足，则对L_i进行转换计算，令L_i＝L_i-(μ-2×σ)／4×σ，所 得到的L_i为最终映射得到的映射距离值。通过上述流程，即将一个团块内的所有 对象分别映射到[0，1]之间。

对于本发明提供的基于无监督与傅里叶核函数的遥感影像像元纯度识别方 法，计算得到的每一个像元的纯度值在[0，1]之间，其中0表示遥感影像中相对最 纯的像元，1表示遥感影像中相对最不纯的像元，通常为处于最边界位置的混合 像元。也就是说，纯度值越大，代表像元纯度最差；而纯度值越小，代表像元 纯度最好。发明人已通过大量实验，证明了本发明的有效性。

下面介绍一种应用本发明提供的基于无监督与傅里叶核函数的遥感影像像 元纯度识别方法的具体实施例：

采用美国Landsat TM-7卫星拍摄我国东北地区一处湿地的遥感影像，该遥 感影像包含Blue、Green、Red、Near IR、Mid IR、Thermal和Far IR共计7个 波段，通过这7个波段可以从多种空间属性角度描述地面的覆盖情况。如图6 所示，为拍摄得到的遥感影像图；在图6中，局部影像1、局部影像2、局部影 像3、局部影像4、局部影像5、局部影像6分别为区域A对应的7个波段影像 数据。

采用本发明提供的方法，对图6所示遥感影像进行分割和聚类操作后，共 产生5个团块，如图7-11所示，为所产生的5个团块的影像图。对于图10所示 团块，将典型区域B放大，分别计算像元A1、B1、C1、D1、E1的映射距离值； 对于图11所示团块，将典型区域C放大，分别计算像元A2、B2、C2、D2、E2 的映射距离值；结果见下表所示。

  A1 B1 C1 D1 E1 映射距离值 0.13 0.27 0.32 0.47 0.87   A2 B2 C2 D2 E2 映射距离值 0.30 0.38 0.54 0.85 0.93

对于图10所示团块，通过实地考察发现，其土地覆盖类型为草地，其典型 区域的5个代表性像元实际情况情况如下：

像元A1：属于该地区相对最纯的草地地区，适合做定量遥感环境监测样本 采集地。

像元B1、C1和D1，与像元A1相比，混杂了一些其它类型地物，纯度均 低于像元A1，不适合做样本。

像元E1，属于该地区相对最不纯的草地地区，属于刚从湿地蜕变为草地的 地区，适合作为土地利用类型自动分类的样本采集地。

可见，对于图10所示地区，像元A1、B1、C1、D1、E1计算得到的映射距 离值所反应的纯度情况，与实际考察的纯度情况相符，可见，本发明计算得到 的像元映射距离值，定量反应的像元的实际纯度情况。

对于图11所示团块，通过实地考察发现，其土地覆盖类型为湿地，其典型 区域的5个代表性像元情况如下：

像元A2：属于该区域相对最纯的湿地，在这个区域中比较适合做定量遥感 环境监测样本采集地，同时也反映出湿地混杂草地的情况已经比较严重，所以 该区域没有更接近于0的像元了。

像元B2、C2和D2：与像元A2相比，混杂了较多的其他地物类型，这些 位置不适合做样本。

像元E2：属于该地区相对最不纯的地区，该位置正在经历从湿地到草地的 退化过程，适合作为土地利用类型自动分类的样本采集地。

可见，对于图10所示地区，像元A2、B2、C2、D2、E2计算得到的映射距 离值所反应的纯度情况，与实际考察的纯度情况相符，可见，本发明计算得到 的像元映射距离值，定量反应的像元的实际纯度情况。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通 技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰， 这些改进和润饰也应视本发明的保护范围。