表1 研究区的岩性类别
Table 1 Lithological classes of the study area
地层符号岩性描述 Q3-4沙质黏土和亚沙土 Q1洪积层 N2砂岩、粉砂岩、洪积层 C-P砂岩、石灰岩、泥岩 D砂岩和洪积层 S泥质岩、粉砂岩、砂岩、底砂岩、泥岩 Є-O白云岩和灰岩
北京大学学报(自然科学版) 第61卷 第3期 2025年5月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 61, No. 3 (May 2025)
doi: 10.13209/j.0479-8023.2025.013
遥感信息与图像分析技术国家级重点实验室基金(6142A010302)资助
收稿日期: 2024–04–01;
修回日期: 2024–05–06
摘要 针对现有基于多光谱–高光谱数据的岩性识别和分类研究中只利用光谱特征进行聚类, 容易产生类别混淆的问题, 利用光谱–空间特征进行多光谱图像的聚类。采取两种策略加入空间信息, 一是利用离散小波变换生成多尺度图像, 对不同尺度的图像依次进行聚类; 二是在各个尺度图像的聚类结果中, 利用马尔可夫随机场模型计算类别的空间概率, 得到类别的光谱–空间概率。综合多个尺度的光谱–空间聚类结果, 得到最终的聚类结果。利用 ASTER 多光谱数据进行岩性识别与分类, 并与基于高斯混合模型的光谱聚类结果、光谱聚类后滤波的结果和原始图像的光谱–空间聚类结果进行对比, 结果表明, 所采用的多尺度光谱–空间聚类方法可获得比上述 3 种方法更高精度的岩性分类结果, 说明综合利用光谱与空间信息进行多尺度的聚类是一种有效的岩性填图方法, 适用于较难获取地面参考样本地区的岩性填图。
关键词 光谱–空间聚类; 岩性分类; 多光谱图像
岩性识别和分类是高光谱数据和多光谱数据重要的地质应用之一[1]。国内外学者在利用不同的多光谱和高光谱数据进行岩性分类方面开展了很多研 究[2–6]。但是, 现有的研究多采用不同的监督分类方法(如随机森林和支持向量机), 很少利用非监督(聚类)方法。在一些研究中, 只将非监督分类方法(如 K-Means 和 SOM(self-organizing maps))作为对比方法或分类过程的一部分, 或用于辅助训练样本的选择, 很少单独将其作为一种岩性分类和识别的方法。Kuhn 等[7]利用地球物理、地球化学和高程数据, 采用随机森林方法进行岩性分类, 同时与 K-Means 和 SOM 两种非监督分类方法做了比较。目前, 还没有专门针对高光谱–多光谱数据的岩性分类的非监督方法研究。
在利用监督分类方法进行岩性识别时, 要求有充足的代表性训练样本, 但实际上往往难以获得, 某些地区(如偏远的、交通不便的地区)甚至没有可用的训练样本。在这样的情况下, 可将非监督分类方法用于研究区岩性的聚类, 得到初始的光谱类别, 这些光谱类别与地表的岩性类别有直接的关联, 可以辅助训练样本的选择。因此, 探索有效的非监督分类方法对高光谱–多光谱数据的岩性识别或分类具有重要作用。
目前已有一些将聚类方法用于岩性分类的研究[7–8], 但这些研究只利用光谱信息进行聚类, 导致具有相似光谱特征的岩性难以区分, 或者数据的噪声给岩性分类带来误差。相关研究表明, 加入空间信息进行聚类(或非监督分类)可以解决上述两方面的问题[9]。现有研究是在图像聚类过程的不同阶段加入空间信息: 1)在聚类前加入, 如用多层次方法估计聚类的初始参数[9], 或用尺度空间方法得到多尺度图像[10]; 2)在聚类中加入, 如通过调整相似性测度(如光谱和空间的隶属度)得到聚类结果[9–10]; 3)在聚类后加入, 即对光谱聚类结果进行滤波[10]。但是, 现有的基于聚类的岩性分类研究中没有加入空间信息的聚类方法。
鉴于上述背景, 本研究拟在图像聚类中加入空间信息, 综合利用光谱–空间信息进行聚类, 并以岩性填图为实例, 评价和验证光谱–空间聚类方法的有效性。
本文选择新疆塔里木盆地西北边缘的柯坪地区为研究区, 地理坐标范围为 40°00′—40°20′N, 78°00′ —78°10′E。研究区为干旱半干旱气候, 植被覆盖少, 岩性出露齐全, 适合利用遥感图像进行岩性识别与分类。
研究区出露的主要岩性单元是寒武纪、奥陶纪、志留纪以及泥盆纪的沉积岩, 共 7 个岩性类别(表 1)。
选用 2009 年 10 月获取的 ASTER 多光谱数据为实验数据。ASTER 共有可见光与近红外(VNIR)、短波红外(SWIR)和热红外(TIR)波长范围的 14 个多光谱波段, 本文使用其中 9 个多光谱波段数据, 包括 3 个 VNIR 波段(空间分辨率为 15m)和 6 个 SWIR波段(空间分辨率为 30m)。以 VNIR 波段为参考图像, 对 SWIR 波段进行配准, 并重采样为 15m 的像元大小。将 3 个 VNIR 波段和 6 个配准后的 SWIR 波段叠加, 合成 9 波段的多光谱图像(像元大小为 15m), 最终采用的 9 波段图像大小为 2488×1888 像元。将波段 1, 2, 3 分别作为 R, G, B, 构建研究区的假彩色影像(图 1)。
本研究参照研究区的 1:20 万地质图[11], 从每个岩性类别中随机选取 300 个样本, 共 2100 个参考样本用于岩性分类的精度计算。
针对在光谱聚类过程中光谱特征相似的像元易产生类别混淆的问题, 在聚类过程中加入空间信息, 综合利用光谱–空间信息进行聚类[9]。根据文献[9–10], 本研究采用两种策略结合的方法: 聚类前生成多尺度图像, 克服图像中的噪声; 聚类过程中根据光谱–空间概率(或隶属度)判定聚类类别。方法流程图如图 2 所示, 主要步骤概括如下。
1)利用离散小波变换对图像进行多尺度分解, 然后在不同尺度(分辨率)的图像中依次进行聚类。
2)在初始尺度(最低分辨率)的图像上, 利用高斯混合模型(Gaussian mixture model, GMM)进行光谱聚类, 得到光谱概率; 然后, 利用马尔可夫随机场(Markov random field, MRF)计算得到空间类别概率, 并结合光谱和空间概率判定像元的类别, 得到初始尺度下的图像聚类结果。
表1 研究区的岩性类别
Table 1 Lithological classes of the study area
地层符号岩性描述 Q3-4沙质黏土和亚沙土 Q1洪积层 N2砂岩、粉砂岩、洪积层 C-P砂岩、石灰岩、泥岩 D砂岩和洪积层 S泥质岩、粉砂岩、砂岩、底砂岩、泥岩 Є-O白云岩和灰岩
图1 研究区ASTER假彩色图像
Fig. 1 ASTER false color images
图2 本文方法流程
Fig. 2 Flowchart of the adopted method
3)在离散小波变换得到的更高分辨率图像上, 依次计算光谱–空间聚类概率, 得到不同尺度(分辨率)的聚类结果。
4)综合上述结果, 得到最终聚类结果。
光谱聚类是以像元的灰度值(即反射率)作为聚类依据, 每个像元的类别标记是独立的, 不考虑周边像元的空间相关性。因此, 如果只利用光谱特征标记像元类别, 由于不同的像元类别可能具有光谱相似性, 会出现类别混淆。
本研究采用高斯混合模型(GMM)作为光谱聚类方法。GMM 是一种软聚类方法, 使用类别后验概率来表达像元的类别归属。利用 GMM 进行图像聚类的基本原理是, 假定在光谱特征空间中, 图像各类别的灰度值分布均符合高斯分布, 即图像总体分布情况可由 K 个高斯模型来拟合[12]。在 GMM 中, 像元的概率密度函数 p(x)由 K 个高斯模型的线性组合表示:
(1)
其中, x 为像元的光谱值, πk 为第 k 个高斯模型的权重, N 为高斯分布, μk 和 Σk 分别为第 k 个类簇像元的均值和方差。
像元 xi 对类别 ωk 的光谱类别概率 pspect(xi|ωk)表达为
(2)
所有类别的概率和为 1, 即
1。
通常采用最大期望(expectation maximization, EM)算法对各类簇像元的光谱均值 μk 和方差 Σk 进行估计[13]。利用模拟退火(simulated annealing, SA)算法来指导 GMM 参数的估计, 得到全局最优解[14]。对于 SA 中初始温度 T0(无单位标量), 本研究采用0.5 [12]。
本研究采用两种加入空间信息的策略进行聚类: 1)基于离散小波变换(discrete wavelet transform, DWT)生成多尺度(分辨率)图像, 在不同尺度的图像中进行聚类; 2)利用马尔可夫随机场(MRF)模型度量邻域像元对中心像元的空间影响, 得到当前尺度下像元的空间类别概率。
2.2.1 利用离散小波变换(DWT)生成多尺度图像
DWT 是一种常用的信号处理方法, 可以将信号在不同尺度上进行分解和重构。它利用一组基函数的线性组合来表达原始信号中不同频率(高频与低频)的成分, 实现信号特征的提取和压缩[15]。
如图 3 所示, 在图像的频率域对原始图像(二维信号)进行离散小波变换, 常采用二进小波作为小波变换函数, 即使用 2 的整数次幂进行划分, 在一个特定的尺度可得到 4 个分量, 其中高频分量为图像的水平、垂直和对角线细节, 低频分量为降分辨率的图像。在低频分量(原始图像的近似部分)上继续分解, 不断得到 4 个分量, 进而降低原始图像的分辨率。
基于离散小波变换方法, 通过设定对原始图像进行分解的尺度数目, 即可得到多尺度(分辨率)用于聚类的图像。
本文用 MRF 模型来表达图像中相邻像元间类别(如利用 GMM 聚类得到的类别)的空间相关性。依照 MRF 的假设, 像元的空间类别概率 pspat(xi|ωk) (ωk 为类别, xi 为像元)由像元 xi 邻域内的像元 xj 决 定[17]。采用一种改进的势函数 U(x|ωk)来描述邻域N(x)的像元对中心像元的空间影响[18]:
(3)
像元 xi 对类别 ωk 的空间类别概率 pspat(xi|ωk)可表达为
(4)
其中, 正常数 β 为空间影响因子, 取值越大代表邻域像元对中心像元类别的空间影响越大; Z 为归一化常数。像元 xi 对所有类别的空间类别概率和为 1:
对于 MRF 中 β 的取值, 许多学者针对不同的图像分类任务进行了研究。Besag[19]建议 β 取值为 1.5, 但认为该参数的最优取值不是一成不变的, 在不同数据源和不同迭代次数的情况下, 最优取值不尽相同。He 等[20]认为 β 在一个合适的范围内取值都可以获得不错的分类效果。Dubes 等[21]和 Li[17]在基于马尔可夫随机场理论的图像聚类研究中, 均采纳Besag[19]β 取值为 1.5 的建议, 得到不错的聚类效果。因此, 本文中 β 的取值采用 1.5。
本文采用的光谱–空间聚类方法, 综合利用上述加入空间信息的聚类策略和基于 GMM 的光谱聚类方法, 流程可概括如下。首先, 采用 DWT 生成多尺度的图像。然后, 在得到的每个尺度图像中, 1)利用 GMM 方法进行光谱聚类, 得到每个像元的光谱类别概率 pspect(xi|ωk)和光谱聚类结果(类别标识); 2)利用 MRF 方法和 GMM 光谱聚类得到的类别标识, 计算每个像元的空间类别概率 pspat(xi|ωk); 3)结合该像元的光谱类别概率与空间类别概率, 得到对各类别的光谱–空间联合类别概率, 即将 pspect(xi|ωk)与 pspat(xi|ωk)相结合。具体地说, 像元 xi 对类别 ωk的光谱–空间联合类别概率 P*表达式[22]为
L 是低通滤波器, H 是高频滤波器, 数字代表进行离散小波变换的次数
图3 利用DWT得到多分辨率图像的示意图[16]
Fig. 3 Generation of multi-scale images using DWT [16]
(5)
邻域内的像元 xj 对类别 ωk 的光谱–空间联合类别概率 P*越大, 像元 xi 对类别 ωk 的空间类别概率就越大。
光谱–空间聚类结果按照以下顺序得到: 从DWT 得到的最低分辨率图像开始进行光谱–空间聚类, 得到光谱–空间联合类别概率 P*后, 对 P*设置一个较大的阈值(如 0.9), 即隶属某类别的概率较大, 判定类别的可靠性较大。如果一个像元的类别概率大于阈值, 其类别标识就是最终的聚类类别, 这个像元不再参与后续处理; 如果一个像元的类别概率小于阈值, 即不满足类别概率阈值条件, 则在下一尺度(更高分辨率)的图像上继续进行光谱–空间聚类, 直至所有像元满足设定的类别概率, 完成聚类过程[9]。
本研究使用的方法需要确定两个重要的参数值, 即 GMM 中的聚类数目 K 和 DWT 中的尺度参数S。GMM 中的聚类数目 K 由用户自行设定, 但聚类数目选择太多会增加运算量, 从而降低效率。同时, 基于有限的样本数据, 选择过多的聚类数目会导致GMM 中各参数无法被可靠地估计, 降低模型的效能[23]。
DWT 的尺度参数 S 为 0, 代表在原始尺度的图像上进行聚类。S 取值越大, 意味着不满足类别阈值设定的像元可在更多的尺度上进行聚类, 获得较高的类别概率, 提升聚类的可靠性。应综合考量运行效率, 选取合适的 S 值。
本研究采用由小到大的方式, 尝试不同 K 值与S 值的组合对聚类结果的影响, 据此选取合适的参数组合。
本研究采用两种定量评价指标来评价不同参数和不同聚类方法的效果。第一种指标是类簇熵(en-tropy), 是在已知地表类别的情况下对聚类效果进行评价, 反映聚类结果中每个类簇内像元的同质性和空间一致性。理想状况下, 每个类簇只含一类真实地物[24]。类簇熵值由每个类簇熵值的加权和求得[25], 熵值越小, 类簇内像元的同质性就越高, 聚类效果越好。第二种指标是 CHI(Calinski-Harabaz index), 是在地表类别未知的情况下对聚类效果进行评价[26], 通过计算类簇中样本点与类簇中心距离的平方和来度量类簇的紧密程度, 并综合考虑类簇间距离与类簇内距离。CHI 值越大, 说明簇内距离越小、簇间距离越大, 聚类效果越好。
本文采用监督分类的精度来评价聚类结果, 即将聚类的类别标识(光谱类)合并为地面实际类别(即信息类), 利用混淆矩阵得到相关的精度指标, 包括总体精度(overall accuracy, OA)、Kappa 值以及各岩性类别的用户精度(user’s accuracy, UA)和生产者精度(producer accuracy, PA)。
最后, 通过与 3 种方法(基于 GMM 的光谱聚类、对光谱聚类结果进行众数滤波平滑以及原始图像(S=0)的光谱–空间聚类)结果的比较, 验证本文采用的光谱–空间聚类方法的有效性。
首先, 保持 S 为常量, 令 K 从 7 变化到 30, 观察两种评价指标。如图 4 所示, 当 K 值不断增大时, CHI 先降后升, 在 K 为 12 时取得极大值, 随后开始下降并逐渐趋于平缓; 熵值总体上变化不大, 在 K为 12 时取得极小值。因此, 选择 K 为 12 可取得较好的聚类结果。
然后, 保持 K 为常量, 令 S 从 1 变化到 6, 观察两种评价指标。如图 5 所示, 当 S 为 2 时, CHI 最大, 熵值最小。因此, 选择最优尺度参数 S 为 2。
原始的 ASTER 图像经离散小波变换后, 得到不同尺度的图像, 我们选取研究区的局部区域来展示经 3 次离散小波变换后各尺度的假彩色图像 (图 6)。
观察表 2 中两种聚类评价指标可以发现, 只利用光谱特征的 GMM 聚类结果熵值最高(0.5462), CHI 值最低(633912); 光谱聚类结果经过滤波平滑(滤波窗口大小为 7×7)后, 熵值降为 0.5421, CHI 值提升至 654266; 在原始图像上的光谱–空间聚类结果熵值降为 0.5392, CHI 值提升至 667623; 本文方法的光谱–空间聚类结果熵值最低(0.5212), CHI 值最高(763509)。
图4 不同参数K的聚类评价指标数值
Fig. 4 Clustering evaluation indices with different value of K
图5 不同参数S的聚类评价指标数值
Fig. 5 Clustering evaluation indices with different value of S
(a) 原始 ASTER 图像; (b) S=1; (c) S=2; (d) S=3
图6 利用DWT得到的不同尺度图像
Fig. 6 Different scale images obtained with DWT methods
表2 不同方法聚类结果评价
Table 2 Evaluation of clustering results with different methods
评价指标GMM聚类GMM聚类后滤波原始尺度光谱–空间聚类本文的光谱–空间聚类 CHI633912654266667623763509 熵值0.54620.54210.53920.5212
我们选取研究区的局部区域, 更直观地对比不同方法的聚类结果。从图 7 可以看出, 相较于 3 种对比方法, 本文方法的光谱–空间聚类结果有明显的改善, 类簇边界更清晰, 类间混淆较少, 类内一致性较高。
(a)研究区局部区域的假彩色图像; (b)GMM 聚类; (c)众数滤波结果; (d)原始图像聚类; (e)本文方法的光谱–空间聚类。10 种颜色代表局部区域的 10 个光谱类别
图7 不同方法的局部聚类结果
Fig. 7 Clustering results in a selected local area with different methods
表3 不同方法聚类结果的精度(%)
Table 3 Accuracy of clustering results with different methods (%)
岩性类别GMM聚类GMM聚类后滤波原始尺度光谱–空间聚类本文的光谱–空间聚类PAUAPAUAPAUAPAUA Q3-472.9680.1771.3382.2380.6276.0986.3382.04 Q146.3767.6351.2570.3451.3673.3666.0283.92 N216.6210.3415.3612.4234.9639.9142.3340.34 C-P46.3736.5350.2639.3467.5742.2576.0867.67 D83.2180.1685.5680.1181.1692.2082.2093.40 S93.0092.0892.6794.3690.1991.3692.9693.00 Є-O72.3681.3273.5782.4676.9881.6780.1582.30 OA67.2969.3572.9381.62 Kappa54.0656.2760.8271.69
根据聚类中的光谱类别与地表实际信息类别(岩性类别)空间分布的对应关系, 以及光谱类别在信息类别中的占比, 将光谱类别合并至占比高的信息类别中, 即将 12 个聚类类别(K=12)合并为 7 个地表岩性类别。从各岩性类别中随机选取 300 个检验样本, 得到分类结果的混淆矩阵, 计算各种精度指标, 结果列于表 3。
从表 3 看出, 只使用光谱信息进行聚类的方法, 总体精度和 Kappa 值较低, 分别为 67.29%和 0.5406; 经过滤波平滑后, 总体精度和 Kappa 值分别为69.35%和 0.5627; 原始图像基于马尔可夫随机场进行光谱–空间聚类后, 总体精度和 Kappa 值进一步提高, 分别为 72.93%和 0.6082; 使用本文多尺度光谱–空间聚类方法得到的总体精度和 Kappa 值分别达到 81.62%和 0.7169, 可见本文加入空间信息的聚类方法显著地提高了岩性分类精度。
(a)研究区实际岩性分布; (b)GMM聚类; (c)众数滤波结果; (d)原始图像聚类; (e)本文方法的光谱–空间聚类。两种颜色的椭圆标记两个重点对比的局部区域
图8 不同方法的岩性分类结果
Fig. 8 Lithological classification results of different methods
表 3 还显示, 在使用本文的光谱–空间聚类后加入空间信息, 各岩性类别的 PA 和 UA 值均有一定程度的提升, Q1 和 C-P 的提升幅度最大, PA 最多分别提升 19.65%和 29.71%, UA 最多分别提升 16.29%和31.14%。
通过定量地评价各种方法得到的聚类结果可知, 相较于 3 种对比方法, 本文的光谱–空间聚类方法的各个精度指标均显著提升。
对比各种方法的分类结果可以发现, 在各类别的空间一致性以及与地质图(图 8(a))的相似程度方面, 使用本文光谱–空间聚类方法的结果(图 8(e))比光谱聚类结果(图 8(b))、光谱聚类后滤波平滑结果(图 8(c))、原始图像的光谱–空间聚类结果(图 8(d))有明显的改善, 表现为同类像元的聚集程度更高, “椒盐”噪声更小, 类间边界更清晰, 类别的空间分布与实际情况更一致, 有效地改善了光谱特征相近岩性的混淆情况。
图 8 中圈出两个岩性类别分布较多且易混淆的区域, 对照光谱聚类结果(图 8(b))和光谱聚类后滤波平滑的结果(图 8(c)), 可以发现图 8(d)和(e)中, 红色椭圆内 N2 的误分情况得到有效的缓解, Є-O 与 Q3-4 的混淆部分明显减少, 图 8(e)中两个类别的边界较为清晰, 与地质图的空间分布更加一致; 黑色椭圆内 Q1 和 D 被较多地识别出来。
本文采用两种策略进行光谱–空间聚类, 并用于多光谱图像的岩性识别与分类。与只利用光谱特征的高斯混合模型聚类方法、光谱聚类后滤波平滑和原始图像的光谱–空间聚类方法相比, 本文方法能获得更高精度的聚类结果。本文方法考虑局部邻域像元类别的空间相关性以及多尺度聚类, 能有效地改善因光谱特征相似造成的类别混淆, 提高岩性分类的精度。本文方法也适用于高光谱图像的岩性聚类, 在对高光谱图像进行光谱–空间聚类时, 应先进行特征提取或选择, 选择合适的特征波段参与聚类, 可以提高运行效率。因此, 综合利用光谱与空间信息进行聚类, 是一种有效的聚类方法, 适用于较难获取地面参考样本地区的岩性填图。
参考文献
[1] van der Meer F D, van der Werff H M A, van Rui-tenbeek F J A, et al. Multi- and hyperspectral geologic remote sensing: a review. International Journal of Ap-plied Earth Observation and Geoinformation, 2012, 14(1): 112–128
[2] Gomez C, Delacourt C, Allemand P, et al. Using AS-TER remote sensing data set for geological mapping, in Namibia. Physics and Chemistry of the Earth, 2005, 30: 97–108
[3] Li P, Yu H, Cheng T. Lithologic mapping using ASTER imagery and multivariate texture. Canadian Journal of Remote Sensing, 2009, 35(1): 117–125
[4] Murphy R J, Monteiro S T, Schneider S. Evaluating classification techniques for mapping vertical geology using field-based hyperspectral sensors. IEEE Transac-tions on Geoscience & Remote Sensing, 2012, 50(8): 3066–3080
[5] Feng J, Rivard B, Rogge D, et al. The longwave in-frared (3–14 μm) spectral properties of rock encrusting lichens based on laboratory spectra and airborne SE-BASS imagery. Remote Sensing of Environment, 2013, 131: 173–181
[6] 柯元楚, 史忠奎, 李培军, 等. 基于Hyperion高光谱数据和随机森林方法的岩性分类与分析. 岩石学报, 2018, 34(7): 2181–2188
[7] Kuhn S, Cracknell M J, Reading A M. Lithological mapping in the Central African Copper Belt using Ran-dom Forests and clustering: strategies for optimized results. Ore Geology Reviews, 2019, 112: 103015
[8] Bedini E. Mapping lithology of the Sarfartoq carbona-tite complex, southern West Greenland, using HyMap imaging spectrometer data. Remote Sensing of Envi-ronment, 2009,113(6): 1208–1219
[9] Hilger K B. Exploratory analysis of multivariate data [D]. Kongens Lyngby: Technical University of Den-mark, 2001
[10] Krooshof P W T, Postma G J, Melssen W J, et al. Effects of including spatial information in clustering multivariate image data, Trends in Analytical Chemi-stry, 2006, 25(11): 1067–1080
[11] 新疆维吾尔自治区地质矿产局. 赞比勒地质图(1:20万). 北京: 地质出版社, 1965
[12] Canty M J. Image Analysis, Classification and change detection in remote sensing (the fourth edition). Boca Raton: CRC Press, 2019: 329–373
[13] Redner R A, Walker H F. Mixture densities, maximum likelihood and the EM algorithm. SIAM Review, 1984, 26(2): 195–239
[14] Laarhoven P J M, Aarts E H L. Simulated Annealing: theory and Applications. London: Springer, 1987: 7–15
[15] Mallat S G. A theory for multiresolution signal decom-position: the wavelet representation. IEEE Transac-tions on Pattern Analysis and Machine Intelligence, 1989, 11(7): 674–693
[16] Pajares G, de la Cruz J M. A wavelet-based image fusion tutorial. Pattern Recognition, 2004, 37: 1855–1872
[17] Li S Z. Markov random field modeling in image anal-ysis. London: Springer, 2009: 183–215
[18] Wiemker R. Unsupervised fuzzy classification of mul-tispectral imagery using spatial-spectral features // Ma-thar R, Schader M. Classification, data analysis, and data highways. Berlin: Springer, 1998: 101–109
[19] Besag J. On the statistical analysis of dirty pictures. Journal of the Royal Statistical Society, 1986, 48(3): 259–302
[20] He L, Parikh N A. Automated detection of white matter signal abnormality using T2 relaxometry: application to brain segmentation on term MRI in very preterm infants. Neuro Image, 2013, 64: 328–340
[21] Dubes R C, Jain A K. Random field models in image analysis. Journal of Applied Statistics, 1989, 16(2): 131–164
[22] He L, Greenshields I R. An MRF spatial fuzzy clus-tering method for fMRI SPMs. Biomedical Signal Processing and Control, 2008, 3: 327–333
[23] Shende A, Mishra S, Kumar S. Comparison of different parameters used in GMM based automatic speaker recognition. International Journal of Soft Computing and Engineering, 2011, 1(3): 14–18
[24] Halkidi M, Batistakis Y, Vazirgiannis M. On clustering validation techniques. Journal of Intelligent Informa-tion Systems, 2001, 17(2): 107–145
[25] Tang H, Shen L, Qi Y, et al. A multiscale latent diri-chlet allocation model for object-oriented clustering of VHR panchromatic satellite images. IEEE Transac-tions of Geoscience and Remote Sensing, 2013, 51(3): 1680–1692
[26] Caliński T, Harabasz J. A dendrite method for cluster analysis. Communications in Statistics, 1974, 3(1): 1–27
Lithological Mapping from Multispectral Images Using Spectral-Spatial Clustering
Abstract To address the class confusion caused by using only spectral features in image clustering for lithological recognition and classification based on multispectral- hyperspectral data, spectral-spatial features are used in multispectral image clustering. Two ways are adopted to include spatial information in clustering. First, multi-scale images are generated using discrete wavelet transform, and are successively clustered from coarse to fine scale images. Second, the spatial class membership probability is calculated using Markov random field model in different image scales, and the combined spectral-spatial class membership probability is then obtained. Finally, the final clustering results are obtained by combining the spectral-spatial clustering results of multiple scales. ASTER multispectral data were used for lithological mapping, and compared with the spectral clustering method based on Gaussian mixture model clustering (spectral clustering), the filtering after spectral clustering, and the spectral-spatial clustering on original image. The results show that the method of spectral-spatial clustering on multiscale image adopted obtained higher accuracies of lithological mapping than the above three comparative methods. Using spectral and spatial information at multi-scales in clustering is an effective method for lithological mapping. The method is applicable to lithological mapping in areas where training samples are difficult to obtain.
Key words spectral-spatial clustering; lithological mapping; multispectral images