SAR Image Lake Shoreline Detection Method and System Based on MRSF Model

Abstract

The present invention discloses a method and system for detecting lake shorelines in SAR images based on the MRSF model. The method of the present invention includes the following steps: Step 1), suppressing speckle noise in the registered SAR image using the BM3D algorithm; Step 2), coarsely segmenting the SAR image after noise reduction using the OSTU algorithm to obtain the initial boundary between water and land; Step 3), constructing an energy function of the MRSF model based on a new boundary energy term and an introduced global energy term; Step 4), evolving the result of the coarse segmentation as the initial contour of the MRSF model to obtain the final fine lake shoreline. The present invention is more effective and accurate in extracting lake shorelines from SAR images.

Description

Technical Field

[0001] The present invention belongs to the technical field of image detection and target recognition, and relates to the shoreline detection technology of synthetic aperture radar (SAR) images. Specifically, it is a shoreline detection method and system based on an improved RSF (MRSF, modified RSF) model. In the present invention, the SAR image is first filtered by BM3D (block-matching 3D) to suppress speckle noise. In the coarse segmentation stage, an initial contour is obtained by combining the classic OSTU algorithm with morphological processing. In the fine segmentation stage, a novel energy function is constructed based on the edge energy term and the global energy term, where the edge energy function is jointly detected by the LoG operator and the ROEWA operator. Background Art

[0002] Since SAR sensors can provide all-weather and all-time comprehensive surface observations, the automatic interpretation of SAR images has wide applications in environmental remote sensing. The segmentation algorithm is an effective processing method for the automatic interpretation of SAR images, such as the detection of shorelines, the segmentation of water areas, and the segmentation of vegetation textures. The detected shoreline can reflect land subsidence, water texture structure, erosion and deposition, etc., which puts forward higher requirements for the accuracy of shoreline detection. The accurate detection of shorelines has also become an important topic in current research.

[0003] In view of the important role of SAR image segmentation in shoreline detection, researchers have proposed segmentation algorithms such as threshold segmentation, watershed, clustering, and active contour models. Traditional algorithms have good effects in spectral images, but have poor scalability, large local limitations, and discontinuous boundaries in the detection of SAR image shorelines. The active contour model based on curve evolution can effectively solve the application limitations of traditional algorithms. The basic idea of the active contour model is to first establish a curve contour, and then construct an energy function. During the process of minimizing the energy function, the contour line continuously evolves along the normal direction and finally stops evolving at the edge of the target. The active contour model is generally divided into two major types: the edge-based model and the region-based model. The edge-based model uses the gradient as the driving force for the evolution of the edge curve and has good effects for images with obvious gradients, but the effects for weak boundary images are not satisfactory. The region-based model uses regional statistical information as the driving force for the evolution of the edge curve. This model is insensitive to noise and also has good performance on weak boundary images. In a complex water environment, due to the speckle noise and intensity non-uniformity of SAR images, directly applying the active contour model to shoreline detection has poor effects. Therefore, it is necessary to construct a more robust energy function. Summary of the Invention

[0004] To solve the above problems, the present invention proposes a method and system for detecting the shoreline of SAR images based on the MRSF model. In the present invention, after the SAR image is filtered by the BM3D algorithm, the OSTU algorithm is used in combination with morphological processing to obtain an initial contour, and then the MRSF model is used for fine segmentation based on this initial contour. Based on the new edge energy term and the introduced global energy term, a robust energy function is constructed in the MRSF model of the present invention. The traditional RSF model only uses a single regional energy as the energy for curve evolution, losing a large amount of image information and making it difficult to accurately detect the shoreline. In contrast, the present invention makes full use of the edge features and global features of the SAR image, and uses the result of the OSTU algorithm as the initial contour, effectively avoiding the wrong segmentation and falling into local minima caused by inappropriate initial contour settings in the traditional RSF model.

[0005] The novel energy function constructed by the present invention includes both a new edge energy term and a global energy term. Among them, the new edge energy term is jointly detected by the LoG operator and the ROEWA operator, which can better locate the weak edges in the SAR image and reduce false edges. Combining the concept of relative entropy, the global energy term fits the global information of the SAR image, makes full use of the global statistical features of the image, and effectively overcomes the shortcoming of the traditional RSF model being limited to local energy.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for detecting the shoreline of SAR images based on the MRSF model, which includes the following steps:

[0008] Step 1): Use the BM3D algorithm to suppress the speckle noise of the registered SAR image;

[0009] Step 2): Use the OSTU algorithm to perform rough segmentation on the SAR image after noise reduction to obtain the initial boundary between water and land;

[0010] Step 3): Construct the energy function of the MRSF model based on the new boundary energy term and the introduced global energy term;

[0011] Step 4): Use the result of the rough segmentation as the initial contour of the MRSF model to evolve and obtain the final fine shoreline.

[0012] Preferably, step 1) is specifically implemented as follows:

[0013] The classical BM3D filtering algorithm is proposed to suppress additive white Gaussian noise (AWGN). Considering that the speckle noise in the SAR image belongs to multiplicative noise, the intensity of the observed signal I(s) can be expressed as:

[0014] I(s) = Q(s)u(s) (1)

[0015] Wherein, Q(s) represents the signal strength without noise, and u(s) represents the speckle noise conforming to the Gamma distribution. Performing a logarithmic transformation on the above formula, the multiplicative noise model is transformed into the following additive noise model:

[0016] logI = logQ + logu (2)

[0017] Therefore, the speckle noise can be regarded as AWGN, and logQ can be estimated by using the BM3D algorithm through logI.

[0018] The BM3D algorithm mainly includes two steps. In the first step, the original noisy image is first filtered to obtain a basic estimate. Then, the actual filtering is performed in the second step based on the estimate generated in the first step. Three operations are performed in both steps:

[0019] 1) Grouping;

[0020] 2) Collaborative filtering;

[0021] 3) Aggregation.

[0022] In the grouping operation, considering the statistical characteristics of the noise, similar target blocks are grouped into a 3D stack based on the non-local principle. In the collaborative filtering operation, a wavelet transform is performed on the 3D stack. Finally, all the target blocks return to their original positions through weighted averaging in the aggregation operation.

[0023] Preferably, step 2) is specifically implemented as follows:

[0024] In the traditional RSF model, if the initial contour is set unreasonably, it will not only lead to a sharp increase in the amount of computation but also cause the evolution of the curve to fall into a local minimum. To solve this limitation, the present invention uses the OSTU algorithm combined with morphological processing to obtain an initial rough segmentation contour. Considering that variance can be used as a measure of the gray-scale distribution uniformity of the SAR image, the present invention divides the SAR image into two parts, namely, lakes and land, by maximizing the between-class variance, thereby obtaining the initial contour of the lakes.

[0025] Preferably, step 3) is specifically implemented as follows:

[0026] To make full use of the global information of the image and enhance the weak edge preservation ability, the present invention constructs a novel energy function on the basis of the RSF model:

[0027] E MRSF = E RSF + ωE edge + ξE global (3)

[0028] Among them, ω and ξ represent weight coefficients. The constructed energy function mainly includes three parts: the traditional RSF energy term, the edge energy term, and the global energy term. The representation forms of each part are described as follows.

[0029] (1) RSF model

[0030] To solve the problem that the CV model cannot handle images with intensity inhomogeneity, researchers proposed the RSF model based on the region variable fitting energy. In the RSF model, the local energy function is defined as:

[0031]

[0032] Among them, λ 1 , λ 2 , μ, ν represent positive constants, K σ (x - y) represents a Gaussian kernel with a standard deviation of σ, and

[0033]

[0034] Among them, H ε represents the regularized Heaviside function, and its definition is:

[0035]

[0036] The smoothing function f 1 (·) and f 2 (·) approximately represent the image intensities inside and outside the contour C in a certain region respectively.

[0037] In addition, this energy function simultaneously considers the distance regularization term P(φ) and the length term L(φ), and their definitions are respectively:

[0038]

[0039]

[0040] Among them, δ ε (φ) represents the regularized Dirac function, and its definition is:

[0041]

[0042] P(φ) can maintain the horizontality and regularity of the function, and L(φ) helps to smooth the curve C.

[0043] By minimizing the level set function φ with respect to E RSF through the steepest descent method, the following variational equation is obtained:

[0044]

[0045] (2) New edge energy term

[0046] Considering that there are a large number of weak edges and false edges in SAR images, the present invention constructs a new edge energy term based on the LoG operator and the ROEWA operator:

[0047]

[0048] And

[0049] g(x) = exp{-α|x|} (12)

[0050] Z = mM LoG + nM ROEWA (13)

[0051] Where α, β, m, and n all represent weighting values, and K σ represents a Gaussian kernel with a standard deviation of σ. g(x) is an edge indicator, whose value is approximately 1 in homogeneous regions and approximately 0 near the water-land boundary. Z represents the fusion result of LoG and ROEWA weighted by m and n respectively, and M LoG and M ROEWA represent the detection results of the LoG operator and the ROEWA operator respectively.

[0052] 1) LoG operator

[0053] The LoG operator is a well-known second-order difference edge detection operator, mainly composed of a Gaussian operator and a Laplace operator, and its definition is:

[0054] M LoG (x, y) = Δ[K σ * I(x, y)] = I(x, y) * LoG (14)

[0055] Where K σ represents a Gaussian kernel with a standard deviation of σ, Δ represents the Laplace operator, and M LoG (x, y) represents the image detected by the LoG operator. The LoG operator is positive in the dark part of the image edge and negative in the bright part of the edge. Slow-changing edges and rapidly changing edges can be detected through the zero-crossing points of positive and negative values. Therefore, the LoG operator can detect the water-land edges of SAR images mainly in the form of steps.

[0056] 2) ROEWA operator

[0057] The ROEWA operator is a multi-edge model based on linear minimum mean square error and is suitable for speckle noise in SAR images. In the two-dimensional case, the definition of the ROEWA operator is:

[0058] h2-D (x,y) = h(x)h(y) (15)

[0059] Where h(x) and h(y) both represent one - dimensional filters, and h(x) can be given by the following formula:

[0060] h(x) = Ce -k|x| (16)

[0061] Where k represents the filtering coefficient and C represents the regularization constant.

[0062] The edge intensity detected by the ROEWA operator is jointly given by the horizontal and vertical components:

[0063]

[0064] Where M X (x,y) is the horizontal edge intensity component, and M Y (x,y) is the vertical edge intensity component.

[0065] To calculate the horizontal edge intensity component of the image I(x,y), first smooth each column of I(x,y) using the filter h(y), and then use the causal filter h 1 (·) and the anti - causal filter h 2 (·) to convolve with each row of the smoothed image to obtain the following exponentially weighted averages μ X1 and μ X2 :

[0066]

[0067] Where * represents the convolution operation. When μ X1 and μ X2 are normalized to be unbiased, the horizontal edge intensity component M X of the image can be given by the following formula:

[0068]

[0069] In addition, the vertical edge intensity component M Y (x,y) can be calculated using a similar process.

[0070] (3) Global energy term

[0071] Considering the limitation that the traditional RSF model only uses local region fitting, and combining the concept of relative entropy, the present invention introduces a global energy function that can fully utilize the global statistical features of SAR images:

[0072]

[0073] Among them, I G (x) is used to fit the global image. For the two-phase case, I G (x) is expressed as:

[0074] I G (x) = c 1 H ε (φ) + c 2 (1 - H ε (φ))(21)

[0075] Among them, c 1 and c 2 represent the average intensities of the two regions inside and outside the curve C.

[0076] Preferably, step 4) is specifically implemented as follows:

[0077] The present invention evolves the curve C by the gradient descent flow method. According to the variational principle, the Euler-Lagrange equation for minimizing the energy functional E can be obtained:

[0078]

[0079] Here, the energy functional E is with respect to the level set function φ.

[0080] Based on the above analysis, the energy functional E constructed by the present invention MRSF can be expressed as:

[0081]

[0082] According to Equation (22), it can be obtained that:

[0083]

[0084] Combining the above variational equation with the initial contour obtained in step 2) and iteratively solving the minimum value of the energy functional, the final fine lake shoreline can be obtained.

[0085] The present invention also discloses a SAR image lake shoreline detection system based on the MRSF model, which includes the following modules:

[0086] Speckle noise suppression module: The registered SAR image is processed by the BM3D algorithm to suppress speckle noise;

[0087] Coarse segmentation module: The denoised SAR image is coarsely segmented by the OSTU algorithm to obtain the initial boundary between water and land;

[0088] Energy function module for constructing the MRSF model: Based on the new boundary energy term and the introduced global energy term, the energy function of the MRSF model is constructed;

[0089] Lake shoreline detection result output module: The result of rough segmentation is used as the initial contour of the MRSF model to evolve and obtain the final fine lake shoreline.

[0090] Preferably, in the speckle noise suppression module, the BM3D algorithm includes two steps: First step, the original noisy image is filtered to obtain a basic estimate; Second step, the actual filtering is performed based on the estimate generated in the first step in the second step; Both steps perform the following three operation steps:

[0091] 1) Grouping: In the grouping operation, considering the statistical characteristics of the noise, similar target blocks are grouped into a 3D stack based on the non-local principle;

[0092] 2) Collaborative filtering: In the collaborative filtering operation, a wavelet transform is performed on the 3D stack;

[0093] 3) Aggregation: All target blocks return to their original positions through weighted averaging in the aggregation operation.

[0094] Preferably, in the rough segmentation module, the SAR image is divided into two parts, lake and land, by maximizing the inter-class variance, thereby obtaining the initial contour of the lake.

[0095] Preferably, the energy function module for constructing the MRSF model is specifically as follows:

[0096] Based on the RSF model, a novel energy function is constructed:

[0097] E MRSF =E RSF +ωE edge +ξE global (3)

[0098] Among them, ω and ξ represent weight coefficients; The constructed energy function mainly includes three parts: traditional RSF energy term, edge energy term, and global energy term; The representation forms of each part are as follows:

[0099] (1) RSF model

[0100] In the RSF model, the local energy function is defined as:

[0101]

[0102] Among them, λ 1 ,λ 2 ,μ,ν represent positive constants, K σ (x - y) represents a Gaussian kernel with a standard deviation of σ, and

[0103]

[0104] Among them, Hε Denotes the regularized Heaviside function, defined as:

[0105]

[0106] Smoothing function f 1 (·) and f 2 (·) approximate the image intensities inside and outside the contour C in a certain region respectively;

[0107] In addition, this energy function simultaneously takes into account the distance regularization term P(φ) and the length term L(φ), defined respectively as:

[0108]

[0109]

[0110] Among them, δ ε (φ) denotes the regularized Dirac function, defined as:

[0111]

[0112] P(φ) maintains the horizontality and regularity of the function, and L(φ) helps to smooth the curve C;

[0113] Minimize the level set function φ with respect to E RSF by the steepest descent method, and the following variational equation is obtained:

[0114]

[0115] (2) New edge energy term

[0116] Construct a new edge energy term based on the LoG operator and the ROEWA operator:

[0117]

[0118] And

[0119] g(x) = exp{-α|x|} (12)

[0120] Z = mM LoG + nM ROEWA (13)

[0121] Among them, α, β, m, n all represent weighted values, and K σ denotes a Gaussian kernel with a standard deviation of σ; g(x) is an edge indicator, whose value is approximately 1 in the homogeneous region and approximately 0 near the water-land boundary; Z represents the fusion result of LoG and ROEWA weighted by m and n respectively, M LoG and M ROEWArespectively represent the detection results of the LoG operator and the ROEWA operator;

[0122] 1) LoG operator

[0123] The LoG operator is a second-order difference edge detection operator, defined as:

[0124] M LoG (x, y) = Δ[K σ *I(x, y)] = I(x, y)*LoG (14)

[0125] where K σ represents the Gaussian kernel with standard deviation σ, Δ represents the Laplacian operator, and M LoG (x, y) represents the image detected by the LoG operator; the LoG operator is positive in the dark part of the image edge and negative in the bright part of the edge. Both slow-changing edges and sharp edges can be detected through the zero-crossing points of positive and negative values; therefore, the LoG operator is detected for the land-water edge of the SAR image mainly in the form of steps;

[0126] 2) ROEWA operator

[0127] The ROEWA operator is a multi-edge model based on linear minimum mean square error, applicable to the speckle noise in SAR images; in the two-dimensional case, the definition of the ROEWA operator is:

[0128] h 2-D (x, y) = h(x)h(y) (15)

[0129] where h(x) and h(y) both represent one-dimensional filters, and h(x) is given by:

[0130] h(x) = Ce -k|x| (16)

[0131] where k represents the filtering coefficient and C represents the regularization constant;

[0132] The edge intensity detected by the ROEWA operator is jointly given by the horizontal component and the vertical component:

[0133]

[0134] where M X (x, y) is the horizontal edge intensity component, and M Y (x, y) is the vertical edge intensity component;

[0135] To calculate the horizontal edge intensity component of the image I(x, y), first smooth each column of I(x, y) using the filter h(y), and then use the causal filter h 1(·) and the anti-causal filter h 2 (·) is convolved with each row of the smoothed image to obtain the following exponentially weighted average μ X1 and μ X2 :

[0136]

[0137] where * represents the convolution operation; when μ X1 and μ X2 are normalized to be unbiased, the horizontal edge intensity component M X of the image is given by:

[0138]

[0139] In addition, the vertical edge intensity component M Y (x, y) is calculated using a similar process;

[0140] (3) Global energy term

[0141] A global energy function that makes full use of the global statistical features of the SAR image is introduced:

[0142]

[0143] where I G (x) is used to fit the global image; for the two-phase case, the expression of I G (x) is:

[0144] I G (x) = c 1 H ε (φ) + c 2 (1 - H ε (φ)) (21)

[0145] where c 1 and c 2 represent the average intensities of the two regions inside and outside the curve C.

[0146] Preferably, the lake shoreline detection result output module is specifically as follows:

[0147] The curve C is evolved by the gradient descent flow method, and according to the variational principle, the Euler-Lagrange equation for minimizing the energy functional E is obtained:

[0148]

[0149] The energy functional E is with respect to the level set function φ;

[0150] Based on the above, the constructed energy functional E MRSF is expressed as:

[0151]

[0152] According to Equation (22), we get:

[0153]

[0154] Combining the above variational equation with the initial contour obtained by the coarse segmentation module, iteratively solve the minimum value of the energy functional to obtain the final fine lake shoreline.

[0155] Advantages of the present invention:

[0156] First of all, the present invention uses the BM3D algorithm to suppress the speckle noise of SAR images, and can better retain the image details while reducing noise. In the coarse segmentation stage, the OSTU algorithm combined with morphological processing is used to obtain the initial contour of water and land, which can effectively avoid unreasonable initial contour settings and prevent the level set function from falling into local minima during the evolution process. In the fine segmentation stage, the present invention constructs an energy function based on a new type of edge energy term and global energy term, which can effectively locate weak edges and fit global statistical information. Based on the initial segmentation result, the present invention uses the MRSF model to obtain the fine lake shoreline, which can better utilize the global statistical features and edge features of the image compared with the traditional RSF model. The present invention is more effective and accurate in extracting the lake shoreline from SAR images. Description of the Drawings

[0157] Figure 1 It is a flow chart of the method of the present invention.

[0158] Figure 2 It is the process of the BM3D algorithm.

[0159] Figure 3 It is a comparison diagram of the results of the method of the present invention and the manually drawn shoreline.

[0160] Figure 4 It is a system block diagram of the present invention. Detailed Embodiments

[0161] The following will describe in detail the preferred embodiments of the present invention with reference to the accompanying drawings.

[0162] As Figures 1-3 shown, the SAR image lake shoreline detection method of this embodiment is based on an improved region variable fitting energy (RSF) model, and the specific implementation steps are as follows:

[0163] Step 1), use the BM3D algorithm to suppress speckle noise:

[0164] The process of the BM3D algorithm is as Figure 2 shown. Figure 2For the input image z, a basic estimate is first filtered in the first step The filtering operation in the second step is based on the result generated by the first estimate to drive and perform actual filtering in the wavelet domain, thus greatly improving the filtering result. And each step includes three operations: 1) grouping; 2) collaborative filtering; 3) aggregation. In the grouping operation, considering the statistical characteristics of noise, similar target blocks are grouped into a 3D stack based on the non-local principle. In the collaborative filtering operation, a wavelet transform is performed on the 3D stack. Finally, all target blocks return to their original positions through weighted averaging in the aggregation operation.

[0165] After the SAR image is filtered by the above BM3D algorithm, the intensity of speckle noise is greatly reduced, and the image details can be better retained.

[0166] Step 2): Obtain the initial contour using the OSTU algorithm:

[0167] The filtered SAR image obtained in step 1) is combined with morphological processing using the OSTU algorithm to obtain an initial rough segmentation contour. Considering that variance can be used as a measure of the gray-level distribution uniformity of the SAR image, the present invention divides the SAR image into two parts, lakes and land, by maximizing the between-class variance, thereby obtaining the initial contour of the lakes.

[0168] Step 3): Construct the energy function of the MRSF model based on the new boundary energy term and the introduced global energy term:

[0169] Based on the RSF model, a novel energy function is constructed:

[0170] E MRSF = E RSF + ωE edge + ξE global (3)

[0171] where ω and ξ represent weight coefficients. The constructed energy function mainly includes three parts, the traditional RSF energy term, the edge energy term, and the global energy term. The representation forms of each part are as follows.

[0172] (1) RSF model

[0173] To solve the problem that the CV model cannot handle images with intensity inhomogeneity, researchers proposed the RSF model based on the region-variable fitting energy. In the RSF model, the local energy function is defined as:

[0174]

[0175] where λ 1 ,λ 2, μ, ν denote positive constants, K σ (x - y) represents a Gaussian kernel with standard deviation σ, and

[0176]

[0177] where, H ε represents a regularized Heaviside function, which is defined as:

[0178]

[0179] The smoothing function f 1 (·) and f 2 (·) approximate the image intensities inside and outside the contour in a certain region respectively.

[0180] In addition, this energy function simultaneously takes into account the distance regularization term P(φ) and the length term L(φ), and their definitions are respectively:

[0181]

[0182]

[0183] where, δ ε (φ) represents a regularized Dirac function, defined as:

[0184]

[0185] P(φ) can maintain the horizontality and regularity of the function, and L(φ) helps to smooth the curve C.

[0186] By minimizing the level set function φ with respect to E RSF , the following variational equation can be obtained:

[0187]

[0188] (2) New edge energy term

[0189] Considering that there are a large number of weak edges and false edges in SAR images, the present invention constructs a new edge energy term based on the LoG operator and the ROEWA operator:

[0190] E edge =∫∫ Ω g(K σ *I(x, y)) × (Z - 0) 2 +(1 - g(K σ *I(x, y))) × (β × Z) 2 dxdy (11)

[0191] and

[0192] g(x) = exp{-α|x|} (12)

[0193] Z = mM LoG + nM ROEWA (13)

[0194] where α, β, m, n all represent weighted values, and K σ represents a Gaussian kernel with a standard deviation of σ. g(x) is an edge indicator, whose value is approximately 1 in homogeneous regions and approximately 0 near the water - land boundary. Z represents the fusion result of LoG and ROEWA weighted by m and n respectively, and M LoG and M ROEWA represent the detection results of the LoG operator and the ROEWA operator respectively.

[0195] 1) LoG operator

[0196] The LoG operator is a well - known second - order difference edge detection operator, mainly composed of two parts: the Gaussian operator and the Laplacian operator. Its definition is:

[0197] M LoG (x, y) = Δ[K σ * I(x, y)] = I(x, y) * LoG (14)

[0198] where K σ represents a Gaussian kernel with a standard deviation of σ, Δ represents the Laplacian operator, and M LoG (x, y) represents the image detected by the LoG operator. The LoG operator is positive in the dark part of the image edge and negative in the bright part of the edge. Both slowly changing edges and sharply changing edges can be detected through the zero - crossing points of positive and negative values. Therefore, the LoG operator can detect the water - land edge of SAR images mainly in the form of steps.

[0199] 2) ROEWA operator

[0200] The ROEWA operator is a multi - edge model based on linear minimum mean - square error, suitable for speckle noise in SAR images. In two - dimensional cases, the definition of the ROEWA operator is:

[0201] h 2-D (x, y) = h(x)h(y) (15)

[0202] where h(x) and h(y) both represent one - dimensional filters, and h(x) can be given by the following formula:

[0203] h(x) = Ce -k|x| (16)

[0204] where k represents the filtering coefficient and C represents the regularization constant.

[0205] The edge intensity detected by the ROEWA operator is jointly given by the horizontal and vertical components:

[0206]

[0207] where M X (x, y) is the horizontal edge intensity component, and M Y (x, y) is the vertical edge intensity component

[0208] To calculate the horizontal edge intensity component of the image I(x, y), first smooth each column of I(x, y) using the filter h(y), and then perform a convolution operation on each row of the smoothed image with the causal filter h 1 (·) and the anti-causal filter h 2 (·) to obtain the following exponentially weighted averages μ X1 and μ X2 :

[0209]

[0210] where * represents the convolution operation. When μ X1 and μ X2 are normalized to be unbiased, the horizontal edge intensity component M X of the image can be given by the following formula:

[0211]

[0212] In addition, the vertical edge intensity component M Y (x, y) can be calculated using a similar process.

[0213] (3) Global energy term

[0214] Considering the limitation that the traditional RSF model only uses local region fitting, and combining the concept of relative entropy, the present invention introduces a global energy function that can fully utilize the global statistical features of SAR images:

[0215]

[0216] where I G (x) is used to fit the global image. For the two-phase case, the expression of I G (x) is:

[0217] I G (x) = c 1 H ε (φ) + c 2 (1 - H ε (φ)) (21)

[0218] Among them, c 1 and c 2 represent the average intensities of the two regions inside and outside the curve C.

[0219] Combining the above energy terms, the energy function of the MRSF model is obtained.

[0220] Step 4): Use the result of the rough segmentation as the initial contour of the MRSF model to evolve and obtain the final fine lake shoreline:

[0221] Use the initial contour obtained in Step 2 as the initial contour of the MRSF model.

[0222] Evolve the curve C by the gradient descent flow method. According to the variational principle, the Euler - Lagrange equation for minimizing the energy functional E can be obtained:

[0223]

[0224] Here, the energy functional E is with respect to the level set function φ.

[0225] The energy functional E constructed by the MRSF model MRSF can be expressed as:

[0226]

[0227] According to Equation (22), it can be obtained that:

[0228]

[0229] Iteratively solve the minimum value of the energy functional to obtain the final fine lake shoreline.

[0230] Applying the above processing procedure, the fine lake shoreline contour (A) as shown in Figure 3 can be obtained, and it can better match the real lake shoreline (B) depicted manually.

[0231] The data used in this embodiment is the Sentinel - 1A single - look SAR image.

[0232] Table 1 shows the experimental parameters used in this illustrative embodiment.

[0233] Table 1 Experimental Parameters

[0234]

[0235]

[0236] According to the parameters in Table 1, the lake shoreline (A) obtained by using the method of the present invention is as shown in Figure 3As shown, it can be seen that the shoreline detected by the present invention is close to the true shoreline (B) depicted manually, and can better retain details such as some weak edges and concave edges.

[0237] Therefore, it can be concluded that the SAR image shoreline detection algorithm of the present invention based on the novel edge energy and global energy is effective and accurate.

[0238] As Figure 4 shown, this embodiment discloses a SAR image shoreline detection system based on the MRSF model, which includes the following modules:

[0239] Speckle noise suppression module: The registered SAR image is processed by the BM3D algorithm to suppress speckle noise; in the speckle noise suppression module, the BM3D algorithm includes two steps: The first step is that the original noisy image is filtered to obtain a basic estimate; the second step is that the actual filtering is performed based on the estimate generated in the first step in the second step; both steps perform the following three operating steps:

[0240] 1) Grouping: In the grouping operation, considering the statistical characteristics of the noise, similar target blocks are grouped into a 3D stack based on the non-local principle;

[0241] 2) Collaborative filtering: In the collaborative filtering operation, a wavelet transform is performed on the 3D stack;

[0242] 3) Aggregation: All target blocks return to their original positions through weighted averaging in the aggregation operation.

[0243] Coarse segmentation module: The denoised SAR image is coarsely segmented by the OSTU algorithm to obtain the initial boundary between water and land; in the coarse segmentation module, the SAR image is divided into two parts, namely, the lake and the land, by maximizing the between-class variance, so as to obtain the initial contour of the lake.

[0244] Energy function module for constructing the MRSF model: Based on the novel boundary energy term and the introduced global energy term, the energy function of the MRSF model is constructed; the energy function module for constructing the MRSF model is specifically as follows:

[0245] On the basis of the RSF model, a novel energy function is constructed:

[0246] E MRSF =E RSF +ωE edge +ξE global (3)

[0247] Among them, ω and ξ represent weight coefficients; the constructed energy function mainly includes three parts: the traditional RSF energy term, the edge energy term, and the global energy term; the representation forms of each part are as follows:

[0248] (1) RSF model

[0249] In the RSF model, the local energy function is defined as:

[0250]

[0251] where λ 1 , λ 2 , μ, ν represent positive constants, K σ (x - y) represents a Gaussian kernel with standard deviation σ, and

[0252]

[0253] where H ε represents the regularized Heaviside function, defined as:

[0254]

[0255] The smoothing functions f 1 (·) and f 2 (·) approximate the image intensities inside and outside the contour C in a certain region respectively;

[0256] In addition, this energy function also takes into account the distance regularization term P(φ) and the length term L(φ), which are defined respectively as:

[0257]

[0258]

[0259] where δ ε (φ) represents the regularized Dirac function, defined as:

[0260]

[0261] P(φ) maintains the horizontality and regularity of the function, and L(φ) helps to smooth the curve C;

[0262] By minimizing the level set function φ with respect to E RSF using the steepest descent method, the following variational equation is obtained:

[0263]

[0264] (2) New edge energy term

[0265] Based on the LoG operator and the ROEWA operator, a new edge energy term is constructed:

[0266]

[0267] and

[0268] g(x) = exp{-α|x|} (12)

[0269] Z = mM LoG + nM ROEWA (13)

[0270] where α, β, m, and n all represent weighting values, and K σ represents a Gaussian kernel with a standard deviation of σ; g(x) is an edge indicator whose value is approximately 1 in homogeneous regions and approximately 0 near the water - land boundary; Z represents the fusion result of LoG and ROEWA weighted by m and n respectively, M LoG and M ROEWA represent the detection results of the LoG operator and the ROEWA operator respectively;

[0271] 1) LoG operator

[0272] The LoG operator is a second - order difference edge detection operator, defined as:

[0273] M LoG (x, y) = Δ[K σ * I(x, y)] = I(x, y) * LoG (14)

[0274] where K σ represents a Gaussian kernel with a standard deviation of σ, Δ represents the Laplacian operator, and M LoG (x, y) represents the image detected by the LoG operator; the LoG operator is positive in the dark part of the image edge and negative in the bright part of the edge. Both slowly changing edges and sharply changing edges can be detected through the zero - crossing points of positive and negative values; therefore, the LoG operator is used to detect the water - land edge of SAR images mainly in the form of steps;

[0275] 2) ROEWA operator

[0276] The ROEWA operator is a multi - edge model based on linear minimum mean - square error and is applicable to speckle noise in SAR images; in two - dimensional cases, the definition of the ROEWA operator is:

[0277] h 2-D (x, y) = h(x)h(y) (15)

[0278] where h(x) and h(y) both represent one - dimensional filters, and h(x) is given by:

[0279] h(x) = Ce -k|x| (16)

[0280] where k represents the filtering coefficient and C represents the regularization constant;

[0281] The edge strength detected by the ROEWA operator is given by both the horizontal and vertical components:

[0282]

[0283] where M X (x, y) is the horizontal edge strength component, and M Y (x, y) is the vertical edge strength component;

[0284] To calculate the horizontal edge strength component of the image I(x, y), first smooth each column of I(x, y) using the filter h(y), and then convolve each row of the smoothed image with the causal filter h 1 (·) and the anti-causal filter h 2 (·) to obtain the following exponentially weighted averages μ X1 and μ X2 :

[0285]

[0286] where * represents the convolution operation; when μ X1 and μ X2 are normalized to be unbiased, the horizontal edge strength component M X of the image is given by:

[0287]

[0288] In addition, the vertical edge strength component M Y (x, y) is calculated using a similar process;

[0289] (3) Global energy term

[0290] Introduce a global energy function that fully utilizes the global statistical features of the SAR image:

[0291]

[0292] where I G (x) is used to fit the global image; for the two-phase case, the expression of I G (x) is:

[0293] I G (x) = c 1 H ε (φ) + c 2 (1 - H ε (φ)) (21)

[0294] where c 1 and c 2Represents the average intensity of the two regions inside and outside the curve C.

[0295] Lake shoreline detection result output module: Uses the result of rough segmentation as the initial contour of the MRSF model to evolve and obtain the final fine lake shoreline. The lake shoreline detection result output module is as follows:

[0296] Evolve the curve C by the gradient descent flow method. According to the variational principle, obtain the Euler-Lagrange equation for minimizing the energy functional E:

[0297]

[0298] The energy functional E is with respect to the level set function φ;

[0299] Based on the above, the constructed energy functional E MRSF Is expressed as:

[0300]

[0301] According to Equation (22), obtain:

[0302]

[0303] Combine the above variational equation with the initial contour obtained by the rough segmentation module, and iteratively solve the minimum value of the energy functional to obtain the final fine lake shoreline.

[0304] Those of ordinary skill in the art in this technical field should recognize that the above embodiments are only used to illustrate the present invention, rather than as a limitation to the present invention. As long as it is within the scope of the present invention, changes and deformations of the above embodiments will fall within the protection scope of the present invention.

Claims

1. SAR image shoreline detection method based on the MRSF model, characterized in that it includes the following steps: Step 1), suppress the speckle noise of the registered SAR image by using the BM3D algorithm; Step 2), coarsely segment the SAR image obtained by noise reduction by using the OSTU algorithm to obtain the initial boundary between water and land; Step 3), construct the energy function of the MRSF model based on the new boundary energy term and the introduced global energy term; Step 4), use the result of the coarse segmentation as the initial contour of the MRSF model to evolve to obtain the final fine shoreline; Step 3) is specifically as follows: Based on the RSF model, construct a novel energy function: E MRSF = E RSF + ωE edge + ξE global (3) where ω and ξ represent weight coefficients; the constructed energy function mainly includes three parts: the traditional RSF energy term, the edge energy term, and the global energy term; the representation forms of each part are as follows: (1) RSF model In the RSF model, the definition of the local energy function is: where λ 1 , λ 2 , μ, ν denote positive constants, K σ (x - y) represents a Gaussian kernel with standard deviation σ, I(y) is the input image, φ is the zero level set function, and where H ε represents the regularized Heaviside function with a smoothing factor of ε and is defined as: Smoothing function f 1 (·) and f 2 (·) approximate the image intensities inside and outside the contour C in a certain region, respectively; In addition, this energy function simultaneously considers the distance regularization term P(φ) and the length term L(φ), and the definitions are respectively: Among them, represents the differential operator, δ ε (φ) represents the regularized Dirac function, which is defined as: P(φ) maintains the horizontality and regularity of the function, and L(φ) helps to smooth the curve C; Minimize the level set function φ with respect to E by the steepest descent method RSF to obtain the following variational equation: (2) Novel edge energy term Construct a new edge energy term based on the LoG operator and the ROEWA operator: and g(x) = exp{-α|x|} (12) Z = mM LoG + nM ROEWA (13) where α, β, m, and n all represent weighting values, and K σ represents a Gaussian kernel with a standard deviation of σ; g(x) is an edge indicator, whose value is approximately 1 in homogeneous regions and approximately 0 near the water-land boundary; Z represents the fusion result of LoG and ROEWA weighted by m and n respectively, M LoG and M ROEWA represent the detection results of the LoG operator and the ROEWA operator respectively; 1) LoG operator The LoG operator is a second-order difference edge detection operator, defined as: M LoG (x,y) = Δ[K σ *I(x,y)] = I(x,y) * LoG (14) Among them, K σ represents a Gaussian kernel with a standard deviation of σ, Δ represents the Laplacian operator, and M LoG (x, y) represents the image detected by the LoG operator; the LoG operator is positive in the dark part of the image edge and negative in the bright part of the edge, and both slowly changing edges and sharply changing edges can be detected through the zero-crossing points of positive and negative values; therefore, the LoG operator is detected for the land-water edge of the SAR image mainly in the form of steps; 2) ROEWA operator The ROEWA operator is a multi-edge model based on linear minimum mean square error, suitable for the speckle noise in SAR images; in the two-dimensional case, the definition of the ROEWA operator is: h 2-D (x,y) = h(x)h(y) (15) where h(x) and h(y) both represent one-dimensional filters, and h(x) is given by the following formula: h(x) = Ce -k|x| (16) where k represents the filtering coefficient and C represents the regularization constant; The edge intensity detected by the ROEWA operator is jointly given by the horizontal component and the vertical component: where M X (x, y) is the horizontal edge intensity component, and M Y (x, y) is the vertical edge intensity component; To calculate the edge intensity component in the horizontal direction of the image I(x, y), first, each column of I(x, y) is smoothed using the filter h(y), and then the causal filter h 1 (·) and the anti-causal filter h 2 (·) are convolved with each row of the smoothed image to obtain the following exponentially weighted averages μ X1 and μ X2 : where * represents the convolution operation; when μ X1 and μ X2 are normalized to be unbiased, the horizontal edge intensity component M X of the image is given by: In addition, calculate the vertical edge intensity component M Y (x, y); (3) Global energy term Introduce a global energy function that makes full use of the global statistical characteristics of the SAR image: Among them, I G (x) is used to fit the global image; for the two-phase case, I G (x) is expressed as: I G f(x) = c 1 H ε (φ) + c 2 (1 - H ε (φ))(21) where c 1 and c 2 represent the average intensities of the two regions inside and outside the curve C.

2. The SAR image shoreline detection method according to claim 1, characterized in that, in Step 1), the BM3D algorithm includes two steps: the first step, the original noisy image is filtered to obtain a basic estimated value; the second step, the actual filtering is carried out based on the estimated value generated in the first step in the second step; both steps are carried out as follows in three operation steps: 1) Grouping: In the grouping operation, considering the statistical characteristics of the noise, similar target blocks are concentrated in a 3D stack based on the non-local principle; 2) Collaborative filtering: In the collaborative filtering operation, perform wavelet transform on the 3D stack; 3) Aggregation: All target blocks return to their original positions through weighted average in the aggregation operation.

3. The SAR image shoreline detection method according to claim 1, characterized in that, in Step 2), the SAR image is divided into two parts, namely the lake and the land, by maximizing the between-class variance, so as to obtain the initial contour of the lake.

4. The SAR image shoreline detection method according to claim 1, characterized in that, Step 4) is specifically as follows: The curve C is evolved by the gradient descent method, and according to the variational principle, the Euler-Lagrange equation that minimizes the energy functional E is obtained: The energy functional E is related to the level set function φ; Based on the above, the constructed energy functional E MRSF is expressed as: According to formula (22), we get: The above variational equation is combined with the initial contour obtained in step 2), and the minimum value of the energy functional is iteratively solved to obtain the final refined lake shoreline.

5. SAR image lake shoreline detection system based on MRSF model, Its characteristics are Includes the following modules: Speckle noise suppression module: The BM3D algorithm is used to suppress the speckle noise of the registered SAR image; Rough segmentation module: The de-noised SAR image is roughly segmented using the OSTU algorithm to obtain the initial boundary of land and water; Construct the energy function module of the MRSF model: construct the energy function of the MRSF model based on the new boundary energy term and the introduced global energy term; Lake shoreline detection result output module: the coarse segmentation result is used as the initial contour evolution of the MRSF model to obtain the final fine lake shoreline; The energy function module for constructing the MRSF model is as follows: Based on the RSF model, a novel energy function is constructed: E MRSF = E RSF + ωE edge + ξE global (3) Among them, ω and ξ represent weight coefficients; the constructed energy function mainly includes three parts: traditional RSF energy term, edge energy term and global energy term; the representation of each part is as follows: (1) RSF model In the RSF model, the local energy function is defined as: where λ 1 , λ 2 , μ, ν denote positive constants, K σ (x - y) represents a Gaussian kernel with standard deviation σ, I(y) is the input image, φ is the zero level set function, and where, H ε represents the regularized Heaviside function with a smoothing factor of ε and is defined as: Smoothing function f 1 (·) and f 2 (·) are respectively approximated as the image intensities inside and outside the contour C in a certain region; In addition, the energy function takes into account the distance regularization term P(φ) and the length term L(φ), which are defined as: Among them, represents a differential operator, δ ε (φ) represents a regularized Dirac function, defined as: P(φ) maintains the regularity of the level and function, and L(φ) helps to smooth the curve C; Minimize the level set function φ with respect to E by the steepest descent method RSF to obtain the following variational equation: (2) New edge energy terms A new edge energy term is constructed based on the LoG operator and the ROEWA operator: and g(x)=exp{-α|x|} (12) Z = mM LoG + nM ROEWA (13) where α, β, m, and n all represent weighting values, and K σ represents a Gaussian kernel with a standard deviation of σ; g(x) is an edge indicator, whose value is approximately 1 in homogeneous regions and approximately 0 near the water-land boundary; Z represents the fusion result of LoG and ROEWA weighted by m and n respectively, M LoG and M ROEWA represent the detection results of the LoG operator and the ROEWA operator respectively; 1) LoG operator The LoG operator is a second-order difference edge detection operator, defined as: M LoG (x,y) = Δ[K σ *I(x,y)] = I(x,y) * LoG (14) Among them, K σ represents a Gaussian kernel with a standard deviation of σ, Δ represents the Laplacian operator, and M LoG (x, y) represents the image detected by the LoG operator; the LoG operator is positive in the dark part of the image edge and negative in the bright part of the edge. Both slowly changing edges and sharply changing edges can be detected through the zero-crossing points of positive and negative values; therefore, the LoG operator is detected for the land-water edge of the SAR image mainly in the form of steps; 2) ROEWA operator The ROEWA operator is a multi-edge model based on the linear minimum mean square error, which is suitable for coherent speckle noise in SAR images. In the two-dimensional case, the ROEWA operator is defined as: h 2-D (x,y) = h(x)h(y) (15) Where h(x) and h(y) both represent one-dimensional filters, and h(x) is given by: h(x) = Ce -k|x| (16) Among them, k represents the filter coefficient, and C represents the regularization constant; The edge strength detected by the ROEWA operator is given by the horizontal component and the vertical component: Among them, M X (x, y) is the horizontal edge intensity component, and M Y (x, y) is the vertical edge intensity component; To calculate the edge intensity component in the horizontal direction of the image I(x, y), first, each column of I(x, y) is smoothed using the filter h(y), and then the causal filter h 1 (·) and the anti-causal filter h 2 (·) are convolved with each row of the smoothed image to obtain the following exponentially weighted averages μ X1 and μ X2 : Among them, * represents the convolution operation; when μ X1 and μ X2 are normalized to be unbiased, the horizontal edge intensity component M X of the image is given by the following formula: In addition, the vertical edge strength component M is calculated with a similar process Y (x, y); (3) Global energy term A global energy function that fully utilizes the global statistical characteristics of SAR images is introduced: Among them, I G (x) is used to fit the global image; for the two-phase case, I G (x) is expressed as: I G (x) = c 1 H ε (φ) + c 2 (1 - H ε (φ)) (21) where c 1 and c 2 represent the average intensities of two regions inside and outside the curve C.

6. The SAR image lake shoreline detection system as claimed in claim 5, Its characteristics are: In the Speckle Suppression module, the BM3D algorithm consists of two steps: in the first step, the original noisy image is filtered to obtain a basic estimate; in the second step, the actual filtering is performed in the second step based on the estimate generated in the first step; both steps perform the following three steps: 1) Grouping: In the grouping operation, similar target blocks are concentrated in a 3D stack based on the non-local principle, considering the statistical characteristics of noise; 2) Collaborative filtering: In the collaborative filtering operation, the 3D stack is subjected to wavelet transform; 3) Aggregation: All target blocks return to their original positions through weighted averaging in the aggregation operation.

7. The SAR image shoreline detection system according to claim 5, characterized in that, in the rough segmentation module, the SAR image is divided into two parts, namely the lake and the land, by maximizing the between-class variance, so as to obtain the initial contour of the lake.

8. The SAR image shoreline detection system according to claim 5, characterized in that, the shoreline detection result output module is specifically as follows: evolve the curve C by the gradient descent flow method, and according to the variational principle, obtain the Euler-Lagrange equation for minimizing the energy functional E: the energy functional E is with respect to the level set function φ; Based on the above, the constructed energy functional E MRSF is expressed as: according to Equation (22), obtain: combine the above variational equation with the initial contour obtained by the rough segmentation module, and iteratively solve the minimum value of the energy functional to obtain the final fine shoreline.

Images