CN101853509B - 基于Treelets和模糊C-均值聚类的SAR图像分割方法 - Google Patents

Abstract

本发明公开了一种基于Treelets和模糊C-均值聚类的SAR图像分割方法，主要解决现有技术在人工情况下无法对大量图像进行目标识别的问题。其实现步骤是：(1)提取输入的待分割SAR图像纹理特征和小波特征；(2)将纹理特征和小波特征组成特征矩阵；(3)对此特征矩阵进行Treelets变换，得到基矩阵；(4)将特征矩阵与基矩阵的尺度函数相乘得到一个结构向量；(5)对此结构向量采用模糊C-均值聚类，得到SAR图像的分割结果。本发明利用Treelets变换可对聚类前的高维数据降维，能有效地抑制噪声并提高SAR图像的分割速度，用于对SAR图像的自动分割以进行目标识别。

Description

基于Treelets和模糊C-均值聚类的SAR图像分割方法

技术领域

本发明属于图像处理领域，涉及SAR图像的分割方法，具体地说是一种基于Treelets和模糊C-均值聚类的SAR图像分割方法，可用于目标识别。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)作为主动式雷达，具有全天时、全天候、多极化、多视角等优点，在军事、遥感等领域得到了广泛的应用。由于SAR图像是由微波回波成像形成的，复杂的地物背景以及由此形成的散射相干斑噪声，使得SAR图像比自然图像和医学影像的分割复杂。目前SAR图像的分割方法主要分为不依赖先验知识的数据驱动的分割方法和建立在先验知识基础上的模型驱动的分割方法两类。模型驱动的SAR图像分割方法需要对SAR图像分割问题构建数学模型，主要包括基于组合优化模型的分割方法、基于Markov随机场模型的分割方法和基于多尺度模型的分割，这些方法依赖于SAR图像的先验知识，需要针对具体的问题假设精确的模型和复杂的算法来实现，通用性不足。数据驱动的SAR图像分割方法主要包括基于阈值的分割方法、基于边缘的分割方法和基于区域的的分割方法。其中基于阈值的分割方法速度快但不精确，基于边缘的分割方法由于相干斑噪声的存在使得边缘定位精度难以提高，基于区域的分割方法将均匀同质的像素划分或聚集为若干区域，分割方法相对较为简单精确，但其相似性检验准则的最优确定、迭代收敛速度的提高等都是难点。

模糊C-均值(Fuzzy C-Means，FCM)聚类算法是一种被广泛应用的聚类算法，该算法是基于数据样本的隶属度划分的，即在样本空间中隶属度高的样本被聚在一类，理想的聚类过程样本间被认为是相互独立的。传统FCM聚类算法用于图像分割时将图像的灰度特征作为样本，而真实图像在像素邻域间通常具有较强的相关性，在一个目标中的相邻像素间通常不是相互独立的。当噪声较大时，FCM聚类过程往往难以得到理想的分割结果，尤其是对SAR图像的相干斑噪声不能很好的抑制。如果考虑每个像素点及其周围的邻域信息，利用多种特征则可以有效抑制相干斑噪声对聚类的影响，但高维的冗余多特征增加了聚类所需的时间。因此就需要将冗余的数据特征向量从高维空间映射到低维空间中，使得降维后的数据能够最大程度的保留最重要的或者最能体现个体差异的特性。

一种能够对高维、无序、含噪数据的自适应多尺度表示和分析的方法——Treelets变换特别适合用于降维和特征选择。其思想来源于PCA、小波和层次聚类树，它通过构建层级树和标准正交基，可以反映出数据的内在结构，实现数据的自适应稀疏表示。因此将其作为聚类前的降维工具可有效地降低聚类时间且对噪声具有一定的鲁棒性。

目前将Treelets用于图像处理中，仅有黄珊等针对高分辨全色遥感图像的桥梁检测方法。该方法首先对全色高分辨率遥感图像根据桥梁先验知识建立桥梁知识库，然后计算像素邻域内的多个特征，利用Treelets变换进行特征降维，最后分割水域。该方法选择了8种不同窗口下的11种特征，即最大值、最小值、中值、对比度和前四阶累积量组成88维的特征矩阵。这样构建的特征数据，虽然维数较高但特征种类单一，特征间冗余较大，使得该方法在特征提取上需要消耗大量的时间，这一定程度上影响了该方法应用的实时性。同时该方法用于分割需要先人为选择一定的训练样本和测试样本来确定Treelets变换的最佳基矢量，最后将图像投影到最佳基矢量上，采用KNN分类器确定分割结果。然而Treelets本身是一种数据驱动的分析方法，经变换后得到的多尺度正交基本身就是对原始数据内部结构的反应，该方法没有充分利用Treelets本身数据驱动的特性，所以上述方法无论是在Treelets最佳基的确定还是最后分类器的训练上都需要通过人为选择的样本来确定，这种加入人工辅助判别图像的方法，对于图像处理的自动化是不利的，无法对大量图像进行目标识别，因此降低了它的实用性。

发明内容

本发明的目的在于避免上述已有技术的不足，提供一种基于Treelets和模糊C-均值聚类的SAR图像分割方法，以实现在无需加入人工辅助判别的情况下，自动获得SAR图像的分割结果，同时利用Treelets对聚类前的高维数据降维，降低运行时间。

为实现上述目的，本发明包括如下步骤：

(1)输入待分割SAR图像，在大小为17×17像素的滑窗内，分别提取每个中心像素点的灰度值z、均值m、标准s、平滑度R、三阶矩μ₃、一致性U、熵S、能量EN和同质性H九种纹理特征；

(2)将输入的SAR图像进行3层非下采样离散小波变换，并在大小为16×16像素的滑窗内，提取3层中每个中心像素点的小波特征：e_LL-1、e_LH-1、e_HL-1、e_HH-1、e_LH-2、e_HL-2、e_HH-2、e_LH-3、e_HL-3和e_HH-3；

(3)将步骤(1)得到的纹理特征与步骤(2)得到的小波特征合并，得到19维的待分割SAR图像的特征矩阵X：

X＝[z，m，σ，R，μ₃，U，S，EN，H，e_LL-1，e_LH-1，e_HL-1，e_HH-1，e_LH-2，e_HL-2，e_HH-2，e_LH-3，e_HL-3，e_HH-3]；

(4)将特征矩阵X进行Treelets变换，得到一个正交基矩阵B；

            B＝[Φ ψ₁…ψ_l-1]^T

其中，Φ和ψ分别是正交基矩阵的尺度函数和细节函数，l是分解层数。

(5)将特征矩阵X与尺度函数Φ相乘，得到一个结构向量FX＝[f₁，f₂，…，f_n]^T；其中，f分别对应于SAR图像每个像素的数据点；

(6)对结构向量FX采用模糊C-均值聚类，得到SAR图像的分割结果。

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

a、纹理特征在一定程度上反映了SAR图像中各灰度级在空间上的分布特性，而小波特征能反映SAR图像本身的几何特性，本发明采用将纹理特征和小波特征两类特征的组合，能够弥补单类特征的缺陷，提高SAR图像分割的效果。

b、本发明利用Treelets数据驱动的特性，无需人工选择样本确定最佳基，而是利用样本在尺度函数上的投影来反映样本本身的结构，可实现对SAR图像的自动分割，解决了人工情况下无法对大量图像进行目标识别的问题。

c、本发明利用Treelets对特征矩阵变换，降低了数据维度，避免了对大规模数据的处理，降低了运行时间，可实现对SAR图像的快速分割。

d、由于Treelets变换对噪声有一定的抑制作用，变换后能减小噪声对于分割的影响，因此本发明对于噪声具有较好的鲁棒性。

仿真实验结果表明，本发明能够有效地应用于SAR图像分割，并可以实现对高维数据的快速处理。

附图说明

图1是本发明的实现流程图；

图2是本发明的第一组实验数据和结果图；

图3是本发明的第二组实验数据和结果图。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1，输入待分割SAR图像，提取每个中心像素点的纹理特征。

(1a)在大小为17×17像素的滑窗内，计算每个中心像素点的均值m、标准s、平滑度R、三阶矩μ₃、一致性U、熵S、能量EN和同质性H纹理特征，具体计算公式如下：

$m=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{z}_{i}p\left(z_i\right)$

$\mathrm{\σ}=\sqrt{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{(z_i-m)}^{2}p\left(z_i\right)}$

R＝1-1/(1+σ²)

${\mathrm{\μ}}_3=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{(z_i-m)}^{3}p\left(z_i\right)$

$U=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{p}^2\left(z_i\right)$

$S=-\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}p\left(z_i\right){\log }_{2}p\left(z_i\right)$

$\mathrm{EN}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{\left(p\left(z_i\right)\right)}^2$

$H=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}\frac{p(z_i,z_j)}{1+{(z_i-z_j)}^2}$

其中，z_i是表示像素第i个灰度级的灰度值，p(z_i)是滑窗内的灰度级的直方图；

(1b)将每个像素点的灰度值z与计算所得的8个特征一起共9个特征作为一个特征行向量，提取出来，图像中所有n个像素点的特征行向量组成n×9的纹理特征矩阵X1。

步骤2，对输入的SAR图像，提取每个中心像素点的小波特征。

(2a)将输入的SAR图像进行3层非下采样离散小波变换，并在大小为16×16像素的滑窗内，计算3层中每个中心像素点的小波特征，在每一层都得到HH、HL、LH和LL四个子图像，在各个子图像中，以大小为me×ne像素的滑窗内，计算每个中心像素点的小波分解能量值，共得到10个能量特征，分别记为e_LL-1、e_LH-1、e_HL-1、e_HH-1、e_LH-2、e_HL-2、e_HH-2、e_LH-3、e_HL-3、e_HH-3，其中HL-1、LH-2和HH-3分别表示第1、2、3层分解的HL、LH和HH子图像，这里的能量特征采用l₁-范数计算，即：

$e=\frac{1}{\mathrm{me}\·\mathrm{ne}}\underset{\mathrm{ie}=1}{\overset{\mathrm{me}}{\mathrm{\Σ}}}\underset{\mathrm{je}=1}{\overset{\mathrm{ne}}{\mathrm{\Σ}}}\left|C(\mathrm{ie},\mathrm{je})\right|$

其中，C为该邻域窗内的小波系数，ie和je表示邻域窗内小波系数的索引，邻域窗的大小me＝ne＝16；

(2b)将每个像素点的10个小波特征作为一个特征行向量提取出来，则图像中所有n个像素点的特征行向量组成一个n×10的小波特征矩阵X2。

步骤3，将步骤1得到的纹理特征X1与步骤2得到的小波特征X2合并，得到19维的待分割SAR图像的特征矩阵X：

X＝[z，m，σ，R，μ₃，U，S，EN，H，e_LL-1，e_LH-1，e_HL-1，e_HH-1，e_LH-2，e_HL-2，e_HH-2，e_LH-3，e_HL-3，e_HH-3]；

步骤4，将步骤3得到的特征矩阵X进行Treelets变换，得到一个正交基矩阵B。

(4a)在Treelets变换的最底层，即聚类层级l＝1层，输入n×19的特征矩阵X，并计算X的协方差矩阵

和相似度矩阵

${\hat{S}}^{\left(0\right)}=\left[\begin{array}{cccc}{\mathrm{\σ}}_{11}& {\mathrm{\σ}}_{12}& \·\·\·& {\mathrm{\σ}}_{1q}\\ {\mathrm{\σ}}_{21}& {\mathrm{\σ}}_{22}& \·\·\·& {\mathrm{\σ}}_{2q}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ {\mathrm{\σ}}_{q1}& {\mathrm{\σ}}_{q2}& \·\·\·& {\mathrm{\σ}}_{\mathrm{qq}}\end{array}\right]$

${\hat{M}}^{\left(0\right)}=\left[\begin{array}{cccc}{\mathrm{\ρ}}_{11}& {\mathrm{\ρ}}_{12}& \·\·\·& {\mathrm{\ρ}}_{1q}\\ {\mathrm{\ρ}}_{21}& {\mathrm{\ρ}}_{22}& \·\·\·& {\mathrm{\ρ}}_{2q}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ {\mathrm{\ρ}}_{q1}& {\mathrm{\ρ}}_{q2}& \·\·\·& {\mathrm{\ρ}}_{\mathrm{qq}}\end{array}\right]$

其中

u＝1，2，…，q，v＝1，2，…，q，q是X的总列数，X_u和X_v分别是特征矩阵的第u列和第v列向量，σ_uv是X_u和X_v两向量的方差，ρ_uv是X_u和X_v两向量的相关系数；

(4b)初始化B₀为q×q的单位矩阵：

并且初始化和变量的下标集为δ＝{1，2，…，q}；

(4c)由相似度矩阵

找出最相似的两个变量α和β：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{x,y\∈\mathrm{\δ}}{\max {{\hat{M}}_{\mathrm{xy}}}^{(l-1)}}$

其中x和y是相似度矩阵的位置索引，且必须属于和变量下标集δ，并且x＜y；

(4d)在变量对α和β上进行局部主成分分析变换，其第一和第二主成分分别为变量s_l和差变量d_l，且使得

并得到雅克比旋转矩阵J：

其中c＝cos(θ_l)，s＝sin(θ_l)，|θ_l|≤π/4；

(4e)根据雅克比旋转矩阵J计算当前聚类层级的基矩阵B_l＝B_l-1J，以及相似度矩阵

和协方差矩阵

尺度函数φ_l和细节函数ψ_l分别为基矩阵B_l的第α和β列，其中，当前层级的尺度函数Φ是尺度函数φ_l和上一层的尺度向量集{φ_l-1，jt}_jt≠α，β的合集，将差变量的下标β从和变量的下标集δ中去除，即δ＝δ\{β}；

(4f)重复步骤(4c)至步骤(4e)直至到l＝18层，得到基矩阵为：

            B＝[Φ ψ₁…ψ_l-1]^T

其中，Φ和ψ分别是正交基矩阵的尺度函数和细节函数，l是分解层数。

步骤5，将特征矩阵X与尺度函数Φ相乘，得到一个结构向量FX＝[f₁，f₂，…，f_n]^T，其中，f分别对应于SAR图像每个像素的数据点。

步骤6，对结构向量FX采用模糊C-均值聚类，得到SAR图像的分割结果。

本发明的效果可以通过以下内容进行说明：

1.实验数据

为了便于定量评价本发明所提出图像分割方法的性能，首先选择带有分割标准参考图的人工合成纹理图像及对其加入乘性噪声的图像。人工合成纹理图像大小为256×256，灰度级为256，如图2(a)所示，图2(b)为图2(a)的标准参考图。将图2(a)加入均值为零、噪声方差分别为0.01、0.02、0.03、0.04和0.05的乘性噪声后的加噪图像分别如图2(c)、图2(d)、图2(e)、图2(f)和图2(g)所示。

四幅真实的SAR图像如图3所示，它们的灰度级均为256，其中图3(a)是一幅ERS-SAR图像，大小为256×256，它包括庄稼和森林两类；图3(b)是一幅包括山地和平原两类的SAR图像，大小为256×256；图3(c)是美国新墨西哥州Albuquerque地区附近的格兰德河的Ku波段SAR图像，大小为256×256，该图分为农田、植被和河流三类。图3(d)是一幅来自美国Sandia国家实验室的Washington地区的Ku波段SAR图像，大小为512×512，其中包含3类地物：河流、城区和其它。

2.实验评价指标

实验评价指标包括聚类错误率、Kappa系数和运行时间。聚类错误率是分割结果图与标准参考图中不为同一类的像素的个数除以图像的总像素数。Kappa系数是检验图像分类精度的一种常用指标，通常，当Kappa＝0.75时类别间变化频率小；当0.4＝Kappa＝0.75时，类别间变化频率一般；当Kappa＝0.4时，类别间变化频率较大。计算公式如下：

$\mathrm{Kappa}=\frac{\mathrm{KN}\underset{\mathrm{io}=1}{\overset{r}{\mathrm{\Σ}}}x{f}_{\mathrm{io},\mathrm{io}}-\underset{\mathrm{io}=1}{\overset{r}{\mathrm{\Σ}}}\left({\mathrm{xf}}_{\mathrm{ip}}{\mathrm{xf}}_{\mathrm{pi}}\right)}{{\mathrm{KN}}^2-\underset{\mathrm{io}=1}{\overset{r}{\mathrm{\Σ}}}\left({\mathrm{xf}}_{\mathrm{ip}}{\mathrm{xf}}_{\mathrm{pi}}\right)}$

式中，r是混淆矩阵中总列数，即总的类别数，xf_io，io是误差矩阵中第io行第io列上像素数量，即正确聚类的数目，xf_ip、xf_pi分别是第ip行和第pi列的总像素数量，KN是总的用于精度评估的像素数量。其中混淆矩阵是通过将像素真实的类别和聚类后像素的相应类别比较计算的。混淆矩阵的每一列中的数值等于真实的类别在聚类后对应于相应类别的数量；混淆矩阵的每一行中的数值等于聚类后相应类别在真实类别中的数量。

3.实验结果与分析

为评价本发明与其它FCM方法的性能，本发明将所提出的基于Treelets的FCM聚类方法(简记为Treelets-FCM方法)与直接FCM聚类(简记为FCM方法)和用PCA降维后FCM聚类(简记为PCA-FCM方法)两种方法，对人工合成纹理及其加噪图像与真实的SAR图像的分割结果进行了比较。

实验硬件环境为Intel(R)酷睿2 6300 1.86GHz，内存2G，实验平台为MATLAB 7.0。

表1对人工合成纹理图像及其加噪图像采用三种分割方法的实验结果对比

对于人工合成纹理及其加噪图像，采用FCM方法、PCA-FCM方法和本发明方法实验所得分类错误率、Kappa系数和运行时间结果分别列于表1中。从表1可以看到在噪声方差为0.01时的类别错误率FCM方法为11.7142％，PCA-FCM方法为3.1998％，而本发明则降为2.4658％，比FCM方法下降了近10％。本发明的kappa系数为0.9233，高于FCM方法的0.6961和PCA-FCM方法的0.9019。其它噪声情况下类似。因此本发明在分类正确性和kappa系数上要明显优于其他两种方法。在噪声方差为0.04时FCM方法为12.3906秒，PCA-FCM方法为1.0469秒，而Treelets-FCM方法则降为0.6875秒，使运行时间降到了1秒以内。所以在运行时间上本发明也要低于另外两种方法。综合三种评价指标，本发明方法都要优于FCM和PCA-FCM两种方法。另外，在不同乘性噪声下三种方法的分割结果的类别错误率和Kappa系数也表明本方法对噪声有一定的鲁棒性。

表2四幅真实SAR图像采用三种分割方法的实验运行时间(s)的结果对比

对于真实的SAR图像，由于没有标准参考图，所以只能给出图像的分割结果进行主观评价。对应图3(a)、图3(b)和图3(c)两分类真实SAR图像，采用FCM方法的分割结果分别如图3(a1)、图3(b1)和图3(c1)所示；采用PCA-FCM方法的分割结果分别如图3(a2)、图3(b2)和图3(c2)所示；采用本发明的Treelets-FCM方法的分割结果分别如图3(a3)、图3(b3)和图3(c3)所示。

从图3三种方法对四幅真实SAR图像的实验结果来看，在分割直观视觉效果上三种方法基本相同，没有大的区别。表2为四幅真实SAR图像运行时间的对比，如第二幅真实SAR图像本发明运行时间为0.6563秒，低于FCM的运行时间2.3906秒和PCA-FCM的运行时间0.8125秒。在图像较大且分类数较多的情况下，优势更加明显。如第四幅真实SAR图像本发明的运行时间为8.4844秒，比FCM的运行时间104.4063秒降低了近百秒。因此本发明在时间的耗费上要低于FCM和PCA-FCM两种方法，可完成对真实SAR图像的快速分割。

从以上两类实验可以看出，本发明在分类错误率和运行时间上都要优于对比的两种方法，对于真实SAR图像的分割也更加快速。

Claims (4)

1.一种基于Treelets和模糊C-均值聚类的SAR图像分割方法，包括如下步骤：

(1)输入待分割SAR图像，在大小为17×17像素的滑窗内，分别提取每个中心像素点的灰度值z、均值m、标准差σ、平滑度R、三阶矩μ₃、一致性U、熵S、能量EN和同质性H九种纹理特征；

(2)将输入的SAR图像进行3层非下采样离散小波变换，并在大小为16×16像素的滑窗内，提取3层中每个中心像素点的小波特征：e_LL-1、e_LH-1、e_HL-1、e_HH-1、e_LH-2、e_HL-2、e_HH-2、e_LH-3、e_HL-3和e_HH-3；

(3)将步骤(1)得到的纹理特征与步骤(2)得到的小波特征合并，得到19维的待分割SAR图像的特征矩阵X：

X＝[z，m，σ，R，μ₃，U，S，EN，H，e_LL-1，e_LH-1，e_HL-1，e_HH-1，e_LH-2，e_HL-2，e_HH-2，e_LH-3，e_HL-3，e_HH-3]；

(4)将特征矩阵X进行Treelets变换，得到一个正交基矩阵B；

B＝[Фψ₁…ψ_l-1]^T

其中，Ф和ψ分别是正交基矩阵的尺度函数和细节函数，l是分解层数；

(5)将特征矩阵X与尺度函数Ф相乘，得到一个结构向量FX＝[f₁，f₂，…，f_n]^T，其中，f分别对应于SAR图像每个像素的数据点；

(6)对结构向量FX采用模糊C-均值聚类，得到SAR图像的分割结果。

2.根据权利要求1所述SAR图像分割方法，其中步骤(1)所述的提取每个中心像素点的均值m、标准差σ、平滑度R、三阶矩μ₃、一致性U、熵S、能量EN和同质性H纹理特征，是先按如下公式先进行计算，再将其计算结果提取出来，即

$m=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{z}_{i}p\left(z_i\right)$

$\mathrm{\σ}=\sqrt{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{(z_i-m)}^{2}p\left(z_i\right)}$

R＝1-1/(1+σ²)

${\mathrm{\μ}}_3=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{(z_i-m)}^{3}p\left(z_i\right)$

$U=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{p}^2\left(z_i\right)$

$S=-\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}p\left(z_i\right){\log }_{2}p\left(z_i\right)$

$\mathrm{EN}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{\left(p\left(z_i\right)\right)}^2$

$H=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}\frac{p(z_i,z_j)}{1+{(z_i-z_j)}^2}$

其中，z_i是表示像素第i个灰度级的灰度值，p(z_i)是滑窗内的灰度级的直方图。

3.根据权利要求1所述SAR图像分割方法，其中步骤(2)所述的提取3层中每个中心像素点的小波特征，是在每一层都得到HH、HL、LH和LL四个子图像，在各个子图像中，以大小为me×ne像素的滑窗内，计算每个中心像素点的小波分解能量值，共得到10个能量特征，分别记为e_LL-1、e_LH-1、e_HL-1、e_HH-1、e_LH-2、e_HL-2、e_HH-2、e_LH-3、e_HL-3、e_HH-3，其中HL-1、LH-2和HH-3分别表示第1、2、3层分解的HL、LH和HH子图像，这里的能量特征采用l₁-范数计算，即：

$e=\frac{1}{\mathrm{me}\·\mathrm{ne}}\underset{\mathrm{ie}=1}{\overset{\mathrm{me}}{\mathrm{\Σ}}}\underset{\mathrm{je}=1}{\overset{\mathrm{ne}}{\mathrm{\Σ}}}\left|C(\mathrm{ie},\mathrm{je})\right|$

其中，C为该邻域窗内的小波系数，ie和je表示邻域窗内小波系数的索引，邻域窗的大小me＝ne＝16。

4.根据权利要求1所述SAR图像分割方法，其中步骤(4)所述的将特征矩阵X进行Treelets变换，按照以下步骤进行：

(4a)在Treelets变换的最底层，即聚类层级l＝1层，输入n×19的特征矩阵X，并计算X的协方差矩阵

和相似度矩阵

${\widehat{\mathrm{\Σ}}}^{\left(0\right)}=\left[\begin{array}{cccc}{\mathrm{\σ}}_{11}& {\mathrm{\σ}}_{12}& ...& {\mathrm{\σ}}_{1q}\\ {\mathrm{\σ}}_{21}& {\mathrm{\σ}}_{22}& ...& {\mathrm{\σ}}_{2q}\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ {\mathrm{\σ}}_{q1}& {\mathrm{\σ}}_{q2}& ...& {\mathrm{\σ}}_{\mathrm{qq}}\end{array}\right]$

${\hat{M}}^{\left(0\right)}=\left[\begin{array}{cccc}{\mathrm{\ρ}}_{11}& {\mathrm{\ρ}}_{12}& ...& {\mathrm{\ρ}}_{1q}\\ {\mathrm{\ρ}}_{21}& {\mathrm{\ρ}}_{22}& ...& {\mathrm{\ρ}}_{2q}\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ {\mathrm{\ρ}}_{q1}& {\mathrm{\ρ}}_{q2}& ...& {\mathrm{\ρ}}_{\mathrm{qq}}\end{array}\right]$

其中

${\mathrm{\ρ}}_{\mathrm{uv}}=\frac{{\mathrm{\σ}}_{\mathrm{uv}}}{\sqrt{{\mathrm{\σ}}_{\mathrm{uu}}{\mathrm{\σ}}_{\mathrm{vv}}}},$ u＝1，2，…，q，v＝1，2，…，q，q是X的总列数，X_u和X_v分别是特征矩阵的第u列和第v列向量，σ_uv是X_u和X_v两向量的方差，ρ_uv是X_u和X_v两向量的相关系数；

(4b)初始化B₀为q×q的单位矩阵：

并且初始化和变量的下标集为δ＝{1，2，…，q}；

(4c)由相似度矩阵

找出最相似的两个变量α和β：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{x,y\∈\mathrm{\δ}}{\max {{\hat{M}}_{\mathrm{xy}}}^{(l-1)}}$

其中x和y是相似度矩阵的位置索引，且必须属于和变量下标集δ，并且x＜y；

(4d)在变量对α和β上进行局部主成分分析变换，其第一和第二主成分分别为变量s_l和差变量d_l，且使得

并得到雅克比旋转矩阵J：

其中c＝cos(θ_l)，s＝sin(θ_l)，|θ_l|≤π/4；

(4e)根据雅克比旋转矩阵J计算当前聚类层级的基矩阵B_l＝B_l-1J，以及相似度矩阵

和协方差矩阵尺度函数φ_l和细节函数ψ_l分别为基矩阵B_l的第α和β列，其中，当前层级的尺度函数Ф是尺度函数φ_l和上一层的尺度向量集{φ_l-1，jt}_jt≠α，β的合集，将差变量的下标β从和变量的下标集δ中去除，即δ＝δ\{β}；

(4f)重复步骤(4c)至步骤(4e)直至到l＝18层，得到基矩阵为：

B＝[Фψ₁…ψ_l-1]^T。

Images