CN110751652A - 基于巴氏距离和纹理模式度量的sar图像分割方法 - Google Patents

Abstract

本发明提出了一种基于巴氏距离及纹理模式度量的SAR图像分割方法，用于解决现有技术中存在的分割精度较低的技术问题，实现步骤为对SAR图像进行预处理；计算double格式SAR图像中每个像素点的巴氏距离；对double格式SAR图像进行初始分割；获取double格式SAR图像的中间分割图Y；计算中间分割图Y中相邻区域的统计相似性度量SSM(Y_r,Y_s)；计算中间分割图Y中相邻区域的纹理模式度量TPS(Y_r,Y_s)；计算Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)；获取double格式SAR图像的最终分割图。可用于灾害监测、地球变化研究、地图绘制、农作物生长评估、冰川分布观测等领域。

Description

基于巴氏距离和纹理模式度量的SAR图像分割方法

技术领域

本发明属于图像处理技术领域，涉及一种SAR图像分割方法，具体涉及一种基于巴氏距离和纹理模式度量的SAR图像分割方法，可用于灾害监测、地球变化研究、地图绘制、农作物生长评估、冰川分布观测等领域。

背景技术

合成孔径雷达SAR(Synthetic Aperture Radar)是一种相干成像雷达，具有分辨率高、全天候工作和穿透掩盖物等特点成为当前遥感观测的重要手段。SAR图像中包含着丰富的目标类别，随着SAR的应用领域的不断扩大，对SAR图像自动解译的要求也越来越高，而SAR图像分割一直是自动解译的核心问题和难点所在，SAR图像分割技术是将一幅SAR图像划分为可数的互不重叠的连通区域，它提供SAR图像中区域的封闭轮廓，为SAR图像分类识别和自动解译提供可靠的信息。然而SAR图像中存在的大量相干斑噪声，使每个像素与其真实值往往差距较大，因此常规的光学图像分割算法在SAR图像中的分割精度并不理想。

在SAR图像分割方法中，主要可分为基于模型优化的SAR图像分割方法和基于图像特征的SAR图像分割方法，在基于图像特征的SAR图像分割方法中，需要利用图像的边缘信息和区域信息，该方法通过边缘信息得到高质量的初始分割图，再利用区域信息构造相邻区域间的相似性度量，通过相似性度量及边缘惩罚项构造合并代价函数，最后利用区域合并迭代技术合并相似区域得到最终分割图，准确率、召回率以及区域覆盖率是衡量分割方法精度的直观指标，而合并代价函数包含的内容是影响最终分割精度的关键所在。

Z.J.Zhang,X.Pan,L.Cheng于2019年发表在期刊“IEEE Access”上，名称为“SARimage segmentation using hierarchical region merging with orientated edgestrength weighted Kuiper’s distance,”提出了一种利用定向边缘强度加权的Kuiper距离进行分层区域合并SAR图像分割方法(OESM)，该方法利用巴氏距离得到图像的边缘强度映射图，对边缘强度映射图进行了分水岭变换得到初始分割图，采取定向边缘强度加权的Kuiper距离构造统计相似性度量，利用公共边界长度构造边缘惩罚项，两者构造了新的代价函数，利用合并代价函数，随着逐渐增加的边缘惩罚强度更新层次区域合并准则，进行迭代合并得到最终分割图。该方法利用了定向边缘强度加权的Kuiper距离构造了统计相似性度量，相较于以往基于均值的统计相似度量，能够反映更多信息，一定程度上提升了SAR图像的分割精度。但该方法并不能表征图像的纹理模式，如两个具有相似的统计度量但纹理模式不同的区域，该方法并不能有效分割开，导致此方法对SAR图像分割的精度较低，无论是分割结果的准确率、召回率还是区域覆盖率都仍有改善的空间。

发明内容

本发明的目的在于克服上述现有技术存在的缺陷，提出了一种基于巴氏距离和纹理模式度量的SAR图像分割方法，用于解决现有技术中存在的分割精度较低的技术问题。

为实现上述目的，本发明采取的技术方案包括如下步骤：

(1)对SAR图像进行预处理：

读取像素值大小为M×N的待分割SAR图像，并将其转换为double格式，得到double格式SAR图像，其中，M≥250,N≥250；

(2)计算double格式SAR图像中每个像素点的巴氏距离：

(2a)对所提取的double格式SAR图像的直方图进行量化，得到量化后的包含Q个灰度级的SAR图像,其中Q≥10；

(2b)计算量化后的SAR图像中每个像素点在多尺度旋转双窗下的巴氏距离BD_θ,s(x_m,y_m)，得到M×N巴氏距离BD_θ,s(x,y)，BD_θ,s(x,y)＝{BD_θ,s(x₁,y₁),BD_θ,s(x₂,y₂)...,BD_θ,s(x_m,y_m),...,BD_θ,s(x_M×N,y_M×N)}，其中，θ,s分别表示多尺度旋转双窗的旋转角度和尺度，θ∈[0,π),θ旋转间隔为π/8，(x_m,y_m)表示第m个像素点的坐标，(x,y)＝{(x₁,y₁),(x₂,y₂)...,(x_m,y_m),..,(x_M×N,y_M×N)}；

(3)对double格式SAR图像进行初始分割：

(3a)根据每个像素点在多尺度旋转双窗下的巴氏距离BD_θ,s(x_m,y_m)，计算double格式SAR图像中每个像素点的边缘映射强度ESM_m，并对ESM_m进行阈值化处理，得到过滤后的边缘映射强度ESM'_m，M×N像素点的过滤后的边缘映射强度表示为ESM'，ESM'＝{ESM'₁,ESM'₂...,ESM'_m,...,ESM'_M×N}；

(3b)对过滤后的边缘映射强度ESM'进行分水岭变换，得到double格式SAR图像的包括n个区域的初始分割图X，X＝{X₁,X₂...,X_j,...,X_n}，其中n≥50，X_j表示X中的第j个区域；

(4)获取double格式SAR图像的中间分割图Y：

(4a)采用均值区域合并方法计算与区域X_j相邻的l个区域的代价

其中，

表示与区域X_j相邻的区域X_p的代价；

(4b)设立阈值H，并将X中像素点个数小于H的区域与

中最小代价对应的区域进行合并，得到double格式SAR图像的包括t个区域的中间分割图Y，Y＝{Y₁,Y₂...,Y_r,Y_s,...,Y_t},其中，Y_r和Y_s表示两个相邻区域，H≤80，t＜n；

(5)计算中间分割图Y中相邻区域的统计相似性度量SSM(Y_r,Y_s)：

(5a)对中间分割图Y的直方图进行量化，得到包含Q个灰度级的量化后的中间分割图Y'，Y'＝{Y'₁,Y'₂...,Y'_r,Y'_s,...,Y'_t}，Y'_r和Y'_s分别表示相邻区域Y_r和Y_s量化后的区域；

(5b)计算量化后的中间分割图Y'中相邻区域Y_r'和Y_s'第i个灰度级的频率并通过

计算Y_r和Y_s之间的巴氏距离BD(Y_r,Y_s)：

(5c)通过BD(Y_r,Y_s)计算Y_r和Y_s的统计相似性度量SSM(Y_r,Y_s)：

其中，#Y_r,#Y_s表示Y_r和Y_s包含的像素个数；

(6)计算中间分割图Y中相邻区域的纹理模式度量TPS(Y_r,Y_s)：

(6a)计算相邻区域Y_r和Y_s的灰度平均值

和

并利用a×a的矩形结构对Y_r和Y_s进行腐蚀，得到区域Y_r ^*和Y_s ^*，其中a≥3；

(6b)计算区域Y_r的空间协方差矩阵

和区域Y_s的空间协方差矩阵

其中，(x_r,y_r)表示Y_r ^*内像素点的坐标，(x_s,y_s)表示Y_s ^*内像素点的坐标，τ₀为空间协方差矩阵

和

的最大位移量，

τ_x为

和的横向位移量，τ_y为

和

的纵向位移量-τ₀≤τ_x,τ_y≤τ₀；

表示Y_r ^*所包含的像素个数，#Y_s ^*表示Y_s ^*所包含的像素个数；

(6c)计算两相邻区域Y_r和Y_s的纹理模式度量TPS(Y_r,Y_s)：

(7)计算Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)：

基于中间分割图Y中相邻区域Y_r和Y_s的统计相似性度量SSM(Y_r,Y_s)和纹理模式度量TPS(Y_r,Y_s)，计算Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)：

其中，

表示Y_r和Y_s之间的边缘惩罚项，

表示Y_r和Y_s的公共边界长度，η表示TPS(Y_r,Y_s)项的权系数，λ₂表示边缘惩罚项的权系数，

表示两区域的公共边界长度；

(8)获取double格式SAR图像的最终分割图：

(8a)设立合并门限Th；

(8b)利用Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)计算中间分割图Y中相邻区域Y_r和Y_s的代价v(Y_r,Y_s)，并将所有相邻区域的代价按照由小到大的顺序进行排列，得到代价集合v_order，v_order＝{v₁,v₂,v₃,...}；

(8c)判断v₁＜Th是否成立，若是，对v₁所对应的相邻区域进行合并，得到结果图Y_mid，并令Y_mid＝Y，执行步骤(8b)，否则，得到SAR图像的最终分割图。

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

第一，本发明最终分割结果所采用的合并代价函数是通过纹理模式度量和巴氏距离获取的，在代价函数中添加了纹理模式度量项，该项通过空间协方差矩阵的计算能够反映区域内像素点的空间排列特性，可以有效分割开具有相似的灰度分布但纹理模式不同的区域，与现有技术相比，有效地提高了分割结果的精度。

第二，本发明的合并代价函数中的统计相似性度量应用了巴氏距离，能很好的反映区域间的统计相似性，同时对巴氏距离进行了区域个数加权，可以有效合并小区域，防止区域间的过度分割，进一步提高了分割结果的精度。

附图说明

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

图2是本发明和现有技术分割精度的仿真结果对比图。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步的详细描述。

参照图1，一种基于巴氏距离和纹理模式度量的SAR图像分割方法，包括如下步骤：

步骤1)对SAR图像进行预处理：

读取像素值大小为M×N的待分割SAR图像，并将其转换为double格式，得到double格式SAR图像，其中，M≥250,N≥250。

在本实施例中，M＝501,N＝501，共M×N＝251001个像素点。转换成double格式可使巴氏距离的计算精度更高。

步骤2)计算double格式SAR图像中每个像素点的巴氏距离：

步骤2a)对double格式SAR图像的直方图进行量化，得到量化后的包含Q个灰度级的SAR图像,其中Q≥10；具体实施步骤为：

(2a1)计算所提取的double格式SAR图像的直方图中每个灰度级的频率

并计算前k个直方图所对应的灰度级频率的和其中1≤k≤256，q表示直方图中第q个灰度级，sum(q)表示直方图中灰度级为q的像素点的数量；

(2a2)计算满足

条件的最小灰度级k_min和最大灰度级k_max，并利用量化后的SAR图像中的第i个灰度级对区间[k_min,k_max]内的全部灰度级进行替换，得到量化后的包含Q个灰度级的SAR图像，其中，i＝{1,2...,Q}。

在本实施例中，Q＝10，量化后的SAR图像仅有10个灰度级。

步骤2b)计算量化后的SAR图像中每个像素点在多尺度旋转双窗下的巴氏距离BD_θ,s(x_m,y_m)，得到M×N巴氏距离BD_θ,s(x,y)，BD_θ,s(x,y)＝{BD_θ,s(x₁,y₁),BD_θ,s(x₂,y₂)...,BD_θ,s(x_m,y_m),...,BD_θ,s(x_M×N,y_M×N)}，其中，θ,s分别表示多尺度旋转双窗的旋转角度和尺度，θ∈[0,π),θ旋转间隔为π/8，(x_m,y_m)表示第m个像素点的坐标，(x,y)＝{(x₁,y₁),(x₂,y₂)...,(x_m,y_m),..,(x_M×N,y_M×N)}；

巴氏距离计算公式如下：

其中，h_u(x,y|s,θ)(i)和h_l(x,y|s,θ)(i)分别表示多尺度旋转双窗在量化后的SAR图像中第i个灰度级所对的频率。

在本实施例中，利用量化后的图像计算巴氏距离，可使巴氏距离反映的信息更加全面，多尺度旋转双窗形状均为矩形，间距为1个像素宽度，尺度数量为3，三种尺度的双窗长宽比例为17:9，面积比为1:4:16，围绕像素点(x_m,y_m)建立，即(x_m,y_m)是多尺度旋转双窗的中心，通过双窗的旋转以及尺度的变化计算(x_m,y_m)在每个角度和每个尺度下的巴氏距离BD_θ,s(x_m,y_m)。

步骤3)对double格式SAR图像进行初始分割：

步骤3a)计算double格式SAR图像中每个像素点的边缘映射强度ESM_m，并对ESM_m进行阈值化处理，得到过滤后的边缘映射强度ESM'_m，M×N像素点的过滤后的边缘映射强度表示为ESM'，ESM'＝{ESM'₁,ESM'₂...,ESM'_m,...,ESM'_M×N}；

每个像素点的边缘映射强度ESM_m计算公式如下：

其中，S,M,L表示多尺度旋转双窗的三种尺度，α(S)表示S尺度的权值，α(M)表示M尺度的权值，α(L)表示L尺度的权值且α(S)+α(M)+α(L)＝1。

在本实施例中，阈值化处理采用Savitzky-Golay滤波器，该滤波器可有效去除ESM_m中的伪极大值，使初始分割图更精确，α(S)＝0.5，α(M)＝0.3，α(L)＝0.2。

步骤3b)对过滤后的边缘映射强度ESM'进行分水岭变换，得到double格式SAR图像的包括n个区域的初始分割图X，X＝{X₁,X₂...,X_j,...,X_n}，其中n≥50，X_j表示初始分割图中的第j个区域；

在本实施例中，n＝1283，即初始分割图中共包含了1283个区域。

步骤4)获取double格式SAR图像的中间分割图Y：

步骤4a)采用均值区域合并方法计算与区域X_j相邻的l个区域的代价

其中，

表示与区域X_j相邻的区域X_p的代价，l≥1，区域X_j相邻的区域X_p的代价

的计算公式为：

其中，#X_j和#X_p分别表示区域X_j和X_p所包含的像素点个数，

和

分别表示区域X_j和X_p所包含像素点的均值。

步骤4b)设立阈值H，并将X中像素点个数小于H的区域与

中最小代价对应的区域进行合并，得到double格式SAR图像的包括t个区域的中间分割图Y，Y＝{Y₁,Y₂...,Y_r,Y_s,...,Y_t},其中，Y_r和Y_s表示两个相邻区域，H≤80，t＜n；

在本实施例中，H＝80,t＝687，由于微小区域包含的像素个数过少，反映的信息有限，采用均值区域合并方法可准确地合并，同时，强制合并可提高方法的整体运算效率。

步骤5)计算中间分割图Y中相邻区域的统计相似性度量SSM(Y_r,Y_s)：

步骤5a)对中间分割图Y的直方图进行量化，得到包含Q个灰度级的量化后的中间分割图Y'，Y'＝{Y'₁,Y'₂...,Y'_r,Y'_s,...,Y'_t}，Y'_r和Y'_s分别表示相邻区域Y_r和Y_s量化后的区域；

步骤5b)计算量化后的中间分割图Y'中相邻区域Y_r'和Y_s'第i个灰度级的频率

并通过

计算Y_r和Y_σ之间的巴氏距离BD(Y_r,Y_s)：

步骤5c)通过BD(Y_r,Y_s)计算Y_r和Y_s的统计相似性度量SSM(Y_r,Y_s)：

其中，#Y_r,#Y_s表示Y_r和Y_s包含的像素个数；

在本实施例中，SSM(Y_r,Y_s)是由BD(Y_r,Y_s)进行区域个数加权得到的，该加权可以促进小区域之间的合并，防止区域的过度分割。

步骤6)计算中间分割图Y中相邻区域的纹理模式度量TPS(Y_r,Y_s)：

步骤6a)计算相邻区域Y_r和Y_s的灰度平均值

和

并利用a×a的矩形结构对Y_r和Y_s进行腐蚀，得到区域Y_r ^*和Y_s ^*，其中a≥3；

在本实施例中，a＝5，对区域Y_r和Y_s进行腐蚀是为了保证在计算空间协方差矩阵时，所有位移后的像素点仍分别在区域Y_r和Y_s中。

步骤6b)计算区域Y_r的空间协方差矩阵

和区域Y_s的空间协方差矩阵

其中，(x_r,y_r)表示Y_r ^*内像素点的坐标，(x_s,y_s)表示Y_s ^*内像素点的坐标，τ₀为空间协方差矩阵

和

的最大位移量，

τ_x为

和

的横向位移量，τ_y为

和

的纵向位移量-τ₀≤τ_x,τ_y≤τ₀；

表示Y_r ^*所包含的像素个数，

表示Y_s ^*所包含的像素个数；

在本实施例中，τ₀＝2，

和

均为5×5的矩阵。

步骤6c)计算两相邻区域Y_r和Y_s的纹理模式度量TPS(Y_r,Y_s)：

在本实施例中，两相邻区域Y_r和Y_s的纹理模式度量TPS(Y_r,Y_s)是利用矩阵的二范数对空间协方差矩阵和

进行归一化得到的，TPS(Y_r,Y_s)∈[0,1)。

步骤7)计算Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)：

基于中间分割图Y中相邻区域Y_r和Y_s的统计相似性度量SSM(Y_r,Y_s)和纹理模式度量TPS(Y_r,Y_s)，计算Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)：

其中，表示Y_r和Y_s之间的边缘惩罚项，

表示Y_r和Y_s的公共边界长度，η表示TPS(Y_r,Y_s)项的权系数，λ₂表示边缘惩罚项的权系数，

表示两区域的公共边界长度。

在本实施例中，使用的最佳权系数为η＝2.5,λ₂＝0.9，代价函数v(Y_r,Y_s)融合了统计相似性度量SSM(Y_r,Y_s)，边缘惩罚项以及纹理模式度量TPS(Y_r,Y_s)。

利用区域像素个数加权的巴氏距离计算SSM(Y_r,Y_s)，区域间的巴氏距离是相较于均值能更好地反映区域灰度分布的统计量，同时进行了区域像素个数加权，防止了区域的过度分割。

利用空间协方差矩阵计算TPS(Y_r,Y_s)，纹理区域的特性是不能由一阶统计量表征的，而空间协方差矩阵在二维空域通过计算像素点之间的相似性，可以反映区域内像素的空间排列特性，因此利用空间协方差矩阵计算TPS(Y_r,Y_s)能很好地分割纹理区域。同时在构造v(Y_r,Y_s)时，对TPS(Y_r,Y_s)进行了非线性映射，可以在纹理区域放大TPS(Y_r,Y_s)的影响，进一步提高了分割结果的精度。

步骤8)获取double格式SAR图像的最终分割图：

步骤8a)设立合并门限Th，Th∈[1.5,4]；

步骤8b)利用Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)计算中间分割图Y中相邻区域Y_r和Y_s的代价v(Y_r,Y_s)，并将所有相邻区域的代价按照由小到大的顺序进行排列，得到代价集合v_order，v_order＝{v₁,v₂,v₃,...}；

步骤8c)判断v₁＜Th是否成立，若是，对v₁所对应的相邻区域进行合并，得到结果图Y_mid，并令Y_mid＝Y，执行步骤(8b)，否则，得到SAR图像的最终分割图。

在本实施例中，Th＝2.0，最终分割图包含了120个区域。

下面结合仿真实验对本发明的技术效果作进一步的说明：

1.仿真条件和内容：

本发明的仿真实验是在计算机配置为core i7-8700 3.20GHZ处理器，16G运行内存，WINDOWS10系统和计算机软件配置为MATLAB R2017a环境下进行的，利用本发明方法和现有技术在上述仿真条件下进行实验，即运用本发明和现有技术对一幅德国terrasar数据库中的SAR图像进行分割，分割结果如图2所示，并得到了如表1所示的仿真精度对比。

2.仿真结果分析：

参照图2，图2(a)是仿真实验所使用的德国terrasar数据库中的一幅SAR图像，该图像大小为501×501像素。

图2(b)是仿真实验所使用的SAR图像人工标记的标准区域边缘图的二值图像。

图2(c)是现有技术对SAR图像进行分割的结果图。

图2(d)是本发明对SAR图像进行分割的结果图。

表1.本发明和现有技术分割结果精度对比表：

方法	准确率	召回率	区域覆盖率
				现有技术	0.81	0.74	0.76
本发明	0.86	0.79	0.78

表1中，方法一栏中现有技术表示现有的OESW分割技术，利用了准确率，召回率以及区域覆盖率三种指标评价两种技术的分割精度，三种评价指标的值越大表示分割结果精度的越高。

从表格中可以看出，本发明相较于现有技术OESW，各个评价指标均有提升，其中，准确率提高了0.05，召回率提高了0.05，区域覆盖率提升了0.02，分割精度得到了全面的提升。

Claims (5)

1.一种基于巴氏距离和纹理模式度量的SAR图像分割方法，其特征在于，包括如下步骤：

(1)对SAR图像进行预处理：

读取像素值大小为M×N的待分割SAR图像，并将其转换为double格式，得到double格式SAR图像，其中，M≥250,N≥250；

(2)计算double格式SAR图像中每个像素点的巴氏距离：

(2a)对所提取的double格式SAR图像的直方图进行量化，得到量化后的包含Q个灰度级的SAR图像,其中Q≥10；

(2b)计算量化后的SAR图像中每个像素点在多尺度旋转双窗下的巴氏距离BD_θ,s(x_m,y_m)，得到M×N巴氏距离BD_θ,s(x,y)，BD_θ,s(x,y)＝{BD_θ,s(x₁,y₁),BD_θ,s(x₂,y₂)...,BD_θ,s(x_m,y_m),...,BD_θ,s(x_M×N,y_M×N)}，其中，θ,s分别表示多尺度旋转双窗的旋转角度和尺度，θ∈[0,π),θ旋转间隔为π/8，(x_m,y_m)表示第m个像素点的坐标，(x,y)＝{(x₁,y₁),(x₂,y₂)...,(x_m,y_m),..,(x_M×N,y_M×N)}；

(3)对double格式SAR图像进行初始分割：

(3a)根据每个像素点在多尺度旋转双窗下的巴氏距离BD_θ,s(x_m,y_m)，计算double格式SAR图像中每个像素点的边缘映射强度ESM_m，并对ESM_m进行阈值化处理，得到过滤后的边缘映射强度ESM'_m，M×N像素点的过滤后的边缘映射强度表示为ESM'，ESM'＝{ESM₁',ESM'₂...,ESM'_m,...,ESM'_M×N}；

(3b)对过滤后的边缘映射强度ESM'进行分水岭变换，得到double格式SAR图像的包括n个区域的初始分割图X，X＝{X₁,X₂...,X_j,...,X_n}，其中n≥50，X_j表示X中的第j个区域；

(4)获取double格式SAR图像的中间分割图Y：

(4a)采用均值区域合并方法计算与区域X_j相邻的l个区域的代价

其中，

表示与区域X_j相邻的区域X_p的代价；

(4b)设立阈值H，并将X中像素点个数小于H的区域与中最小代价对应的区域进行合并，得到double格式SAR图像的包括t个区域的中间分割图Y，Y＝{Y₁,Y₂...,Y_r,Y_s,...,Y_t},其中，Y_r和Y_s表示两个相邻区域，H≤80，t＜n；

(5)计算中间分割图Y中相邻区域的统计相似性度量SSM(Y_r,Y_s)：

(5a)对中间分割图Y的直方图进行量化，得到包含Q个灰度级的量化后的中间分割图Y'，Y'＝{Y₁',Y₂'...,Y_r',Y_s',...,Y_t'}，Y_r'和Y_s'分别表示相邻区域Y_r和Y_s量化后的区域；

(5b)计算量化后的中间分割图Y'中相邻区域Y_r'和Y_s'第i个灰度级的频率

并通过

计算Y_r和Y_s之间的巴氏距离BD(Y_r,Y_s)：

(5c)通过BD(Y_r,Y_s)计算Y_r和Y_s的统计相似性度量SSM(Y_r,Y_s)：

其中，#Y_r,#Y_s表示Y_r和Y_s包含的像素个数；

(6)计算中间分割图Y中相邻区域的纹理模式度量TPS(Y_r,Y_s)：

(6a)计算相邻区域Y_r和Y_s的灰度平均值

和

并利用a×a的矩形结构对Y_r和Y_s进行腐蚀，得到区域Y_r ^*和Y_s ^*，其中a≥3；

(6b)计算区域Y_r的空间协方差矩阵

和区域Y_s的空间协方差矩阵

其中，(x_r,y_r)表示Y_r ^*内像素点的坐标，(x_s,y_s)表示Y_s ^*内像素点的坐标，τ₀为空间协方差矩阵

和

的最大位移量，

τ_x为

和

的横向位移量，τ_y为

和

的纵向位移量-τ₀≤τ_x,τ_y≤τ₀；

表示Y_r ^*所包含的像素个数，#Y_s ^*表示Y_s ^*所包含的像素个数；

(6c)计算两相邻区域Y_r和Y_s的纹理模式度量TPS(Y_r,Y_s)：

(7)计算Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)：

基于中间分割图Y中相邻区域Y_r和Y_s的统计相似性度量SSM(Y_r,Y_s)和纹理模式度量TPS(Y_r,Y_s)，计算Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)：

其中，

表示Y_r和Y_s之间的边缘惩罚项，表示Y_r和Y_s的公共边界长度，η表示TPS(Y_r,Y_s)项的权系数，λ₂表示边缘惩罚项的权系数，

表示两区域的公共边界长度；

(8)获取double格式SAR图像的最终分割图：

(8a)设立合并门限Th；

(8b)利用Y中相邻区域Y_r和Y_s的代价函数v(Y_r,Y_s)计算中间分割图Y中相邻区域Y_r和Y_s的代价v(Y_r,Y_s)，并将所有相邻区域的代价按照由小到大的顺序进行排列，得到代价集合v_order，v_order＝{v₁,v₂,v₃,...}；

(8c)判断v₁＜Th是否成立，若是，对v₁所对应的相邻区域进行合并，得到结果图Y_mid，并令Y_mid＝Y，执行步骤(8b)，否则，得到SAR图像的最终分割图。

2.根据权利要求1所述的基于巴氏距离和纹理模式度量的SAR图像分割方法，其特征在于，步骤(2a)中所述的对所提取的double格式SAR图像的直方图进行量化，实现步骤为：

(2a1)计算所提取的double格式SAR图像的直方图中每个灰度级的频率

并计算前k个直方图所对应的灰度级频率的和其中1≤k≤256，q表示直方图中第q个灰度级，sum(q)表示直方图中灰度级为q的像素点的数量；

(2a2)计算满足

条件的最小灰度级k_min和最大灰度级k_max，并利用量化后的SAR图像中的第i个灰度级对区间[k_min,k_max]内的全部灰度级进行替换，得到量化后的包含Q个灰度级的SAR图像，其中，i＝{1,2...,Q}。

3.根据权利要求1所述的基于巴氏距离和纹理模式度量的SAR图像分割方法，其特征在于，步骤(2b)中所述的计算量化后的SAR图像中每个像素点在多尺度旋转双窗下的巴氏距离BD_θ,s(x_m,y_m)，计算公式为：

其中，h_u(x,y|s,θ)(i)和h_l(x,y|s,θ)(i)分别表示多尺度旋转双窗在量化后的SAR图像中第i个灰度级所对的频率。

4.根据权利要求1所述的基于巴氏距离和纹理模式度量的SAR图像分割方法，其特征在于，步骤(3a)中所述的基于巴氏距离计算double格式SAR图像的边缘映射强度ESM_m，计算公式为：

其中，α(S)、α(M)和α(L)分别表示多尺度旋转双窗中S尺度的权值、M尺度的权值和L尺度的权值，且α(S)+α(M)+α(L)＝1。

5.根据权利要求1所述的基于巴氏距离和纹理模式度量的SAR图像分割方法，其特征在于，步骤(4a)中所述的采用均值区域合并方法计算与区域X_j相邻的l个区域的代价

其中，与区域X_j相邻的区域X_p的代价的计算公式为：

其中，#X_j和#X_p分别表示区域X_j和X_p所包含的像素点个数，

和

分别表示区域X_j和X_p所包含像素点的均值。

Images