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

Abstract

The invention discloses an SAR (synthetic aperture radar) image change detection method based on low-rank matrix factorization, and mainly solves the problem that the existing method cannot detect an SAR image change region accurately. The implementation steps of the method include: firstly, performing speckle reduction pretreatment to two SAR images to be detected to obtain smooth SAR images; secondly, constructing logarithm ratio of the two images after speckle reduction; thirdly, performing low-rank sparse factorization of the logarithm ration to obtain a low-rank part and a sparse part of the logarithm ratio; fourthly, transforming the sparse part into a sparse matrix by column; and fifthly, performing clustering to the obtained sparse matrix by the K-means algorithm to obtain the final change detection results. The method has the advantage of accurate change detection region, and can be used in the fields of public security and video monitoring.

Description

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

Technical field

The invention belongs to radar image processing technology field, particularly SAR image change detection method, can be used for solving the not high problem of change detected precision in SAR Image Change Detection.

Background technology

Image change detection method is a kind of important technology of analyzing and understanding remote sensing images of many times, has caused in recent years research widely.This comes from change detecting method application background widely, for example agricultural investigation, forest monitoring, Natural calamity monitoring, city mutation analysis, the aspects such as battle damage assessment.

Image Change Detection be to obtaining from areal different time many times remote Sensing Image Analysis a kind of method.What it stressed is the variation of atural object in identification two width remote sensing images.Existing change detecting method mainly can be divided into two large classes: have measure of supervision and without measure of supervision.

What is called has measure of supervision, is based on supervised classification, need to obtain the training field in image change region, thereby changes detection; And without supervision law without any need for extra information, the Data Detection of phase when directly different to two.Although there is supervision law accurately to determine that region of variation compares without supervision law and have obvious advantage, owing to obtaining, the real information on ground is difficult especially, so non-supervision variation detection method is conventional change detecting method.

Without supervision, change detection method, mainly comprise Principal Component Analysis Method, Wavelet Fusion method.Principal Component Analysis Method is by logarithm ratioing technigue structural differences figure, then on the basis of disparity map, adopts principal component analysis (PCA) dimensionality reduction to extract disparity map feature, and then by k-average to these feature clusterings.Wavelet Fusion method is to treat detected image to carry out respectively average ratio and logarithm ratio operation, with wavelet transform, extracts respectively low-and high-frequency information, and it is carried out to sparse fusion.These two kinds of shortcomings without measure of supervision are that accuracy of detection is lower.

Summary of the invention

The object of the invention is to overcome the deficiencies in the prior art, propose a kind of new SAR image change detection method based on low-rank matrix decomposition, to obtain comparatively accurate difference in change component, improve the accuracy of detection of SAR image change.

Realizing the object of the invention technical scheme is:

One. know-why

In compressed sensing field, if a matrix has unique " low-rank and sparse " structure, this matrix just can be recovered accurately under suitable condition.In the application of " low-rank and sparse " structure of matrix, a principal component analysis RPCA that outstanding work is robust, low-rank after " low-rank and sparse " of image array being decomposed with RPCA is partly equivalent to the principal ingredient of image array, and sparse part is equivalent to noise or other content.Another representative work is the matrix decomposition of stochastic approximation, and it has proved that the Main Ingredients and Appearance of image array can approach by accidental projection for a given image array.At computer vision field, these algorithms have very large application background, such as, background modeling, the illumination of facial image or shadow removal, image rectification etc.

The present invention, be inspired in above-mentioned thought, " low-rank and sparse " of matrix decomposed in the variation detection that is applied to SAR image, what low-rank wherein partly represented is unchanged part in image, what sparse part represented is changing unit, and this changes the needed result of Detection task just.

Two. implementation

The SAR image change detection method that the present invention is based on low-rank matrix decomposition, comprises the steps:

(1) two width the same area different times to input, equal-sized SAR image, uses the denoise algorithm PPB based on probability piece respectively spot to be fallen in this two width image, and the image array falling after spot is designated as respectively X ₁and X ₂:

(2) two width are fallen to for the image array X after spot ₁and X ₂by following formula, ask for the logarithm ratio matrix D of these two image arrays _l:

$D_l=\left|\log \frac{X_2+1}{X_1+1}\right|=|\log (X_2+1)-\log (X_1+1)|;$

(3) ask for logarithm ratio matrix D _llow-rank part and sparse part:

If 3a) logarithm ratio matrix D _lcolumns be odd number, leave out D _llast row, if D _lline number be odd number, leave out D _llast column, the D after changing _lbe designated as D _l':

3b) matrix D _l' be divided into 2 * 2 minor matrix, and each minor matrix is become to a column vector, all column vectors are laterally merged and become a transformation matrices successively, be designated as D _l":

3c) make the initial low-rank part of the input L in GoDec algorithm ₀=D _l", initial sparse part S ₀=0, use GoDec algorithm transformation matrices D _l" resolve into low-rank part L and sparse part S:

(4) sparse part S step (3) being obtained is by being column-transformed into and initial sparse part S ₀the sparse matrix that dimension is consistent, is designated as S ':

(5) the sparse matrix S ' that uses K mean algorithm that step (4) is obtained is polymerized to 2 classes, and the result obtaining is required variation detection figure.

The present invention has used the denoise algorithm PPB based on probability piece to fall spot to two width input pictures, with two, fall spot result and ask for logarithm ratio matrix, then with GoDec algorithm, logarithm ratio matrix is carried out to low-rank and Its Sparse Decomposition, obtained more accurate difference in change component, finally by K mean algorithm, difference in change component is carried out to cluster, the variation accuracy of detection of SAR image is improved than existing accuracy of detection.

Accompanying drawing explanation

Fig. 1 is realization flow figure of the present invention;

Fig. 2 is two width input SAR images used during to Bern area image simulation with the present invention;

Fig. 3 is the reference diagram in Bern area;

Fig. 4 is two width input SAR images used during to Ottawa area image simulation with the present invention;

Fig. 5 is the reference diagram in Ottawa area;

Fig. 6 is the simulation result figure to Fig. 2;

Fig. 7 is the simulation result figure to Fig. 4.

Embodiment

With reference to Fig. 1, specific embodiment of the invention step is as follows:

Step 1, input two width the same areas, different time, equal-sized SAR image, and fall spot and process:

(1a) with arbitrary pixel c in the first width input SAR image _scentered by, choose the neighborhood of N * N size as the Search Area of this pixel, wherein N=21:

(1b) with pixel c _scentered by, get the piece of M * M size, wherein M=7, forms matrix v with the gray-scale value of all pixels in piece _s:

(1c) to remove central pixel point c in Search Area _seach outer pixel f _tcentered by, get the piece of M * M size, in piece, the gray-scale value of all pixels forms matrix v _t:

(1d) according to the f after i-1 denoising of weight formula calculating below _tto c _sweight w _{s, t} ^i-1:

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

${w_{s,t}}^{i-1}=\exp [-\underset{k}{\mathrm{\Σ}}\left(\frac{1}{h}\right[(2L-1)\log (v_{s,k}+v_{t,k})-(L-1)\log \left(v_{s,k}{v}_{t,k}\right)]+\frac{L}{T}\frac{|{{\hat{u}}_{s,k}}^{i-1}-{{\hat{u}}_{t,k}}^{i-1}}{{{\hat{u}}_{s,k}}^{i-1}{{\hat{u}}_{t,k}}^{i-1}})],$

If the first width input SAR image is magnitude image, use weight formula:

${w_{s,t}}^{i-1}=\exp [-\underset{k}{\mathrm{\Σ}}(\frac{1}{\tilde{h}}\log (\frac{v_{s,k}}{v_{t,k}}+\frac{v_{t,k}}{v_{s,k}})+\frac{L}{T}\frac{|{{\hat{u}}_{s,k}}^{i-1}-{{\hat{u}}_{t,k}}^{i-1}}{{{\hat{u}}_{s,k}}^{i-1}{{\hat{u}}_{t,k}}^{i-1}})]$

Wherein ${\hat{u}}_s^{i-1}=\frac{{\mathrm{\Σ}}_t{w}_{s,t}{\left(v_t^{i-1}\right)}^2}{{\mathrm{\Σ}}_t{w}_{s,t}},$ ${\hat{u}}_t^{i-1}=V_t^{i-1},$ $\tilde{h}=\frac{h}{2L-1},$ V _{s, k}and v _{t, k}represent respectively v _sand v _tthe value of k pixel, with be respectively with the value of k pixel, exp be take the exponential function that natural logarithm e is the truth of a matter, h=11, L=3, T=0.2:

(1e) according to formula calculating pixel point c below _svalue after i denoising

$c_s^i=\frac{{\mathrm{\Σ}}_t{w_{s,t}}^{i-1}{f}_t}{{\mathrm{\Σ}}_t{w_{s,t}}^{i-1}}$

(1f) for each pixel, repetition above-mentioned steps (1a) ~ (1e), obtain the image after the i time denoising

(1g) make i=i+1, with image after denoising as new the first width input SAR image, repeat the described denoising process of step (1a) ~ (1f) 15 times, thereby obtain the image array X after final decline spot ₁:

(1h) by step (1a) ~ (1g) described spot process of falling, spot processing is fallen in the second width input SAR image, obtain falling the image array X after spot ₂.

Step 2, falls the image array X after spot by two width ₁and X ₂, by following formula, ask for the logarithm ratio matrix D of these two image arrays _l:

$D_l=\left|\log \frac{X_2+1}{X_1+1}\right|=|\log (X_2+1)-\log (X_1+1)|.$

Step 3, asks for logarithm ratio matrix D _llow-rank part and sparse part:

If 3a) logarithm ratio matrix D _lcolumns be odd number, leave out D _llast row, if D _lline number be odd number, leave out D _llast column, the logarithm ratio matrix after changing is designated as to D _l':

3b) by the logarithm ratio matrix D after changing _l' be divided into 2 * 2 minor matrix, and each minor matrix is become to a column vector, all column vectors are laterally merged and become a transformation matrices successively, be designated as D _l":

3c) to transformation matrices D _l" carry out low-rank Its Sparse Decomposition:

(3c1) make initial low-rank part L ₀=D _l", initial sparse part S ₀=0, order threshold value r=2, t=0; (3c2) the low-rank part L after t decomposition of checking _tsparse part S after decomposing with t time _twhether meet the condition of convergence ε=10 wherein ^-3if, L _tand S _tmeet the condition of convergence, L _tand S _tbe exactly desired low-rank part L and sparse part S, otherwise continue step below:

(3c3) make t=t+1:

(3c4) generate at random one and initial sparse part S ₀the random Gaussian matrix of dimension coupling, is designated as A ₁:

(3c5) according to random Gaussian matrix A ₁with bilateral projection calculate according to the following formula basic Gauss's matrix A ₂ ^t:

${A_2}^t={\tilde{L}}^t{A}_1$

Wherein, ${\tilde{L}}^t={\left[(L_0-S_{t-1}){(L_0-S_{t-1})}^T\right]}^q(L_0-S_{t-1}),q=3:$

(3c6) according to basic Gauss's matrix A ₂ ^t, by following formula, ask for the second random observation sample Y ₂ ^twith the first random observation sample

${Y_2}^t={\left({\tilde{L}}^t\right)}^T{A_2}^t,$

$Y_1^t={\tilde{L}}^t{Y}_2^t,$

Wherein be transposition:

If (3c7) $\mathrm{rank}\left({\left({A_2}^t\right)}^T{Y}_1^t\right)<r,$ Order $r=\mathrm{rank}\left({\left({A_2}^t\right)}^T{Y}_1^t\right),$ Turn back to step (3c2), otherwise continue step below:

(A wherein ₂ ^t) ^ta ₂ ^ttransposition, be order:

(3c8) press following formula to the first random observation sample with the second random observation sample Y ₂ ^tcarry out QR decomposition:

$Y_1^t=Q_1^t{R}_1^t$

Y ₂ ^t＝Q ₂ ^tR ₂ ^t，

Q wherein ₁ ^tand R ₁ ^trespectively to the first random observation sample carry out Q part and the R part of QR decomposition, Q ₂ ^tand R ₂ ^trespectively to the second random observation sample Y ₂ ^tcarry out Q part and R part after QR decomposition:

(3c9) according to formula below, ask for the low-rank part L after t iteration _t:

$L_t=Q_1^t{\left\{R_1^t{\left[{\left({A_2}^t\right)}^T{Y}_1^t\right]}^{-1}{\left(R_2^t\right)}^T\right\}}^{\frac{1}{2q+1}}{\left(Q_2^t\right)}^T,$

Wherein be transposition, be transposition:

(3c10) according to formula below initial low-rank part L ₀low-rank part L after decomposing with t time _tdifference (L ₀-L _t) project in Ω direction, obtain the sparse part S after t iteration _t:

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

Wherein Ω is | L ₀-L _t| front k " nonzero coefficient of individual maximum, k "=9 * 10 ³, P is projection operator:

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

Step 4, the sparse part S that step 3 is obtained is by being column-transformed into and initial sparse part S ₀the sparse matrix that dimension is consistent, is designated as S '.

Step 5, the sparse matrix S ' that uses K mean algorithm that step (4) is obtained is polymerized to 2 classes, and the result obtaining is required variation detection figure.

Effect of the present invention can further illustrate by following emulation experiment:

1. simulated conditions:

At CPU, be Intel Core (TM) 2Duo, dominant frequency 2.33GHz, inside save as in the WINDOWS XP system of 2G, with MATLAB 7.0.1 software, respectively the two width input SAR images in Ottawa area in the two width input SAR images in Bern area in Fig. 2 and Fig. 4 are carried out to emulation.

Emulation content:

(1) the two width input SAR images in the Bern area described in Fig. 2 are carried out to emulation

(1a) the SAR image before Fig. 2 (a) Bern area changes is as the first width input SAR image, SAR image after Fig. 2 (b) Bern area changes is as the second width input SAR image, order threshold value r is adjusted into 1, by method of the present invention, carries out emulation, obtain simulation result as shown in Figure 6.

(1b) with the reference diagram in Bern area in Fig. 3, the simulation result figure obtaining in (1a) is verified, the false-alarm number of known simulation result of the present invention, undetected number and total wrong number are as shown in table 1.

(1c) use existing several method: the method Wavelet Fusion based on Wavelet Fusion, the method L-N based on logarithm ratio and the method PCA based on principal component analysis, respectively the two width input SAR images in Fig. 2 Bern area are carried out to emulation, with the reference diagram in Bern area in Fig. 3, simulation result is verified, obtained the false-alarm number of this several method simulation result separately, undetected number and total wrong number as shown in table 1.

Table 1: distinct methods is in the experimental result of Bern plat picture

Method	False-alarm number	Undetected number	Total wrong number
				L-N	88	226	314
PCA	247	123	370
				Wavelet Fusion	503	77	580
The present invention	100	179	279

(2) the two width input SAR images in the Ottawa area described in Fig. 4 are carried out to emulation

(2a) the SAR image before Fig. 4 (a) Ottawa area changes is as the first width input SAR image, and the SAR image after Fig. 4 (b) Ottawa area changes, as the second width input SAR image, " is adjusted into 9 * 10 by the number k of maximum nonzero coefficient ⁴, by method of the present invention, carry out emulation, obtain simulation result as shown in Figure 7.

(2b) with the reference diagram in Ottawa area in Fig. 5, the simulation result figure obtaining in (2a) is verified, known false-alarm number, undetected number and total wrong number with simulation result of the present invention is as shown in table 2.

(2c) use existing several method: the method Wavelet Fusion based on Wavelet Fusion, the method L-N based on logarithm ratio and the method PCA based on principal component analysis carry out emulation to the two width input SAR images in Fig. 4 Ottawa area respectively, with the reference diagram in Ottawa area in Fig. 5, simulation result is verified, obtained the false-alarm number of this several method simulation result separately, undetected number and total wrong number as shown in table 2.

Table 2: distinct methods is in the experimental result of Ottawa plat picture

Method	False-alarm number	Undetected number	Total wrong number
				L-N	1674	583	2257
PCA	955	1515	2470
				Wavelet Fusion	896	1073	1969
The present invention	334	1204	1538

2. the simulation experiment result analysis:

As can be drawn from Figure 6, the simulation result figure that the present invention obtains on Bern plat picture has good edge, and speckle noise is considerably less.

As can be drawn from Figure 7, the assorted point of simulation result figure that the present invention obtains on Ottawa plat picture has obtained comparatively significantly eliminating.

As can be drawn from Table 1, the total wrong number that the present invention's emulation on Bern plat picture obtains result figure than Wavelet Fusion few 301, than L-N few 35, than PCA few 91.The total wrong number that the present invention obtains is minimum, and the present invention is the highest to the correct verification and measurement ratio of Bern plat picture as can be seen here.

As can be drawn from Table 2, the total wrong number that the present invention's emulation on Ottawa plat picture obtains result figure than Wavelet Fusion few 431, than L-N few 719, than PCA few 932.The total wrong number that the present invention obtains is minimum, and the present invention is the highest to the correct verification and measurement ratio of Ottawa plat picture as can be seen here.

To sum up, the present invention changes detection with low-rank matrix disassembling method to SAR image on the basis of PPB algorithm, obtains higher correct verification and measurement ratio, and accuracy of detection has obvious raising compared with the existing methods.

Claims (2)

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

(1) two width the same area different times to input, equal-sized SAR image, uses the denoise algorithm PPB based on probability piece respectively spot to be fallen in this two width image, and the image array falling after spot is designated as respectively X ₁and X ₂;

(1a) with arbitrary pixel c in input picture _scentered by, choose the neighborhood of N * N size as the Search Area of this pixel;

(1b) with pixel c _scentered by, get the piece of M * M size, in piece, the gray-scale value of all pixels forms matrix v _s;

(1c) to remove central pixel point c in Search Area _seach outer pixel f _tcentered by, get the piece of M * M size, in piece, the gray-scale value of all pixels forms matrix v _t;

(1d) according to the f after i-1 denoising of weight formula calculating below _tto c _sweight w _s,t ^i-1:

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

${w_{s,t}}^{i-1}=\exp [-\underset{k}{\mathrm{\&Sigma;}}\left(\frac{1}{h}\right[(2z-1)\log (v_{s,k}+v_{t,k})-(z-1)\log \left(v_{s,k}{v}_{t,k}\right)]+\frac{z}{T}\frac{|{{\hat{u}}_{s,k}}^{i-1}-{{\hat{u}}_{t,k}}^{i-1}|}{{{\hat{u}}_{s,k}}^{i-1}{{\hat{u}}_{t,k}}^{-1}})],$

If input SAR image is magnitude image, use weight formula:

${w_{s,t}}^{i-1}=\exp [-\underset{k}{\mathrm{\&Sigma;}}(\frac{1}{\hat{h}}\log (\frac{v_{s,k}}{v_{t,k}}+\frac{v_{t,k}}{v_{s,k}})+\frac{z}{T}\frac{|{{\hat{u}}_{s,k}}^{i-1}-{{\hat{u}}_{t,k}}^{i-1}|}{{{\hat{u}}_{s,k}}^{i-1}{{\hat{u}}_{t,k}}^{-1}})]$

Wherein v _s,kand v _t,krepresent respectively v _sand v _tthe value of k pixel, with be respectively with the value of k pixel, exp be take the exponential function that natural logarithm e is the truth of a matter, h represents smoothing parameter, z is equivalent number, T is local auto-adaptive parameter;

(1e) according to formula calculating pixel point c below _svalue c after i denoising _s ⁱ:

${c_s}^i=\frac{{\mathrm{\&Sigma;}}_t{w_{s,t}}^{i-1}{f}_t}{{\mathrm{\&Sigma;}}_t{w_{s,t}}^{i-1}};$

(1f) for each pixel, repetition above-mentioned steps (1a)～(1e), obtain the image after the i time denoising

(1g) make i=i+1, with image after denoising as new input picture, repeat the described denoising process of step (1a)～(1f) 15 times, thereby obtain the image after final decline spot

(2) two width are fallen to for the image array X after spot ₁and X ₂by following formula, ask for the logarithm ratio matrix D of these two image arrays _l:

$D_1=\left|\log \frac{X_2+1}{X_1+1}\right|=|\log (X_2+1)-\log (X_1+1)|;$

(3) ask for logarithm ratio matrix D _llow-rank part and sparse part:

If 3a) logarithm ratio matrix D _lcolumns be odd number, leave out D _llast row, if D _lline number be odd number, leave out D _llast column, the D after changing _lbe designated as D _l';

3b) matrix D _l' be divided into 2 * 2 minor matrix, and each minor matrix is become to a column vector, all column vectors are laterally merged and become a transformation matrices successively, be designated as D _l";

3c) make the initial low-rank part of the input L in GoDec algorithm ₀=D _l", initial sparse part S ₀=0, use GoDec algorithm transformation matrices D _l" resolve into low-rank part L and sparse part S;

(4) sparse part S step (3) being obtained is by being column-transformed into and initial sparse part S ₀the sparse matrix that dimension is consistent, is designated as S ';

(5) the sparse matrix S ' that uses K mean algorithm that step (4) is obtained is polymerized to 2 classes, and the result obtaining is required variation detection figure.

2. a kind of SAR image change detection method based on low-rank matrix decomposition according to claim 1, wherein said step 3c) with GoDec algorithm transformation matrices D _l" resolve into low-rank part L and sparse part S, according to the following steps:

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

(3c2) checking L _tand S _twhether meet the condition of convergence l wherein _tthe low-rank part after t iteration, S _tthe sparse part after t iteration, ε=10 ^-3if, L _tand S _tmeet the condition of convergence, L _tand S _tbe exactly desired low-rank part L and sparse part S, otherwise execution step (3c3)-(3c7):

(3c3) make t=t+1;

(3c4) generate at random one and initial sparse part S ₀the random Gaussian matrix of dimension coupling, is designated as A ₁;

(3c5) according to formula compute matrix A below ₂ ^t:

${A_2}^t={\tilde{L}}^t{A}_1$

Wherein, for intermediate variable, ${\tilde{L}}^t={\left[(L_0-S_{t-1}){(L_0-S_{t-1})}^T\right]}^q(L_0-S_{t-1}),q=3;$

(3c6) by following formula, ask for the second random observation sample Y ₂ ^twith the first random observation sample Y ₁ ^t:

${Y_2}^T={\left({\tilde{L}}^t\right)}^T{A_2}^t$

${Y_1}^t={\tilde{L}}^t{Y_2}^t,$

Wherein, T is transposed operator;

If (3c7) (A ₂ ^t) ^ty ₁ ^torder be less than r, make r=rank ((A ₂ ^t) ^ty ₁ ^t), turn back to step (3c2), wherein rank ((A ₂ ^t) ^ty ₁ ^t) be (A ₂ ^t) ^ty ₁ ^torder, otherwise continue step below:

(3c8) to the first random observation sample Y ₁ ^twith the second random observation sample Y ₂ ^tcarry out QR decomposition

Y ₁ ^t＝Q ₁ ^tR ₁ ^t

Y ₂ ^t＝Q ₂ ^tR ₂ ^t

Q wherein ₁ ^tand R ₁ ^trespectively to the first random observation sample Y ₁ ^tcarry out Q part and R part after QR decomposition, Q ₂ ^tand R ₂ ^trespectively to the second random observation sample Y ₂ ^tcarry out Q part and R part after QR decomposition;

(3c9) according to formula below, ask for the low-rank part L after t iteration _t:

$L_t={Q_1}^t{\left\{{R_1}^t{\left[{\left({A_2}^t\right)}^T{Y_1}^t\right]}^{-1}{\left(R_2^t\right)}^T\right\}}^{\frac{1}{2q+1}}{\left(Q_2^t\right)}^T,$

(3c10) according to formula below (L ₀-L _t) project in Ω direction, obtain the sparse part S after t iteration _t:

S _t＝P _Ω(L ₀-L _t)，

Wherein, Ω is | L ₀-L _t| front k " nonzero coefficient of individual maximum, k "=9 * 10 ³, P is projection operator;

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