CN104020449A - 一种合成孔径雷达干涉相位图滤波方法和设备 - Google Patents

Abstract

本发明实施例公开了一种SAR干涉相位图滤波方法和设备，该方法包括：对SAR复干涉图进行重叠分块，得到至少两个复数据块；对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；将所述最终更新后的复数据块拼接成更新后的复干涉图，并从所述更新后的复干涉图中提取更新后的干涉相位图。

Description

一种合成孔径雷达干涉相位图滤波方法和设备

技术领域

本发明涉及雷达信号处理领域，尤其涉及一种合成孔径雷达(SAR，Synthetic Aperture Radar)干涉相位图滤波方法和设备。

背景技术

基于马尔可夫随机场模型的图像滤波技术是目前图像滤波技术之一，它利用了图像的局部相关性，可以在有效滤除噪声的同时很好的保持图像的边缘特性，在图像恢复和图像分类等方面有广泛的应用。

而现有的SAR干涉相位图滤波方法通常是对图像进行多视处理，具体采用均匀的多视核进行处理，也就是给予每个像素相同的权重。这样没有考虑到局部的像素特性，会造成相位图的过度平滑，从而不能有效的保持边缘特性，导致接下来的相位解缠产生错误。

发明内容

有鉴于此，本发明实施例期望提供一种SAR干涉相位图滤波方法和设备，在保持良好的滤除噪声效果的同时，还能够有效地消除SAR复干涉图的残差点，从而使得滤波后的复干涉图很好地保持了原SAR复干涉图的边缘特性。

为达到上述目的，本发明的技术方案是这样实现的：

第一方面，本发明实施例提供了一种SAR干涉相位图滤波方法，所述方法包括：

对SAR复干涉图进行重叠分块，得到至少两个复数据块；

对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；

将所述最终更新后的复数据块拼接成更新后的复干涉图，并从所述更新后的复干涉图中提取更新后的干涉相位图。

根据第一种可能的实现方式，结合第一方面，所述对SAR复干涉图进行重叠分块处理，具体包括：

按照预先确定的复数据块尺寸和重复率对复干涉图进行重叠分块处理，并截取出相应的复数据块。

根据第二种可能的实现方式，结合第一方面，所述对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块，具体包括：

步骤a：通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素；

步骤b：根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数；

步骤c：根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，得到初步更新后的复数据块；

步骤d：将所述初步更新后的复数据块对应的相位图中的残差点个数与原复数据块对应的相位图中的残差点个数进行比较：

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数减少的时候，对所述更新后的复数据块重复执行步骤b和步骤c；

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数没有减少的时候，所述初步更新后的复数据块为所述最终更新后的复数据块。

根据第三种可能的实现方式，结合第二种可能的实现方式，通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素，具体包括：

提取所述复数据块对应的相位图；

根据残差点的定义从所述相位图中查找到残差点；

对所述复数据块中处于残差点相应位置的像素及以所述像素为中心的确定尺寸范围的像素进行标记，得到所述复数据块的标记像素。

根据第四种可能的实现方式，结合第二种或第三种可能的实现方式，根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数，具体包括：

根据所述复干涉图的特征以及确定的规则选择合适的马尔可夫随机场模型，并得到所述标记像素对应的能量函数

其中，σ表示估计值的标准差；z_p表示标识为p的标记像素的复数值；N_p表示以所述标记像素p为中心的固定阶数的邻近像素集合；z_t表示标记像素p的邻近像素t的复数值；表示标记像素p对邻近像素t的权值；Σ表示求和符号，|*|表示取模运算；

通过计算E(z_p)对z_p求偏导数为零的驻点得到所述标记像素对应的CMRF模型参数 ${{\widehat{\mathrm{\Λ}}}_P}^{*}={\left[\underset{s\∈I}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{sp}}{z}_s{Q}_s^H\right]\left[\underset{s\∈I}{\mathrm{\Σ}}{Q}_s{Q}_s^H\right]}^{-1},$

其中，I表示所述复数据块内所有像素点集合，Q_s＝[z_s1,z_s2,...,z_s24]^T为N_s中所有像素点复数值组成的列向量，与前述的N_p类似的，N_s为以像素s的为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；表示对Q_s进行共轭转置运算；ω_sp表示像素s相对于标记像素p的权重，具体的，ω_sp＝ω_sρ×ω_sc×ω_sd；其中，ω_sρ由所述复干涉图相关系数图来表示；ω_sc表示像素s和标记像素p的复数值所服从的统计特性之间的差异；ω_sd表示像素s和标记像素p之间的空间距离的关系。

根据第五种可能的实现方式，结合第四种可能的实现方式，所述根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，具体包括：

通过更新所述标记像素的复数值，其中，k+1表示当前的迭代计算，k表示上一次的迭代计算；μ为控制收敛速度的常数，通常选取为大于0小于1的自然数。

第二方面，本发明实施例提供了一种SAR干涉相位图滤波设备，所述设备包括：分块单元、更新单元、拼接单元和提取单元，其中，

所述分块单元用于，对SAR复干涉图进行重叠分块，得到至少两个复数据块；

所述更新单元用于，对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；

所述拼接单元用于，将所述最终更新后的复数据块拼接成更新后的复干涉图；

所述提取单元用于，从所述更新后的复干涉图中提取更新后的干涉相位图。

根据第一种可能的实现方式，结合第二方面，所述分块单元具体用于，按照预先确定的复数据块尺寸和重复率对复干涉图进行重叠分块处理，并截取出相应的复数据块。

根据第二种可能的实现方式，结合第二方面，所述更新单元包括：第一获取子单元、第二获取子单元、更新子单元和比较子单元，其中，

所述第一获取子单元用于，通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素；

所述第二获取子单元用于，根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数；

所述更新子单元用于，根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，得到初步更新后的复数据块；

所述比较子单元用于，将所述初步更新后的复数据块对应的相位图中的残差点个数与原复数据块对应的相位图中的残差点个数进行比较：

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数减少的时候，对所述更新后的复数据块触发第二获取子单元和更新子单元；

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数没有减少的时候，所述初步更新后的复数据块为所述最终更新后的复数据块。

根据第三种可能的实现方式，结合第二种可能的实现方式，所述第一获取子单元包括：提取模块、查找模块和标记模块，其中，

所述提取模块用于，提取所述复数据块对应的相位图；

所述查找模块用于，根据残差点的定义从所述相位图中查找到残差点；

所述标记模块用于，对所述复数据块中处于残差点相应位置的像素及以所述像素为中心的确定尺寸范围的像素进行标记，得到所述复数据块的标记像素。

根据第四种可能的实现方式，结合第二种或第三种可能的实现方式，所述第二获取单元具体用于，

根据所述复干涉图的特征以及确定的规则选择合适的马尔可夫随机场模型，并得到所述标记像素对应的能量函数

其中，σ表示估计值的标准差；z_p表示标识为p的标记像素的复数值；N_p表示以所述标记像素p为中心的固定阶数的邻近像素集合；z_t表示标记像素p的邻近像素t的复数值；表示标记像素p对邻近像素t的权值；Σ表示求和符号，|*|表示取模运算；

通过计算E(z_p)对z_p求偏导数为零的驻点得到所述标记像素对应的CMRF模型参数 ${{\widehat{\mathrm{\Λ}}}_P}^{*}={\left[\underset{s\∈I}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{sp}}{z}_s{Q}_s^H\right]\left[\underset{s\∈I}{\mathrm{\Σ}}{Q}_s{Q}_s^H\right]}^{-1},$

其中，I表示所述复数据块内所有像素点集合，Q_s＝[z_s1,z_s2,...,z_s24]^T为N_s中所有像素点复数值组成的列向量，与前述的N_p类似的，N_s为以像素s的为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；表示对Q_s进行共轭转置运算；ω_sp表示像素s相对于标记像素p的权重，具体的，ω_sp＝ω_sρ×ω_sc×ω_sd；其中，ω_sρ由所述复干涉图相关系数图来表示；ω_sc表示像素s和标记像素p的复数值所服从的统计特性之间的差异；ω_sd表示像素s和标记像素p之间的空间距离的关系。

根据第五种可能的实现方式，结合第四种可能的实现方式，所述更新子单元具体用于，通过更新所述标记像素的复数值，其中，k+1表示当前的迭代计算，k表示上一次的迭代计算；μ为控制收敛速度的常数，通常选取为大于0小于1的自然数。

本发明实施例提供了一种SAR干涉相位图滤波方法和设备，通过对SAR复干涉图进行重叠分块以及在模型参数估计时进行加权参数加权估计，从而能够在保证了良好的滤除噪声效果的同时，还能够有效地消除SAR复干涉图的残差点，从而使得滤波后的复干涉图很好地保持了原SAR复干涉图的边缘特性。

附图说明

图1为本发明实施例提供的一种SAR干涉相位图滤波方法的流程示意图；

图2为本发明实施例提供的一种对每一个复数据块中的标记像素进行更新的流程示意图；

图3为本发明实施例的SAR复干涉图；(a)为无噪声环境下的SAR复干涉图的相位图，(b)复高斯白噪声环境下的SAR复干涉图的幅度图，(c)为复高斯白噪声环境下的SAR复干涉图的相位图，(d)为复高斯白噪声环境下的SAR复干涉图的残差点分布图；

图4为采用无重叠分块对图3(b)和图3(c)所示的高斯白噪声环境下的SAR复干涉图进行滤波后效果图：(a)为采用无重叠分块进行滤波后的相位图，(b)为采用无重叠分块进行滤波后的残差点分布图；

图5为采用在模型参数估计时无加权参数加权对图3(b)和图3(c)所示的高斯白噪声环境下的SAR复干涉图进行滤波后效果图：(a)为采用在模型参数估计时无加权参数加权进行滤波后的相位图，(b)为采用在模型参数估计时无加权参数加权进行滤波后的残差点分布图；

图6为采用本发明实施例的技术方案对图3(b)和图3(c)所示的高斯白噪声环境下的SAR复干涉图进行滤波后效果图：(a)为采用本发明实施例的技术方案进行滤波后的相位图，(b)为采用本发明实施例的技术方案进行滤波后的残差点分布图；

图7为本发明实施例提供的一种SAR干涉相位图滤波设备的结构示意图；

图8为本发明实施例提供的另一种SAR干涉相位图滤波设备的结构示意图；

图9为本发明实施例提供的一种SAR干涉相位图滤波设备的硬件结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

由于本发明实施例是针对SAR干涉相位图进行处理，所以在没有特殊说明的情况下，图像的尺寸单位用像素表示，以下不再赘述。

如图1所示，为本发明实施例提供的一种SAR干涉相位图滤波方法的流程示意图，该方法可以包括：

S101：对SAR复干涉图进行重叠分块处理，得到至少两个复数据块；

可以理解的，由于在一幅复干涉图中，满足相同统计分布的像素点的比例比较小，因此需要将复干涉图划分为若干个对复数据块，以使得每个复数据块内的大多数像素满足相同的统计分布；

但是由于复数据块边缘处的像素点与复数据块内其他像素点的平均距离较大，并且在大多数情况下，复数据块边缘处的像素点与复数据块内其他大部分像素服从不同的统计分布，所以，会造成对复数据块边缘像素点的处理常常不正确；

因此，在本发明实施例中，首先对复干涉图进行重叠分块处理，需要说明的是，具体的复数据块的大小既要足够小，以使得复数据块中的大部分像素点满足同一种统计分布；又要足够大，以使得有充足像素点作为样本进行统计估计的运算。结合以上两点考虑，本实施例中，经验性地将每个复数据块的大小设置为64×64，相邻两个复数据块的重复率设置为50％，也就是说，相邻两个复数据块之间有32×64大小的区域是相同的。

需要说明的是，重复率的选择既要足够大以消除复数据块边缘的估计不准，又不能过大而影响对图像的处理效率，所以通常会按照预先确定的复数据块尺寸和重复率对复干涉图进行重叠分块处理，截取出相应的复数据块用来进行后续处理。

S102：对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；

需要说明的是，本发明实施例中，对每一个复数据块通过加权参数估计的方法进行更新的都是一致的，因此，为了简要清楚的对技术方案进行说明，本实施例仅对一个复数据块的处理过程进行详细说明，本领域技术人员可以很容易的将其应用在对其他复数据块的处理过程中，本实施例在此不再赘述；而且具体对复数据块的处理过程可以是同时对多个复数据块进行并行处理，也可以是对所有复数据块依次进行的串行处理，本实施例对此也不做具体的限定。

示例性的，如图2所示：对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块，具体可以包括：

S1021：通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素；

进一步的，S1021的实现过程可以是：

首先提取所述复数据块对应的相位图；接着根据残差点的定义从所述相位图中查找到残差点；最后对所述复数据块中处于残差点相应位置的像素及以所述像素为中心的确定尺寸范围的像素进行标记，得到所述复数据块的标记像素。

需要说明的是，残差点的定义为本领域技术人员的公知常识，在此不再赘述，在本实施例中，当相位图中任意一点A及其相邻的右侧点、下侧点以及右下点按照顺时针方向依次进行相位相减所得到的四个差值的和不为零的时候，该点A即为一个残差点；

还需要说明的是，尽管标记像素的个数可能不止一个，但是与复数据块类似的，本实施例对所有标记像素的处理过程都是一致的，因此本实施例仅对一个标记像素进行处理的过程进行详细描述，本领域技术人员可以很容易的将其应用在对其他标记像素的处理过程中，本实施例在此不再赘述；而且具体对标记像素的处理过程可以是同时对多个标记像素进行并行处理，也可以是对所有标记像素依次进行的串行处理，本实施例对此也不做具体的限定。

S1022：根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场(CMRF，Complex-valued Markov Random Field)模型参数；

具体的，S1022的详细过程如下：

首先，可以根据所述复干涉图的特征以及预先确定的规则选择合适的马尔可夫随机场模型，并得到所述标记像素对应的能量函数E(z_p)，具体的：

$E\left(z_p\right)=\frac{1}{2{\mathrm{\σ}}^2}\left({|z_p-\underset{t\∈N_p}{\mathrm{\Σ}}{\mathrm{\Λ}}_{\mathrm{pt}}^{*}{z}_t|}^2\right)$

其中，σ表示估计值的标准差，在E(z_p)中，σ相对于z_p为常数；z_p表示标识为p的标记像素的复数值；N_p表示以所述标记像素p为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；z_t表示标记像素p的邻近像素t的复数值；表示标记像素p对邻近像素t的权值；Σ表示求和符号，|*|表示取模运算；

其次，通过计算E(z_p)对z_p求偏导数为零的驻点得到所述标记像素对应的CMRF模型参数，对应的计算过程如下：

通过计算并且将计算结果用向量的形式表示，可以得到 ${{\widehat{\mathrm{\Λ}}}_P}^{*}={\left[\underset{s\∈I}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{sp}}{z}_s{Q}_s^H\right]\left[\underset{s\∈I}{\mathrm{\Σ}}{Q}_s{Q}_s^H\right]}^{-1};$

其中，I表示所述复数据块内所有像素点集合，Q_s＝[z_s1,z_s2,...,z_s24]^T为N_s中所有像素点复数值组成的列向量，与前述的N_p类似的，N_s为以像素s的为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；表示对Q_s进行共轭转置运算；ω_sp表示像素s相对于标记像素p的权重，具体的，ω_sp＝ω_sρ×ω_sc×ω_sd；

其中，ω_sρ由所述复干涉图质量图确定，优选的，本实施例采用所述复干涉图的相关系数图来表示所述复干涉图的质量图，可以理解的，由于质量图由相关系数图表示，因此ω_sρ的取值范围为0到1之间的自然数，当没有合适的质量图的时候，也可令ω_sρ为1；ω_sc表示像素s和标记像素p的复数值所服从的统计特性之间的差异，具体在本实施例中，当|z_p+z_s|≤|z_p|时，可以认为像素s和标记像素p的复数值服从相同统计分布，此时ω_sc＝1，当|z_p+z_s|≥|z_p|时，可以认为像素s和标记像素p的复数值服从不同的统计分布，此时ω_sc＝0；ω_sd表示像素s和标记像素p之间的空间距离的关系，具体的，像素s和标记像素p之间的空间距离越大，ω_sd越小；反之，像素s和标记像素p之间的空间距离越小ω_sd越大，优选的，本实施例中，优选的ω_sd与像素s和标记像素p之间的空间距离成反比关系，且ω_sd的取值范围也在0和1之间。

S1023：根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，得到初步更新后的复数据块；

具体的，S1023可以通过 $z_p^{(k+1)}=z_p^{\left(k\right)}-\mathrm{\μ}\·{\mathrm{\Δz}}_p^{\left(k\right)}$ 来实现，其中， ${\mathrm{\Δz}}_p^{\left(k\right)}=z_p^{\left(k\right)}-{\widehat{\mathrm{\Λ}}}_p^{*\left(k\right)}{Q}_p^{\left(k\right)},$ k+1表示当前的迭代计算，k表示上一次的迭代计算；μ为控制收敛速度的常数，通常选取为大于0小于1的自然数。

通过S1022和S1023对复数据块中每一个标记像素的复数值进行更新之后，就可以得到初步更新后的复数据块。

S1024：将所述初步更新后的复数据块对应的相位图中的残差点个数与原复数据块对应的相位图中的残差点个数进行比较；

具体的，可以根据S1021中的具体描述来查找所述初步更新后的复数据块对应的相位图中的残差点，在此不再赘述，从而就能够得到所述初步更新后的复数据块对应的相位图中的残差点个数。而且根据比较结果的不同，具体的处理过程也相应地会有区别，即：

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数减少的时候，对所述更新后的复数据块重复执行S1021至S1023；可以理解的，在此种情况下，说明所述初步更新后的复数据块对应的相位图中的残差点个数还没有稳定下来，也就是所述初步更新后的复数据块还可以继续的按照S1021至S1023的步骤进行更新；

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数没有减少的时候，S1025：所述初步更新后的复数据块为所述最终更新后的复数据块；可以理解的，在此种情况下，说明所述初步更新后的复数据块对应的相位图中的残差点个数已经稳定下来，此时的所述初步更新后的复数据块也就是所述最终更新后的复数据块。

S103：将所述最终更新后的复数据块拼接成更新后的复干涉图，并从所述更新后的复干涉图中提取更新后的干涉相位图；

示例性的，在得到所有的最终更新后的复数据块之后，可以将所有最终更新后的复数据块进行拼接，从而得到更新后的复干涉图，然后从所述更新后的复干涉图中提取更新后的干涉相位图；具体的实现过程为本领域技术人员的常规实验手段，在此不再赘述。

由上述S101至S103的实现过程可以看出，本实施例所采用的技术方案与传统基于MRF模型的滤波方法的主要区别在于：本实施例所采用的技术方案通过对复数据进行有重叠的分块处理将块边缘处像素值估计得更加准确，并且在模型参数估计中引入了权重，增加了参数估计的准确性，从而能够在有效的滤除相位噪声的同时很好的保持边缘特性。

为了体现本实施例所采用的技术方案的有效性，本发明实施例还按照S101至S103的实现过程对具体的SAR复干涉图进行处理；同时，为了说明本实施例所采用的技术方案的优势，本发明实施例还将无重叠分块以及在模型参数估计时无加权参数加权来对如图3所示的具体的SAR复干涉图进行滤波，并将上述两种滤波方式的结果与本实施例所采用的技术方案进行滤波的结果作为比对；

需要说明的是，无重叠分块以及在模型参数估计时无加权参数加权这两种滤波方式为本领域的公知常识，具体的实现过程可参考文献1(AB.Suksmono,A.Hirose.Adaptive Noise Reduction of InSAR Images Based on a Complex-ValuedMRF Model and Its Application to Phase Unwrapping Problem.IEEE Trans,Geosciand Remote Sensing,vol.40,no.3,2002:pp699-709)；还需要说明的是，由于SAR复干涉图的各像素均为复数值，因此，可以通过SAR复干涉图的幅度图和相位图来表示所述SAR复干涉图，在图3中，(a)为无噪声环境下的SAR复干涉图的相位图，(b)复高斯白噪声环境下的SAR复干涉图的幅度图，(c)为复高斯白噪声环境下的SAR复干涉图的相位图，(d)为复高斯白噪声环境下的SAR复干涉图的残差点分布图；

滤波效果的比对可以直观地如图4至图6所示，图4为采用无重叠分块对图3(b)和图3(c)所示的高斯白噪声环境下的SAR复干涉图进行滤波后效果图，其中，图4(a)为采用无重叠分块进行滤波后的相位图，图4(b)为采用无重叠分块进行滤波后的残差点分布图；图5为采用在模型参数估计时无加权参数加权对图3(b)和图3(c)所示的高斯白噪声环境下的SAR复干涉图进行滤波后效果图，其中，图5(a)为采用在模型参数估计时无加权参数加权进行滤波后的相位图，图5(b)为采用在模型参数估计时无加权参数加权进行滤波后的残差点分布图。图6为采用本发明实施例的技术方案对图3(b)和图3(c)所示的高斯白噪声环境下的SAR复干涉图进行滤波后效果图，其中，图6(a)为采用本发明实施例的技术方案进行滤波后的相位图，图6(b)为采用本发明实施例的技术方案进行滤波后的残差点分布图；从图4至图6的直观对比可以看出，图6(a)所示的采用本发明实施例的技术方案进行滤波后的相位图中的相位偏差现象要少于图4(a)和图5(a)中所示的相位图中的相位偏差现象说明了本发明实施例的技术方案能够较好的保持边缘特性；与此同时，6(b)所示的采用本发明实施例的技术方案进行滤波后的残差点的数量也要明显少于图4(b)和图5(b)中所示的残差点数量，因此，可以认为本发明实施例的技术方案相比于无重叠分块以及在模型参数估计时无加权参数加权两个技术方案来说，更有优势。

从具体的数据上来看，滤波结果中的残差点数统计、残差点数减少的百分比以及均方根误差(RMS，Root Mean Square)的比对如表1所示：

表1

从表1中可以明显的看出本发明实施例的技术方案相比于无重叠分块以及在模型参数估计时无加权参数加权两个技术方案来说，在滤除噪声的效果接近的情况下，更好的对残差点进行消除，从而也就更好的保持了原SAR复干涉图的边缘特性。

本实施例提供了一种SAR干涉相位图滤波方法，通过对SAR复干涉图进行重叠分块以及在模型参数估计时进行加权参数加权估计，从而能够在保证了良好的滤除噪声效果的同时，还能够有效地消除SAR复干涉图的残差点，从而使得滤波后的复干涉图很好地保持了原SAR复干涉图的边缘特性。

参见图7，为本发明实施例提供的一种SAR干涉相位图滤波设备70，该设备70可以包括：分块单元701、更新单元702、拼接单元703和提取单元704，其中，

所述分块单元701用于，对SAR复干涉图进行重叠分块，得到至少两个复数据块；

所述更新单元702用于，对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；

所述拼接单元703用于，将所述最终更新后的复数据块拼接成更新后的复干涉图；

所述提取单元704用于，从所述更新后的复干涉图中提取更新后的干涉相位图。

在本发明实施例中，分块单元701对复干涉图进行重叠分块处理，每个复数据块的大小为64×64，相邻两个复数据块的重复率为50％，也就是说，相邻两个复数据块之间有32×64大小的区域是相同的。

需要说明的是，重复率的选择既要足够大以消除复数据块边缘的估计不准，又不能过大而影响对图像的处理效率，所以分块单元701具体用于，按照预先确定的复数据块尺寸和重复率对复干涉图进行重叠分块处理，并截取出相应的复数据块。

需要说明的是，本发明实施例的SAR干涉相位图滤波设备70对每一个复数据块通过加权参数估计的方法进行更新的都是一致的，因此，为了简要清楚的对技术方案进行说明，本实施例仅对一个复数据块的处理过程进行详细说明，本领域技术人员可以很容易的将其应用在对其他复数据块的处理过程中，本实施例在此不再赘述；而且具体对复数据块的处理过程可以是同时对多个复数据块进行并行处理，也可以是对所有复数据块依次进行的串行处理，本实施例对此也不做具体的限定。

示例性的，参见图8，所述更新单元702包括：第一获取子单元7021、第二获取子单元7022、更新子单元7023和比较子单元7024，其中，

所述第一获取子单元7021用于，通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素；

进一步的，参见图8，所述第一获取子单元7021包括：提取模块70211、查找模块70212和标记模块70213，其中，

所述提取模块70211用于，提取所述复数据块对应的相位图；

所述查找模块70212用于，根据残差点的定义从所述相位图中查找到残差点；

所述标记模块70213用于，对所述复数据块中处于残差点相应位置的像素及以所述像素为中心的确定尺寸范围的像素进行标记，得到所述复数据块的标记像素。

需要说明的是，残差点的定义为本领域技术人员的公知常识，在此不再赘述，在本实施例中，当相位图中任意一点A及其相邻的右侧点、下侧点以及右下点按照顺时针方向依次进行相位相减所得到的四个差值的和不为零的时候，该点A即为一个残差点；

还需要说明的是，尽管标记像素的个数可能不止一个，但是与复数据块类似的，本实施例的SAR干涉相位图滤波设备70对所有标记像素的处理过程都是一致的，因此本实施例仅对一个标记像素进行处理的过程进行详细描述，本领域技术人员可以很容易的将其应用在对其他标记像素的处理过程中，本实施例在此不再赘述；而且具体对标记像素的处理过程可以是同时对多个标记像素进行并行处理，也可以是对所有标记像素依次进行的串行处理，本实施例对此也不做具体的限定。

所述第二获取子单元7022用于，根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数；

进一步的，第二获取单元7022具体用于，

根据所述复干涉图的特征以及确定的规则选择合适的马尔可夫随机场模型，并得到所述标记像素对应的能量函数

其中，σ表示估计值的标准差，在E(z_p)中，σ相对于z_p为常数；z_p表示标识为p的标记像素的复数值；N_p表示以所述标记像素p为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；z_t表示标记像素p的邻近像素t的复数值；表示标记像素p对邻近像素t的权值；Σ表示求和符号，|*|表示取模运算；

通过计算E(z_p)对z_p求偏导数为零的驻点得到所述标记像素对应的CMRF模型参数 ${{\widehat{\mathrm{\Λ}}}_P}^{*}={\left[\underset{s\∈I}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{sp}}{z}_s{Q}_s^H\right]\left[\underset{s\∈I}{\mathrm{\Σ}}{Q}_s{Q}_s^H\right]}^{-1},$

其中，I表示所述复数据块内所有像素点集合，Q_s＝[z_s1,z_s2,...,z_s24]^T为N_s中所有像素点复数值组成的列向量，与前述的N_p类似的，N_s为以像素s的为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；表示对Q_s进行共轭转置运算；ω_sp表示像素s相对于标记像素p的权重，具体的，ω_sp＝ω_sρ×ω_sc×ω_sd；

其中，ω_sρ由所述复干涉图质量图确定，优选的，本实施例采用所述复干涉图的相关系数图来表示所述复干涉图的质量图，可以理解的，由于质量图由相关系数图表示，因此ω_sρ的取值范围为0到1之间的自然数，当没有合适的质量图的时候，也可令ω_sρ为1；ω_sc表示像素s和标记像素p的复数值所服从的统计特性之间的差异，具体在本实施例中，当|z_p+z_s|≤|z_p|时，可以认为像素s和标记像素p的复数值服从相同统计分布，此时ω_sc＝1，当|z_p+z_s|≥|z_p|时，可以认为像素s和标记像素p的复数值服从不同的统计分布，此时ω_sc＝0；ω_sd表示像素s和标记像素p之间的空间距离的关系，具体的，像素s和标记像素p之间的空间距离越大，ω_sd越小；反之，像素s和标记像素p之间的空间距离越小ω_sd越大，优选的，本实施例中，优选的ω_sd与像素s和标记像素p之间的空间距离成反比关系，且ω_sd的取值范围也在0和1之间。

所述更新子单元7023用于，根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，得到初步更新后的复数据块；

所述比较子单元7024用于，将所述初步更新后的复数据块对应的相位图中的残差点个数与原复数据块对应的相位图中的残差点个数进行比较；具体的，比较子单元7024可以通过触发第一获取子单元7021来查找所述初步更新后的复数据块对应的相位图中的残差点，在此不再赘述，从而就能够得到所述初步更新后的复数据块对应的相位图中的残差点个数。而且根据比较结果的不同，具体的处理过程也相应地会有区别，即：

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数减少的时候，对所述更新后的复数据块触发第二获取子单元7022和更新子单元7023；可以理解的，在此种情况下，说明所述初步更新后的复数据块对应的相位图中的残差点个数还没有稳定下来，也就是所述初步更新后的复数据块还可以通过继续触发第二获取子单元7022和更新子单元7023进行更新；

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数没有减少的时候，所述初步更新后的复数据块为所述最终更新后的复数据块。可以理解的，在此种情况下，说明所述初步更新后的复数据块对应的相位图中的残差点个数已经稳定下来，此时的所述初步更新后的复数据块也就是所述最终更新后的复数据块。

示例性的，在得到所有的最终更新后的复数据块之后，拼接单元703可以将所有最终更新后的复数据块进行拼接，从而得到更新后的复干涉图，然后提取单元704从所述更新后的复干涉图中提取更新后的干涉相位图；具体的实现过程为本领域技术人员的常规实验手段，在此不再赘述。

为了体现本实施例所采用的技术方案的有效性，本发明实施例通过SAR干涉相位图滤波设备70对如图3所示的具体的SAR复干涉图进行处理；同时，为了说明本实施例所采用的技术方案的优势，本发明实施例还将无重叠分块以及在模型参数估计时无加权参数加权来对如图3所示的具体的SAR复干涉图进行滤波，并将上述两种滤波方式的结果与本实施例所采用的技术方案进行滤波的结果作为比对，比对的结果如前述实施例中所述，在此不再赘述。

值得注意的是，如图9所示，本发明实施例提供的一种SAR干涉相位图滤波设备70，可以包括至少一个处理器901、至少一个存储器902和总线903，其中，处理器901、存储器902通过总线903连接并完成相互间的通信。

该总线903可以是工业标准体系结构(Industry Standard Architecture，简称为ISA)总线、外部设备互连(Peripheral Component，简称为PCI)总线或扩展工业标准体系结构(Extended Industry Standard Architecture，简称为EISA)总线等。该总线903可以分为地址总线、数据总线、控制总线等。为便于表示，图9中仅用一条粗线表示，但并不表示仅有一根总线或一种类型的总线。

存储器902用于存储可执行程序代码，该程序代码包括计算机操作指令。存储器902可能包含高速RAM存储器，也可能还包括非易失性存储器(non-volatile memory)，例如至少一个磁盘存储器。存储器902中存储：操作系统以及用于实现前述实施例的方法或功能模块的程序代码和应用程序。操作系统用于控制和实现处理单元执行的处理功能。应用程序包含程序代码，如字处理软件、email软件。

处理器901可能是一个中央处理器(Central Processing Unit，简称为CPU)，或者是特定集成电路(Application Specific Integrated Circuit，简称为ASIC)，或者是被配置成实施本发明实施例的一个或多个集成电路。

处理器901用于调用存储器902中的实现前述实施例的方法或功能模块的程序代码和应用程序，以执行以下操作：

对SAR复干涉图进行重叠分块，得到至少两个复数据块；

对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；

将所述最终更新后的复数据块拼接成更新后的复干涉图，并从所述更新后的复干涉图中提取更新后的干涉相位图。

进一步的，处理器901还用于调用存储器902中的实现前述实施例的方法或功能模块的程序代码和应用程序，以执行以下操作：按照预先确定的复数据块尺寸和重复率对复干涉图进行重叠分块处理，并截取出相应的复数据块。

进一步的，处理器901还用于调用存储器902中的实现前述实施例的方法或功能模块的程序代码和应用程序，以执行以下操作：

步骤a：通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素；

步骤b：根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数；

步骤c：根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，得到初步更新后的复数据块；

步骤d：将所述初步更新后的复数据块对应的相位图中的残差点个数与原复数据块对应的相位图中的残差点个数进行比较：

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数减少的时候，对所述更新后的复数据块重复执行步骤b和步骤c；

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数没有减少的时候，所述初步更新后的复数据块为所述最终更新后的复数据块。

进一步的，处理器901还用于调用存储器902中的实现前述实施例的方法或功能模块的程序代码和应用程序，以执行以下操作：

提取所述复数据块对应的相位图；

根据残差点的定义从所述相位图中查找到残差点；

对所述复数据块中处于残差点相应位置的像素及以所述像素为中心的确定尺寸范围的像素进行标记，得到所述复数据块的标记像素。

进一步的，处理器901还用于调用存储器902中的实现前述实施例的方法或功能模块的程序代码和应用程序，以执行以下操作：

根据所述复干涉图的特征以及确定的规则选择合适的马尔可夫随机场模型，并得到所述标记像素对应的能量函数

其中，σ表示估计值的标准差；z_p表示标识为p的标记像素的复数值；N_p表示以所述标记像素p为中心的固定阶数的邻近像素集合；z_t表示标记像素p的邻近像素t的复数值；表示标记像素p对邻近像素t的权值；Σ表示求和符号，|*|表示取模运算；

通过计算E(z_p)对z_p求偏导数为零的驻点得到所述标记像素对应的CMRF模型参数 ${{\widehat{\mathrm{\Λ}}}_P}^{*}={\left[\underset{s\∈I}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{sp}}{z}_s{Q}_s^H\right]\left[\underset{s\∈I}{\mathrm{\Σ}}{Q}_s{Q}_s^H\right]}^{-1},$

其中，I表示所述复数据块内所有像素点集合，Q_s＝[z_s1,z_s2,...,z_s24]^T为N_s中所有像素点复数值组成的列向量，与前述的N_p类似的，N_s为以像素s的为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；表示对Q_s进行共轭转置运算；ω_sp表示像素s相对于标记像素p的权重，具体的，ω_sp＝ω_sρ×ω_sc×ω_sd；其中，ω_sρ由所述复干涉图相关系数图来表示；ω_sc表示像素s和标记像素p的复数值所服从的统计特性之间的差异；ω_sd表示像素s和标记像素p之间的空间距离的关系。

进一步的，处理器901还用于调用存储器902中的实现前述实施例的方法或功能模块的程序代码和应用程序，以执行以下操作：

通过更新所述标记像素的复数值，其中，k+1表示当前的迭代计算，k表示上一次的迭代计算；μ为控制收敛速度的常数，通常选取为大于0小于1的自然数。

本实施例提供了一种SAR干涉相位图滤波设备70，通过分块单元701对SAR复干涉图进行重叠分块以及更新单元702在模型参数估计时进行加权参数加权估计，从而能够在保证了良好的滤除噪声效果的同时，还能够有效地消除SAR复干涉图的残差点，从而使得滤波后的复干涉图很好地保持了原SAR复干涉图的边缘特性。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (12)

1.一种SAR干涉相位图滤波方法，其特征在于，所述方法包括：

对SAR复干涉图进行重叠分块，得到至少两个复数据块；

对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；

将所述最终更新后的复数据块拼接成更新后的复干涉图，并从所述更新后的复干涉图中提取更新后的干涉相位图。

2.根据权利要求1所述的方法，其特征在于，所述对SAR复干涉图进行重叠分块处理，具体包括：

按照预先确定的复数据块尺寸和重复率对复干涉图进行重叠分块处理，并截取出相应的复数据块。

3.根据权利要求1所述的方法，其特征在于，所述对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块，具体包括：

步骤a：通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素；

步骤b：根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数；

步骤c：根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，得到初步更新后的复数据块；

步骤d：将所述初步更新后的复数据块对应的相位图中的残差点个数与原复数据块对应的相位图中的残差点个数进行比较：

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数减少的时候，对所述更新后的复数据块重复执行步骤b和步骤c；

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数没有减少的时候，所述初步更新后的复数据块为所述最终更新后的复数据块。

4.根据权利要求3所述的方法，其特征在于，通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素，具体包括：

提取所述复数据块对应的相位图；

根据残差点的定义从所述相位图中查找到残差点；

对所述复数据块中处于残差点相应位置的像素及以所述像素为中心的确定尺寸范围的像素进行标记，得到所述复数据块的标记像素。

5.根据权利要求3或4所述的方法，其特征在于，根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数，具体包括：

根据所述复干涉图的特征以及确定的规则选择合适的马尔可夫随机场模型，并得到所述标记像素对应的能量函数

其中，σ表示估计值的标准差；z_p表示标识为p的标记像素的复数值；N_p表示以所述标记像素p为中心的固定阶数的邻近像素集合；z_t表示标记像素p的邻近像素t的复数值；表示标记像素p对邻近像素t的权值；Σ表示求和符号，|*|表示取模运算；

通过计算E(z_p)对z_p求偏导数为零的驻点得到所述标记像素对应的CMRF模型参数 ${{\widehat{\mathrm{\Λ}}}_P}^{*}={\left[\underset{s\∈I}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{sp}}{z}_s{Q}_s^H\right]\left[\underset{s\∈I}{\mathrm{\Σ}}{Q}_s{Q}_s^H\right]}^{-1},$

其中，I表示所述复数据块内所有像素点集合，Q_s＝[z_s1,z_s2,...,z_s24]^T为N_s中所有像素点复数值组成的列向量，与前述的N_p类似的，N_s为以像素s的为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；表示对Q_s进行共轭转置运算；ω_sp表示像素s相对于标记像素p的权重，具体的，ω_sp＝ω_sρ×ω_sc×ω_sd；其中，ω_sρ由所述复干涉图相关系数图来表示；ω_sc表示像素s和标记像素p的复数值所服从的统计特性之间的差异；ω_sd表示像素s和标记像素p之间的空间距离的关系。

6.根据权利要求5所述的方法，其特征在于，所述根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，具体包括：

通过更新所述标记像素的复数值，其中，k+1表示当前的迭代计算，k表示上一次的迭代计算；μ为控制收敛速度的常数，通常选取为大于0小于1的自然数。

7.一种SAR干涉相位图滤波设备，其特征在于，所述设备包括：分块单元、更新单元、拼接单元和提取单元，其中，

所述分块单元用于，对SAR复干涉图进行重叠分块，得到至少两个复数据块；

所述更新单元用于，对每一个复数据块中的标记像素进行更新，对应地得到最终更新后的复数据块；

所述拼接单元用于，将所述最终更新后的复数据块拼接成更新后的复干涉图；

所述提取单元用于，从所述更新后的复干涉图中提取更新后的干涉相位图。

8.根据权利要求7所述的设备，其特征在于，所述分块单元具体用于，按照预先确定的复数据块尺寸和重复率对复干涉图进行重叠分块处理，并截取出相应的复数据块。

9.根据权利要求7所述的设备，其特征在于，所述更新单元包括：第一获取子单元、第二获取子单元、更新子单元和比较子单元，其中，

所述第一获取子单元用于，通过查找所述复数据块对应的相位图中的残差点获取所述复数据块的标记像素；

所述第二获取子单元用于，根据最小二乘法获取所述标记像素对应的复数马尔可夫随机场CMRF模型参数；

所述更新子单元用于，根据最速下降算法，利用所述标记像素对应的CMRF模型参数和所述标记像素设置范围内的邻近像素点，更新所述标记像素的复数值，得到初步更新后的复数据块；

所述比较子单元用于，将所述初步更新后的复数据块对应的相位图中的残差点个数与原复数据块对应的相位图中的残差点个数进行比较：

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数减少的时候，对所述更新后的复数据块触发第二获取子单元和更新子单元；

当所述初步更新后的复数据块对应的相位图中的残差点个数比所述原复数据块对应的相位图中的残差点个数没有减少的时候，所述初步更新后的复数据块为所述最终更新后的复数据块。

10.根据权利要求9所述的设备，其特征在于，所述第一获取子单元包括：提取模块、查找模块和标记模块，其中，

所述提取模块用于，提取所述复数据块对应的相位图；

所述查找模块用于，根据残差点的定义从所述相位图中查找到残差点；

所述标记模块用于，对所述复数据块中处于残差点相应位置的像素及以所述像素为中心的确定尺寸范围的像素进行标记，得到所述复数据块的标记像素。

11.根据权利要求9或10所述的设备，其特征在于，所述第二获取单元具体用于，

根据所述复干涉图的特征以及确定的规则选择合适的马尔可夫随机场模型，并得到所述标记像素对应的能量函数

其中，σ表示估计值的标准差；z_p表示标识为p的标记像素的复数值；N_p表示以所述标记像素p为中心的固定阶数的邻近像素集合；z_t表示标记像素p的邻近像素t的复数值；表示标记像素p对邻近像素t的权值；Σ表示求和符号，|*|表示取模运算；

通过计算E(z_p)对z_p求偏导数为零的驻点得到所述标记像素对应的CMRF模型参数 ${{\widehat{\mathrm{\Λ}}}_P}^{*}={\left[\underset{s\∈I}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{sp}}{z}_s{Q}_s^H\right]\left[\underset{s\∈I}{\mathrm{\Σ}}{Q}_s{Q}_s^H\right]}^{-1},$

其中，I表示所述复数据块内所有像素点集合，Q_s＝[z_s1,z_s2,...,z_s24]^T为N_s中所有像素点复数值组成的列向量，与前述的N_p类似的，N_s为以像素s的为中心的固定阶数的邻近像素集合，在本实施例中，所述固定阶数优选为5；表示对Q_s进行共轭转置运算；ω_sp表示像素s相对于标记像素p的权重，具体的，ω_sp＝ω_sρ×ω_sc×ω_sd；其中，ω_sρ由所述复干涉图相关系数图来表示；ω_sc表示像素s和标记像素p的复数值所服从的统计特性之间的差异；ω_sd表示像素s和标记像素p之间的空间距离的关系。

12.根据权利要求11所述的设备，其特征在于，所述更新子单元具体用于，通过 $z_p^{(k+1)}=z_p^{\left(k\right)}-\mathrm{\μ}\·{\mathrm{\Δz}}_p^{\left(k\right)}$ 更新所述标记像素的复数值，其中， ${\mathrm{\Δz}}_p^{\left(k\right)}=z_p^{\left(k\right)}-{\widehat{\mathrm{\Λ}}}_p^{*\left(k\right)}{Q}_p^{\left(k\right)},$ k+1表示当前的迭代计算，k表示上一次的迭代计算；μ为控制收敛速度的常数，通常选取为大于0小于1的自然数。