SAR (Synthetic Aperture Radar) image change detection method fusing gravitational field and difference chart

Abstract

The invention discloses an SAR image change detection method fusing a gravitational field and a difference chart, and the method comprises the following steps: obtaining a dual-time-phase SAR image of the same region, and carrying out the bilateral filtering preprocessing, and obtaining a denoised image; the method comprises the following steps of: establishing mapping between image gray and physical quality, extracting a structural candidate region by using a high-order gravitational field model of cubic-level distance attenuation, introducing geometric anisotropy to screen and eliminate false angular point noise, and generating a physical structure mask through adaptive matching and multi-scale morphological expansion; calculating a logarithm ratio difference chart and generating an enhanced difference chart and a statistical significance mask through local texture statistical factor modulation; performing nonlinear multiplicative fusion on the difference chart by using physical and statistical dual masks as gating factors; and finally, obtaining a final result by adopting a self-adaptive threshold segmentation algorithm with closed-loop feedback. According to the method, the dependence on a difference chart is effectively solved, the false alarm rate is remarkably reduced, meanwhile, detection holes are effectively filled, and the detection precision and integrity are improved.

Description

Technical Field

[0001] This invention relates to the field of radar remote sensing image processing technology, and specifically to a method for detecting changes in SAR images by fusing gravitational field and difference maps. Background Technology

[0002] With the rapid development of remote sensing technology, Synthetic Aperture Radar (SAR) plays an irreplaceable role in fields such as natural disaster assessment, environmental monitoring, military reconnaissance, and urban planning due to its all-weather, all-day imaging capabilities and strong penetration into the physical characteristics of the Earth's surface. SAR image change detection aims to identify and extract areas of change on the Earth's surface by analyzing two or more SAR images of the same geographical area acquired at different times, and is one of the core technologies of remote sensing applications. However, due to the influence of the SAR imaging mechanism, multiplicative speckle noise is unavoidable in the images, and the SAR is quite sensitive to imaging geometric errors, making it difficult to achieve high-precision change detection.

[0003] Despite the numerous existing SAR image change detection methods, they still have significant limitations when dealing with complex scenes and strong noise interference: First, pixel-level difference map-based methods mainly rely on grayscale differences, making it difficult to physically distinguish high-intensity speckle noise from real targets, resulting in extremely high false alarm rates. Second, feature extraction methods based on difference maps generally follow a "single-stream serial dependency" architecture, where feature extraction depends entirely on the quality of the initial difference map. Once the initial stage is contaminated by noise, the error will be cascaded and amplified irreversibly in subsequent steps, leading to serious missed detections or false alarms. In addition, existing technologies are sensitive to registration errors and are prone to edge false alarms. Furthermore, due to the lack of holistic perception, "detection holes" often appear inside complete targets, seriously affecting the completeness of detection. Summary of the Invention

[0004] The technical problem to be solved by this invention is how to provide a SAR image change detection method that can significantly improve the robustness of the algorithm, greatly reduce the false alarm rate of change detection, and ensure that the change area is full and the outline is complete.

[0005] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is: a SAR image change detection method that fuses gravitational field and difference map, the method comprising the following steps:

[0006] S1: Acquire dual-temporal SAR images and Bilateral filtering is used to preprocess the image to obtain the processed SAR image. and ;

[0007] S2: Map SAR image pixel grayscale to physical mass, and calculate the gravitational field magnitude based on cubic range attenuation factors. Extract candidate regions of the structure;

[0008] S3: Calculate the local gradient direction consistency weight and correlate it with the gravitational field magnitude. Nonlinear coupling utilizes geometric anisotropy to eliminate spurious corner noise in the candidate structural regions, resulting in a high-confidence structural feature map. ;

[0009] S4: High-confidence structural feature map from the two-phase description Key points are extracted, and matching is performed based on adaptive matching radius to identify variant structural feature points. A physical structure mask is then generated using multi-scale morphological operations. ;

[0010] S5: Utilizing preprocessed SAR images and Calculate the log-ratio difference map, and generate an enhanced difference map by modulating it with a local texture statistical difference factor. And construct a statistical significance mask. ;

[0011] S6: Utilizing the physical structure mask and statistical significance mask As a dual-gating factor, the enhanced difference map is multiplicatively fused to obtain a fused difference map. ;

[0012] S7: An adaptive threshold segmentation algorithm with closed-loop feedback is used to segment the fused difference map to obtain change detection results.

[0013] The beneficial effects of adopting the above technical solution are as follows: 1) In the existing technology, feature extraction is heavily dependent on the initial difference map. Once the source is contaminated by noise, subsequent errors will be amplified and irreversible. The "physical-statistical" dual-stream parallel decoupled architecture constructed in this invention extracts features directly from the original image through the physical structure stream, which is independent of the difference map generation process. By utilizing the heterogeneous complementary characteristics of physical features and statistical features, the cascading propagation path of noise is cut off from the source. Thus, even when the signal-to-noise ratio of the initial difference map is extremely low, the detection reliability can still be maintained, which significantly improves the robustness of the algorithm.

[0014] 2) To address the technical challenge of distinguishing strong reflective speckle noise from real building corners in SAR images in terms of grayscale intensity, leading to a high false alarm rate, this invention introduces a modified higher-order gravitational field model and a geometric anisotropy screening mechanism. By modifying the gravitational distance attenuation factor to the cubic level to cut off background interference, and utilizing the statistical difference between the neat edges of real structures and the chaotic directions of noise points for nonlinear coupling screening, high-confidence physical structural features are captured. This effectively eliminates high-intensity pseudo-corner noise at the physical level, significantly reducing the false alarm rate of change detection.

[0015] 3) To address the issues of existing technologies being sensitive to SAR image registration errors, which can easily lead to false edge detections, and the fragmentation of detection results causing holes inside the target, the present invention designs a "point-field" upscaling mapping strategy. This strategy dynamically calculates the adaptive matching radius based on the density of image corner points and uses multi-scale morphological dilation to restore discrete variable feature points into a continuous physical structure mask, forming a "point-to-surface" reconstruction mechanism. This effectively fills in detection holes while adapting to nonlinear registration errors, thereby ensuring that the interior of the variable region is full and the outline is complete. Attached Figure Description

[0016] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0017] Figure 1 This is the main flowchart of the method described in an embodiment of the present invention;

[0018] Figure 2 This is a set of experimental simulation images from an embodiment of the present invention. These are SAR images of the Sulzberger region taken at different times, with a size of 256×256. Figure 2 Part a of the film was shot on March 11, 2021. Figure 2 Part b of the film was shot on March 16, 2021;

[0019] Figure 3 This is a diagram showing the results of SAR image change detection using the method described in this embodiment of the invention on the Sulzberger dataset;

[0020] Figure 4 This is a real result image of SAR image change detection in the Sulzberger dataset in an embodiment of the present invention;

[0021] Figure 5 This is a set of experimental simulation images from an embodiment of the present invention. These are SAR images of the Ottawa region taken at different times, each with a size of 290×350. Figure 5 Part a of the film was shot in May 1997. Figure 5 Part b of the film was shot in August 1997;

[0022] Figure 6 This is a diagram showing the results of SAR image change detection on the Ottawa dataset using the method described in this embodiment of the invention.

[0023] Figure 7 This is a real result image of SAR image change detection in the Ottawa dataset in an embodiment of the present invention. Detailed Implementation

[0024] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0025] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0026] Overall, such as Figure 1 As shown in the figure, this invention discloses a SAR image change detection method that integrates gravitational field and difference map. This method constructs a "physical-statistical" dual-stream parallel architecture, achieving high-precision change detection through the collaborative work of the physical structure stream and the statistical difference stream. Specifically, it includes the following steps:

[0027] S1: Acquire dual-temporal SAR images and Bilateral filtering is used to preprocess the image to obtain the processed SAR image. and ;

[0028] S2: Map pixel grayscale values ​​of dual-temporal SAR images to physical masses, and calculate the gravitational field magnitude based on cubic range attenuation factors. Extract candidate regions of the structure;

[0029] S3: Calculate the local gradient direction consistency weight and correlate it with the gravitational field magnitude. Nonlinear coupling utilizes geometric anisotropy to eliminate spurious corner noise in the candidate structural regions, resulting in a high-confidence structural feature map. ;

[0030] S4: High-confidence structural feature map from the two-phase description Key points are extracted, and matching is performed based on adaptive matching radius to identify variant structural feature points. A physical structure mask is then generated using multi-scale morphological operations. ;

[0031] S5: Utilizing preprocessed dual-temporal SAR images and Calculate the log-ratio difference map, and generate an enhanced difference map by modulating it with a local texture statistical difference factor. And construct a statistical significance mask. ;

[0032] S6: Utilizing the physical structure mask and statistical significance mask As a dual-gating factor, the enhanced difference map is multiplicatively fused to obtain a fused difference map. ;

[0033] S7: An adaptive threshold segmentation algorithm with closed-loop feedback is used to segment the fused difference map to obtain change detection results.

[0034] The above steps will be explained in detail below:

[0035] Acquire the first temporal SAR image of the same geographic area Second-phase SAR images ,like Figure 2 and Figure 5 As shown, bilateral filtering is used to perform the same preprocessing on the two temporal SAR images.

[0036] For images Each pixel in First, calculate the weight components and normalization coefficients required for filtering. The specific process is as follows:

[0037] (1) Calculate the Gaussian kernel for spatial proximity This kernel function is used to assign weights based on the geometric distance between pixels. Any neighboring pixel within the centered local window Calculate its spatial weight The calculation formula is:

[0038] ;

[0039] in, For pixels and The Euclidean distance between them; The standard deviation of the spatial domain controls the decay rate of the spatial weights; the closer the distance, the greater the weight.

[0040] (2) Calculate the Gaussian kernel for pixel value similarity. This kernel function is used to assign weights based on the grayscale differences between pixels and to calculate the neighboring pixels. Relative to the center pixel Gray weight The calculation formula is:

[0041] ;

[0042] in, This represents the absolute value of the grayscale difference between two pixels. The standard deviation of the value range is used to ensure that in the edge region (where the gray level difference is large), the weight rapidly decays to close to 0, thereby avoiding cross-edge smoothing and achieving the edge preservation effect.

[0043] (3) Calculate the normalization coefficient To ensure energy conservation in the filtered image, it is necessary to sum all weights within the window. The calculation formula is:

[0044] ;

[0045] (4) Calculate the final filtered and denoised image. Based on the weighted components and normalization coefficients calculated above, perform a weighted summation on all pixels within the window to obtain the image for each pixel. Filtered image The calculation formula is:

[0046] ;

[0047] in, Indicates a local window All neighboring pixels within Perform iterative summation; through the above processing, the denoised first phase images are obtained respectively. and the denoised second phase image .

[0048] By using the joint weighted calculation of spatial proximity and pixel value similarity, Gaussian smoothing can be used to effectively suppress multiplicative speckle noise in SAR images in flat areas. At the same time, in textured areas such as building edges, the Gaussian kernel of pixel value similarity decays rapidly with increasing grayscale difference, effectively preventing the smoothing process from crossing physical boundaries.

[0049] S2: Construct a high-order gravitational field model, extract candidate regions of the structure, reconstruct the feature extraction process into a resultant force calculation process in the physical field, and use the low-pass characteristic of integral transform to suppress noise.

[0050] (1) Image grayscale-physical mass mapping: First, establish the mapping relationship between digital image space and physical gravitational field space. Each pixel in the SAR image is defined as a particle with physical mass.

[0051] To eliminate the magnitude differences in imaging data from different sensors, this invention employs a linear normalization strategy to normalize the preprocessed image. median coordinate pixel grayscale value at Mapped to relative physical mass in the interval [0,1] The calculation formula is as follows:

[0052] ;

[0053] in, This is the original grayscale value of that pixel in the image; This represents the minimum grayscale value in the entire image. This represents the maximum grayscale value in the entire image.

[0054] Through this mapping, the relative mass of the brightest strong scatterer in the image is determined. Approaching 1, it can generate strong gravitational potential energy; while the darkest background region has a relative mass... When the value approaches 0, its effect on the gravitational field is negligible.

[0055] (2) Constructing a short-range strong interaction gravitational field: In view of the problem that the distance decay factor in the traditional Newton's law of universal gravitation is a "long-range force", which makes it easy to be affected by distant background clutter when calculating the gravity of the central pixel, the gravity model of this invention is modified to a high-order cubic distance decay factor to construct an interference-resistant short-range strong interaction field.

[0056] For any pixel in the image In its defined local neighborhood Calculate the components of the net gravitational force acting on it. and The calculation formula is as follows:

[0057] ;

[0058] ;

[0059] in, The center pixel of the current calculation has the following coordinates: ; Center pixel For any neighboring pixel within the neighborhood, its coordinates are ; The signal strength of the connected region is enhanced by the "cumulative effect" of mass, which is the product of the relative mass of the center pixel and the neighboring pixels. As a cubic distance attenuation term, it directly cuts off the interference of distant background noise on the center pixel at the physical level; It is the gravitational constant, used to adjust the numerical order of magnitude;

[0060] (3) Calculate the gravitational field modulus and extract candidate structural regions: Perform vector synthesis on the gravitational components in the two orthogonal directions mentioned above, and calculate the pixel value. gravitational field modulus at the location The calculation formula is as follows:

[0061] ;

[0062] Gravitational field modulus The physical saliency of pixels was quantitatively characterized, and isolated noise lacking neighborhood support was automatically filtered out using the physical field integral and cubic attenuation characteristics. Simultaneously, the cumulative gain generated by consecutive pixels within the neighborhood significantly enhanced the signal strength of real structural feature points. The final result is a candidate structure image with high signal-to-noise ratio and low background clutter.

[0063] S3: Robust feature selection based on geometric anisotropy introduces non-physical "geometric anisotropy" as a second constraint, and uses the difference between the neatness of the real edge direction and the disorder of the noise direction for selection.

[0064] First, the gradient direction angle of each pixel is calculated using the Sobel operator. Then, the standard deviation of the gradient direction within the local window is calculated. Based on this, a directional consistency weight is constructed. :

[0065] ;

[0066] in, This is the weight for directional consistency; a larger value indicates more uniformity in direction. This represents the statistical standard deviation of the gradient directions within a local window, reflecting the degree of dispersion of the directions; when the gradient directions in the neighborhood are disordered, the standard deviation... A large value leads to the consistency weight approaching 0; conversely, when true edges exist, the weight approaches 1.

[0067] Finally, the physical gravitational magnitude is nonlinearly coupled with the geometric direction weights to generate a high-confidence structural feature map. This feature map is available at any pixel. Response value at The calculation formula is as follows:

[0068] ;

[0069] in, The high-confidence structural feature response value under dual constraints; The gravitational field modulus represents the physical structural strength of the pixel. This is the weight for directional consistency; a larger value indicates more uniformity in direction.

[0070] Based on the statistical difference that real structural edges are anisotropic while speckle noise is isotropic, a direction consistency weight is introduced. As a criterion for distinguishing authenticity: even if the noise intensity is high, as long as its direction is disordered, it can still be identified. It will approach 0, thus in the calculation These are then removed. By traversing the entire image and calculating, a high-confidence structural feature map with pseudo-corner noise removed is finally obtained. .

[0071] S4: Variable structure identification and point-field upscaling mapping, an upscaling mapping strategy from "zero-dimensional feature points" to "two-dimensional structural fields".

[0072] (1) Dynamically calculate the adaptive matching radius: In order to solve the matching problem caused by uneven distribution of feature points, the matching radius is calculated based on the corner density of the image. Dynamically adjust matching radius The calculation formula is as follows:

[0073] ;

[0074] in, It is the ratio of the number of feature points to the total number of pixels in the image; The truncation function restricts the radius to between [3, 8] pixels. An adaptive mechanism is used to reduce the radius in feature-dense areas to prevent mismatches and increase the radius in sparse areas to compensate for registration errors.

[0075] (2) Extracting key points and performing bidirectional non-matching search: First, extract local maxima points from the structural response maps of the first and second time phases respectively, and construct a set of key points. and Subsequently, using A bidirectional nearest neighbor search is performed to identify unmatched points as variant structure points. The calculation formula is as follows:

[0076] ;

[0077] in, This is the final set of identified variant structure points; This is the set of key local maxima extracted from the first phase structure feature map; This is the set of local maxima key points extracted from the second phase structure feature map; Key point With key points The Euclidean geometric distance between them; For point to set The minimum distance among all points; For adaptive matching radius;

[0078] (3) Constructing a multi-scale physical structure mask: targeting Each mutation point in This invention simulates the radiation attenuation characteristics of physical fields, using radius sequences respectively. By performing multi-layer morphological expansion and assigning different weights, a continuous two-dimensional structural field is constructed, the calculation formula of which is as follows:

[0079] ;

[0080] in, Let be any pixel coordinate in the image; Attenuation weights for layers with different radii; It is an indicator function;

[0081] Through the aforementioned bidirectional matching logic, the key points where substantial changes occurred were precisely located; subsequently, through multi-scale weighted expansion, the "zero-dimensional discrete mutation points" were mapped to a "two-dimensional continuous structural field," effectively filling the detection gaps of the target and fully restoring the physical outline of the mutated target.

[0082] S5: Enhanced statistical difference map generation, employing an enhancement strategy of "macro-statistics guiding micro-pixels".

[0083] (1) Calculation of the basic logarithmic ratio difference map: Considering the multiplicative noise characteristics of SAR images, the logarithmic ratio operator is used to convert multiplicative noise into additive noise. To prevent division by zero errors, a small constant is introduced. The calculation formula is as follows:

[0084] ;

[0085] in, and The two phase images are in pixels. grayscale value at that location To prevent tiny constants with a denominator of 0.

[0086] (2) Constructing local statistical difference factors: Calculating pixels using a sliding window Local mean within the neighborhood and local standard deviation To eliminate the influence of the magnitude of brightness, the normalized mean differences were calculated separately. and variance differences The calculation formula is as follows:

[0087] ;

[0088] (3) Generate an enhanced difference map: Multiplicatively enhance the pixel-level logarithmic ratio using a local statistical difference factor. The calculation formula is as follows:

[0089] ;

[0090] in, To enhance the coefficients, local texture changes are used as gain factors. The signal is amplified only when the gray value and texture change drastically at the same time, which significantly enhances the detectability of weak texture areas and improves the algorithm's ability to perceive unstructured changes.

[0091] (4) Generating a statistical significance mask: To convert the unbounded difference graph into a probability mask suitable for gating mechanisms, an adaptive Sigmoid mapping based on global statistical properties is adopted. First, the statistical significance mask is calculated. global mean and standard deviation Determine the adaptive statistical threshold The calculation formula is as follows:

[0092] ;

[0093] in, This is the statistical sensitivity coefficient;

[0094] Using the Sigmoid function to generate a statistical significance mask The calculation formula is as follows:

[0095] ;

[0096] in, The slope factor is used to map the difference map as a significance probability field in the interval [0, 1], providing a statistical flow gating signal for the multiplicative fusion in the subsequent step S6.

[0097] S6: Physical-Statistical Dual-Stream Adaptive Multiplicative Fusion, based on a "Dual-Gated" Nonlinear Multiplicative Modulation Fusion Model.

[0098] Using physical structure mask and statistical significance mask The difference map is modulated, and its calculation formula is as follows:

[0099] ;

[0100] in, This is the final difference diagram after fusion; , The gain coefficient is used for modulation; a dual-gated multiplicative mechanism is adopted: in the background region, the gain term is 1 to maintain a low noise level; in the variable region, the gain term is significantly greater than 1 to amplify the target signal in a directional manner; this fusion strategy effectively solves the contradiction between "edge blurring" and "detection hole", and achieves extreme suppression of the background and high-fidelity enhancement of the target.

[0101] S7: Adaptive segmentation with closed-loop feedback introduces a closed-loop feedback control loop based on change priors to generate the final change detection result.

[0102] (1) Initial threshold calculation First, Otsu's fusion difference map is used. Perform global calculations to obtain the initial segmentation threshold. .

[0103] (2) Closed-loop feedback detection: Set the current threshold The system automatically calculates the percentage of pixels that change at this threshold. The calculation formula is as follows:

[0104] ;

[0105] in, This is the final difference diagram after fusion; and These are the height and width of the image, respectively; This is an indicator function; it is 1 if the condition is met, and 0 otherwise.

[0106] (3) Threshold adaptive iteration: Establish feedback control logic and preset the safety change rate threshold. .

[0107] like Determine if the current threshold is valid, and set the final threshold. ;

[0108] like If a false alarm is detected, meaning a large amount of background data is misinterpreted as change, the algorithm will automatically adjust the step size accordingly. Iteratively increase the threshold, i.e. And re-execute step 2 calculation. The optimal threshold is obtained when the rate of change returns to a reasonable range. .

[0109] (4) Generate the final change detection result: using the finally determined optimal threshold The fused difference map is binarized and segmented, and a 3×3 morphological opening operation is performed to remove isolated noise points, resulting in the final change detection result map. The calculation formula is as follows:

[0110] ;

[0111] in, For pixels The final classification result is 1, which represents change, and 0, which represents no change. This represents the morphological opening operation.

[0112] The effectiveness of the present invention will be further illustrated by the following SAR image experiments.

[0113] 1) Experimental setup:

[0114] To verify the effectiveness of the method of this invention, two representative datasets were selected for experiments:

[0115] Sulzberger dataset: such as Figure 2 As shown, the image, taken in March 2021, reflects the collapse of the Sulzberger ice shelf and contains a wealth of texture details, used to verify the algorithm's ability to suppress strong background clutter.

[0116] Ottawa dataset: such as Figure 5 As shown, the images, taken in May and August 1997, depict the flooding in the Ottawa area against a backdrop of complex urban architecture. They were used to verify the algorithm's ability to remove strong reflection noise. The performance of the proposed method was compared with that of the logarithmic ratio method and the local mean method in the experiment.

[0117] 2) Results Analysis:

[0118] This experiment uses false positives (FP), false negatives (FN), correct classification rate (PCC), and Kappa coefficient as quantitative evaluation indicators. A lower FP indicates fewer false alarms, and a higher Kappa coefficient indicates greater consistency between the detection results and the actual changes.

[0119] The present invention was tested on the two sets of real SAR image datasets mentioned above, and the resulting change detection maps are shown below. Figure 3 and Figure 6 As shown, the obtained change detection and evaluation indicators are shown in Table 1 and Table 2.

[0120] Table 1: Comparison of Change Detection Results in the Sulzberger Dataset

[0121]

[0122] Table 2: Comparison Results of Change Detection in the Ottawa Dataset

[0123]

[0124] As shown in Table 1, in the Sulzberger dataset, the false alarm rate (FP) of this invention is only 529, far lower than the 1439 of the log-ratio method and the 3596 of the local mean method. This indicates that although existing technologies can preserve change information, they introduce extremely high background false alarm noise. In contrast, this invention, through the construction of a high-order gravitational field model, utilizes the short-range strong interactions of physical fields to effectively cut off background clutter interference, thereby significantly reducing the false alarm rate. Furthermore, the Kappa coefficient of this invention reaches 94.12%, the highest among the three, demonstrating the effectiveness of the "physical-statistical dual-stream parallel architecture." Figure 3 The detection results of the method of the present invention are shown, and... Figure 4 The actual change diagram is highly consistent, with clear edges.

[0125] As shown in Table 2, the Ottawa dataset contains complex urban architectural backgrounds, making it highly susceptible to false alarms. The traditional log-ratio method and local mean method both have false alarms (FP) exceeding 2000, indicating that traditional operators cannot effectively distinguish between speckle noise and real flood changes. The method of this invention controls the false alarm number to 951, significantly outperforming the comparison methods and demonstrating the algorithm's excellent false alarm suppression capability in complex urban backgrounds. Although the local mean method has a lower number of missed detections, it comes at the cost of a large number of false alarms. This invention achieves a Kappa coefficient of 90.87% while maintaining a low false alarm rate. Figure 6 The detection results show that this invention, by utilizing the "point-field upscaling mapping" strategy, effectively fills the detection gaps of the target, making the detection results more comprehensive and consistent with... Figure 7 The actual changes in land features are highly consistent.

Claims

1. A method for detecting changes in SAR images by fusing gravitational field and difference maps, characterized in that, The method includes the following steps: S1: Acquire dual-temporal SAR images and Bilateral filtering was used to preprocess the image, resulting in a processed dual-temporal SAR image. and ; S2: Map pixel grayscale values ​​of dual-temporal SAR images to physical masses, and calculate the gravitational field magnitude based on cubic range attenuation factors. Extract candidate regions of the structure; S3: Calculate the local gradient direction consistency weight and correlate it with the gravitational field magnitude. Nonlinear coupling is performed, and geometric anisotropy is used to remove spurious corner noise in the candidate region of the structure, resulting in a high-confidence structural feature map. ; S4: High-confidence structural feature map from the two-phase description Key points are extracted, and matching is performed based on adaptive matching radius to identify variant structural feature points. A physical structure mask is then generated using multi-scale morphological operations. ; S5: Utilizing preprocessed dual-temporal SAR images and Calculate the log-ratio difference map, and generate an enhanced difference map by modulating it with a local texture statistical difference factor. And construct a statistical significance mask. ; S6: Utilizing the physical structure mask and statistical significance mask As a dual-gating factor, the enhanced difference map is multiplicatively fused to obtain a fused difference map. ; S7: An adaptive threshold segmentation algorithm with closed-loop feedback is used to fuse the difference map. The change detection results are obtained by segmentation.

2. The SAR image change detection method fused with gravitational field and difference map according to claim 1, characterized in that, The step described in S2, which maps SAR image pixel grayscale to physical mass and calculates the gravitational field magnitude, specifically includes the following steps: (1) Image grayscale-physical mass mapping: First, establish the mapping relationship between digital image space and physical gravitational field space. Each pixel in the SAR image is defined as a particle with physical mass. A linear normalization strategy is used to preprocess the image. median coordinate pixel grayscale value at Mapped to relative physical mass in the interval [0,1] The calculation formula is as follows: ； in, This is the original grayscale value of that pixel in the image; This represents the minimum grayscale value in the entire image. This represents the maximum grayscale value in the entire image. Through this mapping, the relative mass of the brightest strong scatterer in the image is determined. Approaching 1, while the darkest background areas have the highest relative quality. Approaching 0; (2) Construct a short-range strong interaction gravitational field: The gravitational model is modified to a higher order, and a cubic distance decay factor is used to construct an interference-resistant short-range strong interaction field. For any pixel in the image In its defined local neighborhood Calculate the components of the net gravitational force acting on it. and The calculation formula is as follows: ； ； in, The center pixel of the current calculation has the following coordinates: ; Center pixel For any neighboring pixel within the neighborhood, its coordinates are ; It is the product of the relative quality of the center pixel and its neighboring pixels; This is a cubic distance attenuation term; It is the gravitational constant, used to adjust the numerical order of magnitude; (3) Calculate the gravitational field modulus and extract candidate structural regions: Perform vector synthesis on the gravitational components in the two orthogonal directions mentioned above, and calculate the pixel value. gravitational field modulus at the location The calculation formula is as follows: 。 3. The SAR image change detection method fusion of gravitational field and difference map according to claim 1, characterized in that, The calculation of local gradient direction consistency weights and nonlinear coupling described in S3 specifically includes the following steps: First, the gradient direction angle of each pixel is calculated using the Sobel operator. Then, the standard deviation of the gradient direction within the local window is calculated. Based on this, a directional consistency weight is constructed. : ； in, This is the weight for directional consistency; a larger value indicates more uniformity in direction. This represents the statistical standard deviation of the gradient direction within the local window, reflecting the degree of dispersion of the direction; Finally, the physical gravitational magnitude is nonlinearly coupled with the geometric direction weights to generate a high-confidence structural feature map. This feature map is present in any pixel. Response value at The calculation formula is as follows: ； in, The high-confidence structural feature response value under dual constraints; The gravitational field modulus represents the physical structural strength of the pixel. This is the weight for directional consistency; a larger value indicates more uniformity in direction. Based on the statistical difference that real structural edges are anisotropic while speckle noise is isotropic, a direction consistency weight is introduced. As a criterion for distinguishing authenticity: even if the noise intensity is high, as long as its direction is disordered, it can still be identified. It will approach 0, thus in the calculation These are removed; by traversing the entire image for calculation, a high-confidence structural feature map with pseudo-corner noise removed is finally obtained. .

4. The SAR image change detection method fusion of gravitational field and difference map according to claim 1, characterized in that, S4 specifically includes the following steps: (1) Dynamically calculate the adaptive matching radius: based on the image corner density Dynamically adjust matching radius The calculation formula is as follows: ； in, It is the ratio of the number of feature points to the total number of pixels in the image; This is a truncation function that restricts the radius to between [3, 8] pixels; (2) Extracting key points and performing bidirectional non-matching search: First, extract local maxima points from the structural response maps of the first and second time phases respectively, and construct a set of key points. and Subsequently, using A bidirectional nearest neighbor search is performed to identify unmatched points as variant structure points. The calculation formula is as follows: ； in, This is the final set of identified variant structure points; This is the set of key local maxima extracted from the first phase structure feature map; This is the set of local maxima key points extracted from the second phase structure feature map; Key point With key points The Euclidean geometric distance between them; For point to the set The minimum distance among all points; For adaptive matching radius; (3) Constructing a multi-scale physical structure mask: targeting Each mutation point in The radiation attenuation characteristics of the simulated physical field are represented by radius sequences. By performing multi-layer morphological expansion and assigning different weights, a continuous two-dimensional structural field is constructed, the calculation formula of which is as follows: ； in, Let be any pixel coordinate in the image; Attenuation weights for layers with different radii; This is an indicator function.

5. The SAR image change detection method fusion of gravitational field and difference map according to claim 1, characterized in that, The method described in S5 for generating an enhanced difference map using local texture statistical difference factors specifically includes the following steps: First, the fundamental log-ratio difference plot is calculated, and the multiplicative noise is converted into additive noise using the log-ratio operator; to prevent division by zero errors, a small constant is introduced. The calculation formula is as follows: ； in, and The two phase images are in pixels. grayscale value at that location To prevent small constants with a denominator of 0; Secondly, a local statistical difference factor is constructed, and a sliding window is used to calculate the pixel value. Local mean within the neighborhood and local standard deviation Calculate the normalized mean difference separately. and variance differences The calculation formula is as follows: ； Next, an enhanced difference map is generated, and the pixel-level logarithmic ratio is multiplicatively enhanced using a local statistical difference factor. The calculation formula is as follows: ； in, For enhancement coefficient; Finally, a statistical significance mask is generated using an adaptive Sigmoid mapping based on global statistical properties; first, the statistical significance mask is calculated. global mean and standard deviation Determine the adaptive statistical threshold The calculation formula is as follows: ； in, This is the statistical sensitivity coefficient; Using the Sigmoid function to generate a statistical significance mask The calculation formula is as follows: ； in, is the mapping slope factor; the difference map is mapped to a significance probability field in the interval [0, 1].

6. The SAR image change detection method fusion of gravitational field and difference map according to claim 1, characterized in that, The multiplicative fusion using a dual-gated factor described in S6 specifically includes the following steps: Based on a dual-gated nonlinear multiplicative modulation fusion model, a physical structure mask is used. and statistical significance mask The difference map is modulated, and its calculation formula is as follows: ； in, This is the final difference diagram after fusion; , The modulation gain coefficient is used; a dual-gated multiplicative mechanism is employed: in the background region, the gain term is 1 to maintain a low noise level; in the variable region, the gain term is significantly greater than 1 to amplify the target signal in a directional manner.

7. The SAR image change detection method fusing gravitational field and difference map according to claim 1, characterized in that, The adaptive threshold segmentation algorithm with closed-loop feedback described in S7 specifically includes the following steps: Based on a change-priority closed-loop feedback control loop, the final change detection result is generated. First, calculate the initial threshold and use Otsu's fusion difference map. Perform global calculations to obtain the initial segmentation threshold. ; Secondly, closed-loop feedback detection is used to set the current threshold. The system automatically calculates the percentage of pixels that change at this threshold. The calculation formula is as follows: ； in, This is the final difference diagram after fusion; and These are the height and width of the image, respectively; This is an indicator function; it is 1 if the condition is met, and 0 otherwise. Finally, an adaptive threshold iteration is implemented to establish feedback control logic and preset a safe rate of change threshold. ;like Determine if the current threshold is valid, and set the final threshold. ;like If a false alarm is detected, meaning a large amount of background data is misinterpreted as change, the algorithm will automatically adjust the step size accordingly. Iteratively increase the threshold, i.e. And re-execute step 2 calculation. The optimal threshold is obtained when the rate of change returns to the set range. ; The final change detection result is generated using the finally determined optimal threshold. The fused difference map is binarized and segmented, and a 3×3 morphological opening operation is performed to remove isolated noise points, resulting in the final change detection result map. The calculation formula is as follows: ； in, For pixels The final classification result is 1, which represents change, and 0, which represents no change. This represents the morphological opening operation.

Images