Automatic flood detection method based on coupling multidimensional K-Means and GMM

Abstract

The present invention discloses a flood automatic detection method coupling multidimensional K-Means and GMM, constructs 4D-SAR data through time series polarization SAR data, and uses different distance functions to improve K-Means clustering for polarization information in feature dimensions so that it can better distinguish water body features. At the same time, the clustering results are used to initialize 4D-GMM, and multiple matrix Gaussian distributions are used to better fit the distribution characteristics of 4D-SAR data features, which can effectively solve the problem of unreasonable classification of some features caused by "hard clustering", provide a basis for efficiently obtaining the scope of flood damage, and fundamentally overcome the problems of insufficient information utilization and inability to better fit the distribution of SAR image features in the unsupervised change detection method of SAR images.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing image processing, and in particular to an automatic flood detection method coupling multidimensional K-Means and GMM. Background Art

[0002] Remote sensing change detection is the process of identifying and detecting changes in land cover by analyzing remote sensing images acquired at different times at the same location. It is a key technology for Earth observation applications. Floods are characterized by sudden onset and widespread impact. SAR sensors, unaffected by light and fog, can provide high-precision change detection during flood disasters and have become the most effective means of flood detection. Traditional change detection methods include pixel-based change detection and object extraction-based methods. Pixel-based change detection detects areas of change by comparing pixel values ​​before and after the change. It can detect not only significant changes but also subtle changes, offering high sensitivity and accuracy. However, traditional threshold segmentation methods lack an efficient and automatic method for determining the threshold. Current selection of the change threshold is often driven by human factors and relies on the experience accumulated by researchers over long periods of experimentation, leading to false detections and missed detections.

[0003] While adaptive threshold extraction methods for difference images, such as the constant false alarm rate threshold extraction method, the Otsu maximum inter-class variance threshold method, the Kapur entropy threshold method, and the Kittler and Illingworth threshold method, can effectively perform threshold extraction on difference images, they do not consider the probabilistic statistical characteristics of SAR image difference maps. Object extraction-based methods offer certain accuracy and robustness, but they have limitations when dealing with complex water environments. Therefore, to achieve more accurate flood detection results, a more appropriate change detection model should be used for flood extraction.

[0004] In the related art, when using SAR imagery for change detection, it is necessary to consider the statistical characteristics of multi-temporal SAR imagery and the universality of the model under different polarization data. Therefore, how to develop a highly universal, efficient, and accurate flood and waterlogging extraction method is an urgent problem for those skilled in the art. Summary of the Invention

[0005] To address the aforementioned technical deficiencies, the present invention aims to provide an automated flood detection method that couples multidimensional K-Means and GMM. This method fundamentally overcomes the issues of insufficient information utilization and inability to accurately fit the distribution of SAR image features found in unsupervised change detection methods based on SAR imagery. By constructing temporal, spatial, and polarimetric SAR data from time-series SAR data, the present invention improves K-Means clustering using different distance functions to better distinguish water features. Furthermore, post-processing the clustering results using GMM effectively addresses the issue of irrational feature classification caused by "hard clustering," providing a foundation for efficiently determining the extent of flood damage.

[0006] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0007] The present invention provides a flood automatic detection method coupling multidimensional K-Means and GMM, comprising the following steps:

[0008] The following steps are involved:

[0009] Step 1: Preprocess the time-series SAR images, including: performing orbit correction, thermal noise removal, radiometric calibration, speckle filtering, terrain correction, and decibel conversion on each phase and polarization SAR image; and registering the time-series SAR images using the same reference image to achieve sub-pixel registration accuracy.

[0010] Step 2: construct the pre-processed multi-dimensional SAR image in time dimension, space dimension and feature dimension to generate 4D-SAR data; the data is represented by a four-dimensional matrix P hwft =[p hwft ] H×W×F×T To express it, where T is the number of time phases, F is the length of the characteristic dimension of 4D-SAR data, including frequency and spectrum information;

[0011] p hwft Represents the fth feature of the pixel (h, w) in the tth time phase, t∈[1, T]; any pixel (h, w) is vectorized using a two-dimensional matrix P hwft =[p hwft ] F×T To express;

[0012] Step 3: Improve K-Means in multiple dimensions; input 4D-SAR data directly into the improved K-Means, define the number of clusters K, determine whether the number of iterations required is met, and obtain the cluster center. Among them, k∈[1,K], represents the kth class;

[0013] Step 4: Improve GMM in multiple dimensions to construct 4D-GMM, and use the results obtained by 4D-K-Means method to input different category results and matrix Gaussian distribution into 4D-GMM and initialize it, and calculate U in an iterative way. k 、V k The two covariances are estimated until convergence or the number of iterations is reached and the optimal parameters are obtained: w k , C k , U k , V k , where w k is the weight of the kth sub-matrix, C k ∈R F×T is the mean matrix of the kth submatrix, U k ∈R F×F 、V k ∈R T×T C k is the covariance matrix of the mean;

[0014] Step 5: Use the optimal parameters of step 4 to obtain a two-dimensional matrix Distribution And the final category label L for each pixel hw , where C∈R F×T is the mean matrix, U∈R F×F 、V∈R T×T are the covariance matrices with C as the mean;

[0015] Step 6: Binarize the obtained pixel category labels and extract the newly added flood area change categories to obtain the flood change detection results.

[0016] Preferably, in step 3, the 4D-SAR data constructed by polarimetric SAR image is input into the improved K-Means clustering, the number of clusters K is set, the input polarization mode is determined, and the distance d between the pixel and the cluster center is calculated. hwk :

[0017]

[0018] in, Represents the two-dimensional matrix after vectorization of any pixel (h, w), which is the 4D-SAR data matrix P hwft =[p hwft ] H×H×F×T part of the image, H and W are the length and width of the image, F is the length of the feature dimension of the 4D-SAR data, T is the number of time phases, T' is the transposition operation, p hwft represents the fth (f∈[1,F]) feature of the pixel (h,w) in the tth (t∈[1,T]) phase, represents the fth feature of the kth cluster center in the tth time phase, represents the pixel matrix after probabilization, and its calculation is shown in formula (2):

[0019]

[0020] represents the k-th cluster center matrix after probabilization, and its calculation is shown in formula (3):

[0021]

[0022] Calculate the matrix Q hwk =[q hwk ] H×W×K , and its update criterion is shown in formula (4):

[0023]

[0024] Use formula (5) to update the cluster center until the number of iterations is reached;

[0025]

[0026] The above steps are all processed in multidimensional data, and the above steps are repeated until the number of iterations is reached to obtain the final K cluster centers. ”

[0027] Preferably, in step 4, the matrix Gaussian distribution is used to fit Based on the distribution of , a 4D-GMM is constructed, and the clustering results obtained in step 3 are matched to the matrix Gaussian distribution. The 4D-GMM is initialized and iterated. The iterative steps are as follows:

[0028] Step 1 Initialize parameter w k , C k , U k , V k :

[0029]

[0030] Among them, w k is the weight of the kth sub-matrix, C k ∈R F×T is the mean matrix of the kth submatrix, U k ∈R F×F 、V k ∈R T×T C k is the covariance matrix of the mean, The final cluster center of the extended K-Means clustering;

[0031] Step.2 Calculate each pixel The posterior probability of the k-th submatrix belonging to the Gaussian distribution is normalized to l hwk :

[0032]

[0033] in, The calculation of is shown in formula (11):

[0034]

[0035] And w k satisfy:

[0036]

[0037] Step 3 Update parameters:

[0038]

[0039] Step.4 Repeat Step.2-Step.3 until convergence or the number of iterations is reached and the optimal parameter is obtained: w k , C k , U k , V k .”

[0040] Preferably, in step 5, the optimal parameters of step 4 are used to obtain a two-dimensional matrix Distribution As shown in formula (17):

[0041]

[0042] Where C∈R F×T is the mean matrix, U∈R F×F 、V∈R T×T are the covariance matrices with C as the mean; at the same time, obtain each pixel P hwft The final category label L hw , as shown in formula (18):

[0043] When l hwk Maximum, L hw =k, where k=1, 2, ..., K (36)

[0044] The beneficial effects of the present invention are: fundamentally overcoming the problems existing in the unsupervised change detection method of SAR images, such as insufficient information utilization, error accumulation caused by step-by-step processing, and inability to completely fit the SAR image feature distribution, achieving the distinction between permanent water bodies and changing water bodies with high precision, and at the same time introducing time information to achieve effective extraction of weakly scattering ground objects. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0046] Figure 1 4D-SAR data four-dimensional matrix diagram provided by an embodiment of the present invention;

[0047] Figure 2 This is a flow chart of an unsupervised coupled multidimensional K-Means and GMM flood automatic detection method provided by an embodiment of the present invention;

[0048] Figure 3 This is a true value map of the single-polarization SAR image research area and its application results provided by an embodiment of the present invention;

[0049] Figure 4 This is a true value map of the dual-polarization SAR image research area and its application results provided by an embodiment of the present invention;

[0050] Figure 5 This is a true value map of the fully polarimetric SAR imagery research area and its application results provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0051] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0052] like Figure 2 As shown, a flood automatic detection method coupling multidimensional K-Means and GMM includes the following steps:

[0053] Step 1: Preprocess the time-series SAR images, including: performing orbit correction, thermal noise removal, radiometric calibration, speckle filtering, terrain correction, and decibel conversion on each phase and polarization SAR image; and registering the time-series SAR images using the same reference image to achieve sub-pixel registration accuracy.

[0054] Step 2: construct the pre-processed multi-dimensional SAR image in time dimension, space dimension and feature dimension to generate 4D-SAR data; the data is represented by a four-dimensional matrix P hwft =[p hwft ] H×W×F×T To express, such as Figure 1 As shown in Figure 1, where T is the number of time phases and F is the length of the characteristic dimension of 4D-SAR data, including frequency and spectrum information. hwft Represents the fth feature of the pixel (h, w) in the tth phase, t∈[1, T]; any pixel (h, w) is vectorized using a two-dimensional matrix Taking polarimetric SAR images as an example, for single-polarization and dual-polarization SAR data, their features are composed of eigenvalues ​​corresponding to different polarization modes; for full-polarization SAR data, a coherent matrix that can better reflect the characteristics of full-polarization data is selected;

[0055] Step 3: Improve K-Means in multiple dimensions. Input 4D-SAR data directly into the improved K-Means, define the number of clusters K, determine whether the number of iterations required is met, and obtain the cluster center. Among them, k∈[1,K], represents the kth class;

[0056] Input the 4D-SAR data constructed by polarimetric SAR image into the improved K-Means clustering, set the number of clusters K, determine the input polarization mode, and calculate the distance d between the pixel and the cluster center. hwk :

[0057]

[0058] in, Represents the two-dimensional matrix after vectorization of any pixel (h, w), which is the 4D-SAR data matrix P hwft =[p hwft ] H×W×F×T part of the image, H and W are the length and width of the image, F is the length of the feature dimension of the 4D-SAR data, T is the number of time phases, T' is the transposition operation, phwft represents the fth (f∈[1,F]) feature of the pixel (h, w) in the tth (t∈[1,T]) time phase, represents the fth feature of the kth cluster center in the tth time phase, represents the pixel matrix after probabilization, and its calculation is shown in formula (2):

[0059]

[0060] represents the k-th cluster center matrix after probabilization, and its calculation is shown in formula (3):

[0061]

[0062] Calculate the matrix Q hwk =[q hwk ] H×W×K , and its update criterion is shown in formula (4):

[0063]

[0064] Use formula (5) to update the cluster center until the number of iterations is reached;

[0065]

[0066] The above steps are all processed in multidimensional data, and the above steps are repeated until the number of iterations is reached to obtain the final K cluster centers.

[0067] Step 4: Improve GMM in multiple dimensions to construct 4D-GMM. Using the results obtained by 4D-K-Means method, input the different category results and matrix Gaussian distribution into 4D-GMM and initialize it, and calculate U in an iterative way. k 、V k The two covariances are estimated until convergence or the number of iterations is reached and the optimal parameters are obtained: w k , C k , U k , V k , where w k is the weight of the kth sub-matrix, C k ∈R F×T is the mean matrix of the kth submatrix, U k ∈R F×F 、V k ∈R T×T C k is the covariance matrix of the mean;

[0068] Using matrix Gaussian distribution fitting Based on the distribution of , a 4D-GMM is constructed, and the clustering results obtained in step 3 are matched to the matrix Gaussian distribution. The 4D-GMM is initialized and iterated. The iterative steps are as follows:

[0069] Step 1 Initialize parameter w k , C k , U k , V k :

[0070]

[0071] Among them, w k is the weight of the kth sub-matrix, C k ∈R F×T is the mean matrix of the kth submatrix, U k ∈R F×F 、V k ∈R T×T C k is the covariance matrix of the mean, is the final cluster center of the extended K-Means clustering.

[0072] Step.2 Calculate each pixel The posterior probability of the k-th submatrix belonging to the Gaussian distribution is normalized to l hwk :

[0073]

[0074] in, The calculation of is shown in formula (11):

[0075]

[0076] And w k satisfy:

[0077]

[0078] Step 3 Update parameters:

[0079]

[0080] Step.4 Repeat Step.2-Step.3 until convergence or the number of iterations is reached and the optimal parameter is obtained: w k , C k , U k , V k .

[0081] Step 5: Use the optimal parameters described in step 4 to obtain a two-dimensional matrix Distribution And the final category label L for each pixel hw , where C∈R F×T is the mean matrix, U∈R F×F 、V∈RT×T are the covariance matrices with C as the mean.

[0082] Specifically: Use the optimal parameters described in step 4 to obtain the two-dimensional matrix P hwft Distribution As shown in formula (17):

[0083]

[0084] Where C∈R F×T is the mean matrix, U∈R F×F 、V∈R T×T are the covariance matrices with C as the mean. At the same time, get each pixel The final category label L hw , as shown in formula (18):

[0085] When l hwk Maximum, L hw =k, where k=1, 2, ..., K (54)

[0086] Step 6: Binarize the obtained pixel category labels and extract the newly added flood area change categories to obtain the flood change detection results.

[0087] The effects of the embodiments of the present invention can be further verified by the following experiments:

[0088] 1) Experimental environment

[0089] Python 3.7, i9 processor.

[0090] 2) Experimental content

[0091] The method of the present invention, i.e. the flood automatic detection method coupling multidimensional K-Means and GMM, is used to perform change detection on single, dual and full polarization data. The experimental results corresponding to the three polarization modes are as follows: Figure 3 、 Figure 4 、 Figure 5 As shown, Figure 3 (a) is the ground truth map. Figure 3 (b) is a binary graph of the experimental results of the method of the present invention. Figure 4 、 Figure 5 and Figure 3 Similarly, it corresponds to the dual-polarization and full-polarization datasets.

[0092] 3) Accuracy evaluation

[0093] In order to quantitatively evaluate the effect of the present invention, the present invention selects overall accuracy (OA), false alarm rate (FA) and Kappa coefficient as evaluation indicators.

[0094] Table 1 Accuracy evaluation of the change detection results of the method of the present invention

[0095]

[0096]

[0097] 4) Experimental results analysis

[0098] Figure 3 、 Figure 4 、 Figure 5 The experimental results in Table 1 show that the change detection results obtained by the proposed method are closer to the ground truth map of flooded areas, can better distinguish the differences in changes between different ground objects, significantly reduce false alarms, and can effectively detect small areas of change. The quantitative evaluation indicators also have a clear advantage. Therefore, the proposed method significantly improves the change detection effect in multiple aspects.

[0099] Obviously, those skilled in the art may make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if such changes and modifications fall within the scope of the claims and their equivalents, the present invention is intended to include such changes and modifications.

Claims

1. The automatic flood detection method coupling multidimensional K-Means and GMM is characterized by: The following steps are involved: Step 1: preprocess the time series SAR images, including: orbit correction, thermal noise removal, radiation calibration, speckle filtering, terrain correction, and decibel conversion for each phase and different polarization SAR images; select the same reference image to register the time series SAR images, so that the registration accuracy reaches the sub-pixel level; Step 2: construct the pre-processed multi-dimensional SAR image in the time dimension, space dimension and feature dimension to generate 4D-SAR data; the data is represented by a four-dimensional matrix P hwft =[p hwft ] H×W×F×T To express it, where T is the number of time phases, and F is the length of the characteristic dimension of 4D-SAR data, including frequency and spectrum information; p hwft represents the fth feature of the pixel (h, w) in the tth phase, t∈[1,T]; any pixel (h, w) is vectorized using a two-dimensional matrix To express; Step 3: Improve K-Means in multiple dimensions; directly input 4D-SAR data into the improved K-Means, define the number of clusters K, determine whether the number of iterations required is met, and obtain the cluster center. Among them, k∈[1,K], represents the kth class; Specifically, the 4D-SAR data constructed by polarimetric SAR images are input into the improved K-Means clustering, the number of clusters K is set, the input polarization mode is determined, and the distance d between the pixel and the cluster center is calculated. hwk : in, Represents the two-dimensional matrix after vectorization of any pixel (h, w), which is the 4D-SAR data matrix P hwft =[p hwft ] H×Ww×F×T part of the image, H and W are the length and width of the image, F is the length of the characteristic dimension of the 4D-SAR data, T is the number of time phases, T' is the transposition operation, and p hwft represents the fth (f∈[1,F]) feature of the pixel (h,w) in the tth (t∈[1,T]) phase, represents the fth feature of the kth cluster center in the tth phase, represents the probabilistic pixel matrix, which is calculated as shown in formula (2): represents the k-th cluster center matrix after probabilization, and its calculation is shown in formula (3): Calculate the matrix Q hωk =[q hwk ] H×W×K , and its update criterion is shown in formula (4): Use formula (5) to update the cluster center until the number of iterations is reached; The above steps are all processed in multidimensional data, and the above steps are repeated until the number of iterations is reached to obtain the final K cluster centers. Step 4: Improve GMM in multiple dimensions to construct 4D-GMM, and use the results obtained by 4D-K-Means method to input different category results and matrix Gaussian distribution into 4D-GMM one by one and initialize it, and calculate U in an iterative way. k 、V k The two covariance estimates are repeated until convergence or the number of iterations is reached and the optimal parameters are obtained: w k , C k , U k , V k , where w k is the weight of the kth submatrix, C k ∈R F×T is the mean matrix of the kth submatrix, U k ∈R F×F 、V k ∈R T×T C k is the covariance matrix of the mean; Specifically: Use matrix Gaussian distribution to fit P hwft , construct 4D-GMM, and correspond the clustering results obtained in step 3 to the matrix Gaussian distribution, initialize 4D-GMM and iterate; the iterative steps are as follows: Step.1 Initialize parameter w k , C k , U k , V k : Among them, w k is the weight of the kth submatrix, C k ∈R F×T is the mean matrix of the kth submatrix, U k ∈R F×F 、V k ∈R T×T C k is the covariance matrix of the mean, The final cluster center for the extended K-Means clustering; Step.2 Calculate each pixel P hwft The posterior probability of belonging to the Gaussian distribution of the k-th submatrix and normalized to l hwk : Among them, P(P hwft |C k , U k , V k ) is calculated as shown in formula (11): And w k satisfy: Step.3 Update parameters: Step.4 Repeat Step.2-Step.3 until convergence or the number of iterations is reached, and the optimal parameters are obtained: w k , C k , U k , V k ; Step 5: Use the optimal parameters of step 4 to obtain the two-dimensional matrix P hwft The distribution of P(P hwft |C, U, V) and the final category label L for each pixel hw , where C∈R F×T is the mean matrix, U∈R F×F , V∈R T×T are the covariance matrices with c as the mean; Use the optimal parameters of step 4 to obtain the two-dimensional matrix P hwft The distribution of P(P hwft |C, U, V), as shown in formula (17): Where C∈R F×T is the mean matrix, U∈R F×F , V∈R T×T are the covariance matrices with C as the mean; at the same time, get each pixel P hwft The final category label L hw , as shown in formula (18): When l hwk Maximum, L hw = k, where k = 1, 2, ..., K (18); Step 6: Binarize the obtained pixel category labels, extract the newly added flood area change categories, and obtain the flood change detection results.

Images