CN103824302B - The SAR image change detection merged based on direction wave area image - Google Patents

Abstract

The invention discloses a kind of SAR image change detection merged based on direction wave area image, overcome disparity map in prior art and when merging, can not preferably retain disparity map detailed information, the problem that disparity map uses binarization segmentation method can cause region of variation information dropout.The step that the present invention realizes is: (1) input picture；(2) Image semantic classification；(3) disparity map is generated；(4) tectonic syntaxis disparity map；(5) segmentation associating disparity map；(6) output image.The present invention has the advantage describing region of variation with not changing zone boundary more accurately, loss is relatively low, can be applicable to the natural disaster detection in synthetic aperture radar SAR image field of target recognition.

Description

SAR image change detection method based on directional wave domain image fusion

Technical Field

The invention belongs to the technical field of image processing, and further relates to a Synthetic Aperture Radar (SAR) image change detection method based on directional wave domain image fusion in the technical field of SAR image change detection. The method can be used for natural disaster detection in the field of SAR image target identification.

Background

The research object of change detection is multi-temporal data, and requires to acquire imaging data of an observation target at different moments, and the purpose of change detection is to extract change information of the observation object at different moments by analyzing the multi-temporal data. Synthetic Aperture Radar (SAR) imaging is used as an important ground observation means, particularly an active microwave imaging principle, so that the SAR can be observed all weather and all day long without being influenced by weather and illumination, in addition, the penetration capability of microwave enables the SAR to penetrate through forests, and the SAR has the detection capability on maneuvering targets hidden in forest regions, so that SAR images become important information sources for image change detection.

The patent application of the university of electronic science and technology of Xian in the 'sar image change detection method based on neighborhood logarithm ratio and anisotropic diffusion' (patent application number: CN201110005209, publication number: CN 102096921A) provides an sar image change detection method based on neighborhood logarithm ratio and anisotropic diffusion. The method comprises the following specific steps: firstly, constructing a difference graph of two time phase graphs by adopting a neighborhood logarithm ratio method, and then carrying out self-adaptive anisotropic diffusion filtering and maximum inter-class variance method segmentation on the difference graph to obtain a final change detection result. Although the method can better maintain the edge information of the change area, the method still has the defect that the change area information is lost and the undetected rate of change detection is high due to the fact that the binary segmentation is not suitable for describing the pixel types of the combined difference map.

The Sigan electronic science and technology university provides an optical remote sensing image change detection method based on image fusion in the patent application 'optical remote sensing image change detection based on image fusion' (patent application number: CN2012102340761, publication number: CN 102750705A). The method comprises the following specific steps: firstly, two difference maps are fused in a wavelet domain to obtain a combined difference map, and then the fused difference map is segmented by using a fuzzy local C-means clustering method to obtain a change detection result. Although the method improves the precision of the detection result by improving the performance of the difference map of the optical remote sensing image, the method still has the defect that the detail information of the difference map cannot be well reserved when the difference map is fused because the extraction performance of wavelet transformation on line singular information is poor, so that the boundary detection of a changed area and an unchanged area in the change detection result is inaccurate.

Disclosure of Invention

The invention aims to overcome the defects of the prior art and provides a SAR image change detection method based on directional wave domain image fusion. According to the SAR image change detection method, the accuracy of the SAR image change detection result is improved by improving the performance of the difference map and improving the segmentation strategy of the difference map.

In order to achieve the purpose, the method comprises the following specific steps:

(1) inputting an image:

two SAR images S optionally shot in the same area and at different time₁And S₂As the SAR image to be changed.

(2) Image preprocessing:

respectively processing SAR image S by adopting a probabilistic block matching filtering method₁And S₂Filtering to obtain filtered image X₁And X₂。

(3) Generating a difference map:

(3a) using ratio operator formula to filter image X₁And X₂The pixel values in (1) are operated to obtain a difference map I₁A pixel value of (a);

(3b) using a hybrid operator formula to filter the image X₁And X₂The pixel values in (1) are operated to obtain a difference map I₂The pixel value of (2).

(4) Constructing a joint difference map:

(4a) respectively to difference chart I₁And I₂Performing directional wave transformation to obtain a difference chart I₁And I₂The directional wave low-frequency subband coefficient and the high-frequency subband coefficient;

(4b) using the following formula, for difference plot I₁And I₂The low-frequency subband coefficients of the directional waves are fused to obtain a combined difference chart I₃Directional wave low frequency subband coefficient of (2):

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

wherein,representation of joint disparity map I₃The coefficients of the ith row and the jth column in the directional wave low-frequency sub-band,shows a difference chart I₁The coefficients of the ith row and the jth column of the directional wave low-frequency sub-band of (1),shows a difference chart I₂The coefficients of the ith row and the jth column of the directional wave low-frequency sub-band of (1), M represents a difference map I₁And I₂N denotes the disparity map I₁And I₂The number of columns;

(4c) using the following formula, for difference plot I₁And I₂The high-frequency sub-band coefficients of the directional waves are fused to obtain a combined difference chart I₃Directional wave high frequency subband coefficient of (2):

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

wherein,representation of joint disparity map I₃The coefficient of the ith row and the jth column in the directional wave high-frequency sub-band of (1), max (-) represents the maximum value operation,shows a difference chart I₁The ith row and the jth column of the high-frequency subband of (2),shows a difference chart I₂M represents the difference map I₁And I₂N denotes the disparity map I₁And I₂The number of columns;

(4d) to joint difference chart I₃The low-frequency coefficient and the high-frequency coefficient of the directional wave are subjected to inverse directional wave transformation to obtain a combined difference chart I₃。

(5) Partitioning the joint difference map:

(5a) adopting a fuzzy C-means clustering method, and combining the difference chart I according to the following steps₃All the pixels are divided into strong change pixels, weak change pixels and unchanged pixels:

(5a1) will combine difference maps I₃The gray value of all the pixels is used as the characteristic vector of each pixel, and meanwhile, the gray value of all the pixels is used as the characteristic vector of each pixel, and the gray value is extracted from the joint difference map I₃The gray value of 3 pixels is arbitrarily selected as 3 initial clustering centers of the 3 types of pixels;

(5a2) updating the joint difference map I according to₃Membership of feature 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 is_ikRepresentation of joint disparity map I₃The feature vector of the kth pixel belongs to the ith pixel∑ is the summation operation, x_kRepresentation of joint disparity map I₃Characteristic vector of the k-th pixel, v_iRepresenting the cluster center of the ith class of pixels, v_jDenotes the cluster center of the jth pixel, k denotes the joint difference map I₃The serial number of the feature vector of the middle pixel, I represents the current cluster type, j represents any cluster type, and n represents the joint difference graph I₃The number of feature vectors of the middle pixel, | | | · | |, represents the Euclidean distance solving operation;

(5a3) the cluster center vector is updated as follows:

$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 is_iRepresentation of joint disparity map I₃Cluster center vector for class i pixels, ∑ for summation operation, μ_ikRepresentation of joint disparity map I₃The feature vector of the kth pixel is the membership degree, x, of the ith pixel_kRepresentation of joint disparity map I₃The feature vector of the k-th pixel in (a), I represents the class of the cluster, and k represents the joint disparity map I₃The serial number of the feature vector of the middle pixel, n represents the joint difference map I₃The number of feature vectors of the middle pixel;

(5a4) performing the first step to the third step for 20 times in an iteration way to obtain a joint difference graph I₃A membership matrix of eigenvectors of the middle pixel;

(5a5) taking the row number of the maximum value of the membership degree of each column in the membership matrix as the class of the pixel to obtain a strong variation class pixel, a weak variation class pixel and an unchanged class pixel;

(5b) combining the strong variation pixels and the weak variation pixels into variation pixels;

(5c) using the unchanged pixels in the step (5a) as the change detection result image I_rThe non-change pixels in (5b) are used as the change detection result image I_rTo obtain a change detection result image I_r。

(6) Outputting an image:

outputting a change detection result image I_r。

Compared with the prior art, the method of the invention has the following advantages:

firstly, because the two difference maps are fused by adopting the direction wave domain image fusion method, the defect that the detail information of the difference maps cannot be well reserved during the fusion of the difference maps in the prior art is overcome, and the performance of extracting the detail information by combining the difference maps is improved, so that the method has the advantage of more accurately describing the boundary of the changed region and the unchanged region.

Secondly, because the difference image pixels are divided into the strong variation pixels, the weak variation pixels and the unchanged pixels by adopting the fuzzy C-means method, the defect that the information of the variation area is lost due to the adoption of a binarization segmentation method for the difference image in the prior art is overcome, the information of the variation area in the combined difference image is more accurately reserved, and the method has the advantage of low omission ratio.

Drawings

FIG. 1 is a flow chart of the present invention;

FIG. 2 is a comparison graph of the SAR image change detection effect of Bern (Bern) region in the method of the present invention and the prior art;

fig. 3 is a comparison graph of the detection effect of the method of the present invention on the change of the SAR image in the Ottawa (Ottawa) area.

Detailed Description

The invention is further described below with reference to the accompanying drawings.

Referring to the attached figure 1, the implementation steps of the invention are as follows:

step 1: an image is input.

Two SAR images S optionally shot in the same area and at different time₁And S₂The synthetic aperture radar SAR images used in the embodiment of the present invention are shown in fig. 2(a), fig. 2(b), fig. 3(a) and fig. 3(b), respectively, wherein fig. 2(a) is a synthetic aperture radar SAR image of 4 months 1999 in Bern (Bern) area, 301 × 301 in size, fig. 2(b) is a synthetic aperture radar SAR image of 5 months 1999 in Bern (Bern) area, 301 × 301 in size, fig. 3(a) is a synthetic aperture radar SAR image of 5 months 1997 in Ottawa (Ottawa) area, 290 × 350 in size, and fig. 3(b) is a synthetic aperture radar image of 8 months 1997 in Ottawa (Ottawa) area, 290 SAR 290 × 350 in size.

Step 2: and (5) image preprocessing.

Respectively processing SAR image S by adopting a probabilistic block matching filtering method₁And S₂Filtering to obtain filtered image X₁And X₂. The probabilistic block matching filtering method comprises the following specific steps:

firstly, setting SAR image S₁And S₂The size of the neighborhood window for any pixel point in the array is 7 × 7.

Secondly, respectively calculating SAR images S according to the following formula₁And S₂The weight of any pixel point and other pixel points in the neighborhood window of the pixel point is as follows:

$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}}))$

where w (q, t) denotes the SAR image S₁And S₂The weight values of any pixel point q and other pixel points t in a neighborhood window of the pixel point are exp (·) representing the index taking operation, ∑ representing the summation operation, h representing a smooth parameter for controlling exponential attenuation degree, log (·) representing the logarithm taking operation, p representing the serial number of the pixel point of any pixel point q in a neighborhood with the size of 7 × 7, r representing the serial number of the pixel point of pixel point t in a neighborhood with the size of 7 × 7, and r = p, A_q,pExpressing the gray value A of any pixel point q at the position corresponding to the pixel point number p in the neighborhood_t,rAnd expressing the gray value of the pixel point t at the position corresponding to the pixel point sequence number r in the neighborhood.

Thirdly, respectively calculating the SAR images S according to the following formula₁And S₂The gray scale estimation value of any pixel point:

${\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,representing SAR images S₁And S₂The gray scale estimation value of any pixel point q in the tree, ∑ represents the summation operation, W_qRepresenting SAR images S₁And S₂A neighborhood window of 7 × 7 size of any pixel point q in the space, w (q, t) represents SAR image S₁And S₂The weight value of any pixel point q and other pixel points t in the neighborhood window, A_tAnd representing the gray value of the pixel point t.

Fourthly, the first step, the second step and the third step are iterated for 25 times to obtain an SAR image S₁And S₂Filtered image X of₁And X₂。

And step 3: and generating a difference map.

Using ratio operator formula to filter image X₁And X₂The pixel values in (1) are operated to obtain a difference map I₁The pixel value of (2). The 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 a disparity map I₁∑ denotes a sum operation, min (-) denotes a min operation, max (-) denotes a max operation, N_iRepresenting the filtered image X₁And X₂Of a neighborhood of size 7 × 7 centered on the ith pixel, i representing the filtered image X₁And X₂K denotes the pixel number in the neighborhood Ni, k =1,2, ·,49,representing the filtered image X₁Middle neighborhood N_iThe number k of the pixels of (a),representing the filtered image X₂Middle neighborhood N_iThe kth pixel of (1).

Using a hybrid operator formula to filter the image X₁And X₂The pixel values in (1) are operated to obtain a difference map I₂The pixel value of (2). The 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 a disparity map I₂Log (-) denotes a logarithm operation, ∑ denotes a sum operation, N_iRepresenting the filtered image X₁And X₂Of a neighborhood of size 7 × 7 centered on the ith pixel, i representing the filtered image X₁And X₂Pixel number in (1), k represents neighborhood N_iThe pixel number in (1), k =1,2, ·,49,representing the filtered image X₁Middle neighborhood N_iThe number k of the pixels of (a),representing the filtered image X₂Middle neighborhood N_iThe kth pixel of (1).

And 4, step 4: a joint difference map is constructed.

Adopting a direction wave domain image fusion method to perform difference image I₁And I₂Performing fusion to obtain a combined difference map I₃. The method for fusing the directional wave domain images comprises the following specific steps:

first, to difference chart I₁And I₂Performing directional wave transformation to obtain a difference chart I₁And I₂Low frequency subband coefficients and high frequency subband coefficients.

Second, using the following formula, for the difference plot I₁And I₂The low-frequency subband coefficients of the directional waves are fused to obtain a combined difference chart I₃Directional wave low frequency subband coefficient of (2):

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

wherein,representation of joint disparity map I₃The coefficients of the ith row and the jth column in the directional wave low-frequency sub-band,shows a difference chart I₁The coefficients of the ith row and the jth column of the directional wave low-frequency sub-band of (1),shows a difference chart I₂The coefficients of the ith row and the jth column of the directional wave low-frequency sub-band of (1), M represents a difference map I₁And I₂N denotes the disparity map I₁And I₂The number of columns.

Third, using the following formula, for the difference chart I₁And I₂The high-frequency sub-band coefficients of the directional waves are fused to obtain a combined difference chart I₃Directional wave high frequency subband coefficient of (2):

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

wherein,representation of joint disparity map I₃The coefficient of the ith row and the jth column in the directional wave high-frequency sub-band of (1), max (-) represents the maximum value operation,shows a difference chart I₁The ith row and the jth column of the high-frequency subband of (2),shows a difference chart I₂M represents the difference map I₁And I₂N denotes the disparity map I₁And I₂The number of columns.

The fourth step is to combine the difference map I₃The low-frequency coefficient and the high-frequency coefficient of the directional wave are subjected to inverse directional wave transformation to obtain a combined difference chart I₃。

And 5: and (5) segmenting the joint difference map.

Using mouldsFuzzy C mean clustering method, combining difference chart I₃All the pixels are divided into pixels with strong variation, pixels with weak variation and pixels with unchanged variation. The fuzzy C-means clustering method comprises the following specific steps:

first, combine the difference maps I₃The gray value of each pixel is used as the characteristic vector of the pixel, and meanwhile, the joint difference map I₃The gray value of 3 pixels is arbitrarily selected as 3 initial clustering centers of the 3 types of pixels.

Second, the joint difference map I is updated according to the following formula₃Membership of feature 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 is_ikRepresentation of joint disparity map I₃The membership of the eigenvector of the kth pixel to the i-th pixel, ∑ is the summation operation, x_kRepresentation of joint disparity map I₃Characteristic vector of the k-th pixel, v_iRepresenting the cluster center of the ith class of pixels, v_jDenotes the cluster center of the jth pixel, k denotes the joint difference map I₃The serial number of the feature vector of the middle pixel, I represents the current cluster type, j represents any cluster type, and n represents the joint difference graph I₃The number of feature vectors of the middle pixel, | | | · | |, represents the euclidean distance solving operation.

Thirdly, updating the clustering center vector according to the following formula:

$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 is_iRepresentation of joint disparity map I₃Cluster center vector for class i pixels, ∑ for summation operation, μ_ikRepresentation of joint disparity map I₃The feature vector of the kth pixel is the membership degree, x, of the ith pixel_kRepresentation of joint disparity map I₃The feature vector of the k-th pixel in (a), I represents the class of the cluster, and k represents the joint disparity map I₃The serial number of the feature vector of the middle pixel, n represents the joint difference map I₃The number of feature vectors of the middle pixel.

Fourthly, the first step to the third step are iterated for 20 times to obtain a joint difference graph I₃A membership matrix of eigenvectors of the middle pixel.

And fifthly, taking the row number of the maximum value of the membership degree of each column in the membership matrix as the class of the pixel to obtain the strong variation class pixel, the weak variation class pixel and the unchanged class pixel.

Using unchanged pixels obtained by the fuzzy C-means clustering method as a change detection result image I_rThe strong change pixels and the weak change pixels are combined into change pixels and used as change detection result images I_rTo obtain a change detection result image I_r。

Step 6: outputting a change detection resultImage I_r。

The effect of the present invention will be further explained with reference to the simulation diagrams of fig. 2 and fig. 3.

1, simulation experiment conditions:

the hardware test platform of the invention is: the processor is an Inter Core2Duo CPU E8200, the dominant frequency is 2.67GHz, the memory is 2GB, and the software platform is as follows: windows7 flagship version 32-bit operating system and MatlabR2012 b. The input images of the invention are two synthetic aperture radar SAR images in a Bern (Bern) area and two synthetic aperture radar SAR images in an Ottawa (Ottawa) area respectively, wherein the two synthetic aperture radar SAR images in the Bern (Bern) area are 301 multiplied by 301, the two synthetic aperture radar SAR images in the Ottawa (Ottawa) area are 290 multiplied by 350, and the formats of the two synthetic aperture radar SAR images are BMP.

2, simulation content:

the two methods used in the prior art for comparison are respectively as follows:

the patent technology owned by the university of west' an electronic technology "SAR image change detection method based on neighborhood logarithm ratio and anisotropic diffusion" (patent application number: CN201110005209, publication number: CN 102096921B) provides a synthetic aperture radar SAR image change detection method, which is based on the maximum inter-class variance method for short.

An optical remote sensing image change detection method provided in the patent technology ' optical remote sensing image change detection based on image fusion ' (patent application No. CN201210234076, publication No. CN 102750705) owned by the university of Western ' An electronic technology, is called wavelet domain image fusion method for short.

3. And (3) simulation result analysis:

fig. 2 is a comparison diagram of detecting changes of SAR images of berne (Bern) regions according to the present invention and the prior art, where fig. 2(a) is a synthetic aperture radar SAR image of the berne (Bern) region in 4 months 1999, fig. 2(b) is a synthetic aperture radar SAR image of the berne (Bern) region in 5 months 1999, fig. 2(c) is a reference image of detecting changes of the synthetic aperture radar SAR images of the Bern (Bern) regions, fig. 2(d) is a change detection result based on the maximum inter-class variance method, fig. 2(e) is a change detection result based on the wavelet domain image fusion method, and fig. 2(f) is a change detection result of the present invention.

Fig. 3 is a comparison graph of the change detection of the SAR image in the Ottawa (Ottawa) area in comparison with the prior art, where fig. 3(a) is the synthetic aperture radar SAR image in the Ottawa (Ottawa) area in 5 months 1997, fig. 3(b) is the synthetic aperture radar SAR image in the Ottawa (Ottawa) area in 8 months 1997, fig. 3(c) is the change detection reference image of the synthetic aperture radar SAR image in the Ottawa (Ottawa) area, fig. 3(d) is the change detection result based on the maximum inter-class variance method, fig. 3(e) is the change detection result based on the wavelet domain image fusion method, and fig. 3(f) is the change detection result of the method of the present invention.

The experimental result evaluation parameters are respectively the number of missed detections, the number of false detections, the total number of errors and the total error rate, and table 1 shows statistics of change detection results in a simulation experiment, wherein, OSTU represents the maximum inter-class variance method, WDFCD represents the wavelet domain image fusion method, ME represents the number of missed detections, FA represents the number of false detections, OE represents the total error number, and OER represents the total error rate.

TABLE 1Bern area and Ottawa area test results

As can be seen from fig. 2, fig. 3 and table 1, the detection result is the worst based on the variance method between the maximum classes, which mainly shows that the number of missed detections is too large, and is especially reflected in that many small targets are not detected; the performance of the wavelet domain-based image fusion method is superior to that of the maximum inter-class variance method, the false detection number and the missed detection number are reduced, but the phenomenon that a smaller target cannot be detected still exists; the method is superior to the former two methods in terms of error detection number and omission number, and particularly, a plurality of smaller point targets are detected, which shows that the method is superior to the former two methods in terms of detection precision.

The above simulation experiments show that: the method of the invention has good performance in the aspects of robustness to speckle noise and detection precision, can solve the problems of sensitivity to speckle noise, low detection precision and the like in the existing method, and is a very practical synthetic aperture radar SAR image change detection method.

Claims (1)

1. A SAR image change detection method based on directional wave domain image fusion comprises the following steps:

(1) inputting an image:

two SAR images S optionally shot in the same area and at different time₁And S₂As a SAR image to be changed and detected;

(2) image preprocessing:

respectively processing SAR image S by adopting a probabilistic block matching filtering method₁And S₂Filtering to obtain filtered image X₁And X₂；

The probabilistic block matching filtering method is to respectively carry out the SAR image S according to the following steps₁And S₂Filtering to obtain filtered image X₁And X₂：

Firstly, setting SAR image S₁And S₂The size of a neighborhood window of any pixel point in the image is 7 × 7;

secondly, respectively calculating SAR images S according to the following formula₁And S₂The weight of any pixel point and other pixel points in the neighborhood window of the pixel point is as follows:

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

where w (q, t) denotes the SAR image S₁And S₂The weight of any pixel point q and other pixel points t in a neighborhood window of the pixel point is exp (·) to represent the operation of fetching the exponent, ∑ to represent the operation of summation, h to represent the smooth parameter for controlling the exponential attenuation degree, log (·) to represent the operation of fetching the logarithm, p to represent the serial number of the pixel point of any pixel point q in the neighborhood of 7 × 7 size, r to represent the serial number of the pixel point of pixel point t in the neighborhood of 7 × 7 size, and r ═ p, A_q,pExpressing the gray value A of any pixel point q at the position corresponding to the pixel point number p in the neighborhood_t,rExpressing the gray value of the pixel point t at the position corresponding to the pixel point serial number r in the neighborhood;

thirdly, calculating according to the following formulaSAR image S₁And S₂The gray scale estimation value of any pixel point:

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

wherein,representing SAR images S₁And S₂The gray scale estimation value of any pixel point q in the tree, ∑ represents the summation operation, W_qRepresenting SAR images S₁And S₂A neighborhood window of 7 × 7 size of any pixel point q in the space, w (q, t) represents SAR image S₁And S₂The weight value of any pixel point q and other pixel points t in the neighborhood window, A_tExpressing the gray value of the pixel point t;

fourthly, the first step, the second step and the third step are iterated for 25 times to obtain an SAR image S₁And S₂Filtered image X of₁And X₂；

(3) Generating a difference map:

(3a) using ratio operator formula to filter image X₁And X₂The pixel values in (1) are operated to obtain a difference map I₁A pixel value of (a);

the ratio operator formula is as follows:

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

wherein D represents a disparity map I₁∑ denotes a sum operation, min (-) denotes a min operation, max (-) denotes a max operation, N_iRepresenting the filtered image X₁And X₂Of a neighborhood of size 7 × 7 centered on the ith pixel, i representing the filtered image X₁And X₂Pixel number in (1), k represents neighborhood N_iThe pixel number in (1), 2, …,49,representing the filtered image X₁Middle neighborhood N_iThe number k of the pixels of (a),representing the filtered image X₂Middle neighborhood N_iThe kth pixel of (1);

(3b) using a hybrid operator formula to filter the image X₁And X₂The pixel values in (1) are operated to obtain a difference map I₂A pixel value of (a);

the formula of the mixed operator is as follows:

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

wherein F represents a disparity map I₂Log (-) denotes a logarithm operation, ∑ denotes a sum operation, N_iRepresenting the filtered image X₁And X₂Of a neighborhood of size 7 × 7 centered on the ith pixel, i representing the filtered image X₁And X₂Pixel number in (1), k represents neighborhood N_iThe pixel number in (1), 2, …,49,representing the filtered image X₁Middle neighborhood N_iThe number k of the pixels of (a),representing the filtered image X₂Middle neighborhood N_iThe kth pixel of (1);

(4) constructing a joint difference map:

(4a) respectively to difference chart I₁And I₂Performing directional wave transformation to obtain a difference chart I₁And I₂The directional wave low-frequency subband coefficient and the high-frequency subband coefficient;

(4b) using the following formula, for difference plot I₁And I₂The low-frequency subband coefficients of the directional waves are fused to obtain a combined difference chart I₃Directional wave low frequency subband coefficient of (2):

$L_{hv}^{I_3}=\frac{1}{2}(L_{hv}^{I_1}+L_{hv}^{I_2}),h\∈M,v\∈N$

wherein,representation of joint disparity map I₃The coefficients of the h-th row and the v-th column in the directional wave low-frequency sub-band,shows a difference chart I₁The coefficients of the h row and the v column of the directional wave low-frequency sub-band,shows a difference chart I₂M represents the coefficient of the h-th line and v-th column of the directional wave low-frequency subband, and M represents the difference map I₁And I₂N denotes the disparity map I₁And I₂The number of columns;

(4c) using the following formula, for difference plot I₁And I₂The high-frequency sub-band coefficients of the directional waves are fused to obtain a combined difference chart I₃Directional wave high frequency subband coefficient of (2):

$H_{hv}^{I_3}=max(H_{hv}^{I_1},H_{hv}^{I_2}),h\∈M,v\∈N$

wherein,representation of joint disparity map I₃The coefficients of the h-th row and the v-th column in the directional wave high-frequency sub-band of (1), max (-) represents the maximum value operation,shows a difference chart I₁The h-th row and the v-th column of the high-frequency subband,shows a difference chart I₂M represents the coefficient of the h-th line and the v-th column in the directional wave low-frequency subband, and M represents the difference map I₁And I₂N denotes the disparity map I₁And I₂The number of columns;

(4d) to joint difference chart I₃Low frequency coefficient and high frequency of the directional waveThe coefficients are subjected to inverse directional wave transformation to obtain a joint difference map I₃；

(5) Partitioning the joint difference map:

(5a) adopting a fuzzy C-means clustering method, and combining the difference chart I according to the following steps₃All the pixels are divided into strong change pixels, weak change pixels and unchanged pixels:

(5a1) will combine difference maps I₃The gray value of all the pixels is used as the characteristic vector of each pixel, and meanwhile, the gray value of all the pixels is used as the characteristic vector of each pixel, and the gray value is extracted from the joint difference map I₃The gray value of 3 pixels is arbitrarily selected as 3 initial clustering centers of the 3 types of pixels;

(5a2) updating the joint difference map I according to₃Membership of feature vector of middle pixel:

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

wherein u is_ikRepresentation of joint disparity map I₃The membership of the eigenvector of the kth pixel to the i-th pixel, ∑ is the summation operation, x_kRepresentation of joint disparity map I₃Characteristic vector of the k-th pixel, v_iRepresenting the cluster center of the ith class of pixels, v_jDenotes the cluster center of the jth pixel, k denotes the joint difference map I₃The serial number of the feature vector of the middle pixel, I represents the current cluster type, j represents any cluster type, and n represents the joint difference graph I₃The number of feature vectors of the middle pixel, | | | · | |, represents the Euclidean distance solving operation;

(5a3) the cluster center vector is updated as follows:

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

wherein v is_iRepresentation of joint disparity map I₃Cluster center vector for class i pixels, ∑ for summation operation, μ_ikRepresentation of joint disparity map I₃The feature vector of the kth pixel is the membership degree, x, of the ith pixel_kRepresentation of joint disparity map I₃The feature vector of the k-th pixel in (a), I represents the class of the cluster, and k represents the joint disparity map I₃The serial number of the feature vector of the middle pixel, n represents the joint difference map I₃The number of feature vectors of the middle pixel;

(5a4) iteratively executing (5a1) to (5a3) 20 times to obtain a joint difference map I₃A membership matrix of eigenvectors of the middle pixel;

(5a5) taking the row number of the maximum value of the membership degree of each column in the membership matrix as the class of the pixel to obtain a strong variation class pixel, a weak variation class pixel and an unchanged class pixel;

(5b) combining the strong variation pixels and the weak variation pixels into variation pixels;

(5c) using the unchanged pixels in the step (5a) as the change detection result image I_rThe non-change pixels in (5b) are used as the change detection result image I_rTo obtain a change detection result image I_r；

(6) Outputting an image:

outputting a change detection resultImage I_r。