CN102722892B

CN102722892B - SAR (synthetic aperture radar) image change detection method based on low-rank matrix factorization

Info

Publication number: CN102722892B
Application number: CN201210193347.3A
Authority: CN
Inventors: 张向荣; 焦李成; 张祎勃; 郑耀国; 李阳阳; 侯彪; 白静; 杨阳; 翁鹏
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2012-06-13
Filing date: 2012-06-13
Publication date: 2014-11-12
Anticipated expiration: 2032-06-13
Also published as: CN102722892A

Abstract

The invention discloses a SAR image change detection method based on low-rank matrix decomposition, which mainly solves the problem that the existing method cannot accurately detect the SAR image change area. The implementation steps include: (1) pre-processing the two SAR images to be detected to obtain a smoother SAR image; (2) constructing the logarithmic ratio of the two images after speckle reduction; (3) converting the logarithm The low-rank sparse decomposition of the ratio is performed to obtain the low-rank part and the sparse part of the logarithmic ratio; (4) Transform the sparse part into a sparse matrix by column; (5) Cluster the obtained sparse matrix with the K-means algorithm to obtain the final change detection results. The invention has the advantages of relatively accurate detection of changing areas, and can be used in the fields of public security and video surveillance.

Description

SAR image change detection method based on low-rank matrix factorization

技术领域 technical field

本发明属于雷达图像处理技术领域，特别涉及SAR图像变化检测方法，可用于解决在SAR图像变化检测中检测变化精度不高的问题。The invention belongs to the technical field of radar image processing, in particular to a SAR image change detection method, which can be used to solve the problem of low detection change accuracy in the SAR image change detection.

背景技术 Background technique

图像变化检测方法是一种分析和理解多时遥感图像的重要技术，近些年引起了广泛的研究。这源于变化检测方法广泛的应用背景，例如农业调查，森林监测，自然灾害监测，城市变化分析，战场打击效果评估等方面。Image change detection method is an important technique for analyzing and understanding multi-temporal remote sensing images, which has attracted extensive research in recent years. This stems from the wide application background of change detection methods, such as agricultural survey, forest monitoring, natural disaster monitoring, urban change analysis, battlefield strike effect evaluation and so on.

图像变化检测是对从同一地区不同时间获取的多时遥感图像分析的一种方法。它侧重的是识别两幅遥感图像中地物的变化。现有的变化检测方法主要可以分为两大类：有监督方法和无监督方法。Image change detection is a method for analyzing multi-temporal remote sensing images acquired from the same area at different times. It focuses on identifying changes in ground objects in two remote sensing images. Existing change detection methods can be mainly divided into two categories: supervised methods and unsupervised methods.

所谓有监督方法，是基于监督分类法，需要获取图像变化区域的训练样区，从而进行变化检测；而无监督法不需要任何额外的信息，直接对两个不同时相的数据检测。尽管有监督法准确确定变化区域相比无监督法有明显的优势，但是由于得到地面的真实信息特别困难，因此无监督变化检测方法是常用的变化检测方法。The so-called supervised method is based on the supervised classification method, which needs to obtain the training sample area of the image change area to perform change detection; while the unsupervised method does not require any additional information, and directly detects data of two different phases. Although the supervised method has obvious advantages over the unsupervised method in accurately determining the change area, it is particularly difficult to obtain the ground truth information, so the unsupervised change detection method is a commonly used change detection method.

无监督变化检测法，主要包括主成份分析法、小波融合法。主成份分析法通过对数比值法构造差异图，再在差异图的基础上，采用主成分分析降维来提取差异图特征，然后再用k-均值对这些特征聚类。小波融合法是对待检测图像分别进行均值比和对数比操作，用离散小波变换分别提取高低频信息，对其进行稀疏融合。这两种无监督方法的缺点是检测精度比较低。Unsupervised change detection methods mainly include principal component analysis and wavelet fusion. The principal component analysis method constructs a difference map through the logarithmic ratio method, and then uses principal component analysis to reduce the dimensionality of the difference map to extract the features of the difference map, and then uses k-means to cluster these features. The wavelet fusion method is to perform mean ratio and logarithmic ratio operations on the image to be detected, and use discrete wavelet transform to extract high and low frequency information respectively, and perform sparse fusion on them. The disadvantage of these two unsupervised methods is that the detection accuracy is relatively low.

发明内容 Contents of the invention

本发明的目的在于克服现有技术的不足，提出一种新的基于低秩矩阵分解的SAR图像变化检测方法，以获得较为精确的变化差分图，提高SAR图像变化的检测精度。The purpose of the present invention is to overcome the deficiencies of the prior art and propose a new SAR image change detection method based on low-rank matrix decomposition to obtain a more accurate change difference map and improve the detection accuracy of SAR image changes.

实现本发明目的技术方案是：Realize the technical scheme of the object of the present invention is:

一．技术原理one. Technical principle

在压缩感知领域，如果一个矩阵具有唯一的“低秩和稀疏”结构，则在适当的条件下该矩阵就能够被精确的恢复。在矩阵的“低秩和稀疏”结构的应用中，一个杰出的工作是鲁棒的主分量分析RPCA，用RPCA对图像矩阵的“低秩和稀疏”分解后的低秩部分相当于图像矩阵的主要成分，稀疏部分相当于噪声或其它内容。另一个代表性工作是随机逼近的矩阵分解，它证明了对于一个给定的图像矩阵，图像矩阵的主要成份可以通过随机投影进行逼近。在计算机视觉领域，这些算法具有很大的应用背景，比如，背景建模，人脸图像的光照或阴影去除，图像校正等等。In the field of compressed sensing, if a matrix has a unique "low-rank and sparse" structure, the matrix can be accurately recovered under appropriate conditions. In the application of the "low-rank and sparse" structure of the matrix, an outstanding work is the robust principal component analysis RPCA, and the low-rank part of the "low-rank and sparse" decomposition of the image matrix with RPCA is equivalent to the image matrix. The main component, the sparse part is equivalent to noise or other content. Another representative work is matrix factorization for random approximation, which proves that for a given image matrix, the principal components of the image matrix can be approximated by random projections. In the field of computer vision, these algorithms have a great application background, such as background modeling, illumination or shadow removal of face images, image correction and so on.

本发明，受启发于上述思想，把矩阵的“低秩和稀疏”分解应用到SAR图像的变化检测中，其中的低秩部分代表的是图像中未变化的部分，稀疏部分代表的是变化部分，这正是变化检测任务所需要的结果。Inspired by the above ideas, the present invention applies the "low-rank and sparse" decomposition of the matrix to the change detection of SAR images, wherein the low-rank part represents the unchanged part of the image, and the sparse part represents the changed part , which is exactly the desired result for the change detection task.

二．实现方案two. Implementation plan

本发明基于低秩矩阵分解的SAR图像变化检测方法，包括如下步骤：The present invention is based on the SAR image change detection method of low-rank matrix decomposition, comprises the following steps:

（1）对输入的两幅同一区域不同时间,大小相等的SAR图像，用基于概率块的去噪算法PPB分别对这两幅图像进行降斑，降斑后的图像矩阵分别记为X₁和X₂：(1) For the two input SAR images of the same area at different times and equal in size, use the probabilistic block-based denoising algorithm PPB to reduce speckle on the two images respectively, and the image matrices after speckle reduction are denoted as X ₁ and _X2 :

（2）将两幅降斑后的图像矩阵X₁和X₂按如下公式求取这两个图像矩阵的对数比矩阵D_l：(2) Calculate the logarithmic ratio matrix D _l of the two image matrices X ₁ and X ₂ after speckle reduction according to the following formula:

${D D.}_{l l} = = | | log log \frac{{X x}_{22} + + 11}{{X x}_{11} + + 11} | | = = | | log log (({X x}_{22} + + 11)) - - log log (({X x}_{11} + + 11)) | |;;$

（3）求取对数比矩阵D_l的低秩部分和稀疏部分：(3) Find the low-rank part and sparse part of the logarithmic ratio matrix D _l :

3a)如果对数比矩阵D_l的列数为奇数，则删去D_l的最后一列，如果D_l的行数为奇数，则删去D_l的最后一行，把变化后的D_l记为D_l′：3a) If the number of columns of the logarithmic ratio matrix _D1 is an odd number, then delete the last column of _D1 , if the number of rows of _D1 is an odd number, then delete the last row of _D1 , and record the changed _D1 as D _l ':

3b）把矩阵D_l′分成2×2的小矩阵，并把每一个小矩阵变成一个列向量，把所有列向量依次横向合并成为一个变化矩阵，记为D_l″：3b) Divide the matrix D _l ′ into 2×2 small matrices, and turn each small matrix into a column vector, and merge all the column vectors horizontally in turn to form a change matrix, which is denoted as D _l ″:

3c)令GoDec算法中的输入初始低秩部分L₀＝D_l″，初始稀疏部分S₀＝0，用GoDec算法把变化矩阵D_l″分解成低秩部分L和稀疏部分S：3c) Let the input initial low-rank part L ₀ =D _l ″ in the GoDec algorithm, the initial sparse part S ₀ =0, use the GoDec algorithm to decompose the change matrix D _l ″ into a low-rank part L and a sparse part S:

（4）把步骤（3）得到的稀疏部分S按列变换成和初始稀疏部分S₀维数一致的稀疏矩阵,记为S′：(4) Transform the sparse part S obtained in step (3) into a sparse matrix with the same dimension as the initial sparse part S ₀ by column, denoted as S′:

（5）使用K均值算法把步骤（4）得到的稀疏矩阵S′聚成2类，得到的结果即为所求的变化检测图。(5) Use the K-means algorithm to cluster the sparse matrix S′ obtained in step (4) into two categories, and the obtained result is the desired change detection map.

本发明使用了基于概率块的去噪算法PPB对两幅输入图像进行降斑，用两个降斑结果求取对数比矩阵，然后用GoDec算法对对数比矩阵进行低秩和稀疏分解，获得了更为精确的变化差分图，最后用K均值算法对变化差分图进行聚类，使SAR图像的变化检测精度比现有的检测精度有一定的提高。The present invention uses the probability block-based denoising algorithm PPB to reduce speckles on two input images, uses the two speckle reduction results to obtain a logarithmic ratio matrix, and then uses the GoDec algorithm to perform low-rank and sparse decomposition on the logarithmic ratio matrix, A more accurate change difference map is obtained, and finally the K-means algorithm is used to cluster the change difference map, so that the change detection accuracy of SAR images is improved to a certain extent compared with the existing detection accuracy.

附图说明 Description of drawings

图1是本发明的实现流程图；Fig. 1 is the realization flowchart of the present invention;

图2是用本发明对Bern地区图像仿真时所用的两幅输入SAR图像；Fig. 2 is two pieces of input SAR images used when using the present invention to Bern area image simulation;

图3是Bern地区的参考图；Figure 3 is a reference map of the Bern area;

图4是用本发明对Ottawa地区图像仿真时所用的两幅输入SAR图像；Fig. 4 is two pieces of input SAR images used when using the present invention to Ottawa region image simulation;

图5是Ottawa地区的参考图；Figure 5 is a reference map of the Ottawa area;

图6是对图2的仿真结果图；Fig. 6 is the simulation result figure to Fig. 2;

图7是对图4的仿真结果图。FIG. 7 is a simulation result diagram of FIG. 4 .

具体实施方式 Detailed ways

参照图1，本发明的具体实施步骤如下：With reference to Fig. 1, concrete implementation steps of the present invention are as follows:

步骤1,输入两幅同一区域、不同时间,大小相等的SAR图像，并进行降斑处理：Step 1, input two SAR images of the same area, different time, and equal size, and perform speckle reduction processing:

（1a）以第一幅输入SAR图像中任一像素点c_s为中心，选取N×N大小的邻域作为该像素点的搜寻区域，其中N＝21：(1a) Take any pixel point c _s in the first input SAR image as the center, select a neighborhood of N×N size as the search area of the pixel point, where N=21:

（1b）以像素点c_s为中心，取M×M大小的块，其中M＝7，用块内所有像素点的灰度值构成矩阵v_s：(1b) Take the pixel c _s as the center, take a block of M×M size, where M=7, and use the gray values of all pixels in the block to form a matrix v _s :

（1c）以搜寻区域中除中心像素点c_s外的每一个像素点f_t为中心，取M×M大小的块，块内所有像素点的灰度值构成矩阵v_t：(1c) Taking each pixel f _t in the search area except the central pixel c _s as the center, take a block of M×M size, and the gray values of all pixels in the block form a matrix v _t :

(1d）按照下面权重公式计算经过i-1次去噪后f_t到c_s的权重w_s，t ^i-1：(1d) Calculate the weight w _s , t i-1 from f _t to c _s after i- ¹ times of denoising according to the following weight formula:

如果第一幅输入SAR图像是强度图像，使用权重公式：If the first input SAR image is an intensity image, use the weighting formula:

${w w}_{s the s,, t t}^{i i - - 11} = = exp exp [[- - \underset{k k}{Σ Σ} ((\frac{11}{h h} [[((22 L L - - 11)) log log (({v v}_{s the s,, k k} + + {v v}_{t t,, k k})) - - ((L L - - 11)) log log (({v v}_{s the s,, k k} {v v}_{t t,, k k}))]] + + \frac{L L}{T T} \frac{| | {\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} - - {\overset{^^}{u u}}_{t t,, k k}^{i i - - 11}}{{\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} {\overset{^^}{u u}}_{t t,, k k}^{i i - - 11}}))]],,$

如果第一幅输入SAR图像是幅度图像，则使用权重公式：If the first input SAR image is a magnitude image, the weighting formula is used:

${w w}_{s the s,, t t}^{i i - - 11} = = exp exp [[- - \underset{k k}{Σ Σ} ((\frac{11}{\overset{~ ~}{h h}} log log ((\frac{{v v}_{s the s,, k k}}{{v v}_{t t,, k k}} + + \frac{{v v}_{t t,, k k}}{{v v}_{s the s,, k k}})) + + \frac{L L}{T T} \frac{| | {\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} - - {\overset{^^}{u u}}_{t t,, k k}^{i i - - 11}}{{\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} {\overset{^^}{u u}}_{t t,, k k}^{i i - - 11}}))]]$

其中 ${\hat{u}}_{s}^{i - 1} = \frac{Σ_{t} w_{s, t} {(v_{t}^{i - 1})}^{2}}{Σ_{t} w_{s, t}},$ ${\hat{u}}_{t}^{i - 1} = V_{t}^{i - 1},$ $\tilde{h} = \frac{h}{2 L - 1},$ v_s，k和v_t，k分别表示v_s和v_t的第k个像素点的值，和分别是和的第k个像素点的值，exp是以自然对数e为底数的指数函数，h＝11，L＝3，T＝0.2：in ${\hat{u}}_{the s}^{i - 1} = \frac{Σ_{t} w_{the s, t} {(v_{t}^{i - 1})}^{2}}{Σ_{t} w_{the s, t}},$ ${\hat{u}}_{t}^{i - 1} = V_{t}^{i - 1},$ $\tilde{h} = \frac{h}{2 L - 1},$ v _{s, k} and v _{t, k} respectively represent the value of the kth pixel of v _s and v _t , and respectively and The value of the kth pixel, exp is an exponential function with natural logarithm e as the base, h=11, L=3, T=0.2:

（1e）按照下面公式计算像素点c_s经过i次去噪后的值 (1e) According to the following formula, calculate the value of pixel c _s after i times of denoising

${c c}_{s the s}^{i i} = = \frac{{Σ Σ}_{t t} {w w}_{s the s,, t t}^{i i - - 11} {f f}_{t t}}{{Σ Σ}_{t t} {w w}_{s the s,, t t}^{i i - - 11}}$

（1f）对于每一个像素点，重复上述步骤（1a)~（1e)，得到第i次去噪后的图像 (1f) For each pixel point, repeat the above steps (1a)~(1e) to obtain the i-th denoised image

（1g）令i＝i+1，以去噪后图像作为新的第一幅输入SAR图像，重复进行步骤（1a)~（1f)所述的去噪过程15次，从而得到最终降斑后的图像矩阵X₁：(1g) Let i=i+1, to denoise the image As the new first input SAR image, repeat the denoising process described in steps (1a) to (1f) 15 times, so as to obtain the final speckle-reduced image matrix X ₁ :

（1h）按步骤（1a)~（1g)所述的降斑过程对第二幅输入SAR图像进行降斑处理，得到降斑后的图像矩阵X₂。(1h) Perform speckle reduction processing on the second input SAR image according to the speckle reduction process described in steps (1a) to (1g), to obtain a speckle-reduced image matrix X ₂ .

步骤2,将两幅降斑后的图像矩阵X₁和X₂，按如下公式求取这两个图像矩阵的对数比矩阵D_l：Step 2: Calculate the logarithmic ratio matrix D _l of the two image matrices X ₁ and X ₂ after speckle reduction according to the following formula:

${D D.}_{l l} = = | | log log \frac{{X x}_{22} + + 11}{{X x}_{11} + + 11} | | = = | | log log (({X x}_{22} + + 11)) - - log log (({X x}_{11} + + 11)) | | . .$

步骤3,求取对数比矩阵D_l的低秩部分和稀疏部分：Step 3, find the low-rank part and sparse part of the logarithmic ratio matrix D _l :

3a)如果对数比矩阵D_l的列数为奇数，则删去D_l的最后一列，如果D_l的行数为奇数，则删去D_l的最后一行，把变化后的对数比矩阵记为D_l′：3a) If the number of columns of the logarithmic ratio matrix _D1 is an odd number, then delete the last column of _D1 , if the number of rows of _D1 is an odd number, then delete the last row of _D1 , and change the logarithmic ratio matrix Denote as D _l ′:

3b）将变化后的对数比矩阵D_l′分成2×2的小矩阵，并把每一个小矩阵变成一个列向量，把所有列向量依次横向合并成为一个变化矩阵，记为D_l″：3b) Divide the changed logarithmic ratio matrix D _l ′ into 2×2 small matrices, and turn each small matrix into a column vector, and merge all the column vectors horizontally in turn to form a change matrix, denoted as D _l ″ :

3c)对变化矩阵D_l″进行低秩稀疏分解：3c) Perform low-rank sparse decomposition on the change matrix D _l ″:

（3c1)令初始低秩部分L₀＝D_l″，初始稀疏部分S₀＝0，秩阈值r＝2，t＝0；（3c2)验证t次分解后的低秩部分L_t和t次分解后的稀疏部分S_t是否满足收敛条件其中ε＝10^-3,如果L_t和S_t满足收敛条件，则L_t和S_t就是所要求的低秩部分L和稀疏部分S，否则继续下面步骤：(3c1) Make the initial low-rank part L ₀ =D _l ″, the initial sparse part S ₀ =0, the rank threshold r=2, t=0; (3c2) verify the low-rank part L _t and t times of decomposition after t times Whether the decomposed sparse part S _t satisfies the convergence condition Where ε=10 ^-3 , if L _t and S _t meet the convergence conditions, then L _t and S _t are the required low-rank part L and sparse part S, otherwise proceed to the following steps:

（3c3)令t＝t+1：(3c3) Let t=t+1:

（3c4)随机生成一个与初始稀疏部分S₀维数匹配的随机高斯矩阵，记为A₁：(3c4) Randomly generate a random Gaussian matrix matching the dimension of the initial sparse part _S0 , denoted as _A1 :

（3c5)根据随机高斯矩阵A₁和双边投影按照下式计算基础高斯矩阵A₂ ^t：(3c5) According to random Gaussian matrix A ₁ and bilateral projection The basic Gaussian matrix A ₂ ^t is calculated according to the following formula:

${A A}_{22}^{t t} = = {\overset{~ ~}{L L}}^{t t} {A A}_{11}$

其中， ${\tilde{L}}^{t} = {[(L_{0} - S_{t - 1}) {(L_{0} - S_{t - 1})}^{T}]}^{q} (L_{0} - S_{t - 1}), q = 3 :$ in, ${\tilde{L}}^{t} = {[(L_{0} - S_{t - 1}) {(L_{0} - S_{t - 1})}^{T}]}^{q} (L_{0} - S_{t - 1}), q = 3 :$

（3c6)根据基础高斯矩阵A₂ ^t，按如下公式求取第二随机观测样本Y₂ ^t和第一随机观测样本 (3c6) According to the basic Gaussian matrix A ₂ ^t , the second random observation sample Y ₂ t and the first random observation sample Y 2 ^t are obtained according to the following formula

${Y Y}_{22}^{t t} = = {(({\overset{~ ~}{L L}}^{t t}))}^{T T} {A A}_{22}^{t t},,$

${Y Y}_{11}^{t t} = = {\overset{~ ~}{L L}}^{t t} {Y Y}_{22}^{t t},,$

其中是的转置：in yes The transpose of:

（3c7)如果 $rank ({({A_{2}}^{t})}^{T} Y_{1}^{t}) < r,$ 则令 $r = rank ({({A_{2}}^{t})}^{T} Y_{1}^{t}),$ 返回到步骤（3c2)，否则继续下面步骤：(3c7) if $rank ({({A_{2}}^{t})}^{T} Y_{1}^{t}) < r,$ order $r = rank ({({A_{2}}^{t})}^{T} Y_{1}^{t}),$ Return to step (3c2), otherwise continue the following steps:

其中(A₂ ^t)^T是A₂ ^t的转置，是的秩：where (A ₂ ^t ) ^T is the transpose of A ₂ ^t , yes rank of:

（3c8)按如下公式对第一随机观测样本和第二随机观测样本Y₂ ^t进行QR分解：(3c8) according to the following formula for the first random observation sample and the second random observation sample Y ₂ ^t for QR decomposition:

${Y Y}_{11}^{t t} = = {Q Q}_{11}^{t t} {R R}_{11}^{t t}$

Y₂ ^t＝Q₂ ^tR₂ ^t，Y ₂ ^t = Q ₂ ^t R ₂ ^t ,

其中Q₁ ^t和R₁ ^t分别是对第一随机观测样本进行QR分解的Q部分和R部分，Q₂ ^t和R₂ ^t分别是对第二随机观测样本Y₂ ^t进行QR分解后的Q部分和R部分：where Q ₁ ^t and R ₁ ^t are the first random observation sample The Q part and R part of QR decomposition, Q ₂ ^t and R ₂ ^t are respectively the Q part and R part after QR decomposition of the second random observation sample Y ₂ ^t :

（3c9)按照下面公式求取t次迭代后的低秩部分L_t：(3c9) Calculate the low-rank part L _t after t iterations according to the following formula:

${L L}_{t t} = = {Q Q}_{11}^{t t} {{{R R}_{11}^{t t} {[[{(({A A}_{22}^{t t}))}^{T T} {Y Y}_{11}^{t t}]]}^{- - 11} {(({R R}_{22}^{t t}))}^{T T}}}^{\frac{11}{22 q q + + 11}} {(({Q Q}_{22}^{t t}))}^{T T},,$

其中是的转置，是的转置：in yes the transposition of yes The transpose of:

（3c10)按照下面公式把初始低秩部分L₀与t次分解后的低秩部分L_t的差值(L₀-L_t)投影到Ω方向上，得到t次迭代后的稀疏部分S_t：(3c10) According to the following formula, project the difference (L ₀ -L _t ) between the initial low-rank part L ₀ and the low-rank part L _t after t times of decomposition to the Ω direction, and obtain the sparse part S _t after t times of iterations :

S_t=P_Ω(L₀-L_t)，S _t = P _Ω (L ₀ -L _t ),

其中Ω为|L₀-L_t|的前k″个最大的非零系数，k″＝9×10³，P是投影算子：Where Ω is the first k″ largest non-zero coefficients of |L ₀ -L _t |, k″=9×10 ³ , and P is the projection operator:

（3c11)返回步骤（3c2)。(3c11) returns to step (3c2).

步骤4,把步骤3得到的稀疏部分S按列变换成与初始稀疏部分S₀维数一致的稀疏矩阵,记为S′。Step 4, transform the sparse part S obtained in step 3 into a sparse matrix with the same dimension as the initial sparse part S ₀ by column, denoted as S'.

步骤5,使用K均值算法把步骤（4）得到的稀疏矩阵S′聚成2类，得到的结果即为所求的变化检测图。Step 5, use the K-means algorithm to cluster the sparse matrix S′ obtained in step (4) into two categories, and the obtained result is the requested change detection map.

本发明的效果可以通过以下仿真实验进一步说明：Effect of the present invention can be further illustrated by following simulation experiments:

1.仿真条件：1. Simulation conditions:

在CPU为Intel Core(TM)2Duo、主频2.33GHz，内存为2G的WINDOWS XP系统上，用MATLAB 7.0.1软件分别对图2中Bern地区的两幅输入SAR图像和图4中Ottawa地区的两幅输入SAR图像进行仿真。On a WINDOWS XP system with Intel Core(TM) 2Duo CPU, main frequency 2.33GHz, and 2G memory, use MATLAB 7.0.1 software to input the two SAR images of the Bern area in Figure 2 and the Ottawa area in Figure 4 respectively. Two input SAR images are used for simulation.

仿真内容：Simulation content:

（1）对图2所述的Bern地区的两幅输入SAR图像进行仿真(1) Simulate the two input SAR images of the Bern area described in Figure 2

（1a）把图2（a）Bern地区变化前的SAR图像作为第一幅输入SAR图像，图2（b）Bern地区变化后的SAR图像作为第二幅输入SAR图像，将秩阈值r调整为1，用本发明的方法进行仿真，得到仿真结果如图6所示。(1a) Take the SAR image before the change of the Bern area in Figure 2(a) as the first input SAR image, and the SAR image after the change in the Bern area in Figure 2(b) as the second input SAR image, and adjust the rank threshold r to 1, carry out simulation with the method of the present invention, obtain the simulation result as shown in Figure 6.

（1b）用图3中Bern地区的参考图对（1a）中得到的仿真结果图进行验证，可知本发明仿真结果的虚警数、漏检数和总错误数如表1所示。(1b) Verify the simulation result diagram obtained in (1a) with the reference diagram of the Bern area in Fig. 3, it can be seen that the number of false alarms, missed detections and total errors of the simulation results of the present invention are shown in Table 1.

（1c）用现有几种方法：基于小波融合的方法Wavelet Fusion、基于对数比值的方法L-N和基于主分量分析的方法PCA，分别对图2Bern地区的两幅输入SAR图像进行仿真，用图3中Bern地区的参考图对仿真结果进行验证，得到这几种方法各自的仿真结果的虚警数、漏检数和总错误数如表1所示。(1c) Using several existing methods: Wavelet Fusion based on wavelet fusion, L-N based on logarithmic ratio and PCA based on principal component analysis, respectively simulate the two input SAR images in the Bern area in Figure 2. The reference map of the Bern area in 3 is used to verify the simulation results, and the number of false alarms, missed detections and total errors of the simulation results of these methods are shown in Table 1.

表1：不同方法在Bern地区图像的实验结果Table 1: Experimental results of different methods on Bern region images

方法 method 虚警数 Number of false alarms 漏检数 Number of missing detections 总错误数 total errors L-N L-N 88 88 226 226 314 314 PCA PCA 247 247 123 123 370 370 Wavelet Fusion Wavelet Fusion 503 503 77 77 580 580 本发明 this invention 100 100 179 179 279 279

（2）对图4所述的Ottawa地区的两幅输入SAR图像进行仿真(2) Simulate the two input SAR images of the Ottawa area described in Figure 4

（2a）把图4（a）Ottawa地区变化前的SAR图像作为第一幅输入SAR图像，图4（b）Ottawa地区变化后的SAR图像作为第二幅输入SAR图像，将最大的非零系数的个数k″调整为9×10⁴，用本发明的方法进行仿真，得到仿真结果如图7所示。(2a) Take the SAR image of the Ottawa area in Figure 4(a) before the change as the first input SAR image, and the SAR image in Figure 4(b) of the Ottawa area after the change as the second input SAR image, and set the largest non-zero coefficient The number k" of is adjusted to 9×10 ⁴ , and the method of the present invention is used for simulation, and the simulation result is shown in FIG. 7 .

（2b）用图5中Ottawa地区的参考图对（2a）中得到的仿真结果图进行验证，可知用本发明仿真结果的虚警数、漏检数和总错误数如表2所示。(2b) Verify the simulation result diagram obtained in (2a) with the reference diagram of the Ottawa area in Fig. 5. It can be seen that the number of false alarms, missed detections and total errors of the simulation results of the present invention are shown in Table 2.

（2c）用现有几种方法：基于小波融合的方法Wavelet Fusion、基于对数比值的方法L-N和基于主分量分析的方法PCA分别对图4Ottawa地区的两幅输入SAR图像进行仿真，用图5中Ottawa地区的参考图对仿真结果进行验证，得到这几种方法各自的仿真结果的虚警数、漏检数和总错误数如表2所示。(2c) Use several existing methods: Wavelet Fusion based on wavelet fusion, L-N based on logarithmic ratio and PCA based on principal component analysis to simulate the two input SAR images in the Ottawa area in Figure 4, and use Figure 5 The simulation results are verified by the reference map of the Central Ottawa area, and the number of false alarms, missed detections and total errors of the simulation results of these methods are shown in Table 2.

表2：不同方法在Ottawa地区图像的实验结果Table 2: Experimental results of different methods on images in the Ottawa region

方法 method 虚警数 Number of false alarms 漏检数 Number of missing detections 总错误数 total errors L-N L-N 1674 1674 583 583 2257 2257 PCA PCA 955 955 1515 1515 2470 2470 Wavelet Fusion Wavelet Fusion 896 896 1073 1073 1969 1969 本发明 this invention 334 334 1204 1204 1538 1538

2.仿真实验结果分析：2. Analysis of simulation experiment results:

从图6可以得出，本发明在Bern地区图像上得到的仿真结果图具有较好的边缘，斑点噪声非常少。It can be drawn from Fig. 6 that the simulation result map obtained by the present invention on the Bern area image has better edges and very little speckle noise.

从图7可以得出，本发明在Ottawa地区图像上得到的仿真结果图杂点得到了较为明显的消除。It can be drawn from FIG. 7 that the noise points in the simulation result map obtained by the present invention on the image of the Ottawa area have been more obviously eliminated.

从表1可以得出，本发明在Bern地区图像上仿真得到结果图的总错误数比Wavelet Fusion少了301，比L-N少了35，比PCA少了91。本发明得到的总错误数最少，由此可见本发明对Bern地区图像的正确检测率最高。It can be drawn from Table 1 that the total error number of the result map obtained by simulation on the Bern area image of the present invention is 301 less than Wavelet Fusion, 35 less than L-N, and 91 less than PCA. The total number of errors obtained by the present invention is the least, so it can be seen that the present invention has the highest correct detection rate for images in the Bern area.

从表2可以得出，本发明在Ottawa地区图像上仿真得到结果图的总错误数比Wavelet Fusion少了431，比L-N少了719，比PCA少了932。本发明得到的总错误数最少，由此可见本发明对Ottawa地区图像的正确检测率最高。It can be drawn from Table 2 that the total error number of the result map obtained by the present invention is 431 less than Wavelet Fusion, 719 less than L-N, and 932 less than PCA on the Ottawa region image. The total number of errors obtained by the present invention is the least, so it can be seen that the correct detection rate of the present invention is the highest for images in the Ottawa area.

综上，本发明在PPB算法的基础上用低秩矩阵分解方法对SAR图像进行变化检测，得到较高的正确检测率，与现有的方法相比检测精度有较明显的提高。To sum up, the present invention uses a low-rank matrix decomposition method to detect changes in SAR images on the basis of the PPB algorithm, and obtains a higher correct detection rate. Compared with the existing methods, the detection accuracy is significantly improved.

Claims

1. A SAR image change detection method based on low-rank matrix decomposition, comprising the steps:

(1) For the two input SAR images of the same area at different times and with the same size, use the probabilistic block-based denoising algorithm PPB to reduce speckle on the two images respectively, and the image matrices after speckle reduction are denoted as X ₁ and _X2 ;

(1a) Take any pixel point c _s in the input image as the center, select a neighborhood of N×N size as the search area of the pixel point;

(1b) Taking the pixel point c _s as the center, take a block of M×M size, and the gray values of all pixels in the block form a matrix v _s ;

(1c) Take every pixel f _t except the central pixel c _s in the search area as the center, take a block of M×M size, and the gray values of all pixels in the block form a matrix v _t ;

(1d) Calculate the weight w _s ,t i-1 from f _t to c _s after i- ¹ times of denoising according to the following weight formula:

If the input SAR image is an intensity image, use the weight formula:

{w w}_{s the s,, t t}^{i i - - 11} = = exp exp [[- - \underset{k k}{Σ Σ} ((\frac{11}{h h} [[((22 z z - - 11)) log log (({v v}_{s the s,, k k} + + {v v}_{t t,, k k})) - - ((z z - - 11)) log log (({v v}_{s the s,, k k} {v v}_{t t,, k k}))]] + + \frac{z z}{T T} \frac{| | {\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} - - {\overset{^^}{u u}}_{t t,, k k}^{i i - - 11} | |}{{\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} {\overset{^^}{u u}}_{t t,, k k}^{- - 11}}))]],,

If the input SAR image is a magnitude image, the weighting formula is used:

{w w}_{s the s,, t t}^{i i - - 11} = = exp exp [[- - \underset{k k}{Σ Σ} ((\frac{11}{\overset{^^}{h h}} log log ((\frac{{v v}_{s the s,, k k}}{{v v}_{t t,, k k}} + + \frac{{v v}_{t t,, k k}}{{v v}_{s the s,, k k}})) + + \frac{z z}{T T} \frac{| | {\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} - - {\overset{^^}{u u}}_{t t,, k k}^{i i - - 11} | |}{{\overset{^^}{u u}}_{s the s,, k k}^{i i - - 11} {\overset{^^}{u u}}_{t t,, k k}^{- - 11}}))]]

in v _{s, k} and v _{t, k} respectively represent the value of the kth pixel of v _s and v _t , and respectively and The value of the kth pixel of , exp is an exponential function with natural logarithm e as the base, h represents a smoothing parameter, z is an equivalent visual number, and T is a local adaptive parameter;

(1e) Calculate the value c _s ⁱ of the pixel point c _s after i times of denoising according to the following formula:

{c c}_{s the s}^{i i} = = \frac{{Σ Σ}_{t t} {w w}_{s the s,, t t}^{i i - - 11} {f f}_{t t}}{{Σ Σ}_{t t} {w w}_{s the s,, t t}^{i i - - 11}};;

(1f) For each pixel, repeat the above steps (1a)~(1e) to get the i-th denoised image

(1g) Let i=i+1, the image after denoising As a new input image, repeat the denoising process described in steps (1a) to (1f) 15 times to obtain the final speckle-reduced image

(2) Calculate the logarithmic ratio matrix D _l of the two image matrices X ₁ and X ₂ after speckle reduction according to the following formula:

{D D.}_{11} = = | | log log \frac{{X x}_{22} + + 11}{{X x}_{11} + + 11} | | = = | | log log (({X x}_{22} + + 11)) - - log log (({X x}_{11} + + 11)) | |;;

(3) Find the low-rank part and the sparse part of the logarithmic ratio matrix _D1 :

3a) If the number of columns of the logarithmic ratio matrix _D1 is an odd number, then delete the last column of _D1 , if the number of rows of _D1 is an odd number, then delete the last row of _D1 , and record the changed _D1 as D _l ';

3b) Divide the matrix D _l ' into small matrices of 2×2, and turn each small matrix into a column vector, and merge all the column vectors horizontally in turn into a change matrix, which is denoted as D _l ";

3c) Let the input initial low-rank part L ₀ =D _l ″ in the GoDec algorithm, the initial sparse part S ₀ =0, use the GoDec algorithm to decompose the change matrix D ₁ ″ into a low-rank part L and a sparse part S;

(4) Transform the sparse part S obtained in step (3) into a sparse matrix with the same dimension as the initial sparse part _S0 in columns, denoted as S';

(5) Use the K-means algorithm to cluster the sparse matrix S' obtained in step (4) into two categories, and the obtained result is the requested change detection map.

2. a kind of SAR image change detection method based on low-rank matrix decomposition according to claim 1, wherein said step 3c) change matrix D ₁ " is decomposed into low-rank part L and sparse part S with GoDec algorithm, Follow these steps:

(3c1) let L ₀ =D _l ″, S ₀ =0, r=2, t=0;

(3c2) Verify whether L _t and S _t meet the convergence conditions where L _t is the low-rank part after t iterations, S _t is the sparse part after t iterations, ε=10 ^-3 , if L _t and S _t meet the convergence conditions, then L _t and S _t are the required Low-rank part L and sparse part S, otherwise perform steps (3c3)-(3c7):

(3c3) Let t=t+1;

(3c4) Randomly generate a random Gaussian matrix matching the dimension of the initial sparse part _S0 , denoted as _A1 ;

(3c5) Calculate matrix A ₂ ^t according to the following formula:

{A A}_{22}^{t t} = = {\overset{~ ~}{L L}}^{t t} {A A}_{11}

in, as an intermediate variable,

{\tilde{L}}^{t} = {[(L_{0} - S_{t - 1}) {(L_{0} - S_{t - 1})}^{T}]}^{q} (L_{0} - S_{t - 1}), q = 3;

(3c6) Calculate the second random observation sample Y ₂ ^t and the first random observation sample Y ₁ ^t according to the following formula:

{Y Y}_{22}^{T T} = = {(({\overset{~ ~}{L L}}^{t t}))}^{T T} {A A}_{22}^{t t}

{Y Y}_{11}^{t t} = = {\overset{~ ~}{L L}}^{t t} {Y Y}_{22}^{t t},,

Among them, T is the transpose operator;

(3c7) If the rank of (A ₂ ^t ) ^T Y ₁ ^t is less than r, then set r=rank((A ₂ ^t ) ^T Y ₁ ^t ), return to step (3c2), where rank((A ₂ ^t ) ^T Y ₁ ^t ) is the rank of (A ₂ ^t ) ^T Y ₁ ^t , otherwise proceed to the following steps:

(3c8) Perform QR decomposition on the first random observation sample Y ₁ ^t and the second random observation sample Y ₂ ^t

Y ₁ ^t = Q ₁ ^t R ₁ ^t

Y ₂ ^t = Q ₂ ^t R ₂ ^t

Among them, Q ₁ ^t and R ₁ ^t are the Q part and R part after performing QR decomposition on the first random observation sample Y ₁ ^t respectively, and Q ₂ ^t and R ₂ ^t are the QR decomposition on the second random observation sample Y ₂ ^t respectively. The decomposed Q part and R part;

(3c9) Calculate the low-rank part L _t after t iterations according to the following formula:

{L L}_{t t} = = {Q Q}_{11}^{t t} {{{R R}_{11}^{t t} {[[{(({A A}_{22}^{t t}))}^{T T} {Y Y}_{11}^{t t}]]}^{- - 11} {(({R R}_{22}^{t t}))}^{T T}}}^{\frac{11}{22 q q + + 11}} {(({Q Q}_{22}^{t t}))}^{T T},,

(3c10) Project (L ₀ -L _t ) onto the Ω direction according to the following formula to obtain the sparse part S _t after t iterations:

S _t = P _Ω (L ₀ -L _t ),

Wherein, Ω is the first k" largest non-zero coefficients of |L ₀ -L _t |, k"=9×10 ³ , and P is a projection operator;

(3c11) Return to step (3c2).