A global mesoscale frontal automatic identification method based on sliding window threshold

Abstract

The application discloses a global mesoscale frontal automatic identification method based on a sliding window threshold value, and comprises the following steps: performing pretreatment operation of erosion missing value and noise filtering on input data; calculating frontal intensity according to an improved Sobel operator; marking a candidate frontal area according to a sliding window threshold value; obtaining final frontal area product data by using multiple mathematical morphological operators; forming a candidate frontal line according to multiple direction frontal intensity extreme values and a mask of the frontal area product; and finally obtaining frontal line product data by using multiple mathematical morphological operators. The method provided by the application is not susceptible to noise, is not sensitive to absolute gradient intensity, can identify mesoscale fronts missed by existing algorithms, solves the problem that traditional marine frontal identification methods are sensitive to absolute gradient and cannot identify weak fronts, and the identified frontal has better continuity, and has good applicability to global or regional sea areas.

Description

Technical Field

[0001] This invention belongs to the fields of satellite remote sensing applications and ocean dynamics, specifically relating to an automatic identification algorithm for global ocean mesoscale fronts based on a sliding window threshold. Background Technology

[0002] Oceanic fronts are high-energy oceanic features existing in the upper ocean globally, generally defined as narrow regions with relatively high horizontal gradients of specific water properties (such as temperature, salinity, nutrients, and phytoplankton) compared to the surrounding water. The occurrence of fronts is usually associated with various oceanic forcing mechanisms, such as water mass convergence, upwelling, tidal mixing, river plumes, and sudden changes in water depth. Strong vertical transport caused by intense turbulent mixing and transfrontal secondary circulation can be observed in many fronts, bringing nutrient-rich water into the euphotic zone and promoting primary and secondary production. Due to enhanced turbulence and powerful water mass transport in both vertical and horizontal directions, fronts play a crucial role in ocean dynamics, air-sea interactions, marine ecology, biogeochemistry, fishery prediction, pollutant distribution, and military applications. The occurrence and characteristics of fronts are fundamental information for understanding the ecological and climatic impacts of frontal dynamics; therefore, the identification of oceanic fronts is not only of significant scientific importance but also has great practical application value.

[0003] Satellite observations and numerical simulations provide global and local ocean data, which are the only reasonable methods for acquiring and monitoring long-term frontal activity information. By designing an automatic front identification algorithm, frontal occurrence information can be extracted from massive amounts of satellite observation and numerical simulation data, avoiding the complexity and subjectivity of manual identification. Currently, ocean front identification mainly uses the gradient-based Canny algorithm and the histogram-based algorithm proposed by Cayula and Cornillon (1992).

[0004] However, because the Canny algorithm uses a fixed threshold, while it can effectively identify strong fronts, it struggles to identify relatively weak fronts. These weak fronts should not be ignored, as many are characterized by subsurface intensification, meaning a significant increase in the subsurface of the thermal gradient across the front. This fixed threshold parameter causes the Canny algorithm to miss a large number of fronts in seasons or waters with weak gradients, while performing relatively well in seasons or waters with strong gradients. This limitation may lead to misleading results when applying the Canny algorithm to identify long-term series and large areas, as such seasonal or regional variations in gradients should be considered fluctuations in frontal intensity, not fluctuations in the occurrence of the front. Furthermore, histogram algorithms are insensitive to gradients. While they can effectively avoid missing a large number of fronts in seasons or waters with weak gradients, compared to gradient algorithms, histogram algorithms typically use larger windows to ensure the accuracy of statistical tests. This makes them adept at identifying regular, long fronts, but less effective for complex, dense, and near-shore fronts. Summary of the Invention

[0005] To address the shortcomings of existing front identification algorithms, this invention proposes an automatic global ocean mesoscale front identification algorithm based on a sliding window threshold. This algorithm solves the problem that traditional ocean front identification methods are sensitive to absolute gradients and cannot identify weak fronts, and also improves the continuity of the identified fronts.

[0006] This invention is achieved using the following technical solution:

[0007] This invention is achieved using the following technical solution: an automatic identification method for global ocean mesoscale fronts based on a sliding window threshold, comprising the following steps:

[0008] Step A: Obtain the input data and perform preprocessing operations on the input data to remove missing values ​​and filter noise, thus obtaining the preprocessed mesh data;

[0009] Step B: Obtain frontal zone product data:

[0010] Step B1: Calculate the frontal intensity using the improved Sobel operator;

[0011] Step B2: Mark candidate frontal zones based on the sliding window threshold;

[0012] Step B3: Use multiple mathematical morphology operators to obtain the final frontal zone product data;

[0013] Step C: Obtain frontal line product data:

[0014] Step C1: Form candidate front lines based on the extreme values ​​of frontal intensity in multiple directions and the mask of the frontal zone product obtained in step B;

[0015] Step C2: Use multiple mathematical morphology operators to obtain the final frontal line product data;

[0016] Step D: Combining the frontal zone product obtained in step B and the final frontal line product obtained in step C, the final frontal identification result is obtained.

[0017] Furthermore, in step B2, based on the gradient matrix calculated in step B1, a threshold for each grid point is calculated using a sliding window method. The sliding window threshold is determined by the mean and standard deviation of the gradients around each grid point.

[0018] This method first calculates the standard deviation and mean of the gradient within a 1.5° window around each grid point. Simultaneously, it uses the adjusted standard deviation, calculated by dividing the standard deviation of the gradient at each grid point by the power of the mean of the natural logarithm, to reduce missed fronts. The sum of the mean and adjusted standard deviation is considered a threshold for each grid point; grid points with gradient values ​​greater than this threshold are marked as binary matrices of candidate front regions. This sliding window thresholding method marks regions with locally large gradients as candidate front regions, avoiding the loss of numerous weak fronts due to a fixed threshold across the entire field in traditional gradient algorithms.

[0019] Furthermore, in step B3, based on the binary matrix of candidate frontal zones calculated in step B2, the mathematical morphology operators "CLOSE" and "SPUR" are used to fill gaps with a width (smaller) of less than 3 grid points, while simultaneously smoothing the identified frontal zones. Then, frontal zones smaller than 50 km are removed. 2 The smaller frontal zones form the final frontal zone product.

[0020] The smaller frontal zones to be deleted vary depending on the user's focus; if the area is less than 50km², it is not considered. 2 For frontal zones, areas smaller than 50km² should be deleted. 2 The frontal zone.

[0021] Furthermore, in step B1, the gradient matrix of the grid data processed in step A is calculated according to the improved Sobel operator, and the gradient data is further converted to units of ℃ / 100km, while a natural logarithmic transformation is performed to make it present an approximately normal distribution.

[0022] The gradient calculated by the original Sobel operator contains a large amount of sub-mesoscale information, affecting the identification results of mesoscale fronts. The improved operator introduces a parameter n to adjust the kernel spacing of the Sobel operator, filtering out sub-mesoscale information <50km and other noise interference, thus better identifying mesoscale fronts. The calculation formula for the improved Sobel operator is as follows:

[0023]

[0024]

[0025]

[0026] In the formula, t i,j G represents the data of the grid point in the i-th row and j-th column of the grid data after processing in step A. i,j,n,x and G i,j,n,y Let MG be the gradient vector of the grid point in the i-th row and j-th column. i,j,n The gradient value calculated for the Sobel operator with an interval of n at the i-th row and j-th column.

[0027] Furthermore, the parameter n of the improved Sobel operator is adjusted according to the resolution of the input grid data to filter out sub-mesoscale information <50km and other noise. The parameter n must be an integer and is calculated as half of 50km divided by the grid data resolution. If the resolution of the input grid data is 5km, then the parameter n is 5.

[0028] Furthermore, in step C1, based on the gradient matrix calculated in step B1, grid points whose gradient values ​​are local extrema in any direction (east-west, south-north, southeast-northwest, southwest-northeast) are marked as binary matrices of candidate front lines. At the same time, missing value masking and frontal zone masking are performed on the binary matrix to delete candidate front lines located outside the missing values ​​and frontal zones.

[0029] The grid point with a local extremum in any direction is calculated based on the data in the 8-neighborhood of each grid point. If the gradient value of a grid point is simultaneously greater than the gradient values ​​of the grid points to its left and right, top and bottom, southeast and northwest, or southwest and northeast, then that grid point is a grid point with a local extremum in any direction.

[0030] Furthermore, in step C2, based on the binary matrix of the candidate frontal lines calculated in step C1, the discontinuous frontal lines interrupted by 1–2 grid points are connected using the mathematical morphology operators “CLOSE”, “SPUR”, “FILL”, and “THIN”, with a filling area (smaller) less than 50 km². 2 The holes are created, forming frontal lines with a width of one grid point. The shorter frontal lines are then removed to form the final frontal line product.

[0031] The removal of shorter frontal lines varies depending on the user's focus. If the user does not focus on frontal lines shorter than 50km, then frontal lines shorter than 50km will be deleted.

[0032] Furthermore, in step A, the preprocessing operation is performed in the following manner:

[0033] Inverse distance weighted extrapolation is performed on the input data to assign virtual values ​​to missing values ​​in the land or cloud mask. Filtering is then performed to reduce data noise and sub-mesoscale signals. The calculation formula for the inverse distance weighted interpolation method is as follows:

[0034]

[0035]

[0036] In the formula, x0 and y0 are the latitude and longitude positions of the grid points with missing values, respectively; x i and y i Let T be the latitude and longitude of the i-th non-missing value grid point, and T be the latitude and longitude of the grid point. i and W i , where represents the ocean temperature value and weight of the i-th non-missing value grid point, respectively; N is the total number of non-missing value grid points within the search radius, and T0 is the ocean temperature value of the missing value grid point calculated using inverse distance weight interpolation or the virtual ocean temperature value assigned to the land grid point.

[0037] Furthermore, in step A, the input data is selected from any one of the uniform grid data of ocean temperature and other elements obtained by global or regional satellite remote sensing or ocean dynamics model simulation; ocean temperature and other elements include sea surface temperature, sea surface salinity, sea surface chlorophyll, and suspended matter concentration obtained by satellite remote sensing observation or numerical simulation; when identifying the sea surface chlorophyll and suspended matter concentration fronts, the chlorophyll and suspended matter concentrations need to be logarithmically converted to a normal distribution.

[0038] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0039] This invention identifies frontal zones using a sliding window threshold algorithm. The sliding threshold is determined by the mean and standard deviation of the gradient around each grid point, rather than the fixed double threshold in the Canny algorithm. This allows the frontal identification algorithm proposed in this invention to avoid the problem of traditional gradient algorithms missing a large number of fronts in seasons or waters with weak gradients. Furthermore, the identified fronts are more objective and accurate, and the histogram algorithm also avoids the problem of difficulty in identifying complex, dense, and near-shore fronts.

[0040] In addition, by calculating gradients using the improved Sobel operator, the sensitivity of traditional algorithms to noise and sub-mesoscale signals is avoided, making the recognition results more focused on mesoscale fronts. At the same time, by introducing a set of mathematical morphology operators to generate fronts, the front information of the recognition is smoother and more continuous, and the recognition of fronts has good continuity, avoiding interference from noise and sub-mesoscale signals.

[0041] Moreover, it can be easily, widely, and quickly applied to major water bodies around the world. It is also applicable to satellite observation data and numerical simulation data of various marine parameters such as temperature, salinity, and chlorophyll. It can also be used for the identification of operational frontal forecasts, providing a foundation for further research on the climate and ecological effects of fronts. It also helps in finding fishing grounds, selecting aquaculture sites, and designing marine protected areas. Attached Figure Description

[0042] Figure 1 This is a flowchart of the frontal recognition algorithm described in an embodiment of the present invention;

[0043] Figure 2 This is a schematic diagram of the frontal recognition result in an embodiment of the present invention;

[0044] Figure 3 This is a schematic diagram of the frontal recognition results from the Canny algorithm. Detailed Implementation

[0045] To better understand the above-described objects, features, and advantages of the present invention, the present invention will be further described below in conjunction with the accompanying drawings and embodiments. Many specific details are set forth in the following description to provide a thorough understanding of the present invention; however, the present invention may be practiced in other ways than those described herein, and therefore, the present invention is not limited to the specific embodiments disclosed below.

[0046] This embodiment proposes an automatic identification method for global ocean mesoscale fronts based on a sliding window threshold, such as... Figure 1 As shown, it includes the following steps:

[0047] Step A: Obtain the input data and perform preprocessing operations on the input data to remove missing values ​​and filter noise;

[0048] Step B: Obtain frontal zone product data:

[0049] (1) Calculate the frontal intensity and obtain gradient data based on the improved Sobel operator;

[0050] (2) Determine the sliding window threshold based on the mean and standard deviation of the gradient around each grid point, and mark the candidate frontal zone based on the sliding window threshold;

[0051] (3) Use multiple mathematical morphology operators to obtain the final frontal zone product data;

[0052] Step C: Obtain frontal line product data:

[0053] (1) Candidate front lines are formed based on the extreme values ​​of frontal intensity in multiple directions and the mask of the frontal zone product obtained in step B.

[0054] (2) The final frontal line product data is obtained by using multiple mathematical morphology operators;

[0055] Step D: Combining the frontal zone product obtained in step B and the final frontal line product obtained in step C, the final frontal identification result is obtained.

[0056] Specifically, this embodiment uses SST (sea surface temperature) fusion data obtained by thermal infrared radiation inversion from Advanced High Resolution Radiometer (AVHRR) and Along-the-Road Scanning Radiometer (ATSR) series sensors. This data is provided by the Copernicus Marine Environmental Monitoring Service (CMEMS), with a spatial resolution of 5 km and a temporal resolution of 1 day. The automatic global oceanic mesoscale front identification algorithm based on a sliding window threshold, as described in this invention, is used to identify fronts in the SST data of the Northwest Pacific region on January 1, 2012. All operations in this embodiment are implemented using MATLAB 2020a software programming, and the specific implementation method is as follows:

[0057] (1) Obtain SST fusion data for the region of interest. Download the SST data observed by satellite on January 1, 2012, from the CMEMS website, according to the latitude and longitude range of 140°E–180°E and 25°N–50°N.

[0058] (2) Using MATLAB iterative loops, inverse distance weighted extrapolation is performed on the input data obtained in step (1) to assign virtual values ​​to the missing values ​​on the land or cloud mask. The calculation formula for the inverse distance weighted extrapolation method is as follows:

[0059]

[0060]

[0061] In the formula, x0 and y0 are the latitude and longitude positions of the grid points with missing values, respectively; x i and y i Let T be the latitude and longitude of the i-th non-missing value grid point, and T be the latitude and longitude of the grid point. i and W i These are the ocean temperature value and weight of the i-th non-missing value grid point, respectively; N is the total number of non-missing value grid points within the search radius. In this embodiment, the search radius is set to 10 grid points (approximately 50 km), and T0 is the ocean temperature value of the missing value grid point calculated using inverse distance weight interpolation, or a virtual ocean temperature value assigned to the land grid point.

[0062] (3) The median filter and mean filter of the SST data obtained in step (2) are performed using the medfilt2 and imfilter functions in MATLAB to reduce data noise and sub-mesoscale signals. The filtering uses a 15×15km median filter and a 25×25km mean filter. In this embodiment, the filter size corresponds to 3×3 grid points and 5×5 grid points.

[0063] (4) Using MATLAB iterative loops, the gradient data of the SST data processed in step (3) is calculated based on the improved Sobel operator. The gradient data is multiplied by 2 to convert the unit to ℃ / 100km, and the log function is used to perform a natural logarithmic transformation to make it approximately normally distributed. The calculation formula of the improved Sobel operator is as follows:

[0064]

[0065]

[0066]

[0067] In the formula, t i,j For the data of the grid point in the i-th row and j-th column of the grid data after step (3), G i,j,n,x and G i,j,n,y Let MG be the gradient vector of the grid point in the i-th row and j-th column. i,j,n This calculates the gradient value at the i-th row and j-th column of the grid point for the Sobel operator with an interval of n. The parameter n is set to 5 to filter out interference from sub-mesoscale signals <50km, thus better identifying mesoscale fronts.

[0068] (5) Using the gradient matrix calculated in step (4) of the MATLAB loop iteration, calculate the standard deviation and mean of the gradient within a 1.5° window around each grid point. Use MATLAB array operations to divide the standard deviation of the gradient of each grid point by the mean power of the natural logarithm as the adjusted standard deviation. Then mark the grid points whose gradient values ​​are greater than the sum of the mean and the adjusted standard deviation as the binary matrix of the candidate frontal zone.

[0069] (6) Perform mathematical morphological operations on the binary image of the candidate front regions calculated in step (5) using the bwmorph function in the MATLAB Image Toolbox. Use the "CLOSE" and "SPUR" operations in the bwmorph function to fill small gaps and make the identified front regions smoother. Then use the bwpropfilt function in the MATLAB Image Toolbox to delete those smaller than 50 km. 2 The smaller frontal zones form the final frontal zone product. Figure 2 The gray area shows the frontal zone identified in this embodiment.

[0070] (7) Using MATLAB iterative loops, the grid points of local extrema in any direction (east-west, south-north, southeast-northwest, southwest-northeast) of the gradient data obtained in step (2) are marked as a binary matrix of candidate front lines. At the same time, the grid points in the binary matrix that belong to the land area and the original SST missing value area are set to 0 using matrix operations. Furthermore, the frontal zone products formed in step (6) are treated as a mask using MATLAB array operations, and the grid points in the binary matrix that do not belong to the frontal zone are set to 0, thus forming candidate frontal line data.

[0071] (8) Use the bwmorph function in the MATLAB image toolbox to perform mathematical morphological operations on the binary matrix of the candidate front lines calculated in step (7). Use the "CLOSE" and "SPUR" functions in the bwmorph function to connect the discontinuous front lines interrupted by 1–2 grid points, and use the "FILL" operator to fill the gaps less than 50 km. 2 The holes are filled, and the "THIN" operator is used to form a front line with a width of 1 grid point. Finally, the bwpropfilt function is used to remove front lines with a length of less than 5 grid points to form the final product. Figure 2 The black lines indicate the frontal lines identified in this embodiment.

[0072] To further demonstrate the advantages of this invention, this embodiment compares a front identification experiment using the Canny algorithm with the global ocean mesoscale front automatic identification algorithm based on a sliding window threshold proposed in this invention:

[0073] Control experiment: Frontal identification was performed using SST data from the Northwest Pacific region on January 1, 2012, using the Canny algorithm. All operations in this control experiment were implemented using MATLAB 2020a software. The specific implementation method is as follows:

[0074] (1) Obtain the SST fusion data of the region of interest. The SST data of the algorithm is the same as the SST data in the implementation case.

[0075] (2) Use the edge function in the MATLAB image toolbox to perform image edge recognition on the SST data in step (1). The edge detection method parameters are set to the "Canny" algorithm, and the sensitivity thresholds are set to 0.03 and 0.1. Figure 3 The black lines show the frontal results identified by the Canny algorithm.

[0076] By comparison Figure 2 and Figure 3 It can be seen that: Figure 3The Canny algorithm's identification of fronts exhibits a significant problem in weak gradient waters, missing a large number of fronts, such as those in the southern and northern regions of the Kuroshio and Oyashio extensions. Furthermore, the fronts identified by the Canny algorithm show obvious discontinuities, broken frontal lines, and interference from sub-mesoscale signals and noise. Typically, in statistical analysis, mesoscale fronts formed by the Kuroshio and Oyashio extensions are easily disturbed by information from numerous sub-mesoscale fronts. Figure 2 The fronts identified by this invention are significantly smoother and more continuous. Some short sub-mesoscale signals have been removed, highlighting mesoscale front information more effectively. Furthermore, it can identify not only front lines but also frontal regions. The comparative experiment demonstrates that the proposed global ocean mesoscale front automatic identification algorithm based on a sliding window threshold avoids the problem of traditional gradient algorithms missing a large number of fronts in seasons or waters with weak gradients. This makes the identified fronts more objective and accurate, while also exhibiting good front continuity, avoiding interference from noise and sub-mesoscale signals, and thus better identifying mesoscale fronts.

[0077] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments for application in other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. An automatic identification method for global ocean mesoscale fronts based on a sliding window threshold, characterized in that, Includes the following steps: Step A: Obtain the input data and perform preprocessing operations on the input data to remove missing values ​​and filter noise, thus obtaining the preprocessed mesh data; Step B: Obtain frontal zone product data: Step B1: Calculate the gradient matrix of the grid data processed in Step A using the improved Sobel operator to obtain the gradient data. Further convert the gradient data to °C / 100 km units and perform a natural logarithmic transformation to make it approximately normally distributed. The calculation formula for the improved Sobel operator is as follows: ， , ， ， In the formula, For the grid data processed in step A, the first Line number Data of grid points, and For the first Line number Gradient vector of column grid points, The interval is The Sobel operator computes the first Line number Gradient values ​​of column grid points; Step B2: Determine the sliding window threshold based on the mean and standard deviation of the gradient around each grid point, mark the candidate frontal zones based on the sliding window threshold, and obtain the binary matrix of the candidate frontal zones; Candidate frontal zones are specifically marked in the following manner: (1) First, calculate the standard deviation and mean of the gradient within a 1.5° window around each grid point. At the same time, divide the standard deviation of the gradient of each grid point by the mean power of the natural logarithm to use as the adjustment standard deviation to reduce the omission of fronts. (2) The sum of the mean and the adjusted standard deviation is regarded as the threshold for each grid point, and the grid points with gradient values ​​greater than the threshold are marked as the binary matrix of the candidate frontal zone; Step B3: Based on the binary matrix of candidate frontal zones, use multiple mathematical morphology operators to obtain the final frontal zone product data; Step C: Obtain frontal line product data: Step C1: Form candidate front lines based on the extreme values ​​of frontal intensity in multiple directions and the mask of the frontal zone product obtained in step B; Step C2: Use multiple mathematical morphology operators to obtain the final frontal line product data; Step D: Combining the frontal zone product obtained in step B and the final frontal line product obtained in step C, the final frontal identification result is obtained.

2. The automatic identification method for global ocean mesoscale fronts based on sliding window thresholds according to claim 1, characterized in that, In step B3, based on the binary matrix of the candidate frontal zones calculated in step B2, multiple mathematical morphology operators are used to fill the gaps of small-width grid points, and then small frontal zones are deleted to form the final frontal zone product.

3. The automatic identification method for global ocean mesoscale fronts based on sliding window thresholds according to claim 1, characterized in that, In step C1, based on the gradient matrix calculated in step B1, grid points whose gradient values ​​are local extrema in any direction are marked as binary matrices of candidate front lines. At the same time, the binary matrix is ​​processed by missing value masking and frontal zone masking to delete candidate front lines located outside the missing values ​​and frontal zones.

4. The automatic identification method for global ocean mesoscale fronts based on sliding window thresholds according to claim 1, characterized in that, In step C2, based on the binary matrix of the candidate front lines calculated in step C1, multiple mathematical morphology operators are used to connect the discontinuous front lines that are interrupted by 1–2 grid points, fill the small holes, and form front lines with a width of 1 grid point. Then, the short front lines are removed to form the final front line product.

5. The automatic identification method for global ocean mesoscale fronts based on sliding window thresholds according to claim 1, characterized in that, In step A, during the preprocessing operation, inverse distance weighted extrapolation is performed on the input data. The principle is as follows: In the formula, and These represent the latitude and longitude coordinates of the grid points with missing values; and The first The latitude and longitude positions of the non-missing value grid points, and and The first Ocean temperature values ​​and weights for each non-missing grid point; This represents the total number of non-missing value grid points within the search radius. This is used to calculate ocean temperature values ​​for missing grid points using inverse distance weighted interpolation, or to assign virtual ocean temperature values ​​to land grid points.

6. The automatic identification method for global ocean mesoscale fronts based on sliding window thresholds according to claim 1, characterized in that, In step A, the input data is selected as uniform grid data of any marine environmental element acquired by global or regional satellite remote sensing or simulated by marine dynamics models. The marine environmental element includes sea surface temperature, sea surface salinity, sea surface chlorophyll and suspended matter concentration obtained by satellite remote sensing observation or numerical simulation.

Images