Change detection method based on multi-source remote sensing data

Abstract

The invention relates to a change detection method based on multi-source remote sensing data. The method comprises the following steps: obtaining a plurality of radar images according to a ground-based radar image; based on the radar images before and after the change, registering the radar images with the radar images at the corresponding moments; after registration, the radar image and the optical image are fused to obtain a fusion difference chart, and a change detection result is obtained through segmentation and clustering. According to the embodiment of the invention, effective data reference is provided for comprehensive monitoring and analysis of the target area, so that accurate positioning and qualitative monitoring of the change area are realized.

Description

Technical Field

[0001] The present disclosure relates to the field of remote sensing data for geological disaster monitoring, and specifically to a change detection method based on multi-source remote sensing data. Background Art

[0002] In the prior art, there are deficiencies in the information volume of geological disaster monitoring based on single remote sensing data, resulting in the monitoring results not meeting the actual application requirements. Summary of the Invention

[0003] The present disclosure aims to provide a change detection method based on multi-source remote sensing data, providing effective data reference for the comprehensive monitoring and analysis of the target area, so as to achieve the accurate positioning and qualitative monitoring of the changed area.

[0004] According to one of the solutions of the present disclosure, a change detection method based on multi-source remote sensing data is provided, including:

[0005] Obtain a complex radar image according to the ground-based radar image;

[0006] Based on the radar images before and after the change, register them with the radar images at the corresponding time;

[0007] After registration, fuse the radar image and the optical image to obtain a fused difference map, and obtain the change detection result through segmentation and clustering.

[0008] In some embodiments, obtaining a complex radar image according to the ground-based radar image includes:

[0009] Obtain the original radar data of the monitored slope area;

[0010] In one sub-cycle, the frequency band of the output signal of each channel remains unchanged, and in the next sub-cycle, the frequency bands of the output signals of all channels are rotated;

[0011] Perform two-dimensional imaging processing on the original radar data respectively;

[0012] Perform image registration and atmospheric phase correction in sequence to obtain the radar complex image reflecting the monitored area in the nth observation cycle.

[0013] In some embodiments, registering the radar images before and after the change with the radar images at the corresponding time includes:

[0014] Select any two radar images before the change and the radar image after the change obtained at any time;

[0015] Perform denoising processing on the two radar images;

[0016] Obtain the optical image at the corresponding time;

[0017] Perform radiometric correction on the optical image;

[0018] Register the corrected optical image with the radar image at the corresponding time.

[0019] In some embodiments, registering the corrected optical image with the radar image at the corresponding time includes:

[0020] Select feature points in the optical image and select corresponding feature points in the radar image as a pair of feature points;

[0021] Obtain a transformation matrix for perspective transformation according to the pair of feature points;

[0022] Transform the optical image into the radar image coordinate system to obtain the final image matching result;

[0023] Obtain the interpolated result according to the information of the ground-based radar image pixels and the optical image pixels.

[0024] In some embodiments, after registration, fusing the radar image and the optical image to obtain a fusion difference map includes:

[0025] Obtain a fusion difference map based on the neighborhood mean ratio difference map and the improved relative entropy difference map method;

[0026] Perform filtering processing to obtain the final fusion difference map.

[0027] In some embodiments, obtaining a fusion difference map based on the neighborhood mean ratio difference map and the improved relative entropy difference map method includes:

[0028] Apply a mean logarithmic ratio operator to the fused image to generate a mean logarithmic ratio map;

[0029] Apply the improved relative entropy to the fused image to generate an improved neighborhood mean relative entropy image;

[0030] Generate two fusion difference maps using the NSCT transform;

[0031] Fuse the difference maps by the local energy method to obtain the final fusion difference map.

[0032] In some embodiments, generating two fusion difference maps using the NSCT transform includes:

[0033] Use the generated mean logarithmic ratio map and neighborhood mean relative entropy image as the input images for NSCT transform processing;

[0034] Use the NSCT transform to obtain the low-frequency information and high-frequency information in the mean logarithmic ratio map and neighborhood mean relative entropy image;

[0035] Fuse the transformed low-frequency coefficients to obtain the fused low-frequency coefficients;

[0036] Perform NSCT inverse transformation on the fused low-frequency coefficients and the high-frequency coefficients respectively to obtain the fused difference map.

[0037] In some embodiments, wherein, perform filtering processing to obtain the final fused difference map, including:

[0038] Use bilateral filtering for edge preservation and further denoising to obtain the final difference map.

[0039] In some embodiments, wherein, obtain the change detection result through segmentation and clustering, including:

[0040] Convert the bilaterally filtered image into a column vector;

[0041] Initialize the membership degree to meet the preset conditions;

[0042] Calculate the cluster center;

[0043] Calculate the objective function;

[0044] Update the membership degree matrix;

[0045] Iterate cyclically until the objective function is less than the set threshold;

[0046] Compare the magnitudes of the membership degrees of each pixel, and perform class labeling on the transformed image based on the two types of membership degrees corresponding to each pixel.

[0047] In some embodiments, wherein, obtaining the change detection result through segmentation and clustering further includes:

[0048] Normalize the bilaterally filtered fused difference map;

[0049] Calculate the class mean, class variance, and class mean under the initial label distribution;

[0050] Determine the optimal threshold;

[0051] Compare the neighborhood features with the optimal threshold;

[0052] Calculate the non-local means after introducing the optimal threshold

[0053] Introduce the non-local means of the optimal threshold into the likelihood function;

[0054] Iterate cyclically until the likelihood energy function is minimized, and output the change detection result.

[0055] The change detection method based on multi-source remote sensing data according to various embodiments of the present disclosure obtains complex radar images at least based on ground-based radar images; based on the radar images before and after the change, register with the radar images at the corresponding time; after registration, fuse the radar images with optical images to obtain a fused difference map, and obtain the change detection result through segmentation and clustering. Thus, the research on the change detection method using multi-source remote sensing images can make the optical images and ground-based radar images complement each other well, making up for the deficiencies of single remote sensing images. Further, various embodiments of the present disclosure propose an optimal criterion non-local mean ground-based radar image segmentation method for the multi-source fused difference map. On the basis of effectively protecting image details and reducing speckle noise, an non-local mean algorithm based on an optimal threshold is proposed, which can maintain the fine structure of the image and enhance the resistance to speckle noise, further improving the accuracy of change detection.

[0056] It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the claimed present disclosure. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] In the drawings that are not necessarily drawn to scale, like reference numerals in different views may represent like components. Like reference numerals with letter suffixes or like reference numerals with different letter suffixes may represent different instances of like components. The drawings generally illustrate various embodiments by way of example and not limitation, and are used in conjunction with the description and the claims to explain the disclosed embodiments.

[0058] Figure 1 The flowchart of the change detection method based on multi-source remote sensing data according to an embodiment of the present disclosure is shown;

[0059] Figure 2 The schematic diagram of data acquisition of the ground-based radar system according to an embodiment of the present disclosure is shown;

[0060] Figure 3 The schematic flowchart of preprocessing of radar images and optical images according to an embodiment of the present disclosure is shown;

[0061] Figure 4 The schematic flowchart of change detection of the fused image according to an embodiment of the present disclosure is shown;

[0062] Figure 5 The flowchart of NSCT transform processing according to an embodiment of the present disclosure is shown;

[0063] Figure 6 The flowchart of fuzzy C-means clustering according to an embodiment of the present disclosure is shown;

[0064] Figure 7Fig. 0 shows a schematic flow chart of the ground-based radar image segmentation based on the Fisher optimal criterion non-local mean MAP-MRF according to an embodiment of the present disclosure. Detailed implementation manners

[0065] In order to make the objectives, technical solutions, and advantages of the embodiments of the present disclosure clearer, the technical solutions of the embodiments of the present disclosure will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present disclosure. Apparently, the described embodiments are some, but not all, of the embodiments of the present disclosure. All other embodiments obtained by those of ordinary skill in the art based on the described embodiments of the present disclosure without creative efforts shall fall within the scope of protection of the present disclosure.

[0066] In the case of complex geological structures and large fluctuations in terrain and landforms, there are material conditions that are extremely prone to geological disasters such as landslides, collapses, and debris flows; in addition, production and living activities are also very likely to induce the occurrence of landslide disasters, and the instability or sliding of a large number of rock and soil slopes are common, bringing serious hazards and losses to engineering construction and property safety. The monitoring and analysis of the disaster-stricken areas are of great significance for safe production. Obtaining disaster information from remote sensing data is a very important means of disaster research.

[0067] The ground-based radar technology has characteristics such as minute-level, all-day, all-weather, large range, long distance, and non-contact, and has great advantages in geological disaster monitoring and early warning. Using the principle of interferometric measurement, the ground-based radar can detect the minute deformation amount of the target area, and conduct real-time monitoring of geological disasters based on the cumulative deformation amount, deformation speed, deformation acceleration, etc. However, in the long-term geological disaster monitoring, the deformation of some areas does not cause disasters such as landslides, but due to the influence of human factors, geological factors, meteorological factors, etc., the radar images of the monitoring areas are severely decorrelated, bringing great difficulties to long-term quantitative monitoring. Therefore, on the basis of quantitative monitoring, combining with the application of change detection to provide effective information for comprehensively understanding the dynamic changes of the monitoring area in the long term has become an urgent problem to be solved at present.

[0068] Traditional high-resolution optical remote sensing has a high spatial resolution, and its images can show detailed geometric information of ground objects, making it the main data source for remote sensing image change detection. However, the imaging principle of optical remote sensing images also determines that they are affected by external conditions such as weather, observation time, and atmospheric conditions, restricting the acquisition of image information of the target area on schedule. The ground-based radar can obtain a stable data source within the monitoring range and supplement the optical images, but the coherent speckle noise brought by the ground-based radar imaging method determines that it is difficult to extract the change information of the ground-based radar images. Therefore, optical images and ground-based radar images are complementary to each other, making up for the deficiencies of single remote sensing images.

[0069] As described in the background art section above, the present disclosure exemplarily records corresponding solutions in the form of embodiments to address the deficiencies in the prior art, but this does not limit the scope of patent protection required by the present disclosure.

[0070] As one of the solutions, as Figure 1 shows a schematic diagram of the implementation process of a change detection method based on multi-source remote sensing data according to an embodiment of the present disclosure. The embodiments of the present disclosure provide a change detection method based on multi-source remote sensing data, including:

[0071] Obtain a complex radar image from the ground-based radar image;

[0072] Based on the radar images before and after the change, register them with the radar images at the corresponding times;

[0073] After registration, fuse the radar image with the optical image to obtain a fused difference map, and obtain the change detection result through segmentation and clustering.

[0074] Regarding the above-mentioned content, the embodiments of the present disclosure aim to comprehensively understand the dynamic change information of the monitoring area, propose a multi-source remote sensing image change detection method, obtain change information including the fusion of optics and radar scattering, so as to accurately locate and analyze the changes in the observation area, provide effective data reference for geological disaster monitoring, and effectively avoid misjudgment of the deformation results by personnel. In addition, the multi-source fusion image change detection can accurately detect and precisely locate the change areas caused by geological disasters such as landslides and collapses, providing a reference basis for disaster prevention and mitigation. As an efficient and accurate multi-source fusion image change detection method, the embodiments of the present disclosure provide effective data reference for the comprehensive monitoring and analysis of the target area through the fusion and change detection of multi-source remote sensing data in the monitoring area, thereby realizing the accurate location and qualitative monitoring of the change area.

[0075] The steps of the embodiments of the present disclosure, by way of example, are further described below through steps S1 to S4 as an example for the change detection method based on multi-source remote sensing data of the present disclosure. Figure 1 The step S5 shown in the figure exemplifies the step of the change detection result according to the change detection method.

[0076] In some specific implementation schemes, the present disclosure may be to obtain a complex radar image from the ground-based radar image, including:

[0077] Obtain the original radar data of the monitored slope area;

[0078] Within one sub-cycle, the frequency band of the output signal of each channel remains unchanged, and in the next sub-cycle, the frequency bands of the output signals of all channels are rotated;

[0079] Perform two-dimensional imaging processing on the original radar data respectively;

[0080] Perform image registration and atmospheric phase correction in sequence to obtain a radar complex image reflecting the monitored area in the nth observation period.

[0081] Step S1: Acquisition of ground-based radar images. The acquisition of radar images is achieved through a ground-based radar system, as Figure 2 is a schematic diagram of data acquisition for the ground-based radar system. Specific step S1 may include steps S11 to S12.

[0082] Step S11: Obtain the original radar data of the monitored slope area periodically transmitted by the ground-based radar.

[0083] In this embodiment, at least five corner reflectors are placed at the monitoring site for optical image and radar image registration. The ground-based radar is placed statically, and its transmitting / receiving antenna beam irradiates the monitored area. At a certain moment, it receives the control parameters and start instructions from the remote monitoring center, conducts a radar data acquisition once, and records the initially acquired original data as s0, and uploads it to the remote monitoring center.

[0084] The ground-based radar repeats the above acquisition process with an observation period of Δt, that is, at intervals of time Δt, the ground synthetic aperture radar conducts a radar data acquisition under the control of the remote monitoring center, and uploads the acquired data s n to the remote control center until the remote monitoring center issues a termination instruction for data acquisition, where s n represents the data acquired after the nth observation period.

[0085] Step S12: Perform two-dimensional imaging processing on the original radar data s0 and s n respectively to obtain complex radar images, denoted as h0(x, y) and h n (x, y) respectively, where x and y represent the horizontal and vertical distances of the observation target relative to the ground synthetic aperture radar. Taking h0(x, y) as the main image and h n (x, y) as the sub-image, perform image registration and atmospheric phase correction in sequence to obtain the radar complex image d n (x, y) reflecting the monitored area in the nth observation period.

[0086] In this embodiment, the distance-Doppler algorithm, Chirp-Scaling algorithm or range migration algorithm can be used to perform two-dimensional imaging processing on the original radar data to obtain complex radar images.

[0087] In some specific implementation schemes, the present disclosure may be based on the radar images before and after changes and registration with the radar images at corresponding times, including:

[0088] Select any two pre-change radar images and post-change radar images obtained at any moment;

[0089] Denoise two radar images;

[0090] Obtain the optical images at corresponding times;

[0091] Perform radiometric correction on the optical images;

[0092] Register the corrected optical images with the radar images at corresponding times.

[0093] The embodiments of the present disclosure may further include step S2 on the basis of the above step S1.

[0094] Step S2: Preprocessing of radar images and optical images. Select any two radar images h i and h j obtained in step S12. h i is the radar image before change, h j is the radar image after change. Denoise the two radar images; obtain the optical images p i and p j at corresponding times; then perform radiometric correction on the optical images p i and p j ; then register the corrected optical images with the radar images at corresponding times; as Figure 3 shown is the implementation step of the preprocessing.

[0095] In this embodiment, the acquisition of optical images needs to ensure that the angular positions of the acquired images are consistent as much as possible at each moment.

[0096] Specifically, the preprocessing step S2 may include steps S21 to S23.

[0097] Step S21: Acquisition of optical images. Obtain the optical images p i and p j at corresponding times through an optical remote sensing device.

[0098] In the embodiments of the present disclosure, the optical device may adopt a high-definition optical remote sensing device, and the acquisition interval needs to be determined according to actual applications. The acquisition of optical images needs to ensure that the angular positions of the acquired images are consistent as much as possible at each moment.

[0099] In some specific implementation schemes, the present disclosure may be

[0100] Registering the corrected optical images with the radar images at corresponding times includes:

[0101] Selecting feature points in the optical images and selecting corresponding feature points in the radar images as feature point pairs;

[0102] Obtain a transformation matrix for perspective transformation based on the feature point pairs;

[0103] Convert the optical image to the radar image coordinate system to obtain the final image matching result;

[0104] Obtain the interpolated result based on the information of the ground-based radar image pixels and the optical image pixels.

[0105] Step S22: Perform radiometric correction on the optical images p i and p j to obtain the pair of radiometrically corrected images; among them, the histogram method can be used for atmospheric radiometric correction.

[0106] Step S23: Register the radar image and the optical image, and the registration algorithm uses the perspective transformation matching algorithm.

[0107] As Figure 3 shown in the registration process, the specific algorithm steps S23 may include steps S231 to S235.

[0108] Step S231: Manually select the feature points in the optical images p i and p j , denoted as (x' s , y' s ), where s is the number of feature points, s ∈ {1, 2... k}. The selection of feature points is the corner reflectors arranged on site; similarly, select the corresponding feature points (x i and h j ) in the radar images h s , y s .

[0109] Step S232: According to the feature point pairs (x s , y s ) and (x' s , y' s ), calculate the transformation matrix H for perspective transformation; where H is:

[0110]

[0111] The calculation method of the transformation matrix H is obtained by least squares estimation through the coordinates of the feature point pairs.

[0112] Step S233: The calculation formula for converting the optical image to the radar image coordinate system is:

[0113]

[0114] where [x”y”1] T is the coordinate after the transformation of the optical image, and [x'y'1] Tare the coordinates before optical image transformation. From the above formula, the following transformation relationship of the optical image coordinate system can be obtained:

[0115]

[0116] Step S234: According to the generation method of the feature point descriptor, the final image matching result can be obtained.

[0117] Step S235: Interpolate the optical image pixel information based on the information of the ground-based radar image pixels, that is, resampling processing; among them, the nearest neighbor interpolation method, bilinear interpolation method, and cubic convolution interpolation method can be used for resampling processing to obtain the interpolated result: optical image p i ' and p' j .

[0118] The embodiment of the present disclosure may further include step S3 on the basis of the above step S2.

[0119] Step S3: Fuse the registered radar image and the optical image, that is, the optical image p before transformation i ' and the radar image h i are fused to obtain the pre-change fused image X1; the optical image p' after transformation j and the radar image h j are fused to obtain the post-change fused image X2; among them, the fusion algorithm can adopt the HIS transformation-based fusion algorithm, wavelet transformation-based fusion algorithm, or PCA-based fusion algorithm.

[0120] In some specific implementation manners, the present disclosure may be that the registered radar image and the optical image are fused to obtain a fusion difference map, including:

[0121] Obtain the fusion difference map based on the neighborhood mean ratio difference map and the improved relative entropy difference map method;

[0122] Perform filtering processing to obtain the final fusion difference map.

[0123] The embodiment of the present disclosure may further include step S4 on the basis of the above step S3.

[0124] Step S4: Fusion image change detection, using the neighborhood mean ratio difference map and the improved relative entropy difference map method, through NSCT transform fusion, obtain the fusion difference map, and use bilateral filtering to filter the fusion difference map to obtain the final fusion difference map; then segment and cluster the final difference map into changed and unchanged classes to obtain the final change detection result.

[0125] In some specific implementation manners, the present disclosure may be that the fusion difference map is obtained based on the neighborhood mean ratio difference map and the improved relative entropy difference map method, including:

[0126] Apply the mean logarithmic ratio operator to the fused image to generate a mean logarithmic ratio map;

[0127] Apply the improved relative entropy to the fused image to generate an improved neighborhood mean relative entropy image;

[0128] Generate two fused difference maps using the NSCT transform;

[0129] Fuse the difference maps by the local energy method to obtain the final fused difference map.

[0130] As Figure 4 shown in the change detection flow chart, the specific implementation steps S4 may include steps S41 to S47.

[0131] Step S41: Apply the mean logarithmic ratio operator to the fused images X1 and X2 to generate a mean logarithmic ratio map:

[0132]

[0133] In formula (4), m1(i,j) and m2(i,j) are the neighborhood means respectively.

[0134] Step S42: Apply the improved relative entropy to the fused images X1 and X2 to generate an improved neighborhood mean relative entropy image:

[0135]

[0136] where X n is the nth pixel corresponding in X1 and X2, X k is the kth pixel corresponding in X1 and X2, X η is the pixel in the neighborhood Ω of pixel X k m is the mean of the pixels in the neighborhood Ω, and α is the weighting parameter. When α = 0, it is equivalent to not introducing neighborhood calculation.

[0137] In some specific implementation schemes, the present disclosure may be that two fused difference maps are generated using the NSCT transform, including:

[0138] Use the generated mean logarithmic ratio map and neighborhood mean relative entropy image as the input images for NSCT transform processing;

[0139] Use the NSCT transform to obtain the low-frequency information and high-frequency information in the mean logarithmic ratio map and neighborhood mean relative entropy image;

[0140] Fuse the transformed low-frequency coefficients to obtain the fused low-frequency coefficients;

[0141] Inverse NSCT transform is performed on the fused low-frequency coefficients and the high-frequency coefficients respectively to obtain the fused difference map.

[0142] Step S43: Generate two fused difference maps by NSCT transform. Figure 5 The NSCT transform processing flow chart is shown, and specific step S43 may include step S431 to step S434.

[0143] Step S431: Use the coherence coefficient map X generated in step S41 m and the logarithmic mean ratio map X generated in step S42 s as the input images for NSCT transform processing.

[0144] Step S432: Use NSCT transform to obtain the low-frequency information and high-frequency information in the difference maps X s and X m . This method has good translational invariance and direction selectivity and can capture geometric information in the image more accurately; the difference maps X s and X m corresponding low-frequency coefficients L s and L m , high-frequency coefficients H s and H m are obtained. Among them, the contourlet transform uses a "9-7" type Laplacian pyramid filter, uses "pkva" as the direction group filter, and the decomposition layer of each pyramid level direction filter group is 4.

[0145] Step S433: In a weighted average manner, fuse the low-frequency coefficients L s and L m obtained by NSCT transform of the difference maps X s and X m to obtain the fused low-frequency coefficient L f . The calculation process is as follows:

[0146] L f =αL s +(1-α)L m (6)

[0147] where α (α≥0) is the fusion coefficient of weighted average, and α = 0.3 is taken.

[0148] Step S434: Perform inverse NSCT transform on the fused low-frequency coefficient L f and the high-frequency coefficients of the difference maps X s and X m respectively to obtain the fused difference maps F s and F m .

[0149] Step S44: Use the local energy method to fuse the difference maps F s and F m to obtain the final fused difference map F h . The calculation process is as follows:

[0150]

[0151] Among them, E s (i,j) and E m (i,j) respectively represent the local energies of the difference maps F s and F m . L i,j represents a 3×3 neighborhood centered on the pixel (i,j), and F s (q) and F m (q) represent the q-th pixel value in the neighborhood.

[0152] In some specific implementation schemes, the present disclosure may be to perform a filtering process to obtain the final fused difference map, including:

[0153] Use bilateral filtering to perform edge preservation and further denoising to obtain the final difference map.

[0154] Step S45: Use bilateral filtering to perform edge preservation and further denoising on the fused difference map to obtain the final difference map. The filtered image is as shown in Equation (9):

[0155]

[0156] Among them, D(i,j) is the pixel value at the coordinate (i,j) in the fused difference map, and w(i,j) is the weight coefficient of the bilateral filtering, as shown in Equation (10):

[0157] w(i,j) = w v (i,j)w u (i,j) (10)

[0158] Among them, w v (i,j) is the pixel similarity weight, and w u (i,j) is the spatial distance weight, as shown in Equations (11) and (12)

[0159]

[0160] In Equations (11) and (12), x and y are the horizontal and vertical coordinates within the 3×3 neighborhood of the pixel point at the coordinate (i,j) in the fused difference map D(i,j), and δ v and δ u are adjustment parameters.

[0161] In some specific embodiments, the present disclosure may be that the change detection result is obtained through segmentation and clustering, including: converting the bilaterally filtered image into a column vector; initializing the membership degree to meet the preset conditions; calculating the cluster centers; calculating the objective function; updating the membership degree matrix; performing iterative loops until the objective function is less than the set threshold; comparing the magnitudes of the membership degrees of each pixel, and performing class labeling on the converted image based on the two types of membership degrees corresponding to each pixel.

[0162] Step S46: Perform initial clustering segmentation on the fusion difference map, dividing the difference map into a changed class and an unchanged class, where the Fuzzy C-Means (FCM) algorithm is selected for initial segmentation, and the specific steps are as Figure 6 The FCM initial segmentation flowchart is given, and the detailed steps S46 may include steps S461 to S467.

[0163] Step S461: Convert the bilaterally filtered image R m×n (i,j) into a column vector, where {1≤i≤m, 1≤j≤n}, and the converted image is R' mn×1 (i,j), where {1≤i≤m×n, j = 1}.

[0164] Step S462: Initialize the membership degree u to satisfy:

[0165]

[0166] where c is the number of classifications, and in the present disclosure, c = 2 (changed class and unchanged class).

[0167] Step S463: Calculate the cluster centers C s :

[0168]

[0169] where m is the fuzzy weight (generally 2).

[0170] Step S464: Calculate the objective function J:

[0171]

[0172] Step S465: Update the membership degree matrix u according to the following formula, where u is solved by the Lagrange multiplier method:

[0173]

[0174] Step S466: Perform iterative loops. If the objective function is less than the set threshold ε, stop the iteration; otherwise, return to step S463 and continue to execute.

[0175] Step S467: Compare the membership degree u of each pixelki in size, where each pixel corresponds to two types of membership degrees u k1 and u k2 , and the following decision condition is used to perform class labeling on R':

[0176]

[0177] X k is the classification label of the k-th pixel.

[0178] In some specific embodiments, the present disclosure may be that, after segmentation and clustering to obtain a change detection result, it further includes: normalizing the fused difference map after bilateral filtering; calculating the class mean and class variance under the initial label distribution, and the class mean; determining the optimal threshold; comparing the neighborhood features with the optimal threshold; calculating the non-local mean after introducing the optimal threshold and introducing the non-local mean of the optimal threshold into the likelihood function; performing iterative loop until the likelihood energy function is minimized, and outputting the change detection result.

[0179] Step S47: Divide the image into a changed class and an unchanged class, as Figure 7 shown in the schematic flowchart of the non-local mean MAP-MRF ground-based radar image segmentation based on the Fisher optimal criterion, where the segmentation method proposes a non-local mean MAP-MRF ground-based radar image segmentation method based on the Fisher optimal criterion. The specific step S47 may include steps S471 to S477.

[0180] Step S471: Normalize the fused difference map R(i,j) after bilateral filtering,

[0181]

[0182] Define the normalized image as the feature field Y = {Y k = Y ij}}, and define the labeled field X = {X k = X ij} for the image after initial segmentation. X k ∈I (I = 1, 2,..., c), where c is the number of classes, and c = 2 in the present disclosure.

[0183] Step S472: Calculate the class mean and class variance under the initial label distribution. The class mean is:

[0184]

[0185] where h(i) is the normalized histogram, t is the threshold for dividing the image gray level into two categories of target and background (it is considered that the changed image is composed of dark objects on a bright background), and θ(t) and 1 - θ(t) are the prior probabilities of the target and background respectively. The formula for calculating the class variance is:

[0186]

[0187] Step S473: Determine the optimal threshold through the Fisher criterion. The expression of the Fisher evaluation function is:

[0188]

[0189] When the threshold t is the optimal threshold, J(t) is the maximum value. At this time, the separation degree of the target and the background divided by the threshold t is the best, that is, the optimal threshold t * can be expressed as:

[0190] t * = Argmax[J(t)] (22)

[0191] Step S474: Compare the neighborhood feature with the optimal threshold. For each pixel k in the image, extract its neighborhood feature vector as:

[0192] F k = {F(p) ∈ W(p)} (23)

[0193] where W(p) represents the non-local window centered on p.

[0194] The Fisher threshold t calculated by non-local space * is used as the optimal threshold, that is, each pixel k corresponds to a Fisher optimal threshold in its non-local space For each neighborhood pixel k, label it as the target and background classes according to the following comparison results:

[0195]

[0196] Then the homogeneous feature vector is:

[0197]

[0198] where is a new logical feature vector generated by the comparison result of the neighborhood feature vector, is the logical feature vector corresponding to the central pixel.

[0199] Step S475: Calculate the non-local mean after introducing the Fisher optimal threshold. The weight of the non-local mean (NLM) algorithm after introducing the Fisher optimal threshold is:

[0200]

[0201] Then the non-local mean after introducing the Fisher optimal threshold:

[0202]

[0203] where |X i -X s | represents the Euclidean distance between neighborhood feature vectors, Z is the normalization constant of the weight, and v is a parameter that controls the decay rate of the Gaussian function.

[0204] Step S476: The likelihood function calculates the new label value. Introducing the non-local mean of the Fisher optimal threshold into the likelihood function can be expressed as the following formula:

[0205]

[0206] where, V h is the group potential energy function with the mathematical model being the Potts model, and H is the set of all groups.

[0207] Step S477: Iterate cyclically. When the likelihood energy function U(X|Y) is minimized, stop the iteration and output the change detection result. Otherwise, return to Step S472 and continue to execute.

[0208] The change detection method based on multi-source remote sensing data in each embodiment of the present disclosure has at least the following beneficial effects compared with the prior art:

[0209] In the present disclosure, a change detection method is studied using multi-source remote sensing images, which can make the optical image and the ground-based radar image complement each other well, making up for the deficiencies of a single remote sensing image. A non-local mean MAP-MRF ground-based radar image segmentation method based on the Fisher optimal criterion is proposed for the multi-source fusion difference map. Based on the fact that the MRF model can effectively protect image details and reduce speckle noise, the NLM algorithm and the Fisher evaluation function are introduced, and a non-local mean algorithm based on the optimal threshold is proposed to maintain the fine structure of the image and enhance the resistance to speckle noise, further improving the accuracy of change detection. The proposed method has certain advantages over the general MRF algorithm.

[0210] Change detection of multi-source fusion images can provide data reference for geological disaster deformation monitoring, and can effectively avoid misjudgment of the deformation results by staff. In addition, change detection of multi-source fusion images can accurately detect the changed areas caused by geological disasters such as landslides and collapses and accurately locate them, providing a reference basis for disaster prevention and mitigation. The present disclosure is an efficient and accurate change detection method for multi-source fusion images.

[0211] As one of the solutions, an embodiment of the present disclosure provides a change detection device based on multi-source remote sensing data, including at least one processing module configured to execute the method described in the above Steps S1 to S4.

[0212] Specifically, one of the inventive concepts of the present disclosure, the change detection method based on multi-source remote sensing data of various embodiments of the present disclosure, obtains complex radar images at least according to ground-based radar images; based on the radar images before and after the change, registers with the radar images at the corresponding time; the registered radar images are fused with optical images to obtain a fused difference map, and the change detection result is obtained through segmentation and clustering. Thus, the research on the change detection method using multi-source remote sensing images can make the optical images and ground-based radar images complement each other well, making up for the deficiencies of single remote sensing images. Further, each embodiment of the present disclosure proposes an optimal criterion non-local mean ground-based radar image segmentation method for the multi-source fused difference map, which can, on the basis of effectively protecting image details and reducing speckle noise, propose a non-local mean algorithm based on an optimal threshold, and maintain fine image structures and enhance the resistance to speckle noise, further improving the accuracy of change detection.

[0213] The present disclosure also provides a computer-readable storage medium, on which computer-executable instructions are stored. When the computer-executable instructions are executed by a processor, it mainly implements the change detection method based on multi-source remote sensing data as described above, at least including:

[0214] Obtain complex radar images according to ground-based radar images;

[0215] Based on the radar images before and after the change, register with the radar images at the corresponding time;

[0216] The registered radar images are fused with optical images to obtain a fused difference map, and the change detection result is obtained through segmentation and clustering.

[0217] The above embodiments are only exemplary embodiments of the present disclosure and are not used to limit the present disclosure. The protection scope of the present disclosure is defined by the claims. Those skilled in the art can make various modifications or equivalent replacements within the essence and protection scope of the present disclosure, and such modifications or equivalent replacements should also be regarded as falling within the protection scope of the present disclosure.

Claims

1. Change detection method based on multi-source remote sensing data, including: Obtaining a complex radar image based on a ground-based radar image; Based on the radar images before and after the change, the radar images at the corresponding time are registered; After registration, the radar image and the optical image are fused to obtain a fused difference map, and the change detection result is obtained through segmentation and clustering.

2. The method according to claim 1, wherein: The complex radar images are obtained from the ground-based radar images, including: Obtaining raw radar data of the monitored slope area; In one sub-cycle, the frequency band of the output signal of each channel remains unchanged, and in the next sub-cycle, the frequency bands of the output signals of all channels are rotated; Perform two-dimensional imaging processing on the original radar data respectively; Image registration and atmospheric phase correction are performed in sequence to obtain a radar complex image reflecting the monitored area in the nth observation period.

3. The method according to claim 2, wherein: Based on the radar images before and after the change, the radar images at the corresponding time are registered, including: Select two radar images before and after the change at any time; De-noising the two radar images; Acquire an optical image at a corresponding moment; Performing radiation correction on the optical image; The corrected optical image is registered with the radar image at the corresponding moment.

4. The method according to claim 3, wherein: The corrected optical image is registered with the radar image at the corresponding time, including: Selecting feature points in the optical image and corresponding feature points in the radar image as feature point pairs; According to the feature point pairs, the transformation matrix for perspective transformation is obtained; The optical image is converted into the radar image coordinate system to obtain the final image matching result; The interpolation result is obtained according to the pixel information of the ground-based radar image and the pixel information of the optical image.

5. The method according to claim 4, wherein: After registration, the radar image and the optical image are fused to obtain a fusion difference map, including: Based on the neighborhood mean ratio difference map and the improved relative entropy difference map method, a fused difference map is obtained; Perform filtering to obtain the final fusion difference map.

6. The method according to claim 5, wherein: Based on the neighborhood mean ratio difference map and the improved relative entropy difference map method, a fused difference map is obtained, including: Applying the mean log ratio operator to the fused image generates a mean log ratio map; Applying improved relative entropy to the fused image to generate an improved neighborhood mean relative entropy image; The NSCT transform is used to generate two fused difference images; The difference maps are fused by the local energy method to obtain the final fused difference map.

7. The method according to claim 6, wherein: The NSCT transform is used to generate two fused difference maps, including: The generated mean logarithmic ratio image and neighborhood mean relative entropy image are used as the input image for NSCT transformation processing; Use NSCT transformation to obtain low-frequency and high-frequency information in the mean logarithmic ratio image and the neighborhood mean relative entropy image; The transformed low-frequency coefficients are fused to obtain fused low-frequency coefficients; The fused low-frequency coefficients are used to perform NSCT inverse transformation with the high-frequency coefficients respectively to obtain the fused difference map.

8. The method according to claim 7, wherein: Perform filtering to obtain the final fusion difference map, including: Bilateral filtering is used to perform edge preservation and further denoising to obtain the final difference map.

9. The method according to claim 8, wherein: The change detection results obtained through segmentation and clustering include: Convert the bilaterally filtered image into a column vector; Initialize the membership to meet the preset conditions; Calculate cluster centers; Calculate the objective function; Update the membership matrix; Iterate the loop until the objective function is less than the set threshold; The membership of each pixel is compared, and the converted image is labeled based on the two types of membership corresponding to each pixel.

10. The method according to claim 9, wherein: The change detection results obtained by segmentation and clustering also include: Normalize the fused difference map after bilateral filtering; Calculate the class mean, class variance, and class mean under the initial label distribution; Determine the optimal threshold; Compare neighborhood features with the optimal threshold; Calculate the non-local mean after introducing the optimal threshold Introduce the non-local mean of the optimal threshold into the likelihood function; The loop is iterated until the likelihood energy function is minimized and the change detection result is output.

Images