CN103824302A - SAR (synthetic aperture radar) image change detecting method based on direction wave domain image fusion - Google Patents

Abstract

The invention discloses an SAR (synthetic aperture radar) image change detecting method based on direction wave domain image fusion, which overcomes the problems that difference image detail information cannot be kept well when difference images are fused in the prior art, and information of a change region is lost when difference images are processed by a binarization segmentation method. The SAR image change detecting method provided by the invention comprises the following implementation steps: (1) inputting images; (2) pretreating the images; (3) generating difference images; (4) constructing a combined difference image; (5) segmenting the combined difference image; and (6) outputting images. The SAR image change detecting method based on direction wave domain image fusion has the advantages of describing boundaries of a changed region and an unchanged region accurately, and being low in detection leakage rate, and can be applied to the natural disaster detection in the field of SAR image aim recognition.

Description

The SAR image change detection method merging based on direction wave area image

Technical field

The invention belongs to technical field of image processing, further relate to a kind of synthetic-aperture radar (Synthetic Aperture Radar, the SAR) image change detection method merging based on direction wave area image in diameter radar image change detection techniques field.The present invention can be used for disaster and detects in SAR image object identification field.

Background technology

Changing the research object detecting is multidate data, requires to obtain observed object at imaging data in the same time not, and changing the object detecting is by multidate data analysis, extracts the change information that object of observation does not occur in the same time.Synthetic-aperture radar (Synthetic Aperture Radar, SAR) imaging is as the important Observations Means over the ground of one, particularly its active microwave imaging principle, enable not to be subject to the impact of weather and illumination, round-the-clock, round-the-clock is observed over the ground, and in addition, the penetration capacity of microwave makes SAR can penetrate forest, have detectivity for the maneuvering target being hidden in wood land, therefore SAR image has become the important information source of Image Change Detection.

Xian Electronics Science and Technology University has proposed a kind of sar image change detection method based on neighborhood logarithm ratio and anisotropy diffusion in its patented claim " based on the sar image change detection method of neighborhood logarithm ratio and anisotropy diffusion " (number of patent application: CN201110005209, publication number: CN102096921A).The concrete steps of the method are: first, the disparity map of phasor while adopting neighborhood logarithm ratio method construct two width, then, carries out self-adaptation anisotropic diffusion filtering and maximum variance between clusters is cut apart to disparity map, obtains finally changing testing result.Although the method can keep the marginal information of region of variation preferably, but the deficiency still existing is, because binarization segmentation is described not too applicable to the pixel class of associating disparity map, so adopt binarization segmentation method can cause region of variation information dropout to disparity map, cause changing the loss detecting larger.

Xian Electronics Science and Technology University has proposed a kind of remote sensing image change detecting method based on image co-registration in its patented claim " remote sensing image based on image co-registration changes detection " (number of patent application: CN2012102340761, publication number: CN102750705A).The method concrete steps are: first, obtain associating disparity map by merge two width disparity map in wavelet field, then, use fuzzy Local C means Method to cut apart the disparity map after fusion, obtain changing testing result.Although the method is by improving the disparity map performance of remote sensing image, improve the precision of testing result, but the deficiency still existing is, because wavelet transformation is for the extraction poor-performing of the unusual information of line, while causing disparity map to merge, disparity map detailed information be can not retain preferably, region of variation in testing result and the Boundary Detection out of true of region of variation not made to change.

Summary of the invention

The object of the invention is to overcome above-mentioned the deficiencies in the prior art, proposed a kind of SAR image change detection method merging based on direction wave area image.The present invention, by improving disparity map performance and improving disparity map segmentation strategy, has improved the accuracy rate of SAR Image Change Detection result.

For achieving the above object, concrete steps of the present invention are as follows:

(1) input picture:

The two width SAR image S that optional areal, different time are taken ₁and S ₂as waiting to change detection SAR image.

(2) image pre-service:

Adopt probability piecemeal matched filtering method, respectively to SAR image S ₁and S ₂filtering, obtains image X after filtering ₁and X ₂.

(3) generate disparity map:

(3a) utilize ratio operator formula, to image X after filtering ₁and X ₂in pixel value operate, obtain disparity map I ₁pixel value;

(3b) utilize hybrid operator formula, to image X after filtering ₁and X ₂in pixel value operate, obtain disparity map I ₂pixel value.

(4) tectonic syntaxis disparity map:

(4a) respectively to disparity map I ₁and I ₂travel direction wave conversion, obtains disparity map I ₁and I ₂direction wave low frequency sub-band coefficient and high-frequency sub-band coefficient;

(4b) utilize following formula, to disparity map I ₁and I ₂direction wave low frequency sub-band coefficient merge, obtain associating disparity map I ₃direction wave low frequency sub-band coefficient:

$L_{\mathrm{ij}}^{I_3}=\frac{1}{2}(L_{\mathrm{ij}}^{I_1}+L_{\mathrm{ij}}^{I_2}),i\∈M,j\∈N$

Wherein,

represent associating disparity map I ₃direction wave low frequency sub-band in the coefficient of capable, the j of i row,

represent disparity map I ₁the coefficient of capable, the j of direction wave low frequency sub-band i row,

represent disparity map I ₂the coefficient of capable, the j of direction wave low frequency sub-band i row, M represents disparity map I ₁and I ₂line number, N represents disparity map I ₁and I ₂columns;

(4c) utilize following formula, to disparity map I ₁and I ₂direction wave high-frequency sub-band coefficient merge, obtain associating disparity map I ₃direction wave high-frequency sub-band coefficient:

$H_{\mathrm{ij}}^{I_3}=\max (H_{\mathrm{ij}}^{I_1},H_{\mathrm{ij}}^{I_2}),i\∈M,j\∈N$

Wherein, represent associating disparity map I ₃direction wave high-frequency sub-band in the coefficient of capable, the j of i row, max () represents to get maxima operation, represent disparity map I ₁high-frequency sub-band in the coefficient of capable, the j of i row,

represent disparity map I ₂direction wave low frequency sub-band in the coefficient of capable, the j of i row, M represents disparity map I ₁and I ₂line number, N represents disparity map I ₁and I ₂columns;

(4d) to associating disparity map I ₃direction wave low frequency coefficient and the inverse transformation of high frequency coefficient travel direction ripple, obtain associating disparity map I ₃.

(5) cut apart associating disparity map:

(5a) adopt fuzzy C-means clustering method, as follows, will combine disparity map I ₃in all pixels be divided into strong variation class pixel, weak variation class pixel and do not change class pixel:

(5a1) will combine disparity map I ₃in the gray-scale value of all pixels as the proper vector of each pixel, meanwhile, from associating disparity map I ₃in select arbitrarily the gray-scale value of 3 pixels as 3 initial cluster centers of 3 class pixels;

(5a2) according to the following formula, upgrade associating disparity map I ₃the degree of membership of the proper vector of middle pixel:

$u_{\mathrm{ik}}=\mathrm{\Σ}{\left(\frac{\left|\right|x_k-v_i\left|\right|}{\left|\right|x_k-v_j\left|\right|}\right)}^{-2},j=\mathrm{1,2},...,n$

Wherein, u _ikrepresent associating disparity map I ₃in the proper vector of k pixel be under the jurisdiction of the degree of membership of i class pixel, ∑ is sum operation, x _krepresent associating disparity map I ₃in the proper vector of k pixel, v _irepresent the cluster centre of i class pixel, v _jrepresent the cluster centre of j class pixel, k represents associating disparity map I ₃the sequence number of the proper vector of middle pixel, i represents current cluster classification, and j represents cluster classification arbitrarily, and n represents associating disparity map I ₃the number of the proper vector of middle pixel, || || represent to ask Euclidean distance operation;

(5a3) according to the following formula, upgrade cluster centre vector:

$v_i=\frac{\mathrm{\Σ}{\left(u_{\mathrm{ik}}\right)}^2{x}_k}{\mathrm{\Σ}{\left(u_{\mathrm{ik}}\right)}^2},k=\mathrm{1,2},...,n$

Wherein, v _irepresent associating disparity map I ₃the cluster centre vector of i class pixel, ∑ is sum operation, μ _ikrepresent associating disparity map I ₃in the proper vector of k pixel be under the jurisdiction of the degree of membership of i class pixel, x _krepresent associating disparity map I ₃in the proper vector of k pixel, i represents the classification of cluster, k represents associating disparity map I ₃the sequence number of the proper vector of middle pixel, n represents associating disparity map I ₃the number of the proper vector of middle pixel;

(5a4) first step to the three step iteration are carried out 20 times, obtained associating disparity map I ₃the proper vector of middle pixel be subordinate to matrix;

(5a5) in being subordinate to matrix, the peaked line number of every row degree of membership, as the classification of pixel, obtains changing by force class pixel, weak variation class pixel and does not change class pixel;

(5b) will change by force class pixel and weak variation class pixel and merge into variation class pixel;

(5c) using the not variation class pixel in step (5a) as changing testing result image I _rin not variation class pixel, using the variation class pixel in step (5b) as changing testing result image I _rin variation class pixel, obtain changing testing result image I _r.

(6) output image:

Exporting change testing result image I _r.

The inventive method has the following advantages compared with prior art:

First, because the present invention has adopted direction wave area image fusion method, two width disparity map are merged, overcome the shortcoming that can not retain preferably disparity map detailed information when disparity map merges in prior art, improve associating disparity map and extracted the performance of detailed information, made the present invention there is the region of variation of description and the more accurate advantage in region of variation border not.

Second, because the present invention adopts fuzzy C-mean algorithm method, disparity map pixel is divided into strong variation class pixel, weak variation class pixel and does not change class pixel, having overcome in prior art adopts binarization segmentation method can cause the shortcoming of region of variation information dropout to disparity map, retain more accurately the region of variation information in associating disparity map, made the present invention have advantages of that loss is lower.

Accompanying drawing explanation

Fig. 1 is process flow diagram of the present invention;

Fig. 2 is that the inventive method and prior art are to the regional SAR Image Change Detection of Berne (Bern) effect contrast figure;

Fig. 3 is that the inventive method and prior art are to Ottawa (Ottawa) regional SAR Image Change Detection effect contrast figure.

Embodiment

Below in conjunction with accompanying drawing, the present invention will be further described.

With reference to accompanying drawing 1, performing step of the present invention is as follows:

Step 1: input picture.

The two width SAR image S that optional areal, different time are taken ₁and S ₂as waiting to change detection SAR image.The synthetic-aperture radar SAR image using in the embodiment of the present invention is respectively as shown in Fig. 2 (a), Fig. 2 (b) and Fig. 3 (a), Fig. 3 (b).Wherein, Fig. 2 (a) is Berne (Bern) the area synthetic-aperture radar SAR image in April, 1999, size is 301 × 301, Fig. 2 (b) is Berne (Bern) the area synthetic-aperture radar SAR image in May, 1999, size is 301 × 301, Fig. 3 (a) is area, Ottawa (Ottawa) the synthetic-aperture radar SAR image in May, 1997, size is 290 × 350, Fig. 3 (b) is area, Ottawa (Ottawa) the synthetic-aperture radar SAR image in August, 1997, and size is 290 × 350.

Step 2: image pre-service.

Adopt probability piecemeal matched filtering method, respectively to SAR image S ₁and S ₂filtering, obtains image X after filtering ₁and X ₂.The concrete steps of probability piecemeal matched filtering method are as follows:

The first step, sets SAR image S ₁and S ₂in arbitrarily pixel neighborhood of a point window size be 7 × 7.

Second step, according to the following formula, calculates respectively SAR image S ₁and S ₂in any weights of pixel and interior other pixels of this pixel neighborhood of a point window:

$w(q,t)=\exp (-\underset{p=1}{\overset{49}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{q,p}}{A_{t,r}}+\frac{A_{t,r}}{A_{q,p}}))$

Wherein, w (q, t) represents SAR image S ₁and S ₂in any weights of pixel q and interior other pixels t of this pixel neighborhood of a point window, exp () represents the operation of fetching number, ∑ represents sum operation, h represents the smoothing parameter of control characteristic attenuation degree, log () represents to take the logarithm operation, and p represents any pixel q pixel sequence number in the neighborhood of 7 × 7 sizes, and r represents the sequence number of pixel t pixel in 7 × 7 large small neighbourhoods, and r=p, A _q,prepresent the gray-scale value of any pixel q pixel sequence number p corresponding position in neighborhood, A _t,rrepresent the gray-scale value of pixel t pixel sequence number r corresponding position in neighborhood.

The 3rd step, according to the following formula, calculates respectively SAR image S ₁and S ₂in any gray scale estimated value of pixel:

${\hat{z}}_1\left(q\right)=\frac{\underset{t\∈W_q}{\mathrm{\Σ}}w(q,t)A_t^2}{\underset{t\∈W_q}{\mathrm{\Σ}}w(q,t)}$

Wherein,

represent SAR image S ₁and S ₂in any gray scale estimated value of pixel q, ∑ represents sum operation, W _qrepresent SAR image S ₁and S ₂in any neighborhood window of 7 × 7 sizes of pixel q, w (q, t) represents SAR image S ₁and S ₂in any pixel q and its neighborhood window in the weights of other pixels t, A _trepresent the gray-scale value of pixel t.

The 4th step, carries out the first step, second step, the 3rd step iteration 25 times, obtains SAR image S ₁and S ₂filtering after image X ₁and X ₂.

Step 3: generate disparity map.

Utilize ratio operator formula, to image X after filtering ₁and X ₂in pixel value operate, obtain disparity map I ₁pixel value.Ratio operator formula is as follows:

$D=\frac{\underset{k=1}{\overset{49}{\mathrm{\Σ}}}\min (X_{1,k}^{N_i},X_{2,k}^{N_i})}{\underset{k=1}{\overset{49}{\mathrm{\Σ}}}\max (X_{1,k}^{N_i},X_{2,k}^{N_i})}$

Wherein, D represents disparity map I ₁pixel value, ∑ represents sum operation, min () represents to get minimum value operation, max () represents to get maxima operation, N _iimage X after expression filtering ₁and X ₂in centered by i pixel, the neighborhood of size as 7 × 7, i represents image X after filtering ₁and X ₂in pixel sequence number, k represents the pixel sequence number in neighborhood Ni, k=1,2 ..., 49,

image X after expression filtering ₁middle neighborhood N _ik pixel,

image X after expression filtering ₂middle neighborhood N _ik pixel.

Utilize hybrid operator formula, to image X after filtering ₁and X ₂in pixel value operate, obtain disparity map I ₂pixel value.Hybrid operator formula is as follows:

$F=\frac{1}{2}\{\log \left(\frac{1}{2}(\frac{X_{1,k}^{N_i}}{X_{2,k}^{N_i}}+\frac{X_{2,k}^{N_i}}{X_{1,k}^{N_i}})\right)+\frac{1}{49}\underset{k=1}{\overset{49}{\mathrm{\Σ}}}\log \left(\frac{1}{2}(\frac{X_{1,k}^{N_i}}{X_{2,k}^{N_i}}+\frac{X_{2,k}^{N_i}}{X_{1,k}^{N_i}})\right)\}$

Wherein, F represents disparity map I ₂pixel value, log () represents the operation of taking the logarithm, ∑ represents sum operation, N _iimage X after expression filtering ₁and X ₂in centered by i pixel, the neighborhood of size as 7 × 7, i represents image X after filtering ₁and X ₂in pixel sequence number, k represents neighborhood N _iin pixel sequence number, k=1,2 ..., 49,

image X after expression filtering ₁middle neighborhood N _ik pixel, image X after expression filtering ₂middle neighborhood N _ik pixel.

Step 4: tectonic syntaxis disparity map.

Adopt direction wave area image fusion method, to disparity map I ₁and I ₂merge, obtain associating disparity map I ₃.The concrete steps of direction wave area image fusion method are as follows:

The first step, respectively to disparity map I ₁and I ₂travel direction wave conversion, obtains disparity map I ₁and I ₂direction wave low frequency sub-band coefficient and high-frequency sub-band coefficient.

Second step, utilizes following formula, to disparity map I ₁and I ₂direction wave low frequency sub-band coefficient merge, obtain associating disparity map I ₃direction wave low frequency sub-band coefficient:

$L_{\mathrm{ij}}^{I_3}=\frac{1}{2}(L_{\mathrm{ij}}^{I_1}+L_{\mathrm{ij}}^{I_2}),i\∈M,j\∈N$

Wherein,

represent associating disparity map I ₃direction wave low frequency sub-band in the coefficient of capable, the j of i row,

represent disparity map I ₁the coefficient of capable, the j of direction wave low frequency sub-band i row, represent disparity map I ₂the coefficient of capable, the j of direction wave low frequency sub-band i row, M represents disparity map I ₁and I ₂line number, N represents disparity map I ₁and I ₂columns.

The 3rd step, utilizes following formula, to disparity map I ₁and I ₂direction wave high-frequency sub-band coefficient merge, obtain associating disparity map I ₃direction wave high-frequency sub-band coefficient:

$H_{\mathrm{ij}}^{I_3}=\max (H_{\mathrm{ij}}^{I_1},H_{\mathrm{ij}}^{I_2}),i\∈M,j\∈N$

Wherein,

represent associating disparity map I ₃direction wave high-frequency sub-band in the coefficient of capable, the j of i row, max () represents to get maxima operation,

represent disparity map I ₁high-frequency sub-band in the coefficient of capable, the j of i row,

represent disparity map I ₂direction wave low frequency sub-band in the coefficient of capable, the j of i row, M represents disparity map I ₁and I ₂line number, N represents disparity map I ₁and I ₂columns.

The 4th step, to associating disparity map I ₃direction wave low frequency coefficient and the inverse transformation of high frequency coefficient travel direction ripple, obtain associating disparity map I ₃.

Step 5: cut apart associating disparity map.

Adopt fuzzy C-means clustering method, will combine disparity map I ₃in all pixels be divided into strong variation class pixel, weak variation class pixel and do not change class pixel.The concrete steps of fuzzy C-means clustering method are as follows:

The first step, will combine disparity map I ₃in the gray-scale value of each pixel as the proper vector of this pixel, meanwhile, from associating disparity map I ₃in select arbitrarily the gray-scale value of 3 pixels as 3 initial cluster centers of 3 class pixels.

Second step, according to the following formula, upgrades associating disparity map I ₃the degree of membership of the proper vector of middle pixel:

$u_{\mathrm{ik}}=\mathrm{\Σ}{\left(\frac{\left|\right|x_k-v_i\left|\right|}{\left|\right|x_k-v_j\left|\right|}\right)}^{-2},j=\mathrm{1,2},...,n$

Wherein, u _ikrepresent associating disparity map I ₃in the proper vector of k pixel be under the jurisdiction of the degree of membership of i class pixel, ∑ is sum operation, x _krepresent associating disparity map I ₃in the proper vector of k pixel, v _irepresent the cluster centre of i class pixel, v _jrepresent the cluster centre of j class pixel, k represents associating disparity map I ₃the sequence number of the proper vector of middle pixel, i represents current cluster classification, and j represents cluster classification arbitrarily, and n represents associating disparity map I ₃the number of the proper vector of middle pixel, || || represent to ask Euclidean distance operation.

The 3rd step, according to the following formula, upgrade cluster centre vector:

$v_i=\frac{\mathrm{\Σ}{\left(u_{\mathrm{ik}}\right)}^2{x}_k}{\mathrm{\Σ}{\left(u_{\mathrm{ik}}\right)}^2},k=\mathrm{1,2},...,n$

Wherein, v _irepresent associating disparity map I ₃the cluster centre vector of i class pixel, ∑ is sum operation, μ _ikrepresent associating disparity map I ₃in the proper vector of k pixel be under the jurisdiction of the degree of membership of i class pixel, x _krepresent associating disparity map I ₃in the proper vector of k pixel, i represents the classification of cluster, k represents associating disparity map I ₃the sequence number of the proper vector of middle pixel, n represents associating disparity map I ₃the number of the proper vector of middle pixel.

The 4th step, carries out the first step to the three step iteration 20 times, obtains associating disparity map I ₃the proper vector of middle pixel be subordinate to matrix.

The 5th step, in being subordinate to matrix, the peaked line number of every row degree of membership, as the classification of pixel, obtains changing by force class pixel, weak variation class pixel and does not change class pixel.

The not variation class pixel that fuzzy C-means clustering method is obtained is as changing testing result image I _rin not variation class pixel, will change by force class pixel and weak variation class pixel and merge into and change class pixel, and as changing testing result image I _rin variation class pixel, obtain changing testing result image I _r.

Step 6: exporting change testing result image I _r.

Analogous diagram below in conjunction with accompanying drawing 2 and accompanying drawing 3 is described further effect of the present invention.

1, emulation experiment condition:

Hardware test platform of the present invention is: processor is Inter Core2Duo CPU E8200, and dominant frequency is 2.67GHz, internal memory 2GB, and software platform is: Windows7 Ultimate 32-bit operating system and Matlab R2012b.Input picture of the present invention is respectively the two width synthetic-aperture radar SAR images in Berne (Bern) area and the two width synthetic-aperture radar SAR images in Ottawa (Ottawa) area, wherein the regional two width synthetic-aperture radar SAR image sizes of Berne (Bern) are 301 × 301, the regional two width synthetic-aperture radar SAR image sizes in Ottawa (Ottawa) are 290 × 350, and form is all BMP.

2, emulation content:

Two methods that the prior art that the present invention uses contrasts are as follows respectively:

Patented technology " based on the sar image change detection method of neighborhood logarithm ratio and anisotropy the diffusion " (number of patent application: CN201110005209 that Xian Electronics Science and Technology University has, publication number: CN102096921B) in the synthetic-aperture radar SAR image change detection method of the middle proposition mentioned, be called for short based on maximum variance between clusters.

Patented technology " remote sensing image based on image co-registration changes the detection " (number of patent application: CN201210234076 that Xian Electronics Science and Technology University has, publication number: CN102750705) the middle remote sensing image change detecting method proposing, is called for short based on wavelet field image co-registration method.

3. analysis of simulation result:

Fig. 2 is that the inventive method and prior art are to the regional SAR Image Change Detection of Berne (Bern) comparison diagram, wherein, Fig. 2 (a) is Berne (Bern) the area synthetic-aperture radar SAR image in April, 1999, Fig. 2 (b) is Berne (Bern) the area synthetic-aperture radar SAR image in May, 1999, Fig. 2 (c) is the regional synthetic-aperture radar SAR Image Change Detection of Berne (Bern) reference picture, Fig. 2 (d) is the variation testing result based on maximum variance between clusters, Fig. 2 (e) is the variation testing result based on wavelet field image co-registration method, the variation testing result that Fig. 2 (f) is the inventive method.

Fig. 3 the inventive method and prior art are to Ottawa (Ottawa) regional SAR Image Change Detection comparison diagram, wherein, Fig. 3 (a) is area, Ottawa (Ottawa) the synthetic-aperture radar SAR image in May, 1997, Fig. 3 (b) is area, Ottawa (Ottawa) the synthetic-aperture radar SAR image in August, 1997, Fig. 3 (c) is Ottawa (Ottawa) regional synthetic-aperture radar SAR Image Change Detection reference picture, Fig. 3 (d) is the variation testing result based on maximum variance between clusters, Fig. 3 (e) is the variation testing result based on wavelet field image co-registration method, the variation testing result that Fig. 3 (f) is the inventive method.

Experimental result evaluating is respectively undetected number, false retrieval number, total wrong number and total false rate, table 1 is depicted as and in emulation experiment, changes testing result statistics, wherein, OSTU represents based on maximum variance between clusters, WDFCD represents based on wavelet field image co-registration method, ME represents that undetected number, FA represent that false retrieval number, OE represent total wrong number, and OER represents total false rate.

Table 1Bern area and Ottawa area testing result

Can find out in conjunction with Fig. 2, Fig. 3 and table 1, the poorest based on maximum variance between clusters testing result, main manifestations is that undetected number is too many, is especially embodied in multiple small targets perhaps and does not detect; Performance based on wavelet field image co-registration method is better than based on maximum variance between clusters, and false retrieval number and undetected number have decline, but the phenomenon that still exists less target to detect; The inventive method is no matter at false retrieval number be undetectedly all better than first two method aspect several, and especially many less point targets are detected, and illustrate that the inventive method is better than first two method in accuracy of detection.

Above emulation experiment shows: the inventive method is aspect the robustness to coherent speckle noise or accuracy of detection, all to show good performance, can solve in existing method, exist to the problem such as coherent speckle noise sensitivity, accuracy of detection be low, be a kind of very practical synthetic-aperture radar SAR image change detection method.

Claims (4)

1. the SAR image change detection method merging based on direction wave area image, comprises the steps:

(1) input picture:

The two width SAR image S that optional areal, different time are taken ₁and S ₂as waiting to change detection SAR image;

(2) image pre-service:

Adopt probability piecemeal matched filtering method, respectively to SAR image S ₁and S ₂filtering, obtains image X after filtering ₁and X ₂;

(3) generate disparity map:

(3a) utilize ratio operator formula, to image X after filtering ₁and X ₂in pixel value operate, obtain disparity map I ₁pixel value;

(3b) utilize hybrid operator formula, to image X after filtering ₁and X ₂in pixel value operate, obtain disparity map I ₂pixel value;

(4) tectonic syntaxis disparity map:

(4a) respectively to disparity map I ₁and I ₂travel direction wave conversion, obtains disparity map I ₁and I ₂direction wave low frequency sub-band coefficient and high-frequency sub-band coefficient;

(4b) utilize following formula, to disparity map I ₁and I ₂direction wave low frequency sub-band coefficient merge, obtain associating disparity map I ₃direction wave low frequency sub-band coefficient:

$L_{\mathrm{ij}}^{I_3}=\frac{1}{2}(L_{\mathrm{ij}}^{I_1}+L_{\mathrm{ij}}^{I_2}),i\∈M,j\∈N$

Wherein,

represent associating disparity map I ₃direction wave low frequency sub-band in the coefficient of capable, the j of i row,

represent disparity map I ₁the coefficient of capable, the j of direction wave low frequency sub-band i row,

represent disparity map I ₂the coefficient of capable, the j of direction wave low frequency sub-band i row, M represents disparity map I ₁and I ₂line number, N represents disparity map I ₁and I ₂columns;

(4c) utilize following formula, to disparity map I ₁and I ₂direction wave high-frequency sub-band coefficient merge, obtain associating disparity map I ₃direction wave high-frequency sub-band coefficient:

$H_{\mathrm{ij}}^{I_3}=\max (H_{\mathrm{ij}}^{I_1},H_{\mathrm{ij}}^{I_2}),i\∈M,j\∈N$

Wherein,

represent associating disparity map I ₃direction wave high-frequency sub-band in the coefficient of capable, the j of i row, max () represents to get maxima operation,

represent disparity map I ₁high-frequency sub-band in the coefficient of capable, the j of i row,

represent disparity map I ₂direction wave low frequency sub-band in the coefficient of capable, the j of i row, M represents disparity map I ₁and I ₂line number, N represents disparity map I ₁and I ₂columns;

(4d) to associating disparity map I ₃direction wave low frequency coefficient and the inverse transformation of high frequency coefficient travel direction ripple, obtain associating disparity map I ₃;

(5) cut apart associating disparity map:

(5a) adopt fuzzy C-means clustering method, as follows, will combine disparity map I ₃in all pixels be divided into strong variation class pixel, weak variation class pixel and do not change class pixel:

(5a1) will combine disparity map I ₃in the gray-scale value of all pixels as the proper vector of each pixel, meanwhile, from associating disparity map I ₃in select arbitrarily the gray-scale value of 3 pixels as 3 initial cluster centers of 3 class pixels;

(5a2) according to the following formula, upgrade associating disparity map I ₃the degree of membership of the proper vector of middle pixel:

$u_{\mathrm{ik}}=\mathrm{\Σ}{\left(\frac{\left|\right|x_k-v_i\left|\right|}{\left|\right|x_k-v_j\left|\right|}\right)}^{-2},j=\mathrm{1,2},...,n$

Wherein, u _ikrepresent associating disparity map I ₃in the proper vector of k pixel be under the jurisdiction of the degree of membership of i class pixel, ∑ is sum operation, x _krepresent associating disparity map I ₃in the proper vector of k pixel, v _irepresent the cluster centre of i class pixel, v _jrepresent the cluster centre of j class pixel, k represents associating disparity map I ₃the sequence number of the proper vector of middle pixel, i represents current cluster classification, and j represents cluster classification arbitrarily, and n represents associating disparity map I ₃the number of the proper vector of middle pixel, || || represent to ask Euclidean distance operation;

(5a3) according to the following formula, upgrade cluster centre vector:

$v_i=\frac{\mathrm{\Σ}{\left(u_{\mathrm{ik}}\right)}^2{x}_k}{\mathrm{\Σ}{\left(u_{\mathrm{ik}}\right)}^2},k=\mathrm{1,2},...,n$

Wherein, v _irepresent associating disparity map I ₃the cluster centre vector of i class pixel, ∑ is sum operation, μ _ikrepresent associating disparity map I ₃in the proper vector of k pixel be under the jurisdiction of the degree of membership of i class pixel, x _krepresent associating disparity map I ₃in the proper vector of k pixel, i represents the classification of cluster, k represents associating disparity map I ₃the sequence number of the proper vector of middle pixel, n represents associating disparity map I ₃the number of the proper vector of middle pixel;

(5a4) first step to the three step iteration are carried out 20 times, obtained associating disparity map I ₃the proper vector of middle pixel be subordinate to matrix;

(5a5) in being subordinate to matrix, the peaked line number of every row degree of membership, as the classification of pixel, obtains changing by force class pixel, weak variation class pixel and does not change class pixel;

(5b) will change by force class pixel and weak variation class pixel and merge into variation class pixel;

(5c) using the not variation class pixel in step (5a) as changing testing result image I _rin not variation class pixel, using the variation class pixel in step (5b) as changing testing result image I _rin variation class pixel, obtain changing testing result image I _r;

(6) output image:

Exporting change testing result image I _r.

2. the SAR image change detection method merging based on direction wave area image according to claim 1, is characterized in that: the described probability piecemeal matched filtering method of step (2) refers to, as follows respectively to SAR image S ₁and S ₂filtering, obtains image X after filtering ₁and X ₂:

The first step, sets SAR image S ₁and S ₂in arbitrarily pixel neighborhood of a point window size be 7 × 7;

Second step, according to the following formula, calculates respectively SAR image S ₁and S ₂in any weights of pixel and interior other pixels of this pixel neighborhood of a point window:

$w(q,t)=\exp (-\underset{p=1}{\overset{49}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{q,p}}{A_{t,r}}+\frac{A_{t,r}}{A_{q,p}}))$

Wherein, w (q, t) represents SAR image S ₁and S ₂in any weights of pixel q and interior other pixels t of this pixel neighborhood of a point window, exp () represents the operation of fetching number, ∑ represents sum operation, h represents the smoothing parameter of control characteristic attenuation degree, log () represents to take the logarithm operation, and p represents any pixel q pixel sequence number in the neighborhood of 7 × 7 sizes, and r represents the sequence number of pixel t pixel in 7 × 7 large small neighbourhoods, and r=p, A _q,prepresent the gray-scale value of any pixel q pixel sequence number p corresponding position in neighborhood, A _t,rrepresent the gray-scale value of pixel t pixel sequence number r corresponding position in neighborhood;

The 3rd step, according to the following formula, calculates respectively SAR image S ₁and S ₂in any gray scale estimated value of pixel:

${\hat{z}}_1\left(q\right)=\frac{\underset{t\∈W_q}{\mathrm{\Σ}}w(q,t)A_t^2}{\underset{t\∈W_q}{\mathrm{\Σ}}w(q,t)}$

Wherein, represent SAR image S ₁and S ₂in any gray scale estimated value of pixel q, ∑ represents sum operation, W _qrepresent SAR image S ₁and S ₂in any neighborhood window of 7 × 7 sizes of pixel q, w (q, t) represents SAR image S ₁and S ₂in any pixel q and its neighborhood window in the weights of other pixels t, A _trepresent the gray-scale value of pixel t;

The 4th step, carries out the first step, second step, the 3rd step iteration 25 times, obtains SAR image S ₁and S ₂filtering after image X ₁and X ₂.

3. the SAR image change detection method merging based on direction wave area image according to claim 1, is characterized in that: the ratio operator formula described in step (3a) is as follows:

$D=\frac{\underset{k=1}{\overset{49}{\mathrm{\Σ}}}\min (X_{1,k}^{N_i},X_{2,k}^{N_i})}{\underset{k=1}{\overset{49}{\mathrm{\Σ}}}\max (X_{1,k}^{N_i},X_{2,k}^{N_i})}$

Wherein, D represents disparity map I ₁pixel value, ∑ represents sum operation, min () represents to get minimum value operation, max () represents to get maxima operation, N _iimage X after expression filtering ₁and X ₂in centered by i pixel, the neighborhood of size as 7 × 7, i represents image X after filtering ₁and X ₂in pixel sequence number, k represents neighborhood N _iin pixel sequence number, k=1,2 ..., 49,

image X after expression filtering ₁middle neighborhood N _ik pixel,

image X after expression filtering ₂middle neighborhood N _ik pixel.

4. the SAR image change detection method merging based on direction wave area image according to claim 1, is characterized in that: the hybrid operator formula described in step (3b) is as follows:

$F=\frac{1}{2}\{\log \left(\frac{1}{2}(\frac{X_{1,k}^{N_i}}{X_{2,k}^{N_i}}+\frac{X_{2,k}^{N_i}}{X_{1,k}^{N_i}})\right)+\frac{1}{49}\underset{k=1}{\overset{49}{\mathrm{\Σ}}}\log \left(\frac{1}{2}(\frac{X_{1,k}^{N_i}}{X_{2,k}^{N_i}}+\frac{X_{2,k}^{N_i}}{X_{1,k}^{N_i}})\right)\}$

Wherein, F represents disparity map I ₂pixel value, log () represents the operation of taking the logarithm, ∑ represents sum operation, N _iimage X after expression filtering ₁and X ₂in centered by i pixel, the neighborhood of size as 7 × 7, i represents image X after filtering ₁and X ₂in pixel sequence number, k represents neighborhood N _iin pixel sequence number, k=1,2 ..., 49, image X after expression filtering ₁middle neighborhood N _ik pixel,

image X after expression filtering ₂middle neighborhood N _ik pixel.

Images