(辽宁工程技术大学 测绘与地理科学学院，辽宁 阜新 123000)

0 引言

主动遥感技术中的合成孔径雷达(synthetic aperture radar，SAR)具有全天时、全天候等光学传感器所不具有的图像采集能力，因此作为其信息载体的SAR图像已在环境监测，地球资源探测和灾害评估中得到了广泛的应用[1]。图像分割将图像细分为非重叠的均匀区域，是SAR图像处理和应用的基本步骤[2]。阈值法由于其简单直观的特点，在图像分割中占有重要地位[3]。在实际使用中，阈值的选取有多种方法，如直方图法[4]、最大熵法[5]、Otsu法[6]、模糊聚类法[7]等。其中，Otsu法具有原理简单、稳定有效的优点，是使用最广泛的阈值方法之一。Otsu法由OTSU[8]提出，利用图像的光谱测度直方图，使用穷举法选取使图像类间方差最大的阈值将图像划分为前景和背景2个部分。然而，由于图像中的对象具有各自不同的光谱分布，因此需要较多的阈值才能从其背景中分割出感兴趣的对象。与单阈值Otsu法相比，多阈值Otsu法可以更好地提取图像目标，然而，由于光谱测度直方图的零阶和一阶累积矩的大量重复计算，计算量显著增加[9]，特别是在需较多阈值的图像分割中更为突出。

Otsu多阈值图像分割可视为一个组合优化问题，其基本原理为在图像的光谱测度范围内，根据特定的准则和类属数搜索系列单阈值的最佳组合。目前，研究人员开始尝试将优化方法和Otsu多阈值分割方法相结合，以减少计算量，缩短运算时间，提高算法的鲁棒性。现代优化算法主要包括蚁群算法[10]、粒子群算法[11]、遗传算法(genetic algorithm，GA)[12]和禁忌搜索算法[13]等。其中，GA具有过程简单、寻优能力强、易与其他算法相结合等优点，是一种常用的优化算法。GA是模拟自然界的进化机制，对个体由确定的适应度函数逐代演化，以最终演化的个体作为问题的最优解[14]。近几十年来，研究人员已提出了几种基于GA的图像分割方法，如Lee等[15]提出的一种基于GA的图像分割自动阈值选取，以类间方差为进化方向，利用GA自动演化能量函数参数，使输出分割结果满足定义的感知能量。赵夫群等[16]提出的一种基于GA-Otsu法的图像阈值分割，其以图像类间方差为适应度函数；通过适应度值确定进一步的搜索方向和范围，最终确定最优阈值实现图像分割。但由于GA易陷入局部最优，出现早熟现象，且迭代过程中计算量大，搜索速度较慢，降低了该类方法的实际应用价值。

量子计算[17]可以解决传统计算领域的许多难题。Narayanan和Moore[18]首次将量子力学的概念和原理引入到GA中，提出了量子启发遗传算法，并成功应用于求解TSP(traveling salesman problem)问题。由于这种类型的算法不使用量子比特来编码染色体，因此它仍然属于传统意义上的GA。Han和Kim[19]将量子比特和量子门的概念引入到GA，并提出了量子遗传算法(quantum genetic algorithm，QGA)，其采用量子比特编码染色体增加种群多样性和全局寻优能力；采用量子门演化量子染色体加快算法收敛，有效地避免早熟现象。目前，QGA已成功用于许多领域，如信号处理[20]、机组组合问题[21]等。

然而，由于SAR的特殊成像机制，使SAR图像具有强相干斑噪声，降低了图像质量[1]。使得直接对SAR图像进行Otsu多阈值分割的结果中含有大量的杂散点[22]，这无疑会对后续的信息提取造成不可逆的错误。而提前对图像进行去噪则是一种消除这种现象的有效操作，但传统的去噪方法都存在不同程度的局限性，如Lee滤波[23]和曲波滤波法[24]在有限的迭代次数内很难去除原图中无用的纹理信息，无法起到去除大量杂散点的效果，而增加迭代次数则会导致对象的精细边缘全部被平滑掉，降低分割精度。SRAD滤波则可在精确地保持边缘的前提下以很少的迭代次数平滑掉纹理信息，最大限度地去除杂散点，特别适用于具有复杂精细结构的目标。但传统的SRAD由于无法自动控制迭代次数使得其会过度平滑掉部分精细结构。

综上，为了解决Otsu多阈值法对SAR图像分割具有的噪声敏感且计算量大的问题，本文采用图像间平均结构相似性指标来自适应地控制SRAD的迭代次数，以最小迭代次数达到既保持精细边缘又不过度平滑掉纹理结构的目的。后经自适应QGA优化阈值组合，获得最佳的图像分割效果。

1 算法描述

1.1 改进的SRAD滤波

设I0={I0(x，y)，(x，y)∈Ω}为定义在图像域Ω、大小为M×N的原始SAR图像，其中，(x，y)表示像素位置，I0(x，y)∈{0，…，L}为像素(x，y)强度的光谱测度，L为图像的光谱测度量化级。将I0(x，y)作为图像初始输入I(x，y；0)，通过偏微分方程建立SRAD模型，可表示为：

(1)

(2)

式中：q(x，y；t)为滤波的瞬时扩散系数，定义为：

(3)

(4)

q0(t)≈q0exp(-ρt)

(5)

式中：ρ为指数衰减常数；q0是为初始扩散阈值。

对SRAD模型，使用雅可比迭代法进行求解。假设较小的时间步长Δt和在坐标x和y方向上较小的空间步长h，则时间和空间坐标可被离散为：t=nΔt，n=0，1，2，…；x=ih，i=0，1，2，…，M′-1；y=jh，j=0，1，2，…，N′-1；M′h×N′h=M×N。

(6)

(7)

(8)

(9)

(10)

其中，j∈{0，…，N′-1}，i∈{0，…，M′-1}。

(11)

(12)

其满足对称边界条件，

(13)

其中，j∈{0，1，…，N′-1}，i∈{0，1，…，M′-1}。

图1 扩散示意图

使用前向差分来近似时间导数，则SRAD模型的数值近似为：

(14)

为了在平滑图像的同时保持图像的结构特征，在迭代过程中引入平均结构相似性指标(mean structural similarity index measurement，MSSIM)[25]作为定量的标准来控制迭代过程。表达式为：

(15)

(16)

|MSSIM(n)-MSSIM(1)|/MSSIM(1)≤ε

(17)

其中，ε为迭代终止阈值，一般取ε=0.5。

当满足式(17)时，得到平滑图像z={z(x，y)，(x，y)∈Ω }，其中，z(x，y)为像素(x，y)的光谱测度值，z(x，y)∈{1，…，L}。

1.2 基于自适应QGA的Otsu多阈值图像分割

利用QGA在函数寻优方面具有计算量小，收敛速度快等特点，将其在大小为L的搜索空间中，以图像类间方差最大为准则，优化得到C个最优阈值，并据此阈值将图像域Ω分割为C+1个同质区域，即Ω={Ωk，k=1，…，C+1}，其中，k为同质区域索引，Ωk表示第k个同质区域，满足Ωk∩ Ωk′=∅，∀k，k′∈{1，…， C+1}，k≠ k′， Ω1∪… ∪ ΩC+1=Ω。在寻优过程中，包含以下步骤：量子染色体编码、量子染色体测量和解码、适应度评价、量子染色体进化。

(1)量子染色体编码。QGA中的染色体用量子比特表示，一个量子比特的状态可以表示为：

|ψ〉=α|0〉+β|1〉

(18)

式中：α和β为概率幅，满足以下归一化条件：

|α|2+|β|2=1

(19)

式中：|α|2和|β|2分别表示量子比特坍缩到态|0〉和|1〉的概率。

(20)

(21)

对二进制串进行解码，得到阈值的十进制值，其原理为在第K个二进制串中，第d个二进制值1表示十进制值2d-1。则阈值与二进制串的关系可表示为：

(22)

(23)

式中：Tk是第k个阈值，其满足T1

(24)

整个图像光谱测度级的平均值μ为：

(25)

则基于自适应QGA的多阈值Otsu的适应度函数F(·)为：

(26)

图2 量子旋转门调整量子比特

(27)

式中：φK，v为第K个染色体中第v个量子比特的旋转角，其转换过程如下：

(28)

因此，第t代种群的更新过程可描述为：

Q(t+1)=G(t)×Q(t)

(29)

式中：G(t)为第t代的量子门；Q(t)和Q(t+1)分别为第t代及其进化后第t+1代的种群。

表1 旋转角设置方案

为了使量子染色体更快地向适应度高的个体方向进化和避免局部最优，当量子染色体与历史最优染色体的适应度差异较大时，旋转角将适当地增大，反之，将适当减小。同时，根据相邻两代量子染色体的差异，逐代地调整量子旋转角的大小。具体的调整策略为：

(30)

1.3 算法流程

算法流程图如图3所示。

图3 流程图

综上，算法的基本步骤可总结为：

①对原始SAR图像进行改进的SRAD滤波；

②定义Otsu多阈值分割时，图像的类间方差；

③利用量子比特将阈值组合编码为种群Q(t)；

④对种群Q(t)进行测量和解码，得到阈值的十进制值；

⑤将②的类间方差作为适应度函数，根据适应度大小评价各个个体的优劣，保留适应度最高的个体作为进化目标；

⑥利用量子旋转门G(t)作用于Q(t)，得到新种群Q(t+1)；

⑦对新种群Q(t+1)执行④、⑤和⑥；

⑧将算法执行到最大进化代数或相邻两代的最大类间方差的差值小于0.05作为终止条件，如满足终止条件，则输出分割结果，否则继续计算。

2 实验与分析

为了验证该方法的可行性和有效性，对模拟和真实SAR强度图像进行了实验。

2.1 模拟SAR图像

图4(a)示出了尺度为128像素×128像素的模拟图像模板，其中，数字Ⅰ-Ⅴ分别表示不同的均匀区域。表2示出了对应于模拟SAR图像的均匀区域的Gamma分布形状参数γ和尺度参数λ，得到的模拟SAR图像如图4(b)所示。

图4 模拟图像

表2 模拟图像各同质区域的Gamma分布参数

利用该方法对图4(b)进行实验，分割结果如图5所示，其中，图5(a)和图5(b)分别为SRAD滤波结果和分割结果。为了验证提出方法的优越性，利用文献[16]和文献[22]的方法对图4(b)进行实验，分割结果如图5(c)和图5(d)所示。从图5(b)～图5(d)的比较可以看出，文献[16]的分割结果含有较多的分割噪声，而文献[22]的方法在区域Ⅱ-V内的分割精度有所提高，但在区域Ⅰ内仍较低。另外，3种方法对图4(b)的运算时间分别为6.35 s、22.52 s、35.81 s，因此，本文方法能较快得到结果。

图5 分割结果

对于定性评估，提取结果的轮廓线(图6(a))，将其叠加在图4(b)的模拟图像上，如图6(b)所示。由叠加结果可以看出，所提轮廓线与实际轮廓线非常一致，表明本文方法能较好地处理各均匀区域的边界。

图6 视觉评价

对于定量评估，根据模板图像分别计算分割结果图5(b)～图5(d)的混淆矩阵[26]，如表3所示。在此基础上，计算产品精度、用户精度、总体精度和Kappa值(表4)，其中各项指标越高，表明其分割精度越高。

从表4可以看出，本文方法得到的产品精度，用户精度和总精度均最高，都达到91.1%以上，Kappa值为0.966，而一般分类器Kappa值达到0.8以上就被视为高质量分类器[26]。上述定量分析结果表明了该算法的有效性和准确性。

表3 模拟图像分割结果的混淆矩阵

表4 产品精度、用户精度、总精度和Kappa值

2.2 真实SAR图像

选取4幅真实SAR强度图像如图7(a)～图7(d)所示。图7(a)和图7(b)是分辨率为25 m，HV极化、尺度为128像素×128像素的Radarsat-Ⅱ图像，其中，图7(a)是河口图像，人为判读含有2个同质区域，图7(b)是城市图像，人为判读含有3个同质区域。图7(c)和图7(d)是分辨率为30 m、VV极化、尺度分别为128像素×128像素和256像素×256像素的Radarsat-I海冰图像，人为判读含有3个同质区域。

从图7可以看出，原始图像含有大量斑点噪声，为此，利用改进的SRAD滤波原始SAR图像，其设置参数Δt=0.1；ρ=0.3；h=1，由图8滤波结果可发现，图像斑点噪声被很好的滤除且边缘得到了较好的保留。

图7 真实SAR图像

图9显示了4幅真实SAR强度图像的分割结果。其中图9(a1)～图9(e1)为文献[16]的分割结果；图9(a2)～图9(e2)为文献[22]的分割结果；图9(a3)～图9(e3)为本文方法的分割结果。从图9可以看出，文献[16]的分割结果中含有较多的分割噪声，文献[22]的分割结果的分割噪声有所降低，但在边界处存在误分割现象，而本文方法在分割噪声和边界方面均能得到较好地实现。次外，对图7中4幅图像分割的3种方法的平均运行时间分别为23.22 s、35.42 s及6.03 s，该方法的速度明显高于其他2种方法。

对于定性评估，提取分割结果的同质区域的边界线并将其叠加在原始图像上，如图10所示。由图10中可以看出，文献[16]和文献[22]中获得的轮廓线与实际轮廓线不一致，而本文方法获得的轮廓线可精确地匹配实际轮廓。因此，验证了该方法的可行性和有效性。

图9 分割结果

图10 视觉评价

3 结束语

Otsu法是一个简单而有效的阈值法，具有良好的性能，被广泛应用于图像分割。然而，它对噪声敏感，并且当它从单阈值推广到多阈值时，因为需计算每一类的累积概率和平均值，存在计算量大，耗时过长的缺点，这使其在实际应用中受到限制。为了提高分割精度，同时实现更好的时间性能，本文提出一种结合降斑各向异性扩散和自适应量子遗传算法的Otsu多阈值SAR图像分割方法。该方法克服了传统方法对噪声敏感以及边缘细节丢失等问题，通过在边缘两侧以异性扩散来最大限度地去除同质区域的斑点噪声，并采用平均相似度有效的控制迭代过程，同时兼顾了去除斑点噪声与保持精度。将改进的自适应量子遗传算法作为Otsu多阈值的优化工具，快速精确地确定最佳阈值组合，使Otsu多阈值图像分割更加实用。但本文需人为确定类别数，使其分割自动化降低，因此，在未来工作中，拟对本文方法进行改进，快速地实现SAR图像的可变类分割。