合成孔径雷达(synthetic aperture radar，SAR)经过多年的发展^[1-3]，成像方面很多问题已经得到解决^[4-5].目前的热点和难点主要集中在新体制雷达的信号处理^[6-8]或SAR图像的解译与应用研究上^[9-10].

图像分割能够从图像中提取感兴趣的信息，是从图像处理到图像分析的关键步骤.图像分割的方法有多种，主要有基于阈值、边缘、区域、特定理论的分割方法等^[11-13].

聚类作为一种无监督的分类方法，在众多领域应用广泛^[13]，例如生物学上的基因分类和动植物分类，图像处理中的分割、增强、压缩、检索等.K均值聚类算法^[14-17]是一种典型的基于划分的聚类算法，该算法思想简单、计算速度快，已经成为最常用的聚类算法之一.

本文以毫米波高分辨SAR图像为研究对象，采用K均值聚类方法对SAR图像进行分割，并对文献中的动态K均值聚类算法进行分析与改进，具体改进了每个聚类适应度函数的计算方式，即用聚类中样本数的正比函数对原始的适应度进行平均.毫米波SAR图像城市区域、路和桥梁的分割结果验证了本文改进算法相对已有算法在分割效率上的提升.

1 当前的K均值聚类算法

K均值聚类算法的基本思路是在最小化误差函数的基础上将数据划分成K类.算法的处理过程^[11-12]为：先指定聚类数目K及K个初始聚类中心，然后根据一定的准则将每个数据分配给最近的聚类中心.

将样本空间X={x₁, x₂, …, x_i, …, x_n}的样本分成K类，聚类中心为C={c₁, c₂, …, c_j, …, c_K}，用d_ij(x_i, c_j)表示样本x_i与其对应的中心c_j间的距离，样本空间内所有数据点与所属聚类质心距离的总和用目标函数J来表示，为

$J=\sum\limits_{j=1}^{K}{\sum\limits_{i,i\in {{c}_{j}}}{{{d}_{ij}}\left( {{x}_{i}},{{c}_{j}} \right).}}$

(1)

目标函数J越小，表明聚类越紧凑，聚类越优.当选择欧式距离作为样本x_i与其对应的中心c_j间的距离时^[12]，$x_{i}^{\left( j \right)}$为属于聚类j的数据样本，n_j为聚类j的样本个数，式(1)为

$J=\sum\limits_{j=1}^{K}{{{J}_{i}}}=\sum\limits_{j=1}^{K}{\left[ \sum\limits_{i,i\in {{c}_{j}}}{\left\| x_{i}^{\left( j \right)} \right.-{{\left. {{c}_{j}} \right\|}^{2}}} \right].}$

(2)

为使目标函数最小，各聚类中心为

${{c}_{j}}=\frac{1}{{{n}_{j}}}\sum\limits_{i=1}^{{{n}_{j}}}{x_{i}^{\left( j \right)}}.$

(3)

K均值聚类算法流程如下：

1)初始化，输入样本集X及聚类数K，并在X中随机选取K个样本作为初始聚类中心；

2)初始聚类；

3)按式(3)计算新的聚类中心；

4)重新聚类；

5)反复进行3)、4)直至迭代结束，得到聚类结果.

称上述K均值聚类算法为基本K均值聚类算法，记为KM_Basic.KM_Basic算法理论严密，计算简单.但是聚类结果对初始聚类中心有很强的依赖性，容易收敛于局部极值点.初始聚类中心选取不当会对聚类结果产生很大的负面影响.

文献[13]提出全局K均值聚类算法，文献[14]对全局K均值聚类算法进行研究.全局K均值聚类算法通过迭代的方式来产生初始聚类中心，处理效率不高.文献[12]提出动态K均值聚类算法，动态K均值聚类算法能减小对聚类中心初值的依赖，改善性能，并且运算效率也较高.动态K均值聚类算法的适应度函数为聚类中心与属于该中心区域内的所有像素之间的欧式距离之和：

$f\left( {{c}_{j}} \right)=\sum\limits_{i}{\left\| {{x}_{i}} \right.-}{{\left. {{c}_{j}} \right\|}^{2}}.$

(4)

若f(c_j)越小，则中心c_j的适应度越小，聚类越紧凑.动态K均值算法通过调整聚类来使各中心的适应度函数均衡，当适应度均衡时，认为聚类最优.动态K均值聚类算法流程如下：

1)给定初始聚类中心c_j和权值α₀(α₀为常数)，α_a=α_b=α₀；

2)初始聚类；

3)根据式(4)计算每个聚类中心的适应度函数f(c_j)；

4)设f(*)中，最大的f(*)对应的聚类中心为c_l，最小的f(*)对应的聚类中心为c_s.如果f(c_s)<α_a f(c_l)，重新分配c_l聚类中的数据样本，将其中x_i<c_l的数据样本分配给聚类c_s；

5)根据式(3)计算新的聚类中心；

6)更新阈值α_a=α_a－α_a/K，重复3)至5)，直至f(c_s)≥α_a f(c_l)；

7)按最小距离原则聚类一次；

8)根据式(3)计算新的聚类中心；

9)更新权值α_a=α₀和α_b=α_b－α_b/K，重复3)至8)，直至f(c_s)≥α_b f(c_l).

记上述动态K均值聚类算法为KM_MKM.KM_MKM通过调整具有最大适应度函数值和最小适应度函数值的聚类区域的样本，最终达到各区域适应度函数的均衡.KM_MKM能减少对初始聚类中心的依赖，改善并减少陷入局部极值引起的死区中心和中心冗余问题，但该算法对孤立数据和噪声敏感的问题依然存在^[12].

KM_MKM将最大适应度函数聚类里面的部分样本强制分配给最小适应度函数聚类区域，如果这些样本是噪声，那么具有最小适应度的聚类就将代表噪声数据，这将导致错误的分类^[12].基于此，文献[12]提出一种自适应动态K均值聚类算法，改变KM_MKM算法强制分配样本给具有最小适应度聚类的做法，引入最小距离原则来分配这些样本，即将待分配样本配给与这些样本最近距离的聚类区域，以此来减少噪声数据的错误分类.自适应动态K均值算法记为KM_AMKM，KM_AMKM与KM_MKM流程类似，只是在步骤4)中将其中x_i<c_l的数据样本按照最小距离原则分配给最近的聚类中心，而非统一分配给聚类c_s.

2 当前算法分析与改进

KM_Basic以J最小为准则迭代聚类，文献[14]的KM_MKM及文献[12]的KM_AMKM均把各J_j均衡作为准则进行迭代聚类.不考虑运算效率，仅针对分割质量指标F, F′, Q而言，在文献[12]中，KM_MKM、KM_AMKM均优于KM_Basic，表明“J_j均衡准则”大部分时候还是优于“J最小准则”的.本节仍然遵循“J_j均衡准则”的思路，但是对该准则进行相应的改进.

根据式(4)，J_j为

${{J}_{j}}=\sum\limits_{i}{\left\| {{x}_{i}} \right.-}{{\left. {{c}_{j}} \right\|}^{2}}.$

(5)

J_j为各聚类里面所有样本与聚类中心的欧式距离之和.J_j组成元素里面各个部分‖x_i－c_j‖²均为非负的，一般而言，聚类所含样本个数越多，聚类适应度函数越大.对于一个需要聚类的数据库来说，如果数据库中样本本质上属于几个类，且类的大小均差不多，那么J_j均衡准则假设比较合理.但是如果数据库样本在本质上所属的类的样本数相差悬殊，那么可以认为，样本数较多的类的适应度函数J_j较大是合理的，这时再强制要求各个聚类的J_j均衡反而不合理.因此，改进式(5)中J_j的定义，为

${{J}_{j}}=\frac{\sum\limits_{i}{\left\| {{x}_{i}} \right.-}{{\left. {{c}_{j}} \right\|}^{2}}}{g\left( {{c}_{j}} \right)},$

(6)

其中，g(c_j)为与聚类c_j中样本个数n_j成正比的函数，可取g(c_j)=n_j.根据式(6)计算各聚类中心的适应度函数，然后再利用适应度函数均衡准则进行聚类，其他处理步骤与KM_MKM或KM_AMKM相同.

3 毫米波SAR图像分割实验 3.1 图像分割结果评价指标

除了目视判读，文献[12, 15]给出如下的评价指标：

$F\left( I \right)=\sqrt{R}\sum\limits_{i=1}^{R}{\frac{e_{i}^{2}}{\sqrt{{{A}_{i}}}},}$

(7)

其中，I表示原始图像，R表示分割的区域个数，A_i表示第i个区域的尺寸，e_i表示原始图像与分割图像在第i个区域内每个对应像素的欧几里德距离之和.

文献[16]对式(7)进行修正，给出如下评价指标：

$\begin{align} & F'\left( I \right)=\frac{1}{10\ 000\left( M\times N \right)} \\ & \sqrt{\sum\limits_{A=1}^{\text{Max}}{|R\left( A \right){{|}^{1+1/A}}}}\times \sum\limits_{i=1}^{R}{\frac{e_{i}^{2}}{\sqrt{{{A}_{i}}}}}, \\ \end{align}$

(8)

其中，M×N表示图像I的尺寸，R(A)表示尺寸为A的区域数，Max为最大尺寸的区域，1+1/A为加大了的小区域权值.

在对过分割和分割不足充分考虑的基础上，文献[17]提出如下的分割质量评价指标：

$\begin{align} & Q\left( I \right)=\frac{1}{10\ 000\left( M\times N \right)}\sqrt{R}\times \\ & {{\sum\limits_{i=1}^{R}{\left[ \frac{e_{i}^{2}}{1+\log {{A}_{i}}}+\frac{R\left( {{A}_{i}} \right)}{{{A}_{i}}} \right]}}^{2}}. \\ \end{align}$

(9)

本文对分割结果的评价将采用以上3个指标，指标数据越小，分割效果越好.

分割采用Intel(R) Core(TM) i3 CPU 550 @ 3.20 GHz，3.19 GHz，2.99 GB的内存，Microsoft Windows XP Professional Service Pack3系统，软件采用MATLAB 7.5.0(R2007b)，分割效率用分割所用时间来衡量.

3.2 毫米波SAR图像分割

分割图像为Ka波段SAR数据，分辨率0.3 m×0.3 m.图 1(a)为建筑物场景，图 1(b)为城区场景，图 1(c)为高架桥场景，图 1(d)为平地道路场景，图 1(a)和图 1(d)大小为600像素×600像素，图 1(b)和图 1(c)大小为2 000像素×2 000像素.表 1为4个场景对应的分割指标.

表 1中，改进适应度函数计算方法后，4种场景下KM_MKM、KM_AMKM分割所用的时间均比改进前有一定的减少.改进前后算法分割质量指标没有变化，这是由于改进前后算法基于的准则没有变化，都是“适应度函数均衡”准则，随着迭代的进行，最终聚类的结果是一样的，指标也一样.

4 结论

本文对几种K均值聚类算法进行分析，改进了动态K均值聚类算法中每个聚类适应度函数的计算方法，用样本数的正比函数对适应度函数进行一定的削弱与平均.毫米波SAR图像分割实验表明，改进适应度函数计算方法前后分割指标不变，分割质量相同，但是改进后算法运算时间有所减少，分割效率得到了提高.