CN112733888B - 一种基于多层mrf聚类的深空探测影像形貌分析方法 - Google Patents

Abstract

本发明提出了一种基于多层MRF聚类的深空探测影像形貌分析方法，其步骤如下：深空探测影像数据进行光谱直方图表示，确定初始聚类的类别数；根据初始聚类的类别数，采用像素级MRF模型对影像数据进行聚类获得初始聚类结果图；对初始聚类的每一种类别的像素分别进行光谱直方图表示，确定每一个初始聚类类别的细化类别数目，采用修正MRF模型对每一个初始类别内的像素点进行聚类获得细化聚类结果图；将初始聚类结果和细化聚类结果叠加形成多层聚类图，获得深空探测影像数据中的形貌分析结果。本发明以多层聚类结果来表示深空探测光学数据，可更好地凸显光照过亮或光照不足区域的影像细节信息，提高深空影像中行星或小行星形貌分析的精度。

Description

一种基于多层MRF聚类的深空探测影像形貌分析方法

技术领域

本发明涉及深空探测影像智能判读的技术领域，尤其涉及一种基于多层MRF聚类的深空探测影像形貌分析方法。

背景技术

从古时夜观天象，到当今的“嫦娥”飞天，人类一直没有停止过对深空的探测步伐。近年来，随着全球各国相继开展的深空探测任务，人类能够获得的深空探测数据呈指数级增长。然而，与深空数据获取能力形成鲜明对比的是深空数据的解译判读能力相对低下。特别是深空环境的复杂多样，进一步加剧了数据判读的难度，尤其是深空探测影像数据，在成像过程中往往存在光照过亮或光照不足的区域，极大限制了数据的解译精度与信息的充分提取。因此，深空探测影像数据的自动化智能分析与判读是当前研究的一个重要科学问题。

为了有效地分析和判读深空探测影像中行星或小行星的形貌特点，一些方法被国内外学者相继提出。但是，这些方法的研究往往关注某一种典型形貌，如撞击坑、砾石等。然后，通过分析典型形貌的影像特点，如形状、中间峰、光照变化等，在影像数据中识别出这一典型形貌。单一形貌特点的识别方法难以适用于整幅影像中行星或小行星的形貌识别。为了有效地识别整幅影像的形貌特点，不同的聚类方法被相继引入到深空探测影像的形貌分析，如早期的阈值法、OSTU方法、K均值、FCM(基于划分的聚类算法)等，但是，由于深空探测影像中光照不均、行星或小行星表面光谱变化不大等特点，导致这些聚类结果难以有效地判读行星或小行星的形貌。MRF模型是一种概率图模型，它在经典聚类方法的基础上，不仅考虑了影像的特征，而且考虑了影像相邻像素间的空间关系，能够更好地反映影像中的形貌特点。但是，单一层面的MRF模型聚类方法仍然难以充分地识别深空探测影像中的形貌信息，特别是光照过亮或光照不足的区域。由于这些区域内不同形貌单元间像素的光谱差异非常小，不仅人眼难以判读识别，MRF模型也往往会将其聚为一类。深空探测影像复杂的成像环境及日益增多的数据量，为人工判读带来了挑战。

为了提高光照过亮或光照不足区域的识别精度，需要对这些部分进行影像增强。但是，现有的影像增强技术，如直方图均衡化、对比度拉伸、对数变换等，对深空探测数据光照过亮或不足区域的增强效果非常有效，难以满足后续聚类和识别分析的需求。因此，亟须一种方法能有效地增强深空探测影像光照过亮或不足区域，并基于此进行影像的聚类识别和形貌的智能分析。

发明内容

针对现有深空探测影像形貌判读时，深空复杂成像环境导致影像形貌识别精度有限，尤其是光照过亮或光照不足区域形貌识别精度不高的技术问题，本发明提出一种基于多层MRF聚类的深空探测影像形貌分析方法，通过MRF模型对深空探测数据进行两层次的聚类分析，将深空探测影像的增强和聚类过程合二为一，从多层聚类图的视角增强和凸显了影像中的形貌规律，尤其是光照过亮或光照不足区域的影像形貌信息，以提高深空探测影像的形貌解译精度。

为了达到上述目的，本发明的技术方案是这样实现的：一种基于多层MRF聚类的深空探测影像形貌分析方法，其步骤如下：

步骤一：输入并读取一幅深空探测影像数据，进行光谱直方图的表示，并根据直方图分布确定初始聚类的类别数；根据初始聚类的类别数，采用条件迭代模式下的像素级MRF模型对影像数据进行聚类，获得初始聚类结果图；

步骤二：对初始聚类的每一种类别分别像素进行对应的光谱直方图表示，根据光谱直方图分布确定每一个初始聚类类别的细化类别数目，并采用修正MRF模型对每一个初始类别内的像素点进行第二层聚类，获得细化聚类结果图；

步骤三：将步骤一的初始聚类结果和步骤二的细化聚类结果叠加形成多层聚类图，从聚类的角度进行图像增强，分析行星或小行星地形地貌在多层聚类图中的特点，获得深空探测影像数据中的形貌分析结果。

所述步骤一中深空探测影像数据I的大小为M×N，光谱波段数为p，p∈N⁺，N⁺为正整数集；若p＝1，即影像数据I为全色单波段，直接对影像数据I的像素进行直方图表示；若p＞1，将多波段的影像数据I转化为单个波段的影像数据，然后再对单波段的影像数据的像素进行直方图表示。

所述步骤一根据直方图分布确定初始聚类的类别数和步骤二中确定每一个初始聚类类别的细化类别数目的方法为：在相应光谱直方图上，以类内方差尽可能小、类间方差尽可能大为准则，通过分析直方图波峰波谷实现的；分析局部的波峰和波谷以及两侧尾的轻重，将每一个单峰或每一条重尾作为一种类别，计数获得初始聚类的类别数n。

所述步骤一中的像素级MRF模型和步骤二中的修正MRF模型的聚类方法是：分别是对原始的深空探测影像数据和某一初始聚类的类别中的像素进行，都采用了正态分布描述MRF模型的一阶势函数，实现对特征的建模；都采用多层逻辑模型描述MRF模型二阶势函数中能量函数内两点间的相互作用，实现对空间关系的建模；通过逐像素的生成式概率迭代方法求解MRF模型，得到对应的聚类结果。

所述步骤一中采用像素级MRF模型对影像数据进行聚类获得初始聚类结果图的方法为：

影像数据I＝{I_s|s∈S}，其中，s＝(s_x,s_y)表示影像数据I的第s_x行、第s_y列的位置点，且1≤s_x≤M,1≤s_y≤N，S是所有位置点的集合，I_s是位置点s处的p维光谱特征；即像素级MRF模型在集合S上定义一个表示各像素点类别的标记场X＝{X_s|s∈S}，其中，随机变量X_s表示位置点s的聚类类别，随机变量X_s取值于类别集合Λ＝{1,2,…,n}；

标记场X服从空间马氏性，即：

P(X_s＝x_s|X_t＝x_t,t∈S/s)＝P(X_s＝x_s|X_t＝x_t,t∈N_s)；

其中，x_s和x_t分别表示位置点s和t处随机变量X_s和随机变量X_t的具体类别取值，S/s表示位置点的集合S剔除位置点s后剩余位置点的集合，N_s表示与位置点s空间相邻的位置点所组成的空间邻域集合，P(·)为概率密度函数；

根据Harmmersley-Clifford定理，具备马氏性的标记场联合概率p(x)服从Gibbs分布，即：

其中，能量函数U(x)＝∑_s∈SU(x_s)，每个位置点的能量函数U(x_s)由其不同势团的势函数求和得到，x为标记场X的一个实现；采用二阶势团的势函数和作为能量函数U(x_s)，即：

其中，势函数：

其中，β为势函数参数；

在标记场给定的条件下，设每个位置点的光谱特征分布是相互独立，且服从正态分布，即：

其中，μ_l和∑_l是利用EM算法得到的第l类正态分布的均值和方差估计值，det表示矩阵的行列式，上标T表示矩阵的转置；

根据p(x)和p(I_s|x_s＝l)，位置点s在给定光谱特征条件下的后验概率为：

其中，p(I_s)是位置点s处的光谱特征所对应的概率，p(x_s＝l)表示位置s处的标记为l类的概率；

根据最大后验概率准则，像素级MRF模型在位置点s的聚类结果

为：

其中，argmax函数的结果是使得目标函数取得最大值的点集；

利用K均值初始化标记场X的初始标记结果，逐像素的迭代求解每一个位置点的最优聚类结果

当迭代收敛时，将各个位置点最优结果的整体作为像素级MRF模型的聚类结果，初始聚类结果图

所述步骤二中根据初始聚类的各类别的光谱直方图确定相应初始聚类类别的细化类别数目的方法为：定义初始聚类类别为l的所有像素的位置点集合

l∈Λ＝{1,2,…,n}；对相应像素点集合I^l＝{I_s|s∈S^l}将各像素的光谱特征转化为单波段数据，并绘制相应的单波段光谱直方图；根据相应的光谱直方图，分析直方图分布特点即多峰的曲线变化，确定初始类别为第l类的第二层细化聚类类别数k_l,l∈Λ＝{1,2,…,n}。

所述步骤二中采用修正MRF模型对每一个初始类别内的像素点进行聚类获得细化聚类结果图的方法为：

修正MRF模型在包括所有位置点的集合S上定义一个类别标记场

其中，若s∈S/S^l时，令

若s∈S^l时，

表示位置点s细化聚类类别的随机变量，且

取值于细化类别集合Λ^l＝{1,2,…,k_l}，且S/S^l表示位置集合S剔除子集S^l后剩余位置点的集合；

在修正MRF模型中，对于位置点s∈S^l，设类别标记场

服从空间马氏性，即：

其中，

和

表示位置点s和t处的细化类别标记，即

描述能量函数的二阶势团的势函数

为：

势函数

的能量函数

且标记场的联合概率密度函数

每个位置点的光谱特征在细化类别给定的条件下服从正态分布，即：

其中，

为位置点s的光谱特征的先验概率，μ_h和∑_h是利用EM算法得到的正态分布第h类的均值和方差估计值，h∈Λ^l；

根据概率

位置点s在细化类别为

的后验概率为：

根据最大后验概率准则，修正MRF模型在位置点s的细化聚类结果

为：

根据逐像素的迭代求解初始聚类类别为l的位置点集合S^l中每一个位置点的最优聚类结果

当迭代收敛时，得到初始类别l的细化聚类结果

将所有初始类别的细化聚类结果合并显示在同一张结果图上得到第二层聚类结果，即细化聚类结果图

所述步骤三中根据多层聚类结果进行形貌分析的方法为：将初步聚类结果图和细化聚类结果图叠加合并，形成多层聚类结果图

从聚类的角度实现深空探测影像的增强；在初步聚类结果图中，根据不同类别的空间依存关系及空间分布特点，分析行星或小行星的宏观形貌特点，包括但不限于盆地、峡谷、高原、大型撞击坑、行星特有宏观形貌；在细化聚类结果图中，根据同一初始聚类类别内的细化聚类结果及各细化类别的空间分布与细化类别相互依存关系，分析行星或小行星的微观形貌特点，包括但不限于中小型撞击坑、砾石；根据不同形貌在多层聚类结果图中的特点，识别、判读和分析影像数据中对应的形貌单元。

本发明的有益效果：首先，根据深空影像数据的直方图，确定初始层次的聚类类别，并以此进行像素级马尔科夫随机场模型(Markov random field model，MRF)的影像初次聚类。然后，针对初始聚类中每一类别的影像数据，进一步分析其特征的直方图分布特点，并确定每一类别的细分类别数目，利用MRF模型对每一初始类别分别进行差异化的第二层次聚类。最后，叠加两层聚类结果形成向量化的多层聚类图，并以其为基础进行影像中行星或小行星形貌特征的规律识别与智能判读。

本发明通过初始聚类分析，可以更清晰地刻画深空探测影像中行星或小行星的宏观形貌特点；通过细化聚类分析，可以凸显影像中宏观形貌的细节特点，尤其是光照过亮或光照不足区域的形貌细节信息。通过形成多层次的聚类图，不仅可以从聚类的角度实现深空探测影像的增强效果，而且还可以根据其类别的空间分布规律进行形貌的判读与分析，识别或发现新的形貌信息，为深空探测提供科学支持。本发明以多层聚类结果表示深空探测光学数据，相较于传统的图像增强算法，可更好地凸显光照过亮或光照不足区域的影像细节信息，提高深空影像中行星或小行星形貌分析的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程示意图。

图2为本发明在一幅月球勘测轨道器相机获取的月球BELL E撞击坑LROC(LunarReconnaissance Orbiter Camera)数据的具体实施实例图；其中，(a)表示获取的月球BELL撞击坑的LROC深刻探测影像数据，(b)表示(a)的直方图。

图3为本发明获得的初始聚类结果图，其中初始类别数等于3。

图4为本发明初始类别图中不同类别对应的直方图，其中，(a)表示初始类别为第一类的像素直方图，(b)表示初始类别为第二类的像素直方图，(c)表示初始类别为第三类的像素直方图。

图5为本发明各初始类别细化聚类结果及第二层整体的聚类结果图，其中，(a)表示初始类别为第一类的细化聚类结果，(b)表示初始类别为第二类的细化聚类结果，(c)表示初始类别为第三类的细化聚类结果，(d)表示综合各初始类别细化聚类结果的第二层聚类结果图。

图6为本发明在影像增强效果方面的对比图及形貌识别说明图，其中，(a)表示采用了均衡化直方图的影像增强结果图，(b)表示采用了对比度拉伸的影像增强结果图，(c)表示采用了对数变换的影像增强结果图，(d)表示采用了本发明多层聚类结果的影像增强表示图及各形貌单元在该图中的形貌分布规律。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明提出了一种基于多层MRF聚类的深空探测影像形貌分析方法，运用MRF模型对深空探测影像数据进行多层次聚类分析，从聚类的角度增强影像信息，尤其是光照过亮或光照不足区域的信息增强，并通过多层聚类结果同时实现行星或小行星的形貌识别与分析。本实施例选取了一幅于2020年10月2日发布的月球勘测轨道器相机获取的月球BELL E撞击坑的LROC影像数据，如图2(a)所示，以下是具体的建模步骤。

步骤一：输入并读取一幅深空探测影像数据，进行光谱直方图的表示，并根据直方图分布特点确定初始聚类的类别数；根据确定的类别数，采用条件迭代模式下的像素级MRF模型对影像数据进行聚类，获得初始聚类结果图。

根据直方图确定聚类类别数，并采用MRF模型进行聚类分析的方法为：

步骤1.1：根据直方图分布确定初始聚类的类别数

深空探测影像数据I，其大小为M×N，即该影像数据有M行N列共M·N个像素，其光谱波段数为p，p∈N⁺。若p＝1，即影像为全色单波段时，直接对影像数据的像素进行直方图的表示；若p＞1，则首先将多波段的影像数据I转化为单个波段的影像数据，然后再对单波段影像数据的像素进行直方图的表示。

对于可视化的直方图，分析其分布特点，尤其是局部的波峰和波谷以及两侧尾的轻重，将每一个单峰或每一条重尾(两侧尾的轻重因此有重尾)作为一种类别，计数获得初步的聚类类别数目n。

步骤一和二中的聚类类别数的确定，都是在相应光谱直方图的基础上，以类内方差尽可能小、类间方差尽可能大为准则，通过分析直方图波峰波谷实现的。由于LROC数据是全色单波段数据，不需要波段特征降维，所以首先可直接绘制该数据的直方图，如图2(b)所示。观测该直方图，可以发现其表现出明显的三峰两谷的特点，即在0和255这两个端点值附近有两个峰值，其中0附近的峰值较为明显，255的峰值基本集中在255这一个值点，而在中间100附近有一个峰值，而两个峰谷则处在三个峰之间。因此，初始的聚类类别数目确定为3类。

步骤1.2：根据类别数n，利用像素级MRF模型对LROC数据进行聚类

当确定初始聚类的类别数n＝3之后，采用像素级MRF模型对影像I＝{I_s|s∈S}进行聚类，其中，s＝(s_x,s_y)表示影像数据的第s_x行、第s_y列的位置点，且1≤s_x≤M,1≤s_y≤N，M和N表示图像的大小为M×N，本实施例中1≤s_x≤5000,1≤s_y≤5000，M＝N＝5000，S是所有位置点的集合，I_s是位置点s处的p维光谱特征。具体而言，像素级MRF模型在集合S上定义一个表示各像素点类别的标记场X＝{X_s|s∈S}，其中，随机变量X_s表示位置点s的聚类类别，其取值于类别集合Λ＝{1,2,…,n}。在像素级MRF模型中，假设标记场X服从空间马氏性，即：

P(X_s＝x_s|X_t＝x_t,t∈S/s)＝P(X_s＝x_s|X_t＝x_t,t∈N_s)。

其中，x_s和x_t分别表示位置点s和t处随机变量X_s和随机变量X_t的具体类别取值，S/s表示位置点的集合S剔除位置点s后剩余位置点的集合，N_s表示与位置点s空间相邻的位置点所组成的空间邻域集合，它在像素级MRF模型中通常选取为四邻域或八邻域，P(·)为概率密度函数。马氏性表明，在像素级MRF模型中每一个像素点类别的取值仅受到其空间邻域集合N_s影响，而与其他像素点无关。根据Harmmersley-Clifford定理，具备马氏性的标记场联合概率p(x)服从Gibbs分布，即：

其中，能量函数U(x)＝∑_s∈SU(x_s)，每个位置点的能量函数U(x_s)由其不同势团的势函数求和得到，x为标记场X的一个实现。在本发明实施例中，采用二阶势团的势函数和作为能量函数U(x_s)，即：

其中，势函数：

其中，β为势函数参数，本实例取值为0.5，在大多数相关应用中，β一般取值于[0.1,5]。然后，在标记场给定的条件下，假设每个位置点的光谱特征分布是相互独立，且服从正态分布的，即：

在本实例中，图像的光谱波段p＝1，即本实例中图像为单波段图像，μ_l和∑_l是利用EM算法得到的第l类正态分布的均值和方差估计值，det表示矩阵的行列式，上标T表示矩阵的转置。根据p(x)和p(I_s|x_s＝l)，某一位置点s在给定光谱特征条件下的后验概率为：

其中，p(I_s)是位置点s处的光谱特征所对应的概率，p(x_s＝l)表示位置s处的标记为l的概率。

根据最大后验概率准则，像素级MRF模型在位置点s的聚类结果

为：

其中，argmax函数的结果是使得目标函数取得最大值的点集。

利用K均值初始化标记场X的初始标记结果，然后根据上式逐像素的迭代求解每一个位置点的最优聚类结果

当迭代收敛时，将各个位置点最优结果的整体作为像素级MRF模型的聚类结果，初始聚类结果图

如图3所示。

根据初始聚类结果图

可以计算得到这3种初始类别的均值与方差，如表1所示。

表1初始聚类结果

初始聚类结果	初始类别1	初始类别2	初始类别3
				均值μ<sub>l</sub>	19.2304	99.4147	189.2198
方差∑<sub>l</sub>	18.2903<sup>2</sup>	23.7348<sup>2</sup>	35.2191<sup>2</sup>

由图3所示，这三种初始类别结果在原始影像图中，分别对应光照不足的阴影区域、光照条件适度区域、光照过亮区域。

步骤二：对初始聚类的每一种类别分别像素进行对应的光谱直方图表示，根据光谱直方图特点确定每一个初始类别的细化类别数目，并采用修正MRF模型对每一个初始类别内的像素点进行第二层聚类，获得细化聚类结果图。

根据初始各类别的光谱直方图确定不同类别的细化聚类数，并采用修正MRF模型分别进行聚类分析的方法：

步骤2.1：根据初始聚类的各类别的光谱直方图，确定相应初始类别的细化类别数目；

根据初始聚类的结果，定义初始聚类类别为l的所有像素的位置点集合

l∈Λ＝{1,2,…,n}。然后，对相应像素点集合I^l＝{I_s|s∈S^l}，采用步骤1.1的方法，将各像素的光谱特征转化为单波段数据，并绘制相应的单波段光谱直方图。最后，根据相应的光谱直方图，分析其分布特点，尤其是多峰的曲线变化，确定初始类别为第l类的第二层细化聚类类别数k_l,l∈Λ＝{1,2,…,n}。

绘制相应像素集合的直方图，如图4所示。其中，第一类表示光照不足区域的类别其光谱取值位于区间[0,60]，其直方图分布特点包括了0附近的一个峰值、20附近的一个峰谷和之后逐步上升的趋势，如图4(a)所示，但是考虑到影像中BELL E撞击坑内存在大量光照不足的阴影区域，且其光谱值都分布在0附近的峰值，而这些值因为较小的光谱值差异而导致人类的目视判读几乎无法获得该区域任何有价值的形貌信息。因此，确定初始类别为第1类的第二层细化聚类类别数k₁＝15，以凸显光照不足区域的形貌细节信息。

对于第二类初始聚类结果，其对应的是光照条件较为适度的区域，直方图取值于100附近，如图4(b)所示，呈现出标准的单峰值分布特点。考虑到该区域数据在原始影像中有较好的可判读性，因此，确定初始类别为第2类的第二层细化聚类类别数k₂＝3。

对于第三类初始聚类结果，其直方图主要取值于[150,255]，如图4(c)所示，其主要对应光照充足或光照过亮的影像区域。类似于第一类初始结果区域，在光照过亮区域中，如BELL E撞击坑的迎光面区域，一些形貌细节信息难以直接判读，为了凸显这部分区域的细节特点，确定初始类别为第3类的第二层细化聚类类别数k₃＝3。

步骤2.2：根据各初始类别的细化聚类类别数k_l，对n个初始类别分别进行修正MRF模型的细化聚类分析。

当初始类别l＝3的细化聚类类别k_l分别确定之后，l＝1,2,3，采用一种修正的MRF模型对其像素点集合I^l进行第二层的细化聚类。具体而言，修正MRF模型仍然在包括所有位置点的集合S上定义一个类别标记场

其中，若s∈S/S^l时，令

若s∈S^l时，

表示位置点s细化聚类类别的随机变量，且其取值于细化类别集合Λ^l＝{1,2,…,k_l},l＝1,2,3。S/S^l表示位置集合S剔除子集S^l后剩余位置点的集合。

在修正MRF模型中，对于s∈S^l，仍假设类别标记场

服从空间马氏性，即：

其中，

和

表示位置点s和t处的细化类别标记，即

但是，其描述能量函数的二阶势团的势函数

定义为：

在本实施例中，为了更好地保留第一层聚类结果的影响，第二层MRF聚类的势函数参数β都取值为1。势函数

的能量函数

且标记场的联合概率函数

同时，每个位置点的光谱特征在细化类别给定的条件下服从正态分布，即

其中，

为位置点s处的光谱特征的先验概率。在本实施例中，p＝1,μ_h和∑_h是利用EM算法得到的正态分布第h类的均值和方差估计值，h∈Λ^l。根据概率

和

位置点s在细化类别为

的后验概率为：

根据最大后验概率准则，修正MRF模型在位置点s的细化聚类结果

为

根据上式逐像素的迭代求解3个初始聚类类别为l的位置点集合S^l中每一个位置点的最优聚类结果

当迭代收敛时，得到初始类别l的细化聚类结果

如图5(a)到5(c)所示。将所有初始类别的细化聚类结果合并显示在同一张结果图上得到第二层聚类结果，即细化聚类结果图

如图5(d)所示。

根据第二层聚类结果的细化聚类结果图

可以计算得到3种初始类别的细化类别的均值与方差，如下面表2-4所示。

表2初始类别1的细化类别表

表3初始类别1的细化类别表

表4初始类别3的细化类别表

步骤一和二中的MRF聚类方法，分别是对原始影像数据和某一初始聚类类别中的像素进行，都采用了正态分布描述MRF模型的一阶势函数，实现了对特征的建模，都采用了多层逻辑模型描述MRF模型二阶势函数中能量函数内两点间的相互作用，实现了对空间关系的建模。最后，通过逐像素的生成式概率迭代方法求解模型，得到对应的聚类结果。

步骤三：将初始聚类结果和细化聚类结果叠加形成多层聚类图，从聚类的角度进行图像增强，并以此为基础，分析行星或小行星典型地形地貌在多层聚类图中的特点，获得最终的形貌分析结果。

对于获取的初始和细化聚类结果，将其叠加形成多层聚类图，并归纳行星或小行星的典型地貌(如撞击坑、砾石等)对应的多层聚类结果的特点，分析影像中行星或小行星的形貌特点。根据多层聚类结果进行形貌分析的方法的具体实施过程如下：

首先，将初步聚类结果图和细化聚类结果图叠加合并，形成多层聚类结果图

从聚类的角度实现深空探测影像的增强，如图6(d)所示。其次，在初步聚类结果图中，根据不同类别的空间依存关系及空间分布特点，分析行星或小行星的宏观形貌特点，包括但不限于盆地、峡谷、高原、大型撞击坑、行星特有宏观形貌等；在细化聚类结果图中，根据同一初始类别内的细化聚类结果及各细化类别空间分布与细化类别相互依存关系，分析行星或小行星的微观形貌特点，包括但不限于中小型撞击坑、砾石等。最后，根据不同形貌在多层聚类结果图中的特点，识别、判读、分析影像中对应的形貌单元。针对深空探测获取的遥感影像的特点，本发明提出了一种以MRF模型分类为基础的、以两层聚类图为结果的新型影像表示方法，从多层聚类角度实现了遥感影像的增强，可有效凸显深空遥感影像数据在光照过亮或光照不足区域的形貌特点，提高影像判读精度。

为了验证本发明在影像增强方面的有效性，与三种经典的影像增加方法进行了对比，它们分别是均衡化直方图、对比度拉伸和对数变换，相应的影像增强结果在图6(a)到图6(c)中进行了展示。从图6中可以看出，传统的影像增强方法都是在基本保持影像整体特征不变的前提下，通过适当的函数变换实现影像增强。这些方法虽然在一般影像增强中具有良好的表现，但是在深空探测影像中，由于成像条件的复杂，导致一些局部光照不足或光照过强区域的增强效果并不理想。例如，在BELL E的两个坑唇部分，既撞击坑环形的凸起部分，存在着左侧光照过强区域和右侧光照不足区域，如图6(d)两个白色方框区域，尤其是右侧光照不足区域，由于光谱值的变化非常细微，导致肉眼难以识别该区域的形貌细节。在三种经典图像增强方法的结果中，图6(b)的对比度拉伸结果完全无法增强右侧光照不足区域；图6(a)的均衡化直方图结果只能轻微地增强内外坑唇间阴影区域，大致增强出该阴影区域的两部分；图6-(c)的对数变换结果能较好地凸显右侧内外坑唇间阴影区域的两个大的组成部分，但是却同时加剧了左侧内外坑唇间区域的光照亮度，使得这部分区域的影像更难以判读。相比于这些经典的影像增强方法，本发明从聚类结果图的角度对原始数据进行了表示，利用多层次聚类进行影像的增强。例如，在图6(d)中，前述内外坑唇间的光照过亮和不足区域(即白色方框区域)在多层聚类结果的表示下，更理想地展现了对应的形貌细节，如右侧内外坑唇区域间的阴影区域，不仅可以分为左右两个明显不同的区域，而且左侧区域内的细化类别变化也展示出这部分区域的形貌变化不是平坦而连续的，而右侧区域由单一细化类别表示则意味着这部分区域应该是一个陡峭的内坑沿。这些形貌细节信息的增强表明，相比于经典的影像增强方法，本发明利用多层聚类的思路可以更有效地增强深空探测影像，凸显其光照不足或光照过亮区域的形貌细节。

此外，多层次聚类图不仅可以从聚类的角度增强影像，而且还可以基于该图进行形貌的具体分析。例如，在初步聚类结果图中，根据不同类别的空间依存关系及空间分布特点，分析行星或小行星的宏观形貌特点，包括但不限于盆地、峡谷、高原、大型撞击坑、行星特有宏观形貌等。在本实施例中，由于撞击坑存在迎光面和背光面，所以其往往表现为左侧迎光面为初始聚类第3类、右侧背光面为初始聚类第1类，空间相邻的两者往往构成撞击坑，特别是本实施例中，BELL E存在内外两个同心坑，宏观形貌表现出类似“甜甜圈”的形貌特点，图6(d)左侧图中的第1、3类聚类结果很好地表现出了这一形貌规律。

在细化聚类结果图中，根据同一初始类别内的细化聚类结果及各细化类别空间分布与细化类别相互依存关系，分析行星或小行星的微观形貌特点，包括但不限于中小型撞击坑、砾石等。在本实施例中，BELL E撞击坑的内部存在着大量的中小型撞击坑，如图2(a)所示。第二层细化的聚类结果就可以对这些细节形貌进行识别。具体而言，对于BELL E撞击坑内的中小型撞击坑，为了方便起见，称里面相对直径较大的为“大型撞击坑”、直径较小的为“小型撞击坑”，如图6(d)的右侧结果所示。在细化聚类结果中，这些“大型撞击坑”的形貌往往表现为左侧由若干初始类别为第3类的细化类别构成，代表迎光面；右侧由若干初始类别为第1类的细化类别构成，代表背光面。而且由于这些撞击坑相对较大，在影像中光照的变化更多样，所以其包含的细化类别数目就相对较多，在本实施例中一般不少于6类细化类别。相应地，“小型撞击坑”的形貌虽然也表现为左侧迎光面由第3类细化类别构成、右侧背光面由第1类细化类别构成，但是由于这些撞击坑直径较小、坑的深度较浅，所以在影像中光照的变化就没有“大型撞击坑”那么多，因此其包含的细化类别数目相对较少，在本实施例中一般只有2到4类细化类别构成。根据这些形貌在聚类图中表现出的规律特点，可以多层聚类图中识别、判读和分析影像中各类型的形貌单元。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于多层MRF聚类的深空探测影像形貌分析方法，其特征在于，其步骤如下：

步骤一：输入并读取一幅深空探测影像数据，进行光谱直方图的表示，并根据直方图分布确定初始聚类的类别数；根据初始聚类的类别数，采用条件迭代模式下的像素级MRF模型对影像数据进行聚类，获得初始聚类结果图；

所述步骤一中采用像素级MRF模型对影像数据进行聚类获得初始聚类结果图的方法为：

影像数据I＝{I_s|s∈S}，其中，s＝(s_x,s_y)表示影像数据I的第s_x行、第s_y列的位置点，且1≤s_x≤M，1≤s_y≤N，S是所有位置点的集合，I_s是位置点s处的p维光谱特征；即像素级MRF模型在集合S上定义一个表示各像素点类别的标记场X＝{X_s|s∈S}，其中，随机变量X_s表示位置点s的聚类类别，随机变量X_s取值于类别集合Λ＝{1,2,…,n}；深空探测影像数据的大小为M×N；

标记场X服从空间马氏性，即：

P(X_s＝x_s|X_t＝x_t，t∈S/s)＝P(X_s＝x_s|X_t＝x_t，t∈N_s)；

其中，x_s和x_t分别表示位置点s和t处随机变量X_s和随机变量X_t的具体类别取值，S/s表示位置点的集合S剔除位置点s后剩余位置点的集合，N_s表示与位置点s空间相邻的位置点所组成的空间邻域集合，P(·)为概率密度函数；

根据Harmmersley-Clifford定理，具备马氏性的标记场联合概率p(x)服从Gibbs分布，即：

其中，能量函数U(x)＝∑_s∈SU(x_s)，每个位置点的能量函数U(x_s)由其不同势团的势函数求和得到，x为标记场X的一个实现；采用二阶势团的势函数和作为能量函数U(x_s)，即：

其中，势函数：

其中，β为势函数参数；

在标记场给定的条件下，设每个位置点的光谱特征分布是相互独立，且服从正态分布，即：

其中，μ_l和∑_l是利用EM算法得到的第l类正态分布的均值和方差估计值，det表示矩阵的行列式，上标T表示矩阵的转置；

根据p(x)和p(I_s|x_s＝l)，位置点s在给定光谱特征条件下的后验概率为：

其中，p(I_s)是位置点s处的光谱特征所对应的概率，p(x_s＝l)表示位置s处的标记为l类的概率；

根据最大后验概率准则，像素级MRF模型在位置点s的聚类结果

为：

其中，argmax函数的结果是使得目标函数取得最大值的点集；

利用K均值初始化标记场X的初始标记结果，逐像素的迭代求解每一个位置点的最优聚类结果

当迭代收敛时，将各个位置点最优结果的整体作为像素级MRF模型的聚类结果，初始聚类结果图

步骤二：对初始聚类的每一种类别分别像素进行对应的光谱直方图表示，根据光谱直方图分布确定每一个初始聚类类别的细化类别数目，并采用修正MRF模型对每一个初始类别内的像素点进行第二层聚类，获得细化聚类结果图；

所述步骤二中采用修正MRF模型对每一个初始类别内的像素点进行聚类获得细化聚类结果图的方法为：

修正MRF模型在包括所有位置点的集合S上定义一个类别标记场

其中，若s∈S/S^l时，令

若s∈S^l时，

表示位置点s细化聚类类别的随机变量，且

取值于细化类别集合Λ^l＝{1,2,…,k_l}，且S/S^l表示位置集合S剔除子集S^l后剩余位置点的集合；

为初始聚类类别为l的所有像素的位置点集合；k_l,l∈Λ＝{1,2,…,n}为初始类别为第l类的第二层细化聚类类别数；

在修正MRF模型中，对于位置点s∈S^l，设类别标记场

服从空间马氏性，即：

其中，

和

表示位置点s和t处的细化类别标记，即

描述能量函数的二阶势团的势函数

为：

势函数

的能量函数

且标记场的联合概率密度函数

每个位置点的光谱特征在细化类别给定的条件下服从正态分布，即：

其中，

为位置点s的光谱特征的先验概率，μ_h和∑_h是利用EM算法得到的正态分布第h类的均值和方差估计值，h∈Λ^l；像素点集合I^l＝{I_s|s∈S^l}；

根据概率

和

位置点s在细化类别为

的后验概率为：

根据最大后验概率准则，修正MRF模型在位置点s的细化聚类结果

为：

根据逐像素的迭代求解初始聚类类别为l的位置点集合S^l中每一个位置点的最优聚类结果

当迭代收敛时，得到初始类别l的细化聚类结果

将所有初始类别的细化聚类结果合并显示在同一张结果图上得到第二层聚类结果，即细化聚类结果图

步骤三：将步骤一的初始聚类结果和步骤二的细化聚类结果叠加形成多层聚类图，从聚类的角度进行图像增强，分析行星或小行星地形地貌在多层聚类图中的特点，获得深空探测影像数据中的形貌分析结果。

2.根据权利要求1所述的基于多层MRF聚类的深空探测影像形貌分析方法，其特征在于，所述步骤一中深空探测影像数据I的大小为M×N，光谱波段数为p，p∈N⁺，N⁺为正整数集；若p＝1，即影像数据I为全色单波段，直接对影像数据I的像素进行直方图表示；若p＞1，将多波段的影像数据I转化为单个波段的影像数据，然后再对单波段的影像数据的像素进行直方图表示。

3.根据权利要求1或2所述的基于多层MRF聚类的深空探测影像形貌分析方法，其特征在于，所述步骤一根据直方图分布确定初始聚类的类别数和步骤二中确定每一个初始聚类类别的细化类别数目的方法为：在相应光谱直方图上，以类内方差尽可能小、类间方差尽可能大为准则，通过分析直方图波峰波谷实现的；分析局部的波峰和波谷以及两侧尾的轻重，将每一个单峰或每一条重尾作为一种类别，计数获得初始聚类的类别数n。

4.根据权利要求3所述的基于多层MRF聚类的深空探测影像形貌分析方法，其特征在于，所述步骤一中的像素级MRF模型和步骤二中的修正MRF模型的聚类方法是：分别是对原始的深空探测影像数据和某一初始聚类的类别中的像素进行，都采用了正态分布描述MRF模型的一阶势函数，实现对特征的建模；都采用多层逻辑模型描述MRF模型二阶势函数中能量函数内两点间的相互作用，实现对空间关系的建模；通过逐像素的生成式概率迭代方法求解MRF模型，得到对应的聚类结果。

5.根据权利要求4所述的基于多层MRF聚类的深空探测影像形貌分析方法，其特征在于，所述步骤二中根据初始聚类的各类别的光谱直方图确定相应初始聚类类别的细化类别数目的方法为：定义初始聚类类别为l的所有像素的位置点集合

对相应像素点集合I^l＝{I_s|s∈S^l}将各像素的光谱特征转化为单波段数据，并绘制相应的单波段光谱直方图；根据相应的光谱直方图，分析直方图分布特点即多峰的曲线变化，确定初始类别为第l类的第二层细化聚类类别数k_l,l∈Λ＝{1,2,…,n}。

6.根据权利要求5所述的基于多层MRF聚类的深空探测影像形貌分析方法，其特征在于，所述步骤三中根据多层聚类结果进行形貌分析的方法为：将初步聚类结果图和细化聚类结果图叠加合并，形成多层聚类结果图

从聚类的角度实现深空探测影像的增强；在初步聚类结果图中，根据不同类别的空间依存关系及空间分布特点，分析行星或小行星的宏观形貌特点，包括但不限于盆地、峡谷、高原、大型撞击坑、行星特有宏观形貌；在细化聚类结果图中，根据同一初始聚类类别内的细化聚类结果及各细化类别的空间分布与细化类别相互依存关系，分析行星或小行星的微观形貌特点，包括但不限于中小型撞击坑、砾石；根据不同形貌在多层聚类结果图中的特点，识别、判读和分析影像数据中对应的形貌单元。

Images