CN103473559A - 基于nsct域合成核的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于NSCT域合成核的SAR图像变化检测方法，主要解决现有Gabor域合成核导致的无法用多分辨率分析进行SAR图像变化检测的问题。其实现步骤是：(1)将原始两时相SAR图像分别进行多个尺度、多个方向的NSCT分解；(2)提取分解后两时相图像在每一个尺度上的归一化的强度特征和纹理特征，并构造尺度内的强度纹理差值合成核；(3)将任意一个尺度内的差值合成核输入到支撑矢量机中进行检测，得到该尺度的变化检测结果；(4)对所有尺度的变化检测结果进行决策级尺度间融合得到最终的变化检测结果。本发明与现有的Gabor域合成核方法相比，检测精度高，运算效率高，可用于SAR图像变化检测。

Description

基于NSCT域合成核的SAR图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及图像变化检测，可用于对SAR图像进行地物状态变化监测与评估。

背景技术

SAR图像的变化检测指利用不同时期的SAR图像获取地物变化信息的过程，是一种针对SAR图像的特点而建立的数据分析方法，用于识别地表状态的变化。光学遥感虽然成像分辨率高，但只能在白天成像，且受天气影响较重，而SAR系统具有全天候、全天时、覆盖面积大等优点，同时具有固定的重访周期，分辨率高，能斜距成像等特点，所以SAR图像比光学遥感图像更适合于变化检测。随着SAR技术的不断发展，SAR图像的变化检测在资源和环境监测、军事目标探测、农业作物监测、以及自然灾情监测与评估等领域的应用越来越广泛。

依据SAR图像变化检测方法的特点，SAR图像的变化检测可分为：基于简单代数运算的变化检测方法，如图像差值法，图像比值法；基于图像变换的变化检测方法，如主成分分析法；基于多信息融合的变化检测方法，如融合SAR影像强度信息和纹理信息进行变换检测；基于统计模型的变化检测方法，如学者GabrieleMoser提出了利用马尔科夫随机场MRF模型构建SAR图像像素的邻域相关信息的变化检测方法。近年来研究较多的是基于多尺度分析和核方法的SAR图像变化检测方法，例如学者L.Bruzzone将基于细节保留小波变换的多尺度分析方法引入到SAR变化检测中，该方法具有良好的抗噪性能，但是对于细节和噪声的保留没有一个较优的处理方法；Kai-Kuang Ma提出一种基于双树-复小波变换DT-CWT的变化检测方法，它利用DT-CWT对对数比值图进行多尺度分解，但它没有考虑到图像的纹理信息，且其中阈值的选取也是一个棘手的问题；最近学者Gustavo Camps-Valls又提出了基于Gabor域合成核的变化检测方法，它首先利用Gabor滤波器提取纹理信息，然后利用核方法融合原始图像的强度信息和纹理信息，最后用支持矢量积SVM实现变化检测，由于该Gabor域合成核方法能较好的融合原始图像的强度特征和纹理特征，因此其能产生较优的变化检测结果，但Gabor域合成核方法存在如下缺点：该方法由于不是专门针对SAR图像而建立的方法，因此在处理受相干斑噪声影响的SAR图像时其精度不高；由于Gabor变换是一种具有多尺度、单一分辨率的变换，不能进行多分辨分析。

发明内容

本发明的目的在于克服上述已有问题的缺点，提出一种基于非下采样轮廓波变换NSCT域合成核的SAR图像变化检测方法，以解决Gabor域合成核导致的无法用多分辨率分析对SAR图像进行变化检测的问题。

为实现上述目的，本发明的实现步骤包括如下：

A训练步骤：

A1)使用非下采样轮廓波变换NSCT将原始两时相SAR图像X₁={X₁(p，q)|1≤p≤I，1≤q≤J}和X₂={X₂(p，q)|1≤p≤I，1≤q≤J}进行S个尺度的分解，得到S个尺度上的两时相图像，其中I，J分别为原始图像长度和宽度大小，(p，q)为图像像素点；

A2)对步骤A1)得到的S个尺度上的两时相图像提取强度特征和纹理特征，并进行归一化，得到归一化后的S个尺度上的强度特征

和纹理特征再将归一化后的这两种特征联合，构成S个尺度上的两时相图像特征

${T_i}^m=\{g_{X_i}^m,w_{X_i}^m\},$

其中，m=1，2...S，i＝1，2，i是指示原始两时相SAR图像Ｘ₁或者X₂的下标；

A3)从步骤A1)得到的S个尺度的两时相图像中，手动提取N_s个标签训练样本，将这些训练样本用步骤A2)得到的归一化后的S个尺度上的强度特征和纹理特征进行表示，得到S个尺度上的两时相的标签训练样本特征

${\mathrm{SP}}_i^m={\left\{(x_{\mathrm{ij}}^m,y_{\mathrm{ij}}^m)\right\}}_{j=1}^{N_s},m=\mathrm{1,2}...S,1\≤j\≤N_s,i=\mathrm{1,2},$

其中， $x_{\mathrm{ij}}^m=\{g_{X_i}^m(p_j,q_j),w_{X_i}^m(p_j,q_j)\},$

(p_j，q_j)为标签训练样本点的坐标，表示第m个尺度上的标签训练样本点(p_j，q_j)处的强度特征，

表示第m个尺度上的标签训练样本点(p_j，q_j)处的纹理特征，

表示第m个尺度上的标签训练样本点(p_j，q_j)处的标签；

A4)将步骤A3)得到的两时相标签训练样本的强度特征和纹理特征映射到核空间，得到核空间中的两时相标签训练样本的强度特征和纹理特征，然后利用这两种特征构造S个尺度上的强度核

和纹理核

其中

为第m个尺度上的标签训练样本点(p_k，q_k)处的特征，

表示第m个尺度上的标签训练样本点(p_k，q_k)处的强度特征，

表示第m个尺度上的标签训练样本点(p_k，q_k)处的纹理特征；

A5)将步骤A4)得到的核空间中的两时相标签训练样本的强度特征和纹理特征进行加权求和，得到核空间中的两时相的强度纹理合成特征，利用该合成特征构造S个尺度上的强度纹理合成核

A6)对步骤A5)得到的两时相的标签训练样本的强度纹理合成特征做差，得到标签训练样本的差特征，利用该差特征构造S个尺度上的差值合成核

其中

表示第m个尺度上的标签训练样本(p_j，q_j)处的差特征，

表示第m个尺度上的标签训练样本(p_k，q_k)处的差特征，1≤j，k≤N_s，m=1，...，S；

A7)将步骤A6)得到的S个尺度上的差值合成核矩阵

输入到支撑矢量积SVM中对SVM进行训练，得到S个尺度上的支撑矢量分类器SVC_m，其中m=1，2...S；

B测试步骤：

B1)将步骤A2)得到的两时相图像的强度特征和纹理特征以及步骤A3)得到的两时相标签训练样本的强度特征和纹理特征映射到核空间，得到核空间中的两时相的强度特征和纹理特征，利用这两个特征构造S个尺度上的强度核

和纹理核

其中， $T_i^m(p,q)=\{g_{X_i}^m(p,q),w_{X_i}^m(p,q)\}$ 为像素点(p，q)处的特征，

表示第m个尺度上的像素点(p，q)处的强度特征，

表示第m个尺度上的像素点(p，q)处的纹理特征；

B2)将步骤B1)得到的核空间中的两时相的强度特征和纹理特征进行加权求和，得到核空间中的两时相的强度纹理合成特征，利用该合成特征构造S个尺度上的强度纹理合成核

B3)对步骤B2)得到的两时相的强度纹理合成特征做差，得到差特征，利用该差特征构造S个尺度上的差值合成核

其中

表示在第m个尺度上的标签训练样本(p_j，q_j)处的差特征，T^m(p，q)表示第m个尺度上的像素点(p，q)处的差特征，1≤j≤N_s，m＝1，...，S；

B4)将差值合成核矩阵

输入到步骤A7)得到的支撑矢量分类器SVC_m中，对SVC_m进行测试，得到S个尺度上的二值变化检测结果图BCM_m；

B5)对二值变化检测结果图BCM_m进行决策级尺度间融合，得到最终的变化检测结果图。

本发明具有如下优点：

(1)本发明由于利用NSCT对两时相SAR图像进行分解，使得对原始图像既能进行多分辨分析又能进行时频局部分析；同时由于本发明在特征提取过程中具有多方向性、各向异性和移不变性等特点，因此所提取的强度信息和纹理信息能更充分的反映原始图像在多分辨率、多方向下的信息，进而得到更有效的融合特征。

(2)本发明由于利用了核方法，可以将原始空间中的非线性不可分问题转化为高维空间中的线性可分问题，因此本发明更适合处理具有非线性特征的SAR图像。

仿真结果表明，本发明与现有的Gabor域合成核方法相比，具有更高的计算效率和检测精度。

附图说明

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

图2是本发明应用于两时相洪水灾害Real SAR图像的变化检测结果图；

图3是本发明应用于两时相机场跑道Real SAR图像的变化检测结果图；

图4是本发明应用于两时相农田Real SAR图像的变化检测结果图。

具体实施方式

参照图1，本发明的具体实施包括如下训练步骤和测试步骤：

一.训练步骤：

步骤1，SAR图像分解。

使用非下采样轮廓波变换NSCT将原始两时相SAR图像X₁＝{X₁(p，q)|1≤p≤I，1≤q≤J}和X₂＝{X₂(p，q)|1≤p≤I，1≤q≤J}进行S个尺度的分解，得到S个尺度上的两时相图像，其中I，J分别为原始图像长度和宽度大小，(p，q)为图像像素点，S＝3。

步骤2，提取图像特征。

对步骤1得到的S个尺度上的两时相图像，提取其强度特征和纹理特征，其步骤如下：

2.1)令

表示第m个尺度上的两时相图像的低通子带系数，并用该低通子带系数

作为两时相图像在第m个尺度上的强度特征，其中m＝1，...，S，S为尺度数，i是指示原始两时相SAR图像X₁或者X₂的下标；

2.2)令

表示第m个尺度上的两时相图像的高通子带系数，d表示该高通子带系数的方向，在以(p，q)为中心像素、大小为V的窗口上提取高通子带系数

的均值信息

和方差信息

${\mathrm{\μ}}_{X_i}^{m,d}(p,q)=\frac{1}{N^2}\underset{j=p-N/2}{\overset{p+N/2}{\mathrm{\Σ}}}\underset{l=q-N/2}{\overset{q+N/2}{\mathrm{\Σ}}}{H}_{X_i}^{m,d}(j,l),$

${\mathrm{\σ}}_{X_i}^{m,d}(p,q)=\frac{1}{N^2-1}\underset{j=p-N/2}{\overset{p+N/2}{\mathrm{\Σ}}}\underset{l=q-N/2}{\overset{q+N/2}{\mathrm{\Σ}}}{(H_{X_I}^{m,d}(j,l)-{\mathrm{\μ}}_{X_i}^{m,d}(p,q))}^2,$

其中，d＝1，...，n_m，n_m为第m个尺度上的方向数；

2.3)以第m个尺度上的两时相图像中所有像素点为中心像素点，提取上述均值信息和方差信息，得到均值矢量和方差矢量

2.4)用均值矢量

和方差矢量

联合构成两时相图像在第m个尺度上的纹理特征。

步骤3，特征归一化。

3.1)对步骤2得到的S个尺度上的两时相图像的强度特征和纹理特征进行归一化，得到归一化后的S个尺度上的强度特征

和纹理特征

归一化过程按如下公式进行：

$g_{X_i}^m=\frac{L_{X_i}^m-\min \left(L_{X_i}^m\right)}{\max \left(L_{X_i}^m\right)-\min \left(L_{X_i}^m\right)},$

${\widetilde{\mathrm{\μ}}}_{X_i}^{m,d}=\frac{{\mathrm{\μ}}_{X_i}^{m,d}-\min \left({\mathrm{\μ}}_{X_i}^{m,d}\right)}{\max \left({\mathrm{\μ}}_{X_i}^{m,d}\right)-\min \left({\mathrm{\μ}}_{X_i}^{m,d}\right)},$

${\widetilde{\mathrm{\σ}}}_{X_i}^{m,d}=\frac{{\mathrm{\σ}}_{X_i}^{m,d}-\min \left({\mathrm{\σ}}_{X_i}^{m,d}\right)}{\max \left({\mathrm{\σ}}_{X_i}^{m,d}\right)-\min \left({\mathrm{\σ}}_{X_i}^{m,d}\right)},$

其中，

是第m个尺度上的两时相图像的归一化后的强度特征，是归一化后的均值矢量，

是归一化后的方差矢量，和

联合构成归一化后的第m个尺度上的两时相图像的纹理特征

3.2)将归一化后的强度特征和纹理特征进行联合，构成归一化后的S个尺度上的两时相图像特征

$T_i^m=\{g_{X_i}^m,w_{X_i}^m\},$

其中，m＝1，2...S，i＝1，2。

步骤4，提取训练样本特征。

4.1)从步骤1得到的S个尺度的两时相图像中，手动提取N_s个标签训练样本；

4.2)将训练样本用步骤3得到的归一化后的S个尺度上的强度特征和纹理特征进行表示，得到S个尺度上的两时相的标签训练样本特征

${\mathrm{SP}}_i^m={\left\{(x_{\mathrm{ij}}^m,y_{\mathrm{ij}}^m)\right\}}_{j=1}^{N_s},m=\mathrm{1,2}...S,1\≤j{\≤N}_s,i=\mathrm{1,2},$

其中， $x_{\mathrm{ij}}^m=\{g_{X_i}^m(p_j,q_j),w_{X_i}^m(p_j,q_j)\},$

(p_j，q_j)为标签训练样本点的坐标，

表示第m个尺度上的标签训练样本点(p_j，q_j)处的强度特征，表示第m个尺度上的标签训练样本点(p_j，q_j)处的纹理特征，

表示第m个尺度上的标签训练样本点(p_j，q_j)处的标签。

步骤5，构造强度核和纹理核。

5.1)将步骤4得到的两时相标签训练样本的强度特征和纹理特征映射到核空间，得到核空间中的标签训练样本的强度特征和纹理特征；

5.2)利用步骤5.1)得到的核空间中标签训练样本的强度特征和纹理特征构造S个尺度上的强度核

和纹理核

其公式如下：

$K_g^m(x_{\mathrm{ij}}^m,x_{\mathrm{ik}}^m)=<{\mathrm{\&phi;}}_g\left(x_{\mathrm{ij}}^m\right),{\mathrm{\&phi;}}_g\left(x_{\mathrm{ik}}^m\right)>$

$=K_{\mathrm{gas}}{(g_{X_i}^m(p_j,q_j),g_{X_i}^m(p_k,q_k))}^{,}$

$K_w^m(x_{\mathrm{ij}}^m,x_{\mathrm{ik}}^m)=<{\mathrm{\&phi;}}_w\left(x_{\mathrm{ij}}^m\right),{\mathrm{\&phi;}}_w\left(x_{\mathrm{ik}}^m\right)>$

$=K_{\mathrm{gas}}{(w_{X_i}^m(p_j,q_j),w_{X_i}^m(p_k,q_k))}^{,}$

其中，K_gas为高斯径向基核函数，为核空间中第m个尺度上的标签训练样本(p_j，q_j)处的强度特征，

为核空间中第m个尺度上的标签训练样本(p_j，q_j)处的纹理特征，

为核空间中第m个尺度上的标签训练样本(p_k，q_k)处的强度特征，

为核空间中第m个尺度上的标签训练样本(p_k，q_k)处的纹理特征，为第m个尺度上的标签训练样本点(p_k，q_k)处的特征，为第m个尺度上的标签训练样本点(p_k，q_k)处的强度特征，

为原始空间中第m个尺度上的标签训练样本点(p_k，q_k)处的纹理特征，1≤k，j≤N_s，m=1，...，S。

步骤6，构造强度纹理合成核。

6.1)将步骤5得到的核空间中的两时相标签训练样本的强度特征和纹理特征进行加权求和，得到核空间中两时相标签训练样本的强度纹理合成特征；

6.2)利用步骤6.1)得到的强度纹理合成特征构造S个尺度上的强度纹理合成核

其公式如下：

其中，为核空间中第m个尺度上的标签训练样本(p_j，q_i)处的强度纹理合成特征，

为核空间中第m个尺度上的标签训练样本(p_k，q_k)处的强度纹理合成特征，α为大于0小于1的常数。

步骤7，构造差值合成核。

7.1)对步骤6得到的两时相标签训练样本的强度纹理合成特征做差，得到标签训练样本的差特征；

7.2)利用步骤7.1)得到的差特征构造S个尺度上的差值合成核

其公式如下：

$K_D^m(x_j^m,x_k^m)=\&lang;\mathrm{\&Phi;}\left(x_j^m\right),\mathrm{\&Phi;}\left(x_k^m\right)\&rang;$

$=K_{\mathrm{gw}}^m(x_{1j}^m,x_{1k}^m)+K_{\mathrm{gw}}^m(x_{2j}^m,x_{2k}^m)-K_{\mathrm{gw}}^m(x_{1j}^m,x_{2k}^m)-K_{\mathrm{gw}}^m(x_{2j}^m,x_{1k}^m)$

$=\mathrm{\&alpha;}{K}_g^m(x_{1j}^m,x_{1k}^m)+(1-\mathrm{\&alpha;})K_w^m(x_{1j}^m,x_{1k}^m)+\mathrm{\&alpha;}{K}_g^m(x_{2j}^m,x_{2k}^m)+(1-\mathrm{\&alpha;})K_w^m(x_{2j}^m,x_{2k}^m)$

$=\mathrm{\&alpha;}{K}_g^m(x_{1j}^m,x_{2k}^m)+(1-\mathrm{\&alpha;})K_w^m(x_{1j}^m,x_{2k}^m)+\mathrm{\&alpha;}{K}_g^m(x_{2j}^m,x_{1k}^m)+(1-\mathrm{\&alpha;})K_w^m(x_{2j}^m,x_{1k}^m),$

$=\mathrm{\&alpha;}[K_{\mathrm{gas}}(g_{X_1}^m(p_j,q_j),g_{X_1}^m(p_k,q_k))+K_{\mathrm{gas}}(g_{X_2}^m(p_j,q_j),g_{X_2}^m(p_k,q_k))$

$-K_{\mathrm{gas}}(g_{X_1}^m(p_j,q_j),g_{X_2}^m(p_k,q_k))-K_{\mathrm{gas}}(g_{X_2}^m(p_j,q_j),g_{X_1}^m(p_k,q_k))]$

$+(1-\mathrm{\&alpha;})[K_{\mathrm{gas}}(w_{X_1}^m(p_j,q_j),w_{X_1}^m(p_k,q_k))+K_{\mathrm{gas}}(w_{X_2}^m(p_j,q_j),w_{X_2}^m(p_k,q_k))$

$-K_{\mathrm{gas}}(w_{X_1}^m(p_j,q_j),w_{X_2}^m(p_k,q_k))-K_{\mathrm{gas}}(w_{X_2}^m(p_j,q_j),w_{X_1}^m(p_k,q_k))]$

其中，

为核空间中第m个尺度上的标签训练样本(p_j，q_j)处的差特征，为核空间中第m个尺度上的标签训练样本(p_k，q_k)处的差特征，

表示第m个尺度上的标签训练样本(p_j，q_j)处的差特征，

表示第m个尺度上的标签训练样本(p_k，q_k)处的差特征，K_gas为高斯径向基核函数，α为大于0小于1的常数，1≤j，k≤N_s，m=1，...，S。

步骤8，对支撑矢量积进行训练。

将步骤7得到的S个尺度上的差值合成核矩阵

输入到支撑矢量积SVM中对SVM进行训练，得到S个尺度上的支撑矢量分类器SVC_m，其中m=1，2...S。

二.测试步骤：

步骤A，构造强度核和纹理核。

(A1)将步骤3得到的两时相图像的强度特征和纹理特征以及步骤4得到的两时相标签训练样本的强度特征和纹理特征映射到核空间，得到核空间中的两时相的强度特征和纹理特征；

(A2)利用步骤(A1)得到的两时相的强度特征和纹理特征构造S个尺度上的强度核

和纹理核

$K_g^m(x_{\mathrm{ij}}^m,T_i^m(p,q))=\&lang;{\mathrm{\&phi;}}_g\left(x_{\mathrm{ij}}^m\right),{\mathrm{\&phi;}}_g\left(T_i^m(p,q)\right)\&rang;,$

$=K_{\mathrm{gas}}(g_{X_i}^m(p_j,q_j),g_{X_i}^m(p,q))$

$K_w^m(x_{\mathrm{ij}}^m,T_i^m(p,q))=\&lang;{\mathrm{\&phi;}}_w\left(x_{\mathrm{ij}}^m\right),{\mathrm{\&phi;}}_w\left(T_i^m(p,q)\right)\&rang;$

$=K_{\mathrm{gas}}{(w_{X_i}^m(p_j,q_j),w_{X_i}^m(p,q))}^{,}$

其中，K_gas为高斯径向基核函数，

为像素点(p，q)处的特征，

表示第m个尺度上的像素点(p，q)处的强度特征，

表示第m个尺度上的像素点(p，q)处的纹理特征，

为核空间中第m个尺度上的标签训练样本(p_j，q_j)处的强度特征，为核空间中第m个尺度上的标签训练样本(p_j，q_j)处的纹理特征，φ_g(T_i ^m(p，q))为核空间中第m个尺度上的像素点(p，q)处的强度特征，φ_w(T_i ^m(p，q))为核空间中第m个尺度上的像素点(p，q)处的纹理特征，1≤k，j≤N_s，m=1，...，S。

步骤B，构造强度纹理合成核。

(B1)将步骤(A1)得到的核空间中的两时相的强度特征和纹理特征进行加权求和，得到核空间中的两时相的强度纹理合成特征；

(B2)利用步骤(B1)得到的两时相的强度纹理合成特征构造S个尺度上的强度纹理合成核

其中，

为核空间中第m个尺度上的标签训练样本(p_j，q_j)处的强度纹理合成特征，为核空间中第m个尺度上的像素点(p，q)处的强度纹理合成特征，α为大于0小于1的常数。

步骤C，构造差值合成核。

(C1)对步骤(B1)得到的两时相的强度纹理合成特征做差，得到差特征；

(C2)利用步骤(C1)得到的差特征构造S个尺度上的差值合成核 $K_D^m(x_j^m,T^m(p,q)):$

$K_D^m(x_j^m,T^m(p,q))=\&lang;\mathrm{\&Phi;}\left(x_j^m\right),\mathrm{\&Phi;}\left(T^m(p,q)\right)\&rang;$

$=K_{\mathrm{gw}}^m(x_{1j}^m,T_1^m(p,q))+K_{\mathrm{gw}}^m(x_{2j}^m,T_2^m(p,q))-K_{\mathrm{gw}}^m(x_{1j}^m,T_2^m(p,q))-K_{\mathrm{gw}}^m(x_{2j}^m,T_1^m(p,q))$

$={\mathrm{\&alpha;K}}_g^m(x_{1j}^m,T_1^m(p,q))+(1-\mathrm{\&alpha;})K_w^m(x_{1j}^m,T_1^m(p,q))+{\mathrm{\&alpha;K}}_g^m(x_{2j}^m,T_2^m(p,q))+(1-\mathrm{\&alpha;})K_w^m(x_{2j}^m,T_2^m(p,q))$

$-{\mathrm{\&alpha;K}}_g^m(x_{1j}^m,T_m^2(p,q))-(1-\mathrm{\&alpha;})K_w^m(x_{1j}^m,T_2^m(p,q))-{\mathrm{\&alpha;K}}_g^m(x_{2j}^m,T_1^m(p,q))-(1-\mathrm{\&alpha;})K_w^m(x_{2j}^m,T_1^m(p,q){)}^{,}$

$=\mathrm{\&alpha;}[K_{\mathrm{gas}}(g_{X_1}^m(p_j,q_j),g_{X_1}^m(p,q))+K_{\mathrm{gas}}(g_{X_2}^m(p_j,q_j),g_{X_2}^m(p,q))$

$-K_{\mathrm{gas}}(g_{X_1}^m(p_j,q_j),g_{X_2}^m(p,q))-K_{\mathrm{gas}}(g_{X_2}^m\left(p_{j,}{q}_j\right),g_{X_1}^m(p,q))]$

$+(1-\mathrm{\&alpha;})[K_{\mathrm{gas}}(w_{X_1}^m(p_j,q_j),w_{X_1}^m(p,q))+K_{\mathrm{gas}}(w_{X_2}^m(p_j,q_j),w_{X_2}^m(p,q))$

$-K_{\mathrm{gas}}(w_{X_1}^m(p_j,q_j),w_{X_2}^m(p,q))-K_{\mathrm{gas}}(w_{X_2}^m(p_j,q_j),w_{X_1}^m(p,q))]$

其中，Φ(T^m(p，q))为核空间中第m个尺度上的像素点(p，q)处的差特征，T^m(p，q)表示第m个尺度上的像素点(p，q)处的差特征，K_gas为高斯径向基核函数，α为大于0小于1的常数，1≤j，k≤N_s，m＝1，...，S。

步骤D，对支撑矢量分类器进行测试。

将差值合成核矩阵

输入到步骤8得到的支撑矢量分类器SVC_m中对SVC_m进行测试，得到S个尺度上的二值变化检测结果图BCM_m。

步骤E，尺度间融合。

对二值变化检测结果图BCM_m进行决策级尺度间融合，得到最终的变化检测结果图，融合过程通过如下公式进行：

$\mathrm{CM}=\underset{m=1}{\overset{S}{\&cap;}}{\mathrm{BCM}}_m,$

其中，∩表示二元交运算，CM为最终的变化检测结果图，S为尺度数。

本发明效果可以通过以下实验进一步证实：

1)实验条件

实验仿真环境为：MATLAB R2009b，Intel(R)Pentium(R)2CPU2.7GHz，Window7旗舰版。

2)实验内容

实验1，分别应用Gabor域合成核方法和本发明方法对两时相洪水灾害RealSAR图像进行变化检测，实验结果如图2所示，其中图2(a)是第一时刻洪水灾害JERS SAR channel1图像，图2(b)是第二时刻洪水灾害JERS SAR channel1图像，图2(c)是变化检测结果参考图，图2(d)是Gabor域合成核的变化检测结果，图2(e)是本发明的变化检测结果。

实验2，分别应用Gabor域合成核方法和本发明方法对两时相机场跑道RealSAR图像进行变化检测，实验结果如图3所示，其中图3(a)是第一时刻机场跑道ERS-1SAR图像，图3(b)是第二时刻机场跑道ERS-1SAR图像，图3(c)是变化检测结果参考图，图3(d)是Gabor域合成核变化检测结果，图3(e)是本发明的检测结果。

实验3，分别应用Gabor域合成核方法和本发明方法对两时相农田Real SAR图像进行变化检测，实验结果如图4所示，其中图4(a)是第一时刻农田SAR图像，图4(b)是第二时刻农田SAR图像，图4(c)是变化检测结果参考图，图4(d)是Gabor域合成核变化检测结果，图4(e)是本发明的变化检测结果。

3)实验结果

为验证变化检测结果质量，本发明选取总体精度、Kappa系数和执行时间作为性能指标参数，评价本发明方法的检测精度和效率。该实验结果如表1和表2所示，表1给出Gabor域合成核方法和本发明方法检测精度的比较结果，表2给出Gabor域合成核方法和本发明方法执行效率的比较结果。

表1本发明方法与Gabor域合成核方法检测精度比较

4)实验结果分析

由图2、图3、图4和表1可以看到，本发明方法的检测精度高于Gabor域合成核方法，这是由于本发明利用的NSCT变换相对Gabor变换不仅能对图像进行多分辨分析，也可以进行时频局部分析，并且在特征提取时具有的多方向性、各向异性和移不变性使得本发明能得到更优越的强度特征和纹理特征。

由表2中可以看出，本发明方法在执行效率上优于现有Gabor域合成核方法，这是由于NSCT变换的复杂度低于Gabor滤波器的复杂度。

Claims (7)

1.一种基于NSCT域合成核的SAR图像变化检测方法，包括：

A训练步骤：

A1)使用非下采样轮廓波变换NSCT将原始两时相SAR图像X₁={X₁(p，q)|1≤p≤I，1≤q≤J}和X₂={X₂(p，q)|1≤p≤I，1≤q≤J}进行S个尺度的分解，得到S个尺度上的两时相图像，其中I，J分别为原始图像长度和宽度大小，(p，q)为图像像素点；

A2)对步骤A1)得到的S个尺度上的两时相图像提取强度特征和纹理特征，并进行归一化，得到归一化后的S个尺度上的强度特征

和纹理特征

再将归一化后的这两种特征联合，构成S个尺度上的两时相图像特征

$T_i^m=\{g_{X_i}^m,w_{X_i}^m\},$

其中，m=1，2...S，i＝1，2，i是指示原始两时相SAR图像X₁或者X₂的下标；

A3)从步骤A1)得到的S个尺度的两时相图像中，手动提取N_s个标签训练样本，将这些训练样本用步骤A2)得到的归一化后的S个尺度上的强度特征和纹理特征进行表示，得到S个尺度上的两时相的标签训练样本特征

m=1，2...S，1≤j≤N_s，i=1，2，

其中， $x_{\mathrm{ij}}^m=\{g_{X_i}^m(p_j,q_j),w_{X_i}^m(p_j,q_j)\},$

(p_j，q_j)为标签训练样本点的坐标，

表示第m个尺度上的标签训练样本点(p_j，q_j)处的强度特征，

表示第m个尺度上的标签训练样本点(p_j，q_j)处的纹理特征，

表示第m个尺度上的标签训练样本点(p_j，q_j)处的标签；

A4)将步骤A3)得到的两时相标签训练样本的强度特征和纹理特征映射到核空间，得到核空间中的两时相标签训练样本的强度特征和纹理特征，然后利用这两种特征构造S个尺度上的强度核

和纹理核

其中

为第m个尺度上的标签训练样本点(p_k，q_k)处的特征，表示第m个尺度上的标签训练样本点(p_k，q_k)处的强度特征，

表示第m个尺度上的标签训练样本点(p_k，q_k)处的纹理特征；

A5)将步骤A4)得到的核空间中的两时相标签训练样本的强度特征和纹理特征进行加权求和，得到核空间中的两时相的强度纹理合成特征，利用该合成特征构造S个尺度上的强度纹理合成核

A6)对步骤A5)得到的两时相的标签训练样本的强度纹理合成特征做差，得到标签训练样本的差特征，利用该差特征构造S个尺度上的差值合成核

其中

表示第m个尺度上的标签训练样本(p_j，q_j)处的差特征，

表示第m个尺度上的标签训练样本(p_k，q_k)处的差特征，1≤j，k≤N_s，m=1，...，S；

A7)将步骤A6)得到的S个尺度上的差值合成核矩阵

输入到支撑矢量积SVM中对SVM进行训练，得到S个尺度上的支撑矢量分类器SVC_m，其中m=1，2...S；

B测试步骤：

B1)将步骤A2)得到的两时相图像的强度特征和纹理特征以及步骤A3)得到的两时相标签训练样本的强度特征和纹理特征映射到核空间，得到核空间中的两时相的强度特征和纹理特征，利用这两个特征构造S个尺度上的强度核

和纹理核 $K_w^m(x_{\mathrm{ij}}^m,T_i^m(p,q)),$ 其中， $T_i^m(p,q)=\{g_{X_i}^m(p,q),w_{X_i}^m(p,q)\}$ 为像素点(p，q)处的特征，

表示第m个尺度上的像素点(p，q)处的强度特征，表示第m个尺度上的像素点(p，q)处的纹理特征；

B2)将步骤B1)得到的核空间中的两时相的强度特征和纹理特征进行加权求和，得到核空间中的两时相的强度纹理合成特征，利用该合成特征构造S个尺度上的强度纹理合成核

B3)对步骤B2)得到的两时相的强度纹理合成特征做差，得到差特征，利用该差特征构造S个尺度上的差值合成核

其中表示在第m个尺度上的标签训练样本(p_j，q_j)处的差特征，T^m(p，q)表示第m个尺度上的像素点(p，q)处的差特征，1≤j≤N_s，m=1，...，S；

B4)将差值合成核矩阵输入到步骤A7)得到的支撑矢量分类器SVC_m中，对SVC_m进行测试，得到S个尺度上的二值变化检测结果图BCM_m；

B5)对二值变化检测结果图BCM_m进行决策级尺度间融合，得到最终的变化检测结果图。

2.根据权利要求1所述方法，其中所述步骤A2)中提取S个尺度上的两时相图像强度特征和纹理特征，按如下步骤进行：

A21)令

表示第m个尺度上的两时相图像的低通子带系数，并用该低通子带系数

作为两时相图像在第m个尺度上的强度特征，其中m=1，...，S，S为尺度数，i是指示原始两时相SAR图像X₁或者X₂的下标；

A22)令

表示第m个尺度上的两时相图像的高通子带系数，d表示该高通子带系数的方向，在以(p，q)为中心像素、大小为N的窗口上提取高通子带系数

的均值信息

和方差信息

${\mathrm{\&mu;}}_{X_i}^{m,d}(p,q)=\frac{1}{N^2}\underset{j=p-N/2}{\overset{p+N/2}{\mathrm{\&Sigma;}}}\underset{l=q-N/2}{\overset{q+N/2}{\mathrm{\&Sigma;}}}{H}_{X_i}^{m,d}(j,l),$

${\mathrm{\&sigma;}}_{X_i}^{m,d}(p,q)=\frac{1}{N^2-1}\underset{j=p-N/2}{\overset{p+N/2}{\mathrm{\&Sigma;}}}\underset{l=q-N/2}{\overset{q+N/2}{\mathrm{\&Sigma;}}}{(H_{X_i}^{m,d}(j,l)-{\mathrm{\&mu;}}_{X_i}^{m,d}(p,q))}^2,$

其中，d=1，...，n_m，n_m为尺度m上的方向数；

A23)以第m个尺度上的两时相图像中所有像素点为中心像素点，提取上述均值信息和方差信息，得到均值矢量和方差矢量

A24)用均值矢量

和方差矢量联合构成两时相图像在第m个尺度上的纹理特征。

3.根据权利要求1所述方法，其中所述步骤A2)中对S个尺度上的两时相图像提取的强度特征和纹理特征进行归一化，按如下公式进行：

$g_{X_i}^m=\frac{L_{X_i}^m-\min \left(L_{X_i}^m\right)}{\max \left(L_{X_i}^m\right)-\min \left(L_{X_i}^m\right)},$

${\widetilde{\mathrm{\&mu;}}}_{X_i}^{m,d}=\frac{{\mathrm{\&mu;}}_{X_i}^{m,d}-\min \left({\mathrm{\&mu;}}_{X_i}^{m,d}\right)}{\max \left({\mathrm{\&mu;}}_{X_i}^{m,d}\right)-\min \left({\mathrm{\&mu;}}_{X_i}^{m,d}\right)},$

${\widetilde{\mathrm{\&sigma;}}}_{X_i}^{m,d}=\frac{{\mathrm{\&sigma;}}_{X_i}^{m,d}-\min \left({\mathrm{\&sigma;}}_{X_i}^{m,d}\right)}{\max \left({\mathrm{\&sigma;}}_{X_i}^{m,d}\right)-\min \left({\mathrm{\&sigma;}}_{X_i}^{m,d}\right)},$

其中

是第m个尺度上的两时相图像的归一化后的强度特征，是归一化后的均值矢量，

是归一化后的方差矢量，

和

联合构成归一化后的第m个尺度上的两时相图像的纹理特征

4.根据权利要求1所述方法，其中所述步骤A4)中强度核

和纹理核

其公式如下：

$K_g^m(x_{\mathrm{ij}}^m,x_{\mathrm{ik}}^m)=K_{\mathrm{gas}}(g_{X_i}^m(p_j,q_j),g_{X_i}^m(p_k,q_k)),$

$K_w^m(x_{\mathrm{ij}}^m,x_{\mathrm{ik}}^m)=K_{\mathrm{gas}}(w_{X_i}^m(p_j,q_j),w_{X_i}^m(p_k,q_k)),$

其中，K_gas为高斯径向基核函数，1≤k，j≤N_s，m=1，...，S。

5.根据权利要求1所述方法，其中所述步骤A5)中的强度纹理合成核

其公式如下：

$K_{\mathrm{gw}}^m(x_{\mathrm{ij}}^m,x_{\mathrm{ik}}^m)=\mathrm{\&alpha;}{K}_g^m(x_{\mathrm{ij}}^m,x_{\mathrm{ik}}^m)+(1-\mathrm{\&alpha;})K_w^m(x_{\mathrm{ij}}^m,x_{\mathrm{ik}}^m),$

其中α为大于0小于1的常数。

6.根据权利要求1所述方法，其中所述步骤A6)中差值合成核其公式如下：

$\begin{array}{c}{K}_D^m(x_j^m,x_k^m)=K_{\mathrm{gw}}^m(x_{1j}^m,x_{1k}^m)+K_{\mathrm{gw}}^m(x_{2j}^m,x_{2k}^m)-K_{\mathrm{gw}}^m(x_{1j}^m,x_{2k}^m)-K_{\mathrm{gw}}^m(x_{2j}^m,x_{1k}^m)\\ ={\mathrm{\&alpha;K}}_g^m(x_{1j}^m,x_{1k}^m)+(1-\mathrm{\&alpha;})K_w^m(x_{1j}^m,x_{1k}^m)+{\mathrm{\&alpha;K}}_g^m(x_{2j}^m,x_{2k}^m)+(1-\mathrm{\&alpha;})K_w^m(x_{2j}^m,x_{2k}^m)\\ -{\mathrm{\&alpha;K}}_g^m\left(x_{1j}^m{,x}_{2k}^m\right)-(1-\mathrm{\&alpha;})K_w^m(x_{1j}^m,x_{2k}^m)-{\mathrm{\&alpha;K}}_g^m(x_{2j}^m,x_{1k}^m)-(1-\mathrm{\&alpha;})K_w^m(x_{2j}^m,x_{1k}^m)\\ =\mathrm{\&alpha;}\&lsqb;K_{\mathrm{gas}}(g_{X_1}^m(p_j,q_j),g_{X_1}^m(p_k,q_k))-K_{\mathrm{gas}}(g_{X_2}^m(p_j,q_j),g_{X_2}^m(p_k,q_k))\\ -K_{\mathrm{gas}}(g_{X_1}^m(p_j,q_j),g_{X_2}^m(p_k,q_k))-K_{\mathrm{gas}}(g_{X_2}^m(p_j,q_j),g_{X_1}^m(p_k,q_k))\&rsqb;\\ +(1-\mathrm{\&alpha;})\&lsqb;K_{\mathrm{gas}}(w_{X_1}^m(p_j,q_j),w_{X_1}^m(p_k,q_k))+K_{\mathrm{gas}}(w_{X_2}^m(p_j,q_j),w_{X_2}^m(p_k,q_k))\\ -K_{\mathrm{gas}}(w_{X_1}^m(p_j,q_j),w_{X_2}^m(p_k,q_k))-K_{\mathrm{gas}}(w_{X_2}^m(p_j,q_j),w_{X_1}^m(p_k,q_k))\&rsqb;\end{array},$

其中，K_gas为高斯径向基核函数，α为大于0小于1的常数。

7.根据权利要求1所述方法，其中所述步骤B5)中对二值变化检测结果图BCM_m进行决策级尺度间融合，通过如下公式进行：

$\mathrm{CM}=\underset{m=1}{\overset{S}{\&cap;}}{\mathrm{BCM}}_m,$

其中∩表示二元交运算，CM为最终的变化检测结果图，S为尺度数。

Images