A multi-method fusion SAR image shadow detection method and system

Abstract

The application discloses a kind of multi-method fusion's SAR image shadow detection method and system, including according to wavelet algorithm, Markov random field and region growing algorithm, process original image set, obtain a scene cutting image set and two scene cutting image set;According to amplitude detection method and coherence coefficient detection method, process a scene cutting image set and two scene cutting image set, obtain first shadow area set, first shadow pending area set, second shadow area set and second shadow pending area set;According to four area sets, set operation obtains first partial shadow area set and third shadow pending area set;According to radar line of sight and digital elevation model analysis algorithm, process third shadow pending area set, obtain second partial shadow area set;According to first partial shadow area set and second partial shadow area set, and set obtains final detection result.This application improves the accuracy and reliability of shadow detection by fusing multiple algorithms.

Description

Technical Field

[0001] The present invention belongs to the technical field of computer vision, and in particular relates to a multi-method fusion SAR image shadow detection method and system. Background Art

[0002] Spaceborne synthetic aperture radar interferometry (InSAR) is an advanced space geodetic technique that calculates phase information between multiple coherent radar images to obtain information such as subtle surface deformation and elevation. However, because synthetic aperture radar (SAR) is an active, side-looking imaging technology, it is prone to shadowing and overlapping in areas with large terrain undulations. Shadows can lead to information loss in SAR images and also hinder object recognition and phase unwrapping in InSAR data processing. Therefore, accurately detecting and marking shadow areas in SAR images is crucial to ensuring data quality and the accuracy of subsequent processing.

[0003] However, existing shadow detection is often performed together with overlap detection. The complete single shadow detection technology has not received enough attention in previous research, and each shadow detection algorithm has its own shortcomings. Summary of the Invention

[0004] In response to the shortcomings of the existing technology, this application provides a multi-method fusion SAR image shadow detection method and system, which integrates multiple algorithms to solve the problem of accurate detection of shadow areas in SAR images and improve the accuracy and reliability of shadow detection.

[0005] To achieve the above objectives, in a first aspect, the present invention provides a multi-method fusion SAR image shadow detection method, comprising:

[0006] According to the preset wavelet multi-scale bilateral filtering algorithm, the original two-view image set is processed to obtain the two-view image set;

[0007] Processing the two image sets according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set;

[0008] Processing the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set;

[0009] Processing the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set;

[0010] According to the first shadow area set, the first shadow area set to be determined, the second shadow area set and the second shadow area set to be determined, a set operation is performed to obtain a first partial shadow area set and a third shadow area set to be determined;

[0011] Processing the third set of undetermined shadow areas according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second set of shadow areas;

[0012] The final shadow area detection result is obtained by performing a union calculation based on the first partial shadow area set and the second partial shadow area set.

[0013] Compared with the prior art, the embodiments of the present application have the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, noise can be effectively removed while retaining important features such as edges; bilateral filtering can smooth the image while maintaining edges, thereby improving image quality; a Markov random field model is constructed to fully utilize the spatial dependence between pixels to improve segmentation accuracy; a region growing algorithm is used to more accurately segment the target area by selecting seed points and gradually expanding the region; an amplitude adaptive threshold detection method is adopted to adapt to the image characteristics in different scenarios and improve the robustness of detection; coherence calculation can effectively identify shadow areas caused by physical characteristics, especially those caused by terrain or object occlusion; the results obtained by multiple methods are subjected to set operations to combine the advantages of various methods, reduce false detections and missed detections, and improve the overall detection effect; through radar line of sight and digital elevation model analysis, shadow areas caused by insufficient radar depression angle are further detected, thereby improving the detection accuracy in areas with large terrain undulations; finally, a multi-method fusion SAR image shadow detection result is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0014] In some embodiments of the first aspect of the present application, processing the original two-view image set according to a preset wavelet multiscale bilateral filtering algorithm to obtain the two-view image set includes:

[0015] According to a preset wavelet transform decomposition algorithm, the original two images in the original two-view image set are processed respectively to obtain each first high-frequency sub-band and first low-frequency sub-band;

[0016] Processing each of the first high frequency sub-bands according to a preset bilateral filtering algorithm to obtain each of the second high frequency sub-bands;

[0017] Processing each of the first low-frequency sub-bands according to a preset smoothing algorithm to obtain each of the second low-frequency sub-bands;

[0018] The second high-frequency sub-band and the second low-frequency sub-band are processed according to a preset inverse wavelet transform reconstruction algorithm to obtain two images, that is, a two-image set.

[0019] Compared with the existing technology, the above embodiment has the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, noise can be effectively removed while retaining important features such as edges. The bilateral filtering process can smooth the image while maintaining the edges, thereby improving the image quality.

[0020] In some embodiments of the first aspect of the present application, processing the two-view image sets according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set includes:

[0021] According to the normalization algorithm, the one-view image set and the two-view image set are processed respectively to obtain the one-view normalized image set and the two-view normalized image set;

[0022] Constructing Markov random field models based on the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a model of a label field;

[0023] According to a preset optimization algorithm, each of the Markov random field models is processed respectively to obtain a one-scene preliminary segmentation image set and a two-scene preliminary segmentation image set;

[0024] According to a preset automatic region growing algorithm, the first scene preliminary segmentation image set and the second scene preliminary segmentation image set are processed respectively to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set;

[0025] According to a preset morphological algorithm, the one-scene secondary segmentation image set and the two-scene secondary segmentation image set are processed respectively to obtain a one-scene cut image set and a two-scene cut image set.

[0026] Compared with the existing technology, the above embodiment has the following beneficial effects: the Markov random field model can well capture the spatial correlation between pixels in the image, perform global optimization, effectively reduce noise interference, and provide fine preliminary segmentation results; region growing can further refine the segmentation area, especially in the edge or transition area, to ensure accurate division of shadow and non-shadow areas; morphological processing can eliminate isolated noise points in the segmentation results, smooth the boundaries, and improve segmentation accuracy and completeness.

[0027] In some embodiments of the first aspect of the present application, a Markov random field model is constructed based on the one-scene normalized image set and the two-scene normalized image set, respectively; wherein the Markov random field model includes a model of a label field, including:

[0028] A Markov random field is defined, each pixel of each image in the one-view normalized image set and the two-view normalized image set is regarded as a random variable, the grayscale distribution in the image is modeled according to a finite Gaussian mixture model, and an original Markov random field model is constructed;

[0029] Optimizing the parameters of the original Markov random field model using an expectation-maximization algorithm;

[0030] The Markov random field model is obtained by adopting Gibbs distribution as the prior probability distribution of the marking field and adopting Porter model as the model of the marking field.

[0031] Compared with the existing technology, the above embodiment has the following beneficial effects: defining a Markov random field model, representing the spatial relationship between pixels by establishing a neighborhood system, and constructing a probability distribution using the Gibbs distribution, which is beneficial to capturing the spatial correlation between pixels in the image and facilitating subsequent global optimization.

[0032] In some embodiments of the first aspect of the present application, the Markov random field models are processed separately according to a preset optimization algorithm to obtain a one-scene preliminary segmented image set and a two-scene preliminary segmented image set, including:

[0033] Adaptive weights are introduced, where the weight estimation model is as follows:

[0034]

[0035]

[0036] Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and y i Represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, N i represents the set of neighboring pixels of pixel i;

[0037] Introducing penalty function The definition is as follows:

[0038]

[0039] Marking field energy function U(x i )for: Where β is a custom weighting coefficient ranging from 0 to 1, V2(x i ,x l ) represents the weight of the corresponding edge between pixels, and exp(*) represents the exponential operation;

[0040] The final energy function expression of the posterior probability solved pixel by pixel is:

[0041]

[0042] According to the preset number of classifications and the K-means algorithm, the images in the one-scene normalized image set and the two-scene normalized image set are respectively subjected to initial segmentation to obtain initial labels of each pixel;

[0043] According to the preset number of iterations, each Markov random field model is iteratively optimized to obtain a one-scene preliminary segmentation image set and a two-scene preliminary segmentation image set, and at each iteration, each model parameter of the Markov random field model is estimated according to a preset expectation algorithm, and a maximum likelihood parameter estimation of each category model is calculated according to a preset iteration algorithm, until a preset convergence criterion is met, and the model parameters of each category calculated by the expectation maximum algorithm at the current iteration number are substituted into the final energy function expression to update the segmentation result of each pixel.

[0044] Compared with the prior art, the above embodiment has the following beneficial effects: the Potter model as the model of the label field realizes global optimal segmentation through iterative optimization, effectively reduces noise interference, and provides a fine preliminary segmentation result.

[0045] In some embodiments of the first aspect of the application, the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set are respectively processed according to a preset automatic region growing algorithm to obtain a one-scene secondary segmentation image set and a two-scene secondary segmentation image set, including:

[0046] According to a preset rule, a plurality of seed points are selected from the images in the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set, respectively;

[0047] According to a preset initial window size, the window size is gradually increased, and the texture difference of each pixel point under each window size is calculated in sequence, and the window size with the minimum texture difference is taken as the best window size;

[0048] According to a preset initial window, the probability distribution of each pixel point belonging to each texture category is calculated, and the entropy of the initial window is calculated according to the probability distribution;

[0049] According to the best threshold and the best window size, the texture difference between each pixel point in the best window of each seed point and the corresponding seed point is calculated, the pixel points with the texture difference meeting the best threshold are screened out, and are merged into the initial region of the corresponding seed point until all the pixel points are screened out, to obtain a one-scene secondary segmentation image set and a two-scene secondary segmentation image set, respectively; wherein the texture difference is obtained by calculating the difference of the gray level co-occurrence value of each pixel point in the window, and the best threshold is calculated by the standard deviation of each initial window category.

[0050] Compared with the prior art, the above embodiment has the following beneficial effects: using the automatic region growing method, by selecting seed points and gradually expanding the region, the target area can be segmented more accurately, especially in the edge or transition area, ensuring accurate division of shadow and non-shadow areas.

[0051] In some embodiments of the first aspect of the present application, processing the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set includes:

[0052] Obtaining the original amplitude image corresponding to each image in the scene cut image set;

[0053] compressing the original amplitude image according to a grayscale compression algorithm to obtain an amplitude image;

[0054] Repeatedly dividing the amplitude image into shadow areas and non-shadow areas according to a preset initial threshold, and after each division, calculating the mean and standard deviation of the shadow areas and non-shadow areas, and updating the initial threshold according to the mean and standard deviation until the updated change of the initial threshold meets the preset difference range, thereby obtaining an optimal threshold;

[0055] The amplitude image is processed according to the optimal threshold to obtain a first shadow area set and a first shadow undetermined area set.

[0056] Compared with the existing technology, the above embodiment has the following beneficial effects: by calculating the mean and variance of the shadow area and the background area, using an iterative algorithm to optimize and obtain the optimal threshold, and using the optimal threshold to separate the shadow area and the non-shadow area, when the difference in grayscale values ​​between the shadow and non-shadow areas is small, the false detection rate can be effectively reduced and the accuracy of shadow detection can be improved.

[0057] In some embodiments of the first aspect of the present application, the step of processing the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set includes:

[0058] According to the one-scene cut image set and the two-scene cut image set, the coherence value of each two-scene image after the phase correction factor is introduced is calculated as follows:

[0059] Where z1, z2 represent the complex values ​​of the two images; * represents the conjugate operator; i represents the pixel number; represents the phase correction factor, Indicates the phase caused by subtracting the terrain factor from the original phase δ flatteneddenotes the coherence value after introducing the phase correction factor, LL is the number of pixels in the window, and j is an imaginary unit;

[0060] The coherence value after introducing the phase correction factor is processed according to a preset spatial average algorithm to obtain a final coherence value;

[0061] The first and second sets of segmented images are processed according to the final coherence value to obtain a second set of shadow regions and a second set of shadow regions to be determined.

[0062] Compared with the prior art, the above embodiment has the following beneficial effects: the coherence coefficients of two images are calculated, the phase correction factor and the spatial neighborhood average value are used to improve the unbiasedness of coherence estimation, the coherence threshold is set, and the shadow region is extracted, so that the shadow region with low coherence can be more accurately identified, and the deviation caused by the terrain or noise is reduced.

[0063] In some embodiments of the first aspect of the application, processing the third set of shadow regions to be determined according to the preset radar line of sight and digital elevation model analysis algorithm to obtain a second set of shadow regions, comprises:

[0064] The preset radar line of sight and digital elevation model analysis algorithm comprises a geometric relationship:

[0065] H / M≥ksinβ1; 0<k≤1; wherein M, H, and β1 are the slant distance, the height difference, and the radar depression angle corresponding to two pixel points, respectively;

[0066] According to the geometric relationship, the images in the third set of shadow regions to be determined are processed, and the pixel points meeting the geometric relationship are classified as shadow regions to obtain a second set of shadow regions.

[0067] Compared with the prior art, the above embodiment has the following beneficial effects: the geometric relationship of the ground objects in the image is analyzed in combination with the radar line of sight and digital elevation model data, the shadow region is determined by calculating the relationship between the ground height difference and the slant distance, and the shadow region caused by insufficient radar depression angle is further detected, thereby improving the detection accuracy in the area with large terrain undulations.

[0068] In the second aspect, the embodiments of the application also provide a SAR image shadow detection system based on multi-method fusion, comprising: a preprocessing module, a two-step segmentation module, an amplitude detection module, a coherence detection module, a set operation module, a radar elevation detection module, and a result output module;

[0069] The preprocessing module is configured to process an original two-set of images according to a preset wavelet multi-scale bilateral filtering algorithm to obtain a two-set of images.

[0070] The two-step segmentation module is used to process the two-scene image set according to a preset Markov random field and region growing algorithm to obtain a one-scene cut image set and a two-scene cut image set;

[0071] The amplitude detection module is used to process the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set;

[0072] The optical coherence detection module is used to process the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set;

[0073] The set operation module is configured to perform set operation on the first shadow area set, the first shadow area set to be determined, the second shadow area set, and the second shadow area set to be determined to obtain a first partial shadow area set and a third shadow area set to be determined;

[0074] The radar elevation detection module is used to process the third set of shadow undetermined areas according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second set of shadow areas;

[0075] The result output module is used to calculate the final shadow area detection result based on the union of the first partial shadow area set and the second partial shadow area set.

[0076] Compared with the prior art, the embodiments of the present application have the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, noise can be effectively removed while retaining important features such as edges; bilateral filtering can smooth the image while maintaining edges, thereby improving image quality; a Markov random field model is constructed to fully utilize the spatial dependence between pixels to improve segmentation accuracy; a region growing algorithm is used to more accurately segment the target area by selecting seed points and gradually expanding the region; an amplitude adaptive threshold detection method is adopted to adapt to the image characteristics in different scenarios and improve the robustness of detection; coherence calculation can effectively identify shadow areas caused by physical characteristics, especially those caused by terrain or object occlusion; the results obtained by multiple methods are subjected to set operations to combine the advantages of various methods, reduce false detections and missed detections, and improve the overall detection effect; through radar line of sight and digital elevation model analysis, shadow areas caused by insufficient radar depression angle are further detected, thereby improving the detection accuracy in areas with large terrain undulations; finally, a multi-method fusion SAR image shadow detection result is output, which significantly improves the accuracy and comprehensiveness of shadow detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1: A flowchart of a multi-method fusion SAR image shadow detection method provided in some embodiments of the present invention.

[0078] Figure 2 : A structural diagram of a multi-method fusion SAR image shadow detection system provided in some embodiments of the present invention.

[0079] Figure 3 : A schematic flow chart of a wavelet multi-scale bilateral filtering algorithm provided in some embodiments of the present invention.

[0080] Figure 4 : A flow chart of a Markov random field and region growing algorithm provided in some embodiments of the present invention.

[0081] Figure 5 : A schematic diagram of the geometric relationship between a radar line of sight and a digital elevation model analysis algorithm provided in some embodiments of the present invention. DETAILED DESCRIPTION

[0082] 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. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0083] Example 1:

[0084] Please refer to Figure 1 , a multi-method fusion SAR image shadow detection method provided by an embodiment of the present invention, including steps S1 to S7:

[0085] Step S1: Processing the original two-view image set according to a preset wavelet multi-scale bilateral filtering algorithm to obtain a two-view image set.

[0086] In some embodiments of the present application, Figure 3 The flowchart of a wavelet multi-scale bilateral filtering algorithm is shown in FIG. 1 , wherein step S1 can be implemented by the following preferred implementation:

[0087] Step S11: processing the original two-view image in the original two-view image set according to a preset wavelet transform decomposition algorithm to obtain first high-frequency sub-bands and first low-frequency sub-bands.

[0088] In this step, the wavelet function of the wavelet transform decomposition algorithm The following conditions should be met:

[0089] in, is the basis function or mother wavelet, which can form a cluster of function systems through scaling and translation on the time axis: in, It is called sub-wavelet, a is the scale factor, which reflects the period length of the wavelet, b is the translation factor, which reflects the time translation, R represents a real number, L 2 (R) represents real two-dimensional space.

[0090] In this step, the wavelet transform function may be selected from Haar wavelet, Dobecy wavelet, Mexican hat wavelet, Mollet wavelet and Mayer wavelet, or other similar wavelet functions may be selected for implementation, without further limitation here.

[0091] Preferably, the calculation formula of wavelet transform can be as follows:

[0092] Among them, w f (a,b) represents continuous wavelet transform, f(t) represents the function to be transformed by wavelet, a represents the scale factor, b represents the time translation factor, Represents a wavelet transform function that has undergone time translation and scale expansion, also known as a sub-wavelet.

[0093] Step S12: processing each of the first high frequency sub-bands respectively according to a preset bilateral filtering algorithm to obtain each of the second high frequency sub-bands.

[0094] Bilateral filtering can preserve image edges while smoothing noise. The calculation formula of bilateral filtering can be as follows:

[0095] Where p and q represent the coordinates of two pixels, i.e. p = (x0, y0), q = (x1, y1), BF[I] P represents the bilateral filtering centered at pixel p, I p , I q Represents the pixel value at position p, q, that is, the intensity value of the pixel, represents Gaussian weighting in the spatial domain, Represents Gaussian weighting within the pixel value range, W p is a normalization factor used to ensure that the weighted sum of the pixels is 1, W p The expression is as follows: Q represents the set of all pixels in a local neighborhood centered on pixel p.

[0096] Step S13: processing each of the first low-frequency sub-bands according to a preset smoothing algorithm to obtain each of the second low-frequency sub-bands.

[0097] Among them, the smoothing process can remove noise while retaining the main structure of the image. The smoothing process algorithm can be mean filtering: Among them, M represents a pixel, Num represents the total number of pixels, and f(i,j) represents the pixel value.

[0098] Preferably, when designing a linear smoothing filter, it should be ensured that the filter has a single peak, ie, a main lobe, and is symmetrical in both vertical and horizontal directions. For example, the window of a typical smoothing filter is 3×3.

[0099] Preferably, the smoothing algorithm may further use a Gaussian filter to remove noise that obeys a normal distribution; wherein a first-order zero-mean Gaussian function may be expressed as: The two-dimensional discrete zero-mean Gaussian function can be expressed as: The distribution exponent σ represents the width of the Gaussian distribution.

[0100] Step S14: processing the second high-frequency sub-band and the second low-frequency sub-band according to a preset inverse wavelet transform reconstruction algorithm to obtain two images, that is, a two-image set.

[0101] In step S1 of this embodiment, the advantages of wavelet transform and bilateral filtering are combined, which can remove noise while retaining image details, and is suitable for application scenarios requiring high-quality SAR images.

[0102] Step S2: processing the two-view image sets according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set.

[0103] In some embodiments of the present application, Figure 4 The flowchart of a Markov random field and region growing algorithm is shown in FIG. 1 , wherein step S2 can be implemented by the following preferred implementation:

[0104] Step S21: According to a normalization algorithm, the one-view image set and the two-view image set are processed respectively to obtain a one-view normalized image set and a two-view normalized image set.

[0105] Preferably, the normalization algorithm may be as follows:

[0106] Among them, I represents the pixel value in the original image; Min and Max represent the maximum and minimum values ​​of the pixels in the original image, and I N Indicates the normalized value.

[0107] Step S22: constructing a Markov random field model based on the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a label field model.

[0108] Furthermore, the construction of the Markov random field model can be achieved through the following preferred implementations:

[0109] Define a Markov random field, treat each pixel of each image in the one-view normalized image set and the two-view normalized image set as a random variable, model the grayscale distribution in the image according to the finite Gaussian mixture model, and construct an original Markov random field model, where the probability distribution of each pixel in the image can be expressed as:

[0110]

[0111] Where S represents the set of random fields, i represents the subset of random fields; P(Y|X=x) represents the probability distribution function or likelihood function of the feature field when the label field is known; y i represents the features in the random field subset i; k = 1, 2, 4, ... L represents the classification label; w represents the dimension of the observed feature y; μ k and ∑ k is the expectation and variance of the k-th feature, π k represents the proportion of pixels marked as segmentation k, and T represents vector transpose.

[0112] Preferably, the expectation maximization algorithm can be used to optimize the parameters of the original Markov random field model, which is achieved by the following preferred implementation:

[0113] The energy function corresponding to the likelihood function of a pixel in the image is expressed as follows:

[0114]

[0115] The parameters θ of the likelihood function are estimated according to the maximum pseudo-likelihood method based on expectation maximization, iterating in the following way:

[0116]

[0117]

[0118] in, Indicates the proportion of pixels with k-type features in this iteration; θ k represents the parameters of the likelihood function, θ k ={μ k ,Σ k}; P(k|y i ,θ k ) represents the conditional probability of each pixel; Represents the expectation of k-class features in this iteration; Represents the expectation of k-class features in this iteration; Represents the variance of the feature labeled k in this iteration.

[0119] Preferably, the Gibbs distribution can be used as the prior probability distribution of the label field, which is expressed as follows:

[0120] P(x)=Z -1 e -U(x) ; Z = ∑ x∈Ω e -U(x) ; U(x)=∑ c∈C V c (x); where Z is the normalized value of the partitioning function, and the energy function U(x) is a series of potential functions V defined on the group c. c (x); the complement C represents the total set of groups; P(x) is the probability of a certain configuration x appearing in the random field X, and Ω is the set of all possible state vectors.

[0121] Among them, the energy function U(x) can be calculated by the following formula:

[0122]

[0123] Among them, V1, V2, V3 represent potential functions of different orders, and the order represents the complexity of the Gibbs distribution model; C1, C2, C3 represent the group sets corresponding to different orders; x i ,x j ,x k Represents the random variable in the pixel or group at position i, j, k.

[0124] Preferably, when only the second-order neighborhood system is considered, the energy function U(x) is expressed as:

[0125] Preferably, the Bode model can be used as the model of the marking field to finally obtain the Markov random field model, wherein the potential function is defined as follows:

[0126]

[0127] The local probability of the Porter model is:

[0128] where n i (x i ) is the number of cells in the neighborhood of position i that is equal to x i The number of pixel neighborhood positions, β ′ is the coupling coefficient.

[0129] Step S23: processing each of the Markov random field models according to a preset optimization algorithm to obtain a one-scene preliminary segmented image set and a two-scene preliminary segmented image set.

[0130] The preferred algorithm can be implemented as follows:

[0131] Adaptive weights are introduced, where the weight estimation model can be as follows:

[0132]

[0133] Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and y i Represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, N i Represents the set of neighboring pixels of pixel i.

[0134] Introducing penalty function The definition is as follows:

[0135]

[0136] Marking field energy function U(x i )for: Where β is a custom weighting coefficient ranging from 0 to 1, V2(x i ,x l ) represents the weight of the corresponding edge between pixels, and exp(*) represents the exponential operation;

[0137] The final energy function expression of the posterior probability solved pixel by pixel is:

[0138]

[0139] According to the preset number of classifications and the K-means algorithm, the images in the one-view normalized image set and the two-view normalized image set are respectively initially segmented to obtain an initial label for each pixel;

[0140] According to a preset number of iterations, each of the Markov random field models is iteratively optimized to obtain a set of preliminary segmented images of one scene and a set of preliminary segmented images of two scenes, and in each iteration, according to a preset expectation algorithm, the model parameters of the Markov random field model are estimated and calculated, and according to a preset iterative algorithm, the maximum likelihood parameter estimate of each category model is calculated until the preset convergence criterion is met, and the characteristic field energy weight under the current number of iterations and the model parameters of each category calculated by the expectation maximization algorithm are calculated, substituted into the final energy function expression, and the segmentation result of each pixel is updated; wherein, the convergence criterion can be that the absolute value of the parameter difference before and after the update is sufficiently small.

[0141] Step S24: processing the one-scene preliminary segmented image set and the two-scene preliminary segmented image set respectively according to a preset automatic region growing algorithm to obtain a one-scene secondary segmented image set and a two-scene secondary segmented image set.

[0142] The automatic region growing algorithm can be implemented in the following preferred manner:

[0143] According to a preset rule, a number of seed points are selected from the images in each of the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set;

[0144] According to the preset initial window size, the window size is gradually increased, and the texture difference of each pixel under each window size is calculated in turn, and the window size with the smallest texture difference is taken as the optimal window size; wherein, the initial window size can be 3×3, and the texture difference can be obtained by calculating the joint entropy, which is calculated as follows:

[0145]

[0146] in, represents the joint co-occurrence matrix of adjacent pixels in the same window, Represents the window texture corresponding to adjacent pixels p, q; represents the gray-level co-occurrence value; N c is the number of texture types, E n represents the joint entropy.

[0147] If the textures of the two windows are equally rough or smooth, the joint matrix will have only one single pixel, that is, the value of the joint entropy is zero.

[0148] According to the preset initial window, the probability distribution of each pixel belonging to each texture class is calculated, and the entropy of the initial window is calculated based on the probability distribution, wherein the calculation method is as follows:

[0149] Among them, α represents the set of texture classes, n represents the size of the window group centered on each pixel, and ε represents the entropy value.

[0150] Furthermore, we define the distance between two entropies as

[0151] in, Represents the standard deviation of each initial window class; the expression of the distance between this entropy measures the average similarity of two texture objects and quantifies the effectiveness of the segmentation. Once the optimal window is established, the average value of each initial window class is calculated. and standard deviation Therefore, if the entropy value of a pixel (εp ) is less different from its average value, then the pixel is aggregated into one region.

[0152] According to the optimal threshold and the optimal window size, the texture difference between each pixel point in the optimal window where various sub-points are located and the corresponding seed point is calculated, and the pixel points whose texture difference meets the optimal threshold are screened and merged into the initial area of ​​the corresponding seed point until all pixels are screened, thereby obtaining a one-scene secondary segmentation image set and a two-scene secondary segmentation image set respectively; wherein, the texture difference is obtained by calculating the difference in the grayscale co-occurrence value of each pixel point in the window; and the optimal threshold is calculated by the standard deviation of each initial window class, which can be obtained according to the following preferred embodiment:

[0153] Calculate the optimal threshold as follows: Among them, ∈ α represents the optimal threshold of texture type a, ε′ p represents the joint entropy value of pixel p, μ′ n represents the entropy mean of the initial window of this class, p→Classτ means aggregating pixels p into τ classes, Indicates the standard deviation of the texture of type A for each initial window.

[0154] Step S25: processing the one-scene secondary segmentation image set and the two-scene secondary segmentation image set respectively according to a preset morphological algorithm to obtain a one-scene cut image set and a two-scene cut image set.

[0155] Morphological operations are defined using sets. Generally, two types of pixel sets are used in image processing: target elements and structuring elements. A target element is defined as a set of foreground pixels, while a structuring element can be determined by combining foreground and background pixels. Basic morphological operations include erosion, dilation, opening, and closing. This solution utilizes opening and closing operations to process the results, removing speckle noise and filling holes.

[0156] For example, the morphological algorithm can be completed by opening and closing operations:

[0157] Among them, the opening operation of the structure element B on the set A (expressed as ) is defined as: It means that B first erodes A, and then B expands the erosion result;

[0158] Geometrically, the opening operation of B on A is the union of all translations of B. When the structural element B cannot fit the narrow line segments or figures in A, these small line segments or figures will be deleted. Therefore, the characteristic of the opening operation is to delete areas narrower than the structural element, which can usually smooth the contours of the object, disconnect narrow gaps, and eliminate slender protrusions.

[0159] The closing operation of the structure element B on the set A (expressed as A·B) is defined as: Indicates that B first expands A, and then B corrodes the expansion result. represents the dilation operation, Represents the erosion operation.

[0160] The closing operation is the complement of the union of all translations of B that do not overlap with A. This operation can bridge narrow breaks and elongated gullies, fill gaps in contours, smooth object edges, and fill small holes inside objects.

[0161] In step S2 of this embodiment, by constructing a Markov random field model, the spatial correlation between pixels in the image can be well captured, global optimization can be performed, noise interference can be effectively reduced, and a fine preliminary segmentation result can be provided; region growing can further refine the segmented area, especially in the edge or transition area, to ensure accurate division of shadow and non-shadow areas; morphological processing can eliminate isolated noise points in the segmentation results, smooth the boundaries, and improve segmentation accuracy and integrity.

[0162] Step S3: processing the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set.

[0163] In some embodiments of the present application, step S3 may be implemented by the following preferred implementations:

[0164] S31: Obtaining the original amplitude image corresponding to each image in the scene cut image set;

[0165] S32: compressing the original amplitude image according to a grayscale compression algorithm to obtain an amplitude image;

[0166] S33: According to the preset initial threshold, repeatedly divide the shadow area and the non-shadow area of ​​the amplitude image, and after each division, calculate the mean and standard deviation of the shadow area and the non-shadow area, and update the initial threshold according to the mean and standard deviation until the updated change of the initial threshold meets the preset difference range, thereby obtaining the optimal threshold; wherein, the preset difference range is, for example, 10 -6 .

[0167] In this preferred implementation step, an improved threshold segmentation method is used. Assuming that the amplitude image contains two grayscale areas, z represents the grayscale value in the image, and p(z) represents the probability density function corresponding to the histogram, then the mixed probability density function can be expressed as: p(z)=P1p1(z)+P2p2(z); wherein P1 represents the probability that the pixel belongs to the shadow area, P2 represents the probability that the pixel belongs to the non-shadow area, and P1+P2=1.

[0168] Furthermore, assuming that there is a threshold Γ that minimizes the average error in determining whether a pixel belongs to a shadow area or a non-shadow area, the probability of mistakenly determining a pixel belonging to a non-shadow area as a pixel belonging to a shadow area is: The probability of judging a pixel belonging to the shadow area as a pixel in the non-shadow area is: Then the overall error probability is: E(Γ)=P2E1(Γ)+P1E2(Γ).

[0169] Furthermore, in order to find the threshold that minimizes E(Γ), it is necessary to calculate the first-order derivative of E(Γ) with respect to Γ and make the first-order derivative equal to 0, and we can obtain: P1p1(Γ)=P2p2(Γ).

[0170] Substituting the Gaussian probability density function into the mixed probability density function, we can obtain:

[0171] Among them, μ1 and Represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the shadow area; μ2 and Represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the non-shaded area.

[0172] Furthermore, the mean and standard deviation of the shaded and non-shaded areas can be calculated as follows:

[0173]

[0174] Among them, μ S , represents the mean and standard deviation of the shaded area, ∑ (i，j) f S (i, j) represents the sum of the pixel values ​​in the shadow area, Num S Indicates the total number of pixels in the shadow area, μ B , represents the mean and standard deviation of the non-shaded area, ∑ (i，j) f B (i, j) represents the sum of the pixel values ​​in the non-shadow area, Num B represents the total number of pixels in the non-shadow area, ∑ (i，j) (f S (i, j)-μ S ) 2 represents the sum of squares of the difference between the pixels in the shadow area and the mean, ∑ (i，j) (f B (i, j)-μ B ) 2 Represents the sum of squared differences between the pixels in the non-shadow area and the mean.

[0175] According to P1p1(Γ)=P2p2(Γ), we can get: AΓ 2 +BΓ+C=0;wherein,

[0176] Furthermore, when the variances are the same, We can get: When P1=P2, the threshold Γ is the average of the means of the two distributions.

[0177] When the above improved algorithm is applied to a small area in an image, the threshold Γ corresponding to a pixel can be obtained by calculating the grayscale average of the pixels in a window centered on the pixel, as shown below:

[0178] Among them, z(i, j) represents the pixel value at (i, j), and m and n represent the size of the window.

[0179] S34: Processing the amplitude image according to the optimal threshold to obtain a first shadow area set and a first shadow undetermined area set.

[0180] In step S3 of this embodiment, by calculating the mean and variance of the shadow area and the background area, an iterative algorithm is used to optimize and obtain the optimal threshold, and the optimal threshold is used to separate the shadow area and the non-shadow area. When the difference in grayscale values ​​between the shadow and non-shadow areas is small, the false detection rate can be effectively reduced and the accuracy of shadow detection can be improved.

[0181] Step S4: processing the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set.

[0182] In some embodiments of the present application, step S4 can be implemented by the following preferred implementations:

[0183] S41: Calculate the coherence value of each of the two images after introducing the phase correction factor based on the one-scene cut image set and the two-scene cut image set. The calculation method is as follows:

[0184] Where z1, z2 represent the complex values ​​of the two images; * represents the conjugate operator; i represents the pixel number; represents the phase correction factor, Indicates the phase caused by subtracting the terrain factor from the original phase δ flattened represents the coherence value after the phase correction factor is introduced, LL is the number of pixels in the window, and j is the imaginary unit.

[0185] In step S41, if the phase correction factor is not introduced, the conventional coherence value Δ is directly calculated or the maximum likelihood coherence value δ of the pixels in the fixed window is calculated. L , as shown below:

[0186] The calculation of these two coherence values ​​does not fully take into account the image of terrain factors, and the deviation will be larger. Therefore, this application introduces a phase correction factor for calculation.

[0187] S42: Processing the coherence value after the phase correction factor is introduced according to a preset spatial averaging algorithm to obtain a final coherence value;

[0188] S43: Processing the one-scene cut image set and the two-scene cut image set according to the final coherence value to obtain a second shadow area set and a second shadow undetermined area set.

[0189] In step S4 of this embodiment, the coherence coefficient of the two images is calculated, and the phase correction factor and the spatial neighborhood average are used to improve the unbiasedness of the coherence estimation; setting the coherence threshold and extracting the shadow area can more accurately identify the low-coherence shadow area and reduce the deviation caused by terrain or noise.

[0190] Step S5: performing a set operation based on the first shadow area set, the first shadow area set to be determined, the second shadow area set and the second shadow area set to be determined, to obtain a first partial shadow area set and a third shadow area set to be determined.

[0191] In some embodiments of the present application, step S5 may be implemented by the following preferred implementations:

[0192] S51: Calculating the intersection of the first shadow area set and the second shadow area set to obtain a first partial shadow area set;

[0193] S52: After finding the intersection of the first set of shadow undetermined regions and the second set of shadow undetermined regions, subtract the first set of shadow undetermined regions to obtain a third set of shadow undetermined regions;

[0194] S53: Obtain a third shadow undetermined area set by calculating the intersection based on the third partial shadow undetermined area set, the first shadow undetermined area set, and the second shadow undetermined area set.

[0195] In step S5 of this embodiment, the advantages of different methods are integrated through set operations to improve the accuracy and robustness of shadow detection.

[0196] Step S6: processing the third set of shadow undetermined areas according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second set of shadow areas.

[0197] The pixel value in the SAR image reflects the echo energy value received by the receiver after the electromagnetic wave interacts with the target point. The incident angle has a great influence on the scattering coefficient. Based on the radar line of sight and combined with the digital elevation model data, the shadow area where the third shadow undetermined area is concentrated can be effectively detected.

[0198] In some embodiments of the present application, Figure 5 The schematic diagram of the geometric relationship between radar line of sight and digital elevation model analysis algorithm is shown. Step S6 can be implemented by the following preferred implementation:

[0199] S61: Along the slant range direction, for two points A and B, assuming that the slant range between the two pixel points is M, the elevation difference between the two points is H, and the radar depression angle is set to β1, when these three variables satisfy: H / M=sinβ1, then it can be determined that the point is a shadow area.

[0200] Furthermore, due to the error in the measurement data, the H / M value of the shadow is relatively large, so the threshold condition of H / M should be changed to: H / M ≥ ksinβ1; where 0 <k≤1。

[0201] S62: Processing the third set of shadow undetermined regions according to the above threshold condition to obtain a second set of partial shadow regions.

[0202] In step S6 of this embodiment, the geometric relationship of the objects in the image is analyzed by combining the radar line of sight and the digital elevation model data. The shadow area is determined by calculating the relationship between the ground height difference and the slant range. The shadow area caused by the insufficient radar depression angle is further detected, thereby improving the detection accuracy in areas with large terrain undulations.

[0203] Step S7: obtaining a final shadow area detection result by performing a union calculation based on the first partial shadow area set and the second partial shadow area set.

[0204] In step S7 of this embodiment, the SAR image shadow detection result obtained by fusion of multiple methods is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0205] Compared with the prior art, the embodiments of the application have the following beneficial effects: the original image is processed by the wavelet transform decomposition algorithm, which can effectively remove noise while retaining important features such as edges, and the bilateral filtering processing can smooth the image while maintaining the edges, thereby improving the quality of the image; the Markov random field model is constructed, the spatial dependence between pixels is fully utilized, the accuracy of segmentation is improved, the region growing algorithm is used, the target region is segmented more accurately by selecting a seed point and gradually expanding the region; the amplitude adaptive threshold detection method is used to adapt to the image characteristics in different scenes and improve the robustness of detection; the coherence calculation can effectively identify the shadow area caused by physical characteristics, especially for the shadow caused by terrain or object obstruction; the results obtained by various methods are subjected to set operation, the advantages of various methods are integrated, false detection and missed detection are reduced, and the overall detection effect is improved; through radar line of sight and digital elevation model analysis, the shadow area caused by insufficient radar depression angle is further detected, and the detection accuracy in the area with large terrain undulation is improved; finally, the SAR image shadow detection result of multi-method fusion is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0206] Embodiment two:

[0207] Please refer to Figure 2 , based on the same inventive concept, the embodiment of the application discloses a kind of multi-method fusion SAR image shadow detection system, comprising: preprocessing module M1, two-step segmentation module M2, amplitude detection module M3, coherence detection module M4, set operation module M5, radar elevation detection module M6 and result output module M7;

[0208] Wherein, the preprocessing module M1 is used to process original two-scene image set according to preset wavelet multiscale bilateral filtering algorithm, and obtains two-scene image set.

[0209] The preprocessing module M1 includes: wavelet decomposition unit, bilateral filtering unit, smoothing processing unit and inverse transformation unit;The wavelet decomposition unit processes each original two-scene image in the original two-scene image set according to the preset wavelet transform decomposition algorithm, respectively, to obtain each first high-frequency sub-band and first low-frequency sub-band;The bilateral filtering unit is used to process each first high-frequency sub-band according to the preset bilateral filtering algorithm, respectively, to obtain each second high-frequency sub-band;The smoothing processing unit is used to process each first low-frequency sub-band according to the preset smoothing processing algorithm, respectively, to obtain each second low-frequency sub-band;The inverse transformation unit is used to process each second high-frequency sub-band and second low-frequency sub-band according to the preset wavelet inverse transform reconstruction algorithm, to obtain each two-scene image, i.e. to obtain two-scene image set.

[0210] Preferably, the processing of the original two-view image in the original two-view image set according to a preset wavelet transform decomposition algorithm to obtain the first high-frequency sub-band and the first low-frequency sub-band can be achieved by the following preferred implementation:

[0211] The wavelet function of the wavelet transform decomposition algorithm Should meet the following requirements:

[0212] in, is the basis function or mother wavelet, which can form a cluster of function systems through scaling and translation on the time axis: in, It is called a sub-wavelet, a is the scale factor, which reflects the period length of the wavelet, b is the translation factor, which reflects the time translation, R represents a real number, and L2(R) represents a real two-dimensional space.

[0213] Preferably, the calculation formula of wavelet transform can be as follows:

[0214] Among them, w f (a, b) represents continuous wavelet transform, f(t) represents the function to be transformed by wavelet, a represents the scale factor, b represents the time translation factor, Represents a wavelet transform function that has undergone time translation and scale expansion, also known as a sub-wavelet.

[0215] Preferably, the processing of each of the first high frequency sub-bands according to a preset bilateral filtering algorithm to obtain each of the second high frequency sub-bands can be achieved by the following preferred implementations:

[0216] Among them, the calculation formula of bilateral filtering can be as follows:

[0217] Where p and q represent the coordinates of two pixels, i.e. p = (x0, y0), q = (x1, y1), BF[I] P represents the bilateral filtering centered at pixel p, I p , I q Represents the pixel value at position p, q, that is, the intensity value of the pixel, represents Gaussian weighting in the spatial domain, Represents Gaussian weighting within the pixel value range, W p is a normalization factor used to ensure that the weighted sum of the pixels is 1, W p The expression is as follows: Q represents the set of all pixels in a local neighborhood centered on pixel p.

[0218] Preferably, the processing of each of the first low-frequency sub-bands according to a preset smoothing algorithm to obtain each of the second low-frequency sub-bands can be achieved by the following preferred implementations:

[0219] The smoothing algorithm may be a mean filter: Among them, M represents a pixel, Num represents the total number of pixels, and f(i,j) represents the pixel value.

[0220] Preferably, when designing a linear smoothing filter, it should be ensured that the filter has a single peak, ie, a main lobe, and is symmetrical in both vertical and horizontal directions. For example, the window of a typical smoothing filter is 3×3.

[0221] Preferably, the smoothing algorithm may also use a Gaussian filter, wherein a first-order zero-mean Gaussian function may be expressed as: The two-dimensional discrete zero-mean Gaussian function can be expressed as: The distribution exponent σ represents the width of the Gaussian distribution.

[0222] The preprocessing module M1 of this embodiment processes the original image through the wavelet transform decomposition algorithm, which can effectively remove noise while retaining important features such as edges. The bilateral filtering process can smooth the image while maintaining the edges, thereby improving the image quality.

[0223] The two-step segmentation module M2 is used to process the two-scene image sets according to a preset Markov random field and region growing algorithm to obtain a one-scene cut image set and a two-scene cut image set.

[0224] The two-step segmentation module includes: a normalization unit, a random field construction unit, a first optimization unit, a region growing unit and a morphological processing unit; the normalization unit is used to process the one-scene image set and the two-scene image set respectively according to the normalization algorithm to obtain the one-scene normalized image set and the two-scene normalized image set; the random field construction unit is used to construct Markov random field models according to the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a model of a labeled field; the first optimization unit is used to process each of the Markov random field models according to a preset optimization algorithm to obtain the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set; the region growing unit is used to process the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set respectively according to a preset automatic region growing algorithm to obtain the one-scene secondary segmentation image set and the two-scene secondary segmentation image set; the morphological processing unit is used to process the one-scene secondary segmentation image set and the two-scene secondary segmentation image set respectively according to a preset morphological algorithm to obtain the one-scene cut image set and the two-scene cut image set.

[0225] Preferably, the normalization algorithm may be as follows:

[0226] Among them, I represents the pixel value in the original image; Min and Max represent the maximum and minimum values ​​of the pixels in the original image, and I N Indicates the normalized value.

[0227] Preferably, the construction of the Markov random field model can be achieved through the following preferred implementations:

[0228] Define a Markov random field, treat each pixel of each image in the one-view normalized image set and the two-view normalized image set as a random variable, model the grayscale distribution in the image based on the finite Gaussian mixture model, and construct an original Markov random field model. The probability distribution of each pixel in the image can be expressed as:

[0229]

[0230] Where S represents the set of random fields, i represents the subset of random fields; P(Y|X=x) represents the probability distribution function or likelihood function of the feature field when the label field is known; y i represents the features in the random field subset i; k = 1, 2, 4, ... L represents the classification label; w represents the dimension of the observed feature y; μ k and ∑ k is the expectation and variance of the k-th feature, π k represents the proportion of pixels marked as segmentation k, and T represents vector transpose.

[0231] Preferably, the expectation maximization algorithm can be used to optimize the parameters of the original Markov random field model, which is achieved by the following preferred implementation:

[0232] The energy function corresponding to the likelihood function of a pixel in the image is expressed as follows:

[0233]

[0234] The parameters θ of the likelihood function are estimated according to the maximum pseudo-likelihood method based on expectation maximization, iterating as follows:

[0235]

[0236] in, Indicates the proportion of pixels with k-type features in this iteration; θ k represents the parameters of the likelihood function, θ k ={μ k ,∑ k}; P(k|y i ,θ k) represents the conditional probability of each pixel; Represents the expectation of k-class features in this iteration; Represents the expectation of k-class features in this iteration; Represents the variance of the feature labeled k in this iteration.

[0237] Preferably, the Gibbs distribution can be used as the prior probability distribution of the label field, which is expressed as follows:

[0238] P(x)=Z -1 e -U (x); Z = ∑ x∈Ω e -U (x); U(x)=∑ c∈C V c (x); where Z is the normalized value of the partitioning function, and the energy function U(x) is a series of potential functions V defined on the group c. c (x); the complement C represents the total set of groups; P(x) is the probability of a certain configuration x appearing in the random field X, and Ω is the set of all possible state vectors.

[0239] Among them, the energy function U(x) can be calculated by the following formula:

[0240]

[0241] Among them, V1, V2, V3 represent potential functions of different orders, and the order represents the complexity of the Gibbs distribution model; C1, C2, C3 represent the group sets corresponding to different orders; x i , x j , x k Represents the random variable in the pixel or group at position i, j, k.

[0242] Preferably, when only the second-order neighborhood system is considered, the energy function U(x) is expressed as:

[0243] Preferably, the Bode model can be used as the model of the marking field to finally obtain the Markov random field model, wherein the potential function is defined as follows:

[0244]

[0245] The local probability of the Porter model is:

[0246] where n i (x i ) is the number of cells in the neighborhood of position i that is equal to x i is the number of pixel neighborhood positions, and β′ is the coupling coefficient.

[0247] Preferably, the preferred algorithm can be implemented by the following preferred implementation:

[0248] Adaptive weights are introduced, where the weight estimation model can be as follows:

[0249]

[0250] Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and y i Represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, and Ni represents the set of neighborhood pixels of pixel i.

[0251] The penalty function is introduced and defined as follows:

[0252]

[0253] Marking field energy function U(x i )for:

[0254] The final energy function expression of the posterior probability solved pixel by pixel is:

[0255]

[0256] According to the preset number of classifications and the K-means algorithm, the images in the one-view normalized image set and the two-view normalized image set are respectively initially segmented to obtain an initial label for each pixel;

[0257] According to a preset number of iterations, each of the Markov random field models is iteratively optimized to obtain a set of preliminary segmented images of one scene and a set of preliminary segmented images of two scenes, and in each iteration, according to a preset expectation algorithm, the model parameters of the Markov random field model are estimated and calculated, and according to a preset iterative algorithm, the maximum likelihood parameter estimate of each category model is calculated until the preset convergence criterion is met, and the characteristic field energy weight under the current number of iterations and the model parameters of each category calculated by the expectation maximization algorithm are calculated, substituted into the final energy function expression, and the segmentation result of each pixel is updated; wherein, the convergence criterion can be that the absolute value of the parameter difference before and after the update is sufficiently small.

[0258] Preferably, the automatic region growing algorithm can be implemented in the following preferred manner:

[0259] According to a preset rule, a plurality of seed points are selected from the images in each of the one-scene primary segmentation image set and the two-scene primary segmentation image set;

[0260] According to a preset initial window size, the window size is gradually increased, and the texture difference of each pixel point under each window size is calculated in sequence, and the window size with the minimum texture difference is taken as the best window size; wherein the texture difference can be obtained by calculating the joint entropy, and the calculation method is as follows:

[0261]

[0262] wherein, indicates the joint co-occurrence matrix of adjacent pixels in the same window, indicates the window texture corresponding to the adjacent pixels p, q; indicates the gray co-occurrence value; N c is the number of texture types, E n indicates the joint entropy.

[0263] According to a preset initial window, the probability distribution of each pixel point belonging to each texture class is calculated, and the entropy of the initial window is calculated according to the probability distribution, wherein the calculation method is as follows:

[0264] wherein, α represents the set of texture classes, n represents the window group size with each pixel as the center point, and ε represents the entropy value.

[0265] Further, the distance between two entropies is defined as

[0266] wherein, indicates the standard deviation of each initial window class; the expression of the distance between the two entropies measures the average similarity of the two texture objects, and quantifies the effectiveness of the segmentation. Once the best window is established, the average value and the standard deviation of each initial window class are calculated. p Therefore, if the entropy value (ε p ) of a pixel is relatively small compared to its average value, the pixel is aggregated into a region.

[0267] According to the best threshold and the best window size, the texture difference between each pixel point in the best window of each seed point and the corresponding seed point is calculated, the pixel points with the texture difference meeting the best threshold are screened, and are merged into the initial region of the corresponding seed point, until all pixel points are screened, and one-scene secondary segmentation image set and two-scene secondary segmentation image set are obtained respectively; wherein the texture difference is obtained by calculating the difference of the gray co-occurrence value of each pixel point in the window, and the best threshold is calculated by the standard deviation of each initial window class.

[0268] Furthermore, the optimal threshold can be calculated as follows:

[0269] Calculate the optimal threshold as follows: Among them, ∈ α represents the optimal threshold of texture type a, ε′ p represents the joint entropy value of pixel p, μ′ n represents the entropy mean of the initial window of this class, p→Classτ means aggregating pixels p into τ classes, Indicates the standard deviation of the texture of type A for each initial window.

[0270] Preferably, the morphological algorithm can be completed by opening and closing operations:

[0271] Among them, the opening operation of the structure element B on the set A (expressed as ) is defined as: It means that B first erodes A, and then B expands the erosion result;

[0272] The closing operation of the structure element B on the set A (expressed as A·B) is defined as: Indicates that B first expands A, and then B corrodes the expansion result. represents the dilation operation, Represents the erosion operation.

[0273] The two-step segmentation module M2 of this embodiment can well capture the spatial correlation between pixels in the image by constructing a Markov random field model, perform global optimization, effectively reduce noise interference, and provide fine preliminary segmentation results; region growing can further refine the segmentation area, especially in the edge or transition area, to ensure accurate division of shadow and non-shadow areas; morphological processing can eliminate isolated noise points in the segmentation results, smooth the boundaries, and improve segmentation accuracy and completeness.

[0274] The amplitude detection module M3 is used to process the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set.

[0275] The amplitude detection module M3 includes: an original acquisition unit, a compression unit, an optimal threshold calculation unit, and an amplitude segmentation unit. The original acquisition unit is used to obtain the original amplitude image corresponding to each image in the scene cut image set; the compression unit is used to compress the original amplitude image according to a grayscale compression algorithm to obtain an amplitude image; the optimal threshold calculation unit is used to repeatedly divide the amplitude image into shadow areas and non-shadow areas according to a preset initial threshold, and after each division, calculate the mean and standard deviation of the shadow area and the non-shadow area, and update the initial threshold according to the mean and standard deviation until the updated change of the initial threshold meets the preset difference range to obtain the optimal threshold; the amplitude segmentation unit is used to process the amplitude image according to the optimal threshold to obtain a first shadow area set and a first shadow undetermined area set.

[0276] Preferably, the optimal threshold calculation unit can be implemented by the following preferred implementations:

[0277] Assuming that the amplitude image contains two grayscale areas, z represents the grayscale value in the image, and p(z) represents the probability density function corresponding to the histogram, the mixed probability density function can be expressed as: p(z) = P1p1(z) + P2p2(z); where P1 represents the probability that the pixel belongs to the shadow area, P2 represents the probability that the pixel belongs to the non-shadow area, and P1+P2=1;

[0278] Furthermore, assuming that there is a threshold Γ that minimizes the average error in determining whether a pixel belongs to a shadow area or a non-shadow area, the probability of mistakenly determining a pixel belonging to a non-shadow area as a pixel belonging to a shadow area is: The probability of judging a pixel belonging to the shadow area as a pixel in the non-shadow area is: Then the overall error probability is: E(Γ)=P2E1(Γ)+P1E2(Γ).

[0279] Furthermore, in order to find the threshold that minimizes E(Γ), it is necessary to calculate the first-order derivative of E(Γ) with respect to Γ and make the first-order derivative equal to 0, and we can obtain: P1p1(Γ)=P2p2(Γ).

[0280] Substituting the Gaussian probability density function into the mixed probability density function, we can obtain:

[0281] Among them, μ1 and Represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the shadow area; μ2 and Represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the non-shaded area.

[0282] Furthermore, the mean and standard deviation of the shaded and non-shaded areas can be calculated as follows:

[0283]

[0284] Among them, μ S , represents the mean and standard deviation of the shaded area, ∑ (i，j) f S (i, j) represents the sum of the pixel values ​​in the shadow area, Num S Indicates the total number of pixels in the shadow area, μ B , represents the mean and standard deviation of the non-shaded area, ∑ (i，j) f B (i, j) represents the sum of the pixel values ​​in the non-shadow area, Num B represents the total number of pixels in the non-shadow area, ∑ (i，j) (f S (i, j)-μ S ) 2 represents the sum of squares of the difference between the pixels in the shadow area and the mean, ∑ (i，j) (f B (i, j)-μ B ) 2 Represents the sum of squares of the differences between the pixels in the non-shadow area and the mean;

[0285] According to P1p1(Γ)=P2p2(Γ), we can get: AΓ 2 +BΓ+C=0;wherein,

[0286] Furthermore, when the variances are the same, We can get:

[0287] When P1=P2, the threshold Γ is the average of the means of the two distributions.

[0288] When the above improved algorithm is applied to a small area in an image, the threshold Γ corresponding to a pixel can be obtained by calculating the grayscale average of the pixels in a window centered on the pixel, as shown below:

[0289] Among them, z(i, j) represents the pixel value at (i, j), and m and n represent the size of the window.

[0290] The amplitude detection module M3 of the embodiment obtains the optimal threshold value by calculating the mean and variance of the shadow area and the background area, using an iterative algorithm to optimize, and separating the shadow area and the non-shadow area using the optimal threshold value, and can effectively reduce the false detection rate and improve the accuracy of shadow detection in the case of small difference between the shadow and non-shadow area gray values.

[0291] The coherence detection module M4 is configured to process the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow pending area set.

[0292] The coherence detection module M4 comprises a coherence value calculation unit, a spatial average algorithm processing unit and a coherence cutting unit. The coherence value calculation unit is configured to calculate the coherence value of each two-scene image after introducing a phase correction factor according to the one-scene cut image set and the two-scene cut image set. The spatial average algorithm processing unit is configured to process the coherence value after introducing the phase correction factor according to a preset spatial average algorithm to obtain a final coherence value. The coherence cutting unit is configured to process the one-scene cut image set and the two-scene cut image set according to the final coherence value to obtain the second shadow area set and the second shadow pending area set.

[0293] Preferably, the coherence value of each two-scene image after introducing the phase correction factor is calculated in the following manner:

[0294] Wherein, z1 and z2 represent the complex values of the two-scene images respectively; * represents the conjugate operator; i represents the pixel serial number. The phase correction factor is represented as The phase caused by the terrain factor in the original phase is represented as δ flattened The coherence value after introducing the phase correction factor is represented as

[0295] The coherence detection module M4 of the embodiment calculates the coherence coefficient of the two-scene images, uses the phase correction factor and the spatial neighborhood average value to improve the unbiasedness of the coherence estimation, sets the coherence threshold value, extracts the shadow area, and can more accurately identify the shadow area with low coherence and reduce the deviation caused by the terrain or noise.

[0296] The set operation module M5 is configured to perform set operation on the first shadow area set, the first shadow pending area set, the second shadow area set and the second shadow pending area set to obtain a first partial shadow area set and a third shadow pending area set.

[0297] The set operation module M5 includes: a first operation unit, a second operation unit and a third operation unit; the first operation unit is used to calculate the intersection of the first shadow area set and the second shadow area set to obtain a first partial shadow area set; the second operation unit is used to calculate the intersection of the first shadow undetermined area set and the second shadow undetermined area set, and then subtract the first partial shadow area set to obtain a third partial shadow undetermined area set; the third operation unit is used to calculate the intersection of the third partial shadow undetermined area set, the first shadow undetermined area set and the second shadow undetermined area set to obtain a third shadow undetermined area set.

[0298] The set operation module M5 of this embodiment integrates the advantages of different methods through set operations to improve the accuracy and robustness of shadow detection.

[0299] The radar elevation detection module M6 is used to process the third set of shadow undetermined areas according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second partial shadow area set.

[0300] Preferably, the processing of the third set of shadow undetermined areas according to a preset radar line of sight and digital elevation model analysis algorithm to obtain the second set of partial shadow areas can be achieved by the following preferred implementations:

[0301] Consider two points A and B along the slant range direction. Assume that the slant range between the two pixels is M, the elevation difference between the two points is H, and the radar depression angle is set to β1. When these three variables satisfy: H / M = sinβ1, then the point can be determined to be in the shadow area.

[0302] Furthermore, due to the error in the measurement data, the H / M value of the shadow is relatively large, so the threshold condition of H / M should be changed to: H / M ≥ ksinβ1; where 0 <k≤1;

[0303] According to the above threshold condition, the third set of shadow undetermined regions is processed to obtain a second set of partial shadow regions.

[0304] The radar elevation detection module M6 of this embodiment combines radar line of sight and digital elevation model data to analyze the geometric relationship of the objects in the image, determines the shadow area by calculating the relationship between the ground height difference and the slant range, and further detects the shadow area caused by insufficient radar depression angle, thereby improving the detection accuracy in areas with large terrain undulations.

[0305] The result output module M7 is configured to calculate a final shadow area detection result based on the union of the first partial shadow area set and the second partial shadow area set.

[0306] The result output module M7 of this embodiment significantly improves the accuracy and comprehensiveness of shadow detection by outputting the SAR image shadow detection results fused by multiple methods.

[0307] In summary, compared with the prior art, the embodiments of the present application have the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, noise can be effectively removed while retaining important features such as edges. The bilateral filtering process can smooth the image while maintaining the edges, thereby improving the image quality. A Markov random field model is constructed to fully utilize the spatial dependency between pixels to improve the accuracy of segmentation. The region growing algorithm is used to more accurately segment the target area by selecting seed points and gradually expanding the region. The amplitude adaptive threshold detection method is adopted to adapt to the image characteristics in different scenarios and improve the robustness of detection. The coherence calculation can effectively identify shadow areas caused by physical characteristics, especially those caused by terrain or object occlusion. The results obtained by multiple methods are subjected to set operations to combine the advantages of various methods, reduce false detections and missed detections, and improve the overall detection effect. The shadow area caused by insufficient radar depression angle is further detected through radar line of sight and digital elevation model analysis, thereby improving the detection accuracy in areas with large terrain undulations. Finally, the SAR image shadow detection result fused by multiple methods is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0308] The specific working process of each module described above can refer to the corresponding process in the aforementioned method embodiment and will not be repeated here. The division of the modules is only a logical function division. In actual implementation, there may be other division methods, such as combining multiple modules or integrating them into another system.

[0309] The specific embodiments described above further illustrate the objectives, technical solutions, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A multi-method fusion SAR image shadow detection method, characterized in that: include: According to the preset wavelet multi-scale bilateral filtering algorithm, the original two-view image set is processed to obtain the two-view image set; Processing the two image sets according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set; Processing the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set; Processing the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set; According to the first shadow area set, the first shadow area set to be determined, the second shadow area set and the second shadow area set to be determined, a set operation is performed to obtain a first partial shadow area set and a third shadow area set to be determined; Processing the third set of undetermined shadow areas according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second set of shadow areas; A final shadow area detection result is obtained by performing a union calculation based on the first shadow area set and the second shadow area set; The method of processing the original two-view image set according to the preset wavelet multi-scale bilateral filtering algorithm to obtain the two-view image set includes: According to a preset wavelet transform decomposition algorithm, the original two images in the original two-view image set are processed respectively to obtain each first high-frequency sub-band and first low-frequency sub-band; Processing each of the first high frequency sub-bands according to a preset bilateral filtering algorithm to obtain each of the second high frequency sub-bands; Processing each of the first low-frequency sub-bands according to a preset smoothing algorithm to obtain each of the second low-frequency sub-bands; Processing each of the second high-frequency sub-bands and the second low-frequency sub-bands according to a preset inverse wavelet transform reconstruction algorithm to obtain two images, that is, obtaining a two-image set; The method of processing the two-view image sets according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set includes: According to the normalization algorithm, the one-view image set and the two-view image set are processed respectively to obtain the one-view normalized image set and the two-view normalized image set; Constructing Markov random field models based on the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a model of a label field; According to a preset optimization algorithm, each of the Markov random field models is processed respectively to obtain a one-scene preliminary segmentation image set and a two-scene preliminary segmentation image set; According to a preset automatic region growing algorithm, the first scene preliminary segmentation image set and the second scene preliminary segmentation image set are processed respectively to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set; According to a preset morphological algorithm, the one-scene secondary segmentation image set and the two-scene secondary segmentation image set are processed respectively to obtain a one-scene cut image set and a two-scene cut image set.

2. The multi-method fusion SAR image shadow detection method according to claim 1, characterized in that: The Markov random field model is constructed based on the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a model of a label field, including: A Markov random field is defined, each pixel of each image in the one-view normalized image set and the two-view normalized image set is regarded as a random variable, the grayscale distribution in the image is modeled according to a finite Gaussian mixture model, and an original Markov random field model is constructed; Optimizing the parameters of the original Markov random field model using an expectation-maximization algorithm; The Markov random field model is obtained by adopting Gibbs distribution as the prior probability distribution of the marking field and adopting Porter model as the model of the marking field.

3. The multi-method fusion SAR image shadow detection method according to claim 2, characterized in that: The Markov random field models are processed according to a preset optimization algorithm to obtain a first-scene preliminary segmented image set and a second-scene preliminary segmented image set, including: Adaptive weights are introduced, where the weight estimation model is as follows: Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and y i Represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, N i represents the set of neighboring pixels of pixel i; Introducing penalty function The definition is as follows: Marking field energy function U(x i )for: Where β is a custom weighting coefficient ranging from 0 to 1, V2(x i ,x l ) represents the weight of the corresponding edge between pixels, and exp(*) represents the exponential operation; The final energy function expression of the posterior probability solved pixel by pixel is: According to the preset number of classifications and the K-means algorithm, the images in the one-view normalized image set and the two-view normalized image set are respectively initially segmented to obtain an initial label for each pixel; According to the preset number of iterations, each of the Markov random field models is iteratively optimized to obtain a set of preliminary segmented images of one scene and a set of preliminary segmented images of two scenes. In each iteration, according to the preset expectation algorithm, the model parameters of the Markov random field model are estimated and calculated. According to the preset iterative algorithm, the maximum likelihood parameter estimate of each category model is calculated until the preset convergence standard is met, and the characteristic field energy weight under the current number of iterations and the model parameters of each category calculated by the expectation maximization algorithm are calculated, substituted into the final energy function expression, and the segmentation result of each pixel is updated.

4. The multi-method fusion SAR image shadow detection method according to claim 3, characterized in that: The method of processing the first scene preliminary segmentation image set and the second scene preliminary segmentation image set respectively according to a preset automatic region growing algorithm to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set includes: According to a preset rule, a number of seed points are selected from the images in each of the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set; According to the preset initial window size, the window size is gradually increased, and the texture difference of each pixel under each window size is calculated in turn, and the window size with the smallest texture difference is taken as the optimal window size; Calculating the probability distribution of each pixel belonging to each texture class according to a preset initial window, and calculating the entropy of the initial window according to the probability distribution; According to the optimal threshold and the optimal window size, the texture difference between each pixel point in the optimal window where various sub-points are located and the corresponding seed point is calculated, and the pixel points whose texture difference meets the optimal threshold are screened and merged into the initial area of ​​the corresponding seed point until all pixel points are screened, thereby obtaining a one-scene secondary segmentation image set and a two-scene secondary segmentation image set respectively; wherein, the texture difference is obtained by calculating the difference in the grayscale co-occurrence value of each pixel point in the window, and the optimal threshold is calculated by the standard deviation of each initial window class.

5. The multi-method fusion SAR image shadow detection method according to claim 4, characterized in that: The step of processing the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set includes: Obtaining the original amplitude image corresponding to each image in the scene cut image set; Compress the original amplitude image according to the gray-scale compression algorithm to obtain an amplitude image; According to a preset initial threshold, repeatedly divide the shadow area and non-shadow area of the amplitude image. After each division, calculate the mean and standard deviation of the shadow area and non-shadow area, and update the initial threshold according to the mean and standard deviation until the update change of the initial threshold meets the preset difference range to obtain an optimal threshold; Process the amplitude image according to the optimal threshold to obtain a first set of shadow areas and a first set of shadow pending areas.

6. The multi-method fusion SAR image shadow detection method according to claim 5, characterized in that: The processing of the first set of segmented images and the second set of segmented images according to the preset coherence coefficient detection method to obtain a second set of shadow areas and a second set of shadow pending areas includes: According to the first set of segmented images and the second set of segmented images, calculate the coherence values of each pair of images after introducing a phase correction factor, and the calculation method is as follows: Where z1, z2 represent the complex values ​​of the two images; * represents the conjugate operator; i represents the pixel number; represents the phase correction factor, Indicates the phase caused by subtracting the terrain factor from the original phase δ flattened represents the coherence value after the phase correction factor is introduced, LL is the number of pixels in the window, and j is the imaginary unit; Process the coherence value after introducing the phase correction factor according to the preset spatial averaging algorithm to obtain a final coherence value; Process the first set of segmented images and the second set of segmented images according to the final coherence value to obtain a second set of shadow areas and a second set of shadow pending areas.

7. The multi-method fusion SAR image shadow detection method according to claim 6, characterized in that: The processing of the third set of shadow pending areas according to the preset radar line of sight and digital elevation model analysis algorithm to obtain a second set of partial shadow areas includes: The preset radar line of sight and digital elevation model analysis algorithm includes a geometric relationship formula: H / M≥ksinβ1; 0<k≤1; where M, H, and β1 are the slant range, elevation difference, and radar depression angle corresponding to two pixel points respectively; Process the images in the third set of shadow pending areas according to the geometric relationship formula, and classify the pixel points that meet the geometric relationship formula as shadow areas to obtain a second set of partial shadow areas.

8. A multi-method fusion SAR image shadow detection system, characterized by: Includes: A preprocessing module, a two-step segmentation module, an amplitude detection module, a coherence detection module, a set operation module, a radar elevation detection module, and a result output module; Among them, the preprocessing module is used to process the original two sets of images according to the preset wavelet multi-scale bilateral filtering algorithm to obtain two sets of images; The two-step segmentation module is used to process the two sets of images according to the preset Markov random field and region growing algorithm to obtain a first set of segmented images and a second set of segmented images; The amplitude detection module is used to process the first set of segmented images according to the preset amplitude adaptive threshold detection method to obtain a first set of shadow areas and a first set of shadow pending areas; The coherence detection module is used to process the first set of segmented images and the second set of segmented images according to the preset coherence coefficient detection method to obtain a second set of shadow areas and a second set of shadow pending areas; The set operation module is used to perform set operations on the first set of shadow areas, the first set of shadow pending areas, the second set of shadow areas, and the second set of shadow pending areas to obtain a first set of partial shadow areas and a third set of shadow pending areas; The radar elevation detection module is used to process the third set of shadow pending areas according to the preset radar line of sight and digital elevation model analysis algorithm to obtain a second set of partial shadow areas; The result output module is used to calculate the final shadow area detection result based on the union of the first partial shadow area set and the second partial shadow area set; The preprocessing module includes: a wavelet decomposition unit, a bilateral filtering unit, a smoothing processing unit and an inverse transform unit; the wavelet decomposition unit is used to process each original two-scene image in the original two-scene image set according to a preset wavelet transform decomposition algorithm to obtain each first high-frequency sub-band and each first low-frequency sub-band; the bilateral filtering unit is used to process each first high-frequency sub-band according to a preset bilateral filtering algorithm to obtain each second high-frequency sub-band; the smoothing processing unit is used to process each first low-frequency sub-band according to a preset smoothing processing algorithm to obtain each second low-frequency sub-band; the inverse transform unit is used to process each second high-frequency sub-band and second low-frequency sub-band according to a preset wavelet inverse transform reconstruction algorithm to obtain each two-scene image, that is, to obtain a two-scene image set; The two-step segmentation module includes: a normalization unit, a random field construction unit, a first optimization unit, a region growing unit and a morphological processing unit; the normalization unit is used to process the one-scene image set and the two-scene image set respectively according to the normalization algorithm to obtain the one-scene normalized image set and the two-scene normalized image set; the random field construction unit is used to construct Markov random field models according to the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a model of a labeled field; the first optimization unit is used to process each of the Markov random field models according to a preset optimization algorithm to obtain the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set; the region growing unit is used to process the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set respectively according to a preset automatic region growing algorithm to obtain the one-scene secondary segmentation image set and the two-scene secondary segmentation image set; the morphological processing unit is used to process the one-scene secondary segmentation image set and the two-scene secondary segmentation image set respectively according to a preset morphological algorithm to obtain the one-scene cut image set and the two-scene cut image set.

Images