CN104200471B

CN104200471B - SAR image change detection based on adaptive weight image co-registration

Info

Publication number: CN104200471B
Application number: CN201410437415.5A
Authority: CN
Inventors: 慕彩红; 焦李成; 张文婷; 熊涛; 侯彪; 刘红英; 马文萍; 陈锋; 吴雪文
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2014-08-30
Filing date: 2014-08-30
Publication date: 2017-03-01
Anticipated expiration: 2034-08-30
Also published as: CN104200471A

Abstract

The invention discloses a kind of SAR image change detection based on adaptive weight image co-registration, mainly solve the problems, such as that in prior art, single type disparity map accuracy of detection is low.Implementation step is：1. read in not two width SAR image of areal in the same time, carry out PPB filtering respectively and obtain filtered image X₁And X₂；2., according to filtered image, obtain differential chart D_sWith log ratio figure D_l；3. respectively to D_sAnd D_lCarry out mean filter, obtain the differential chart D after mean filter_s' with medium filtering after ratio figure D_l'；4. calculate fusion parameters η (i, j)；5. according to fusion parameters to differential chart D_s' and ratio figure D_l' carry out merging generation disparity map D；6. disparity map D is polymerized to two different classes, obtains change-detection result.The present invention has the advantages that simple to operate, noise immunity is good and accuracy of detection is high, can be applicable to environmental monitoring, hazard prediction.

Description

SAR image change detection method based on self-adaptive weight image fusion

Technical Field

The invention belongs to the technical field of image processing, and relates to an SAR image change detection method which can be used in the fields of urban planning, natural disaster assessment, disaster prediction, land utilization, dynamic monitoring of land coverage and the like.

Background

Synthetic aperture radar images, also called SAR images, belong to a kind of microwave radar imaging, compared with visible light and infrared remote sensing, microwave remote sensing has incomparable advantages: firstly, the microwave energy penetrates through cloud, fog, rain and snow, and has all-weather and all-day working capability. Secondly, the microwave has certain penetrating power to the ground objects. And thirdly, the imaging is carried out in a side-looking mode, so that the coverage area is large. Due to these advantages of SAR images, SAR images are increasingly becoming an important data source for change detection. The remote sensing image change detection means that two or more images in different periods of the same region are compared and analyzed, and change information of ground objects is obtained according to the difference between the images. The change detection technology can detect the change information of the radiation value and the local texture of the remote sensing image in different periods, and the characteristic enables the remote sensing image to be widely applied to the aspects of resource and environment monitoring, such as land utilization and coverage change, forest and vegetation change, urban expansion and the like, geographic space data updating, agricultural crop growth monitoring, natural disaster monitoring and evaluation, military field research and the like.

In the change detection research, scholars classify the existing multiple change detection methods into different categories, and the most common of them is the change detection method based on the difference image analysis. The method generally comprises 3 key steps: (1) preprocessing an image; (2) construction of a disparity map; (3) and extracting change information. Wherein:

the difference method and the ratio method are two most basic methods, and the difference image is obtained by subtracting the corrected two time phase remote sensing images pixel by pixel and performing division operation. When the difference value method is used for constructing the difference map, the difference value of 0 or close to 0 is regarded as a constant area, and the difference value of 0 is not regarded as a variable area. In the SAR image, due to the influence of speckle noise, a difference map is mostly constructed by adopting a ratio method, and the ratio method is insensitive to calibration errors. At present, the difference map is mainly constructed by a ratio method in the form of logarithmic ratio or mean ratio. Dekker proposes to use a logarithmic ratio method to construct a difference map, wherein the logarithmic ratio method is to take the logarithm of a ratio method, which can convert multiplicative noise into additive noise and compress the variation range of a ratio image, so that the variation of a high-pixel-value area is weakened, the pixel gray value of the outline of a part of real variation area is closer to the gray value of non-variation pixels, and the preservation of the outline of the variation area and the improvement of the variation detection precision are not facilitated. The mean ratio method is to take the neighborhood mean of the corresponding pixel and then calculate the ratio, which has stronger robustness to noise, but brings too much pseudo-variation information to the background region, namely the unchanged region. Therefore, the change detection using a single type of difference map has problems of low detection accuracy and narrow application range.

Therefore, the fusion of the single type difference maps plays a crucial role in change detection, and the existing methods for fusing the single type difference maps include methods such as parameter weight-based difference map fusion and discrete wavelet fusion. Wherein:

the method for fusing discrete wavelet transform comprises the following steps: and respectively carrying out discrete wavelet decomposition on the average ratio graph and the logarithmic ratio graph to obtain a series of subband images, namely low-frequency components, horizontal high-frequency components, vertical high-frequency components and diagonal high-frequency components, carrying out fusion processing on each decomposition layer on different scales to form wavelet decomposition representation of the fused image, and finally obtaining the fused image through inverse discrete wavelet transformation. The method has the disadvantage that the fusion process is relatively complicated.

The difference map fusion method based on the parameter weight comprises the following processes: obtaining their difference graph and log ratio graph for two images to be fused, extracting difference graph D_sIn which the part with better detail information is reservedAnd logarithmic ratio map D_lThe part with better anti-noise performance in the method is combined with the advantages of the two types of difference graphs and fused according to different weight parameters α, the fused difference graph is D, and the fusion formula is D- α × D_s+(1-α)*D_lThe difference graph and the logarithmic ratio graph are fused according to different weights, the value of α is taken to perform multiple tests from the range of 0-1, the value with the best fusion effect is selected as the value of α to obtain the finally fused difference graph, and the method has the disadvantages that (1) the fusion result cannot be judged to be good or bad without a standard reference graph, the application range is narrow, (2) multiple tests are required to compare the fusion effect when the parameter α is selected, the parameter α cannot be self-adaptive, the parameter needs to be manually adjusted, and the operation is complicated.

Disclosure of Invention

The invention aims to provide an SAR image change detection method based on self-adaptive weight image fusion so as to improve the detection precision and solve the problems that parameters are difficult to select and cannot be self-adaptive in image fusion in the prior art.

In order to achieve the above object, the present invention achieves the steps comprising:

(1) reading in two SAR images I acquired in the same area at different time₁And I₂And for the two SAR images I₁And I₂Respectively carrying out block-based weight probability PPB filtering to obtain two filtered images X₁And X₂；

(2) Calculating the two filtered images X₁And X₂The difference of the corresponding pixel gray value is normalized to obtain a difference value graph D_s；

(3) Calculating the two filtered images X₁And X₂Corresponding to the quotient of the gray value of the pixel, and carrying out normalization to obtain a log ratio graph D_l；

(4) Difference map D_sCarry out 11 × 11 window mean filtering, eliminating difference map D_sObtaining a difference graph D after mean value filtering for the isolated pixel points in the image_s'；

(5) Will log ratio chart D_lPerforming median filtering in 3 × 3 window to suppress the log ratio diagram D_lObtaining a logarithmic ratio figure D after median filtering for the isolated pixel points in the middle_l'；

(6) Respectively calculating difference maps D_s' 3 × 3 neighborhood Ω of each pixel point x (i, j)_xMean value mu of pixel_x(i, j) and variance σ_x(i, j), obtaining fusion parameters:

the fusion parameter η (i, j) varies with the pixel variation, and is used to represent the 3 × 3 neighborhood Ω of the measurement pixel point x (i, j)_xThe homogeneous region refers to a smooth region in the image, the heterogeneous region refers to a noise part or an edge part in the image, and the fusion parameter η (i, j) of the pixel point x (i, j) in the heterogeneous region is larger than the fusion parameter η (i, j) of the pixel point x (i, j) in the homogeneous region;

(7) the mean value filtered difference map D is obtained according to the obtained fusion parameters η (i, j)_s' logarithmic ratio map D after median filtering_lPerforming image fusion to obtain a fused difference map D, wherein pixel points of each horizontal and vertical coordinate at i, j in the difference map D are as follows:

D(i,j)＝η(i,j)×D_s′(i,j)+(1-η(i,j))×D_l′(i,j)，

wherein, D'_s(i, j) is a point, D ', of which the horizontal and vertical coordinates are i, j respectively in the mean value filtered difference value diagram'_l(i, j) are points of which the horizontal and vertical coordinates are i, j respectively in the logarithmic ratio value chart after median filtering;

(8) clustering the fused difference graph D into two different classes by using a k-means clustering method, respectively calculating mean values of the two different classes, defining the class with the larger mean value as a variable class and the class with the smaller mean value as a non-variable class, and obtaining a final change detection result.

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

1. difference value map D 'after mean value filtering'_s(i, j) and the median filtered logarithmic ratio map D'_l(i, j) utilizes D 'at the time of fusion'_s(i, j) each pixel 3 x 3 neighborhood Ω_xVariance σ of_x(i, j) and the mean value μ_xThe ratio η (i, j) of (i, j) is used as a fusion parameter, so that the original fusion parameter can be automatically adjusted according to different pixel points, the fusion parameter does not need to be manually adjusted through multiple tests, and the operation is simple.

2. The fusion parameter eta (i, j) adopted by the invention considers the gray information and the spatial information of the neighborhood pixels, improves the anti-noise performance and the capability of retaining the image details, and further improves the change detection precision.

Simulation experiments show that the method can obtain satisfactory change detection results for a plurality of SAR images and has better robustness.

Drawings

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

FIG. 2 is a diagram of a yellow river SAR image dataset used in the simulation of the present invention;

FIG. 3 is a graph of the results of a simulation experiment conducted on FIG. 2 using the present invention and a comparative method;

FIG. 4 is a diagram of a Wttawa region SAR image dataset used in simulation according to the present invention;

FIG. 5 is a graph of the results of a simulation experiment of FIG. 4 with the present invention and comparative method;

FIG. 6 is a plot of a Burnely SAR image dataset used in the simulation of the present invention;

FIG. 7 is a graph of the results of a simulation experiment of FIG. 6 with the present invention and comparative method.

Detailed Description

Referring to fig. 1, the method comprises the following specific steps:

step 1: reading in two SAR images I which are obtained in the same area at different time₁And I₂。

Step 2: for two SAR images I₁And I₂Respectively carrying out PPB filtering to obtain filtered images X₁And X₂。

For the first SAR image I₁Carrying out fundamental patch-based weights, namely PPB filtering, and obtaining a filtered image X₁；

For the second SAR image I₂PPB filtering is carried out to obtain a filtered image X₂；

The PPB filtering is a filtering method based on a noise distribution model and having a more general statistical similarity criterion. When PPB filtering is carried out, the process of removing noise is expressed as a weighted maximum likelihood estimation problem, the weights are updated based on the similarity of a noise block and the similarity of a previously calculated block in an iterative process, and the iterative process improves the performance of noise removal.

And step 3: obtaining a difference image D of two SAR images_s。

3a) For the first SAR image X₁Gray value X at coordinate (i, j)₁(i, j) and a second SAR image X₂Gray value X at coordinate (i, j)₂(i, j) performing difference operation to obtain the difference value of the difference matrix D at the coordinate (i, j): d (i, j) ═ X₁(i,j)-X₂(i, j) |, and then a difference matrix is obtained: d ═ D (i, j) };

3b) normalizing the difference matrix to obtain a difference map D_sValue D at coordinate (i, j)_s(i,j)：

Further obtaining a difference map: d_s＝{D_s(i,j)}。

And 4, step 4: calculating the two filtered images X₁And X₂Corresponding to the quotient of the gray value of the pixel, and carrying out normalization to obtain a log ratio graph D_l。

4a) Calculating and calculating the two filtered images X by the following formula₁And X₂Quotient R (i, j) for the gray value of the pixel:

wherein X₁(i, j) is a first SAR image X₁Pixel points, X, at abscissa and ordinate i, j₂(i, j) is a second SAR image X₂Pixel points at the positions of the horizontal and vertical coordinates i, j;

4b) normalizing the quotient of the pixel gray values by the following formula to obtain a result of the quotient normalization of each pixel gray value:

wherein R ═ { R (i, j) };

4c) obtaining a log ratio map by using the result of quotient normalization of each pixel gray value: d_l＝{D_l(i,j)}。

And 5: difference map D_sCarrying out mean value filtering of 11 × 11 window to eliminate difference value image D_sThe isolated pixel points in the image are more complete in area, the local area keeps better continuity, the image is smoother, and finally the difference image D after mean filtering is obtained_s'。

Step 6: will log ratio chart D_lPerforming median filtering in 3 × 3 window to suppress the log ratio diagram D_lThe isolated pixel points in the middle are used for reserving the edge information to finally obtain a logarithmic ratio chart D after median filtering_l'。

And 7: the fusion parameters η (i, j) are calculated.

7a) Calculating a difference map D using the following formula_s' 3 × 3 neighborhood Ω of each pixel point x (i, j)_xMean value mu of pixel_x(i,j)

Wherein,is 3 x 3 neighborhood omega_xThe sum of the gray values of each pixel x (i, j);

7b) according to the mean value mu_x(i, j), calculating a difference map D using the following formula_s' 3 × 3 neighborhood Ω of each pixel point x (i, j)_xVariance σ of pixel_x(i,j)：

Wherein,is 3 x 3 neighborhood omega_xSubtracting the neighborhood omega from the gray value of each pixel point x (i, j)_xInternal mean value mu_x(i, j) the sum of the squares of the results;

7c) according to the mean value mu_x(i, j) and variance σ_x(i, j) calculating fusion parameters η (i, j):

the fusion parameter η (i, j) varies with the pixel, and its variation range is between 0 and 1, which is used to represent the 3 × 3 neighborhood Ω of the measurement pixel point x (i, j)_xThe homogeneous region refers to a smooth region in the image, the heterogeneous region refers to a noise part or an edge part in the image, and the fusion parameter η (i, j) of the pixel point x (i, j) in the heterogeneous region is larger than the fusion parameter η (i, j) of the pixel point x (i, j) in the homogeneous region.

And 8: fusion mean filtered difference map D_s' logarithmic ratio map D after median filtering_l', obtaining a fused difference map D,

according to the fusion parameters η (i, j) obtained in step 7, the difference map D after mean value filtering_s' logarithmic ratio map D after median filtering_lPerforming image fusion to obtain a fused difference map D, wherein pixel points of each horizontal and vertical coordinate at i, j in the difference map D are as follows:

D(i,j)＝η(i,j)×D_s'(i,j)+(1-η(i,j))×D_l' (i, j), and obtaining a fused difference map D ═ D (i, j) },

wherein, D'_s(i, j) is a point, D ', of which the horizontal and vertical coordinates are i, j respectively in the mean value filtered difference value diagram'_l(i, j) is the point whose horizontal and vertical coordinates are i, j respectively in the logarithmic ratio diagram after median filtering, and since the coefficient η (i, j) is a parameter for measuring the heterogeneity of the neighborhood region x (i, j), when the value η (i, j) is small, x (i, j) is in the homogeneous region, D_l' the difference map is generated to be dominant, so that speckle noise of the image is better suppressed, x (i, j) is in a heterogeneous region when η (i, j) is large, and D_s' dominate when generating disparity maps, better preserving image detail information.

Step 9: and carrying out k-means clustering on the difference graph D to obtain a final change detection result.

9a) Randomly selecting 2 pixel points from the difference graph D as an initial clustering center;

9b) respectively calculating the distance between each pixel point x (i, j) in the difference graph D and two clustering centers, and dividing the corresponding pixel points again according to the minimum distance to classify x (i, j) into the class closest to the pixel points;

9c) recalculating the mean value of each class as a new clustering center;

9d) repeating the steps 9b) -9c) until each cluster is not changed any more, and obtaining the mean value of the two categories;

9e) and converting the pixels with small mean values into 0, namely the unchanged class, converting the pixels with large mean values into 255, namely the changed class, and outputting the obtained result graph to obtain the final change detection result.

The effect of the present invention can be further illustrated by the following simulation results:

1. conditions of the experiment

The experimental environment is as follows: windows XP, SPI, CPU Pentium (R)4, basic frequency 2.4GHZ, and software platform MatlabR2010 a.

The first data set was a partial cut-out of the SAR image data set of the yellow river, as shown in fig. 2(a) and 2(b), acquired in eastern province of the mountains at a resolution of 8m by 8 m. In this case, 2(a) is an image of 2008. 6 month, and 2(b) is an image of 2009. 6 month. The image size is 257 × 289 pixels each, and the image gray scale is 256. Fig. 2(c) is a diagram of the standard variation detection result, which includes 13423 variation pixels and 60841 non-variation pixels.

The second data set is a wotta region SAR image data set, which is composed of two Radarsat SAR images captured at different times, and the variation area is mainly caused by flood disasters, as shown in fig. 4(a) and 4(b), wherein fig. 4(a) is an image of 5 months 1997, fig. 4(b) is an image of 8 months 1997, the image sizes are 290 × 350 pixels, and the gray level is 256. FIG. 4(c) is a standard change detection graph of the data set, comprising 85451 unchanged pixels and 16049 changed pixels.

The third data set is a bernly SAR image data set, which is composed of two Radarsat SAR images taken at different times, as shown in fig. 6(a) and 6(b), where fig. 6(a) is an image of 4 months 1999, fig. 6(b) is an image of 5 months 1999, the image sizes are all 301 × 301 pixels, and the gray level is 256. Fig. 6(c) is a graph of the standard variation detection result of the data set, which includes 1269 changed pixels and 89332 unchanged pixels.

2. Evaluation index of experiment

The evaluation indexes used in the experiment are the number of missed detections FN, the number of false detections FP, the total number of errors OE, the correct detection rate PCC and the Kappa coefficient. The missing detection number is the sum of pixels which are not detected and actually change, the false detection number is the sum of pixels which are not detected and actually change and are detected as changes, the total error number is the sum of the missing detection number and the false detection number, and the correct detection rate is 1-total error number/total pixel number. The total number of pixels is N, the total number of changed points is Nc, the total number of unchanged points is Nu, the correctly detected changed points are TP, TP is total changed points-missing detection number, the correctly detected unchanged points are TN, TN is total unchanged points-false detection number, and the order is thatTo obtainIf the variation graph and the reference graph are completely consistent, the Kappa coefficient is 1, otherwise, the Kappa coefficient is 0. The Kappa coefficient contains more classified detail information and is therefore a more accurate evaluation criterion than PCC.

3. Content of the experiment

The method of the invention and the existing 4 kinds of change detection methods are used for detecting the change of 3 image data sets, and the 4 kinds of comparison methods are respectively as follows: image X filtered with PPB₁And X₂Difference operation is carried out to construct a difference map, and then the difference map is recorded by a K-means clustering methodA difference method; image X filtered with PPB₁And X₂Carrying out logarithmic ratio operation to construct a difference graph, and then marking the difference graph as a logarithmic ratio method by using a K-means clustering method; image X filtered with PPB₁And X₂Carrying out mean ratio operation to construct a difference graph, and then recording the difference graph as a mean ratio method by using a K-means clustering method; the method proposed by Yaoguozheng in the article "Using Combined Difference Image and K-means Cluster for SAR Image Change Detection" is referred to as the CDI-K method.

Experiment 1

The change detection was performed on the first data set of fig. 2(a) and 2(b) using the method of the present invention and 4 conventional change detection methods, 2(c) is a reference graph, and the experimental results are shown in fig. 3, where 3(a) is a graph of the change detection result by the difference method, 3(b) is a graph of the change detection result by the logarithmic ratio method, 3(c) is a graph of the change detection result by the mean ratio method, 3(d) is a graph of the change detection result by the CDI-K method, and 3(e) is a graph of the change detection result by the present invention.

Experiment 2

The second data set, FIGS. 4(a) and 4(b), was subjected to change detection using the method of the present invention and 4 conventional change detection methods, 4(c) is a reference graph, and the experimental results are shown in FIG. 5, where 5(a) is a graph of the change detection results of the difference method, 5(b) is a graph of the change detection results of the logarithmic ratio method, 5(c) is a graph of the change detection results of the mean ratio method, 5(d) is a graph of the change detection results of the CDI-K method, and 5(e) is a graph of the change detection results of the present invention.

Experiment 3

The third data set, FIGS. 6(a) and 6(b), was subjected to change detection using the method of the present invention and 4 conventional change detection methods, 6(c) is a reference graph, and the experimental results are shown in FIG. 7, where 7(a) is a graph of the change detection results of the difference method, 7(b) is a graph of the change detection results of the logarithmic ratio method, 7(c) is a graph of the change detection results of the mean ratio method, 7(d) is a graph of the change detection results of the CDI-K method, and 7(e) is a graph of the change detection results of the present invention.

4. Analysis of Experimental results

As can be seen from fig. 3, 5, and 7, the conventional single-type difference map has poor noise immunity, large total error number, low detection accuracy, and poor visual effect. The CDI-K method utilizes the parameter alpha to fuse the difference map after mean filtering and the logarithm ratio map after median filtering to generate the difference map, but the parameter is not adaptive in the fusion process, the parameter alpha needs to select the optimal result of change detection from 0 to 1 through a plurality of experiments to determine the value of the parameter alpha, the process needs to manually adjust the parameter and is complicated, only the gray information of pixels is considered, and the spatial information of neighborhoods is not considered, so the anti-noise performance is poor, and the detection precision is low. The invention enables the parameters to be self-adaptive when generating the difference map, does not need to manually adjust the parameters, greatly simplifies the experimental steps, ensures that the method is simple and efficient, considers the pixel gray information and the neighborhood space information when calculating the parameters, enhances the anti-noise capability, reduces the total error number and improves the precision of change detection.

The results of statistics of the change detection values of 3 data sets by the method of the present invention and the existing 4 change detection methods are shown in table 1.

TABLE 1 comparison of Change detection values for the present invention and the existing four methods

As can be seen from table 1, the total error number of the present invention is less for all of the three SAR images, and the PCC and Kappa coefficients are better than for the other methods.

Claims

1. A SAR image change detection method based on self-adaptive weight image fusion comprises the following steps:

(1) reading in two SAR images I acquired in the same area at different time₁And I₂And for the two SAR images I₁And I₂Respectively carrying out weight probability filtering based on blocks to obtain two filtered images X₁And X₂；

(2) Calculating the two filtered images X₁And X₂The difference of the corresponding pixel gray values is normalized to obtain a difference valueDrawing D_s；

(4) Difference map D_sCarrying out mean value filtering of 11 × 11 window to eliminate difference value image D_sObtaining a difference image D after mean value filtering for the noise pixel points in the image_s'；

(6) Respectively calculating difference maps D_s' 3 × 3 neighborhood Ω of each pixel point x (i, j)_xMean value mu of pixel gray values_x(i, j) and variance σ_x(i, j), obtaining fusion parameters:

the fusion parameter η (i, j) varies with the pixel, and its variation range is between 0 and 1, which is used to represent the 3 × 3 neighborhood Ω of the measurement pixel point x (i, j)_xThe homogeneous region refers to a smooth region in the image, the heterogeneous region refers to a noise part or an edge part in the image, and the fusion parameter η (i, j) of the pixel point x (i, j) in the heterogeneous region is larger than the fusion parameter η (i, j) of the pixel point x (i, j) in the homogeneous region;

D(i,j)＝η(i,j)×D_s′(i,j)+(1-η(i,j))×D_l′(i,j)，

wherein, D'_s(i, j) is a point whose horizontal and vertical coordinates are i, j respectively in the difference diagram after mean filtering, D_l' (i, j) is the horizontal and vertical coordinate of the logarithmic ratio chart after median filteringPoints other than i, j;

2. The method as set forth in claim 1, wherein said step (3) is performed as follows:

3a) calculating and calculating the two filtered images X by the following formula₁And X₂Quotient R (i, j) for the gray value of the pixel:

wherein X₁(i, j) is image X₁Pixel points, X, at abscissa and ordinate i, j₂(i, j) is image X₂Pixel points at the positions of the horizontal and vertical coordinates i, j;

3b) normalizing the quotient of the pixel gray values by the following formula to obtain a result of the quotient normalization of each pixel gray value:

wherein R ═ { R (i, j) },

3c) the result of the quotient normalization of the gray value of each pixel yields a log ratio map: d_l＝{D_l(i,j)}。

3. The method as claimed in claim 1, wherein said step (6) of calculating a difference map D_s' 3 × 3 neighborhood Ω of each pixel point x (i, j)_xMean value mu of pixel_x(i, j) and variance σ_x(i, j), the calculation formula is as follows:

wherein,is 3 x 3 neighborhood omega_xThe sum of the gray values of each pixel x (i, j),

is 3 x 3 neighborhood omega_xSubtracting the neighborhood omega from the gray value of each pixel point x (i, j)_xInternal mean value mu_x(i, j) sum of squares of the results.

4. The method according to claim 1, wherein the step (8) of clustering the fused difference maps D into two different classes by using a k-means clustering method comprises the following steps:

8a) randomly selecting 2 pixel points from the difference graph D as an initial clustering center;

8b) respectively calculating the distance between each pixel point x (i, j) in the difference graph D and two clustering centers, and dividing the corresponding pixel points again according to the minimum distance to classify x (i, j) into the class closest to the pixel points;

8c) recalculating the mean value of each class as a new clustering center;

8d) repeating steps 8b) -8c) until no more changes occur to each cluster.