CN102930532A - 基于mrf迭代的sar图像非监督变化检测方法和装置 - Google Patents

Abstract

针对现有SAR检测技术检测能力弱、检测精度低和速度慢等不足，本发明提出了一种基于MRF迭代的SAR图像非监督变化检测方法，采用改进的EM算法稳健的估计差异图混合分布模型统计参数，采用MRF模型定义像素间的依赖关系，并利用基于图割的MRF分割产生更可靠更精确的图像变化检测结果；本发明还提供基于上述方法的装置，包括光纤链路输入模块、DSP处理模块、SDRAM存储器、CPCI输出模块和显示终端；本发明有益的效果有：解决了多类、高精度、快速稳健的非监督变化检测的难题；本发明具备海量数据批量处理的工程应用能力，满足复杂条件下的图像分析需求；本发明在处理精度和处理速度方面较传统方法有显著提高。

Description

基于MRF迭代的SAR图像非监督变化检测方法和装置

技术领域

本发明属于遥感图像数据处理技术领域，具体涉及一种雷达图像的非监督变化检测方法。

背景技术

合成孔径雷达SAR(Synthetic Aperture Radar)图像变化检测技术是基于机载或星载SAR获取的不同时相图像，检测地表变化的一种技术。近年来，SAR图像变化检测技术成为国内外研究热点。光学数据受到气候、覆盖等因素的影响，并不能满足所有变化检测的需求。SAR作为一种主动微波传感器，具有全天候、全天时、强穿透的工作能力，利用SAR图像进行变化检测具有重要意义。目前，SAR图像变化检测已经在多个方面取得广泛应用，如土地利用分析、森林采伐监测、灾情估计、军事侦察、打击效果评估等。

在变化检测的处理技术方面，非监督变化检测问题是国内外关注的重点。非监督变化检测不需要先验变化信息的支持，直接从配准好的前后两个时刻SAR图像获取变化信息。这样处理的好处，一方面降低人为误差的影响，另一方面符合实际应用中先验变化信息缺失的现实情况。SAR图像变化检测的流程概括为三个步骤，分别是图像预处理、差异图获取和差异图分割。非监督变化检测方法的核心工作就是针对这三个步骤展开的。

在预处理阶段，图像配准和斑噪滤波是必不可少的。图像配准是变化检测的前提，它保证前后两幅图像的像元尺寸、地理位置一致。配准精度要求在1个像素以内。斑噪滤波能够有效抑制相干斑，提高图像质量。目前很多变化检测算法直接针对配准的图像进行处理，忽略了斑噪滤波，从而导致检测效果不好，算法不稳定。

在差异图获取阶段，针对不同数据源有不同的处理方法。比如，针对强度或幅度数据，常采用对数比值法获取差异图；针对单视复数据，常采用相干系数和相关系数获取差异图；针对多波段或多极化数据，常采用特征矢量的方式获取差异图。不同的差异图获取方法影响着差异图分割技术的选择。通常，差异图包含变化和未变化两类信息，服从混合分布模型。因此，从差异图上进行自动阈值分割是获取变化区域的常用方法。研究者将关注点集中在差异图统计分布模型选择和参数估计上，进而在此基础上选择阈值分割算法，比如基于最小错误准则的阈值分割、K&I阈值分割、CFAR阈值分割和聚类分割等。这些分割算法都存在一个明显的不足，那就是没有考虑像素的邻域特征，从而导致最终的分割效果仍然有限。因为，如果一个像素标记为变化或没有变化的区域，那么它周围的像素极有可能是同样的标记。因此，需要采用一种可以利用邻域信息的技术提高变化检测结果的可靠性与精确性。

此外，现有技术的非监督变化检测方法在同时检测图像的变强和变弱区域、提高检测精度和速度、实时处理等方面不能满足实际应用的要求。

发明内容

本发明的目的是针对现有SAR图像变化检测技术无法满足同时检测变强和变弱区域、较高的检测精度和速度以及实时处理等方面的需求的问题，提出了一种基于马尔可夫随机场MRF(Markov Random Field)迭代的SAR图像非监督变化检测方法，方法中采用改进的EM算法稳健的估计差异图混合分布模型统计参数，采用MRF模型定义像素间的依赖关系，并利用基于图割的MRF分割产生更可靠更精确的图像变化检测结果；本发明方法充分考虑斑噪、图像统计分布等因素，解决多类、高精度、快速稳健的非监督变化检测的难题。

一种基于MRF迭代的SAR图像非监督变化检测方法，按如下的步骤进行：

(1)获取待检测图像并滤波：首先输入已配准的t₁时刻的SAR图像I₁和t₂时刻的SAR图像I₂，并分别对I₁和I₂进行斑噪滤波，其中，滤波算法采用均值法或Lee滤波算法；滤波矩形窗口大小设置为3×3、5×5或7×7中的任一个；所述滤波窗口的单位为像素；

(2)计算对数比值差异图：将滤波后的图像I₁上的像素同位于滤波后图像I₂上的相应位置的像素进行对数比较，并生成对数比值差异图，其中计算对数比值差异图时所采用的计算公式为：

(3)采用改进的EM算法估计差异图混合分布模型统计参数；在高斯混合模型的条件下，采用改进的EM算法迭代估计出差异图的变强区域w₁、变弱区域w₂和未发生变化区域w₃三类的统计分布参数和权重参数。具体为：

EM迭代的各类初始参数估计通过如下阈值获得，保证检测结果的正确收敛。假设对数比值差异图的像素均值为m，标准差为s，在差异图直方图上设定两个阈值T₁和T₂分别作为变弱和变强区域的初始分割：

T₁＝m-a·s；                    (I)

T₂＝m+a·s；                    (II)

其中，a是控制系数，取值在1~2之间，默认为1。在高斯混合模型假设下，对数比值图的概率密度函数表示为：

$f\left(x\right|\mathrm{\θ})=\underset{i=1}{\overset{g}{\mathrm{\Σ}}}{\mathrm{\α}}_i{f}_i\left(x\right|{\mathrm{\θ}}_i)---\left(\mathrm{III}\right)$

其中，g表示类别数，α₁,...α_g是每类密度函数的权重，且有

f_i(x|θ_i)是每类的密度函数， $f_i\left(x\right|{\mathrm{\θ}}_i)={\left(2\mathrm{\π}\right)}^{-d./2}{\left|{\mathrm{\Σ}}_i\right|}^{-1/2}\exp (-\frac{1}{2}{(x-{\mathrm{\μ}}_i)}^T{{\mathrm{\Σ}}_i}^{-1}(x-{\mathrm{\μ}}_i)),$ θ_i＝(μ_i,∑_i)是密度模型参数。在阈值T₁和T₂下，各类的均值和方差计算公式如下：

${\mathrm{\μ}}_i=\frac{1}{n_i}\underset{k=1,x_k\∈w_i}{\overset{n_i}{\mathrm{\Σ}}}{x}_k---\left(\mathrm{IV}\right)$

${\mathrm{\Σ}}_i=\frac{1}{n_i}\underset{k=1,x_k\∈w_i}{\overset{n_i}{\mathrm{\Σ}}}\left[{(x_k-{\mathrm{\μ}}_i)}^2\right]---\left(V\right)$

其中，n_i,i＝1,2,3表示各类的像素个数。

改进的EM算法迭代步骤如下：

E步骤(t+1次迭代)：

$p_{\mathrm{ij}}^{(t+1)}=p_i(x_j;{\mathrm{\θ}}_i^{\left(t\right)})=\frac{{\mathrm{\α}}_i^{\left(t\right)}{f}_i(x_j;{\mathrm{\θ}}_i^{\left(t\right)})}{\underset{k=1}{\overset{g}{\mathrm{\Σ}}}{\mathrm{\α}}_k^{\left(t\right)}{f}_k(x_j;{\mathrm{\θ}}_k^{\left(t\right)})}---\left(\mathrm{VI}\right)$

M步骤(t+1次迭代)：

${\mathrm{\α}}_i^{(t+1)}=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}^{(t+1)},---\left(\mathrm{VII}\right)$

${\mathrm{\μ}}_i^{(t+1)}=\frac{k_i{{\mathrm{\μ}}_i}^{\left(0\right)}+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}^{(t+1)}{x}_j}{k_i+n{\mathrm{\α}}_i^{(t+1)}},---\left(\mathrm{VIII}\right)$

${\mathrm{\Σ}}_i^{(t+1)}=\frac{k_i{{\mathrm{\Σ}}_i}^{\left(0\right)}+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}^{(t+1)}(x_j-{\mathrm{\μ}}_i^{(t+1)})(x_j-{\mathrm{\μ}}_i^{(t+1)})}{k_i+n{\mathrm{\α}}_i^{(t+1)}}.---\left(\mathrm{IX}\right)$

其中，n为所有像素的个数；μ_i ⁽⁰⁾和∑_i ⁽⁰⁾是各类均值和方差的初始估计；k_i是非负整数，越大表示初始估计在最终结果中的比重越大，其作用是防止EM迭代收敛到偏离真实值的位置；

(4)最大后验概率初始分割：对差异图上每个像素，在最大后验概率准则下对其进行类别划分，分别归入变强区域w₁、变弱区域w₂和未发生变化区域w₃，判别准则如下：

$c=\underset{i\∈\left(\mathrm{1,2,3}\right)}{\arg }\max {\mathrm{\α}}_i{f}_i\left(x\right|{\mathrm{\θ}}_i)---\left(X\right)$

(5)基于图割的MRF分割：在前一步分割的基础上，采用基于图割的MRF分割算法提高变化检测的精度。网络流加权图G=[V，E]由节点集合V和边集合E组成(如图3)。s,t∈V是特殊的两个顶点，分别称为源和汇。除s,t外，图中每个节点代表图像中一个像素。图中相邻结点之间的边权值代表它们之间的相似程度。在此加权图中寻找最小割的过程就是实现MRF分割的过程。寻找最小割的算法采用图论中成熟的最大流/最小割算法。基于图割的MRF分割算法比基于ICM、SA的MRF分割算法具有更高的精度和速度，对比结果如表1所示。

(6)在完成基于图割的MRF分割后进行图像的收敛性判断：如果收敛则得出变化检测结果；反之根据当前分割结果，分别采用公式(IV)和(V)重新估计均值和方差，循环步骤(5)和(6)的操作，直到统计分布参数收敛，收敛判断准则如下式：

$\underset{i=1}{\overset{g}{\mathrm{\Σ}}}\left(\right|{\mathrm{\μ}}_i^{t+1}-{\mathrm{\μ}}_i^t|+|{\mathrm{\Σ}}_i^{t+1}-{\mathrm{\Σ}}_i^t\left|\right)\≤\mathrm{\ϵ}---\left(\mathrm{XI}\right)$

其中，ε是一个极小的正数，默认取值1e^-4。

(7)将符合收敛的结果输出。

采用上述一种基于MRF迭代的SAR图像非监督变化检测方法的SAR图像处理装置，包括光纤链路输入模块100、DSP处理模块200、SDRAM存储器300、CPCI输出模块400、显示终端500，其组成关系为：

所述的光纤链路输入模块100由光纤协议转换模块101、输入端数据缓冲时序控制模块102和DSP链路口协议模块103组成，所述的光纤协议转换模块101、输入端数据缓冲时序控制模块102和DSP链路口协议模块103分别由FPGA的内部单元实现；其中，所述的光纤链路输入模块100负责接收通过光纤传输过来的t₁、t₂时刻的图像数据，解译协议并将高速串行数据转化成并行数据，利用FPGA内部的DPRAM实现数据缓冲，通过输入端数据缓冲时序控制模块102的时序控制和链路口协议将t₁、t₂时刻的图像数据经DSP链路口协议模块103发送给DSP处理模块200，实现图像数据的输入；

所述的DSP处理模块200，采用ADSP-TS201处理芯片，并通过链路口实现图像数据的接收，此外通过其自带的SDRAM控制器将从光纤链路输入模块100获取的t₁、t₂时刻的图像数据存储在SDRAM存储器300中；DSP处理模块200中运行设定的算法程序，按照权利要求1所述的方法对存储在SDRAM存储器300中的图像数据进行处理，并将处理结果通过CPCI输出模块400显示在显示终端500上；

CPCI输出模块400由DSP总线接口401、输出端数据缓冲时序控制模块402以及CPCI总线协议模块组成，所述的CPCI输出模块400由FPGA的内部单元实现具体的功能，并采用IP核实现CPCI总线协议；所述的CPCI输出模块400接收自DSP处理模块200发出的数据并通过CPCI总线传递给显示终端500，实现处理结果的实时显示。

本发明的有益效果有：

1）本发明充分考虑乘性噪声、像素特性之间的邻域相关关系、图像幅度的统计分布等因素，突破统计模型参数估计、阈值自动获取、多类变化区域快速收敛等瓶颈，解决多类、高精度、快速稳健的非监督变化检测的难题；

2）本技术发明具备海量数据批量处理的工程应用能力，满足复杂条件下的图像分析需求；

3）本发明在处理精度和处理速度方面，比传统采用ICA和SA的MRF分割有显著提高。

附图说明

图1是本发明处理方法的流程图。

图2是本发明的硬件框图。

图3是网络流和图割示意图。

图4是检测用的t₁时刻的ERS-1SAR幅度图I₁。

图5是检测用的t₂时刻的ERS-2SAR幅度图I₂。

图6是对图5进行仿真处理后的ERS-2SAR图像I₂。

图7是图6与图4的对数比值差异图的直方图。

图8是对数比值差异图的EM参数估计残差收敛曲线。

图9是对数比值差异图的EM估计混合密度函数曲线。

图10是对数比值差异图最大后验概率初始检测结果，其中黑色区域变弱，白色区域变强，灰色区域未变化。

图11是对数比值差异图MRF迭代分割最终检测结果，其中黑色区域变弱，白色区域变强，灰色区域未变化。

具体实施方式

以下结合附图对本发明作进一步实施说明。

图4、5分别是两幅间隔1天的ERS-1SAR图像I₁和ERS-2SAR图像I₂经过精配准之后的实验区，两幅图所显示的实验区的像素均为400×400。除了噪声之外，图像区域的地物类别没有发生变化。

如图6所示，为了证明本技术的有效性，基于I₂进行了变化区域仿真，形成图像I_2c，其中在(101,101)的位置填补了一块40×40的变弱区域；在(251,251)的位置填补了一块50×50的变强区域。本实验过程采用图4和图6进行数据分析比较。

参见图1，按本发明的方法进行如下操作：采用均值滤波法，分别对I₁和I_2c进行滤波，滤波像素窗口为3×3；经过滤波之后，I₁和I_2c图像上的噪声得到很好的抑制；通过计算对数比值，我们得到如图7所示差异图直方图。采用高斯混合模型对直方图进行建模，并假设类别包括变强区域、变弱区域和未发生变化的区域。在进行模型参数估计时，阈值控制系数a取默认值1。图8给出了差异图的EM参数估计的残差收敛曲线。图9给出了三类概率密度函数曲线。

在EM参数估计的基础上，利用最大后验概率准则，对差异图进行三类分类，初始分类结果如图10所示。在经过基于图割的MRF迭代分割处理之后，最终收敛的变化检测结果如图11所示。表1给出了本文方法和在图像分割阶段采用基于ICM、SA的MRF迭代分割方法的检测混淆矩阵，并从检测率和计算耗时两个指标进行了评价估计。从表中可以看出，本发明所提出的方法在精度和效率上都优于传统基于ICM或SA的MRF迭代方法。

表1变化检测混淆矩阵分析(像素个数)

注：计算环境CPU Intel奔腾双核2.0GHz,内存2.0GB。

参见图2，采用上述基于MRF迭代的SAR图像非监督变化检测方法的SAR图像处理装置，包括光纤链路输入模块100、DSP处理模块200、SDRAM存储器300、CPCI输出模块400、显示终端500，其组成关系为：

所述的光纤链路输入模块100由光纤协议转换模块101、输入端数据缓冲时序控制模块102和DSP链路口协议模块103组成，所述的光纤协议转换模块101、输入端数据缓冲时序控制模块102和DSP链路口协议模块103分别由FPGA的内部单元实现；其中，所述的光纤链路输入模块100负责接收通过光纤传输过来的t₁、t₂时刻的图像数据，解译协议并将高速串行数据转化成并行数据，利用FPGA内部的DPRAM实现数据缓冲，通过输入端数据缓冲时序控制模块102的时序控制和链路口协议将t₁、t₂时刻的图像数据经DSP链路口协议模块103发送给DSP处理模块200，实现图像数据的输入；

所述的DSP处理模块200，采用ADSP-TS201处理芯片，并通过链路口实现图像数据的接收，此外通过其自带的SDRAM控制器将从光纤链路输入模块100获取的t₁、t₂时刻的图像数据存储在SDRAM存储器300中；DSP处理模块200中运行设定的算法程序，按照权利要求1所述的方法对存储在SDRAM存储器300中的图像数据进行处理，并将处理结果通过CPCI输出模块400显示在显示终端500上；

CPCI输出模块400由DSP总线接口401、输出端数据缓冲时序控制模块402以及CPCI总线协议模块组成，所述的CPCI输出模块400由FPGA的内部单元实现具体的功能，并采用IP核实现CPCI总线协议；所述的CPCI输出模块400接收自DSP处理模块200发出的数据并通过CPCI总线传递给显示终端500，实现处理结果的实时显示。

Claims (2)

1.基于MRF迭代的SAR图像非监督变化检测方法，其特征在于按如下步骤进行：

(1)获取待检测图像并滤波：首先输入已配准的t₁时刻的SAR图像I₁和t₂时刻的SAR图像I₂，并分别对I₁和I₂进行斑噪滤波处理，其中，滤波算法为均值法或Lee滤波算法；

(2)计算对数比值差异图：对滤波后的图像I₁和滤波后的图像I₂按像素进行对数比较，并生成对数比值差异图；

(3)采用改进的EM算法估计差异图混合分布模型统计参数；在高斯混合模型的条件下，采用改进的EM算法迭代估计出差异图的变强区域w₁、变弱区域w₂和未发生变化区域w₃三类的统计分布参数和权重参数；

(4)最大后验概率初始分割：在最大后验概率准则下，对差异图上每个像素进行类别划分，归类至变强区域w₁、变弱区域w₂和未发生变化区域w₃；

(5)基于图割的MRF分割：在前一步分割的基础上，采用基于图割的MRF分割算法提高变化检测的精度；

(6)在完成基于图割的MRF分割后进行图像的收敛性判断：如果判断为收敛，则得出变化检测结果；反之，根据当前分割结果重新估计各区域的统计分布参数，重复步骤(5)和(6)的操作，直到统计分布参数收敛；

（7）将变化检测结果输出。

2.采用权利要求1所述基于MRF迭代的SAR图像非监督变化检测方法的装置，包括光纤链路输入模块（100）、DSP处理模块（200）、SDRAM存储器（300）、CPCI输出模块（400）、显示终端（500），其特征在于：

所述的光纤链路输入模块（100）由光纤协议转换模块（101）、输入端数据缓冲时序控制模块（102）和DSP链路口协议模块（103）组成，所述的光纤协议转换模块（101）、输入端数据缓冲时序控制模块（102）和DSP链路口协议模块（103）分别由FPGA的内部单元实现；其中，所述的光纤链路输入模块（100）负责接收通过光纤传输过来的t₁、t₂时刻的图像数据，解译协议并将高速串行数据转化成并行数据，利用FPGA内部的DPRAM实现数据缓冲，通过输入端数据缓冲时序控制模块（102）的时序控制和链路口协议将t₁、t₂时刻的图像数据经DSP链路口协议模块（103）发送给DSP处理模块（200），实现图像数据的输入；

所述的DSP处理模块（200），采用ADSP-TS201处理芯片，并通过链路口实现图像数据的接收，此外通过其自带的SDRAM控制器将从光纤链路输入模块（100）获取的t₁、t₂时刻的图像数据存储在SDRAM存储器（300）中；DSP处理模块（200）中运行设定的算法程序，按照权利要求1所述的方法对存储在SDRAM存储器（300）中的图像数据进行处理，并将处理结果通过CPCI输出模块（400）显示在显示终端（500）上；

CPCI输出模块（400）由DSP总线接口（401）、输出端数据缓冲时序控制模块（402）以及CPCI总线协议模块（403）组成，所述的CPCI输出模块（400）由FPGA的内部单元实现具体的功能，并采用IP核实现CPCI总线协议；所述的CPCI输出模块（400）接收自DSP处理模块（200）发出的数据并通过CPCI总线传递给显示终端（500），实现处理结果的实时显示。

Images