Early warning method for sea surface temperature anomaly detection

Abstract

The invention provides an early warning method for sea surface temperature anomaly detection, which belongs to the technical field of ocean forecasting, and comprises the following steps of: constructing a multi-source sea temperature data fusion system and a multi-scale adaptive grid, generating a high-resolution temperature field by adopting ensemble Kalman filtering data assimilation, extracting an abnormal component by utilizing ensemble empirical mode decomposition, and carrying out early warning on the abnormal component. Multi-level anomaly discrimination is performed based on a sparse coding recognition model and fractal dimension mutation detection, anomaly types are distinguished in combination with an atmospheric compulsive event feature library and a random forest classifier, and partial differential equation inverse problem reverse deduction is performed on ocean endogenous anomaly to reconstruct a three-dimensional anomaly structure. And finally, the early warning level is determined through the early warning decision function and the multi-dimensional indexes, and the technical problem that the real-time performance and the accuracy of sea surface temperature anomaly detection are difficult to guarantee at the same time is solved.

Description

Technical Field

[0001] This invention belongs to the field of marine forecasting technology, and more specifically, relates to an early warning method for detecting sea surface temperature anomalies. Background Technology

[0002] Sea surface temperature anomaly detection is crucial for marine disaster early warning and climate monitoring. Traditional methods for sea surface temperature anomaly detection mainly rely on statistical thresholding and numerical model forecasting, identifying anomalous events by setting fixed temperature deviation thresholds or establishing statistical models based on historical climate states. These methods have significant limitations in practical applications. Statistical thresholding is insufficiently sensitive to slowly evolving anomalous events, and numerical model forecasting depends heavily on initial field quality and consumes enormous computational resources. Both methods struggle to achieve high spatiotemporal resolution real-time monitoring in complex marine environments. In current marine observation systems, insufficient fusion of multi-source observation data, limited anomaly feature extraction methods, and insufficient physical mechanism inversion capabilities lead to high false alarm and false negative rates in anomaly detection. In other words, existing technologies face the challenge of simultaneously ensuring real-time performance and accuracy in sea surface temperature anomaly detection. Summary of the Invention

[0003] In view of this, the present invention provides an early warning method for detecting sea surface temperature anomalies, which can solve the technical problem that it is difficult to simultaneously ensure the real-time performance and accuracy of sea surface temperature anomaly detection in the prior art.

[0004] This invention is implemented as follows: It provides an early warning method for detecting sea surface temperature anomalies, comprising the following steps: Constructing a multi-source sea surface temperature observation data fusion system, collecting satellite remote sensing data, buoy measured data, and Argo profile data; performing quality control to remove outliers; and forming an initial observation dataset after unifying the spatiotemporal reference. Establishing a multi-scale adaptive grid system, calculating the temperature gradient field based on the initial observation dataset, setting a fine grid in regions where the temperature gradient is greater than a gradient threshold, and using a reference grid in regions where the temperature gradient is less than a gradient threshold. Employing ensemble Kalman filtering data assimilation technology, the initial observation dataset is integrated into the numerical oceanography model background field; the model's initial field is adjusted using a variational method to generate a physically consistent high-resolution temperature field, while simultaneously outputting a grid-by-grid analysis error covariance matrix. Ensemble empirical mode decomposition is performed on the physically consistent high-resolution temperature field, decomposing it into an intrinsic mode function sequence; instantaneous frequencies are extracted using Hilbert-Huang transform; non-stationary anomaly components are identified and the anomaly temperature field is reconstructed. The anomaly temperature field is then input into the temperature field sparse... The coding and identification model outputs a reconstruction error field and sparse coding coefficients. When the reconstruction error exceeds an adaptive threshold, it is marked as a suspected anomaly region. For suspected anomaly regions, fractal dimension mutation identification is performed. The fractal dimension time series is calculated using the difference box counting method. When the rate of change of the fractal dimension exceeds the fractal control limit, a primary warning is triggered. A feature database of atmospheric forcing events is established by real-time access to meteorological forecast information. A random forest classifier is used to determine the anomaly type corresponding to the primary warning region. When the anomaly duration exceeds the characteristic time scale of atmospheric disturbance and the spatial distribution does not conform to the wind field structure, it is determined to be an endogenous marine anomaly. For the region corresponding to the endogenous marine anomaly, a partial differential equation inverse problem is initiated for inverse deduction. Based on the heat conduction equation, the internal heat source distribution is inverted. Iterative optimization using the adjoint method and quasi-Newton method is employed to reconstruct the three-dimensional anomaly structure and estimate the anomaly intensity. An early warning decision function is established. The warning level index is calculated by comprehensively considering the reconstruction error, the rate of change of the fractal dimension, the anomaly duration, and the anomaly intensity. The warning response level is determined based on the interval to which the warning level index belongs, and the warning information is issued.

[0005] Specifically, the ensemble empirical mode decomposition step involves adding normally distributed white noise to a physically consistent high-resolution temperature field to form ensemble samples, performing empirical mode decomposition on each ensemble sample to obtain single-sample intrinsic mode functions, and averaging the single-sample intrinsic mode functions corresponding to all ensemble samples to eliminate the influence of white noise, thereby obtaining an intrinsic mode function sequence.

[0006] The Hilbert-Huang transform obtains an analytic signal by performing a Hilbert transform on the eigenmode functions in the eigenmode function sequence, and calculates the instantaneous frequency by calculating the derivative of the instantaneous phase of the analytic signal with respect to time. Abnormal fluctuations in the instantaneous frequency indicate the occurrence of non-stationary temperature anomalies.

[0007] The structure of the temperature field sparse coding recognition model is as follows: the input layer receives the spatial distribution data of the abnormal temperature field; the sparse coding layer uses an overcomplete dictionary to perform sparse linear representation of the spatial distribution data; the sparse coding layer solves the sparse coefficient vector through the orthogonal matching pursuit algorithm; the reconstruction layer reconstructs the temperature field using the atomic vectors and sparse coefficient vectors in the overcomplete dictionary; and the output layer calculates the root mean square error between the abnormal temperature field and the temperature field reconstructed by the reconstruction layer as the reconstruction error.

[0008] The overcomplete dictionary contains 256 atomic vectors, and the dimension of each atomic vector is consistent with the number of grid points in the anomalous temperature field.

[0009] Specifically, the steps for establishing the training dataset for the sparse coding recognition model of the temperature field are as follows: collecting normal sea surface temperature observation data from the past 10 years, removing data from the periods when known abnormal events occurred, normalizing the removed normal sea surface temperature observation data so that the temperature values ​​are distributed between 0 and 1, dividing the normalized temperature field into training samples according to spatial regions, with each training sample containing the temperature field evolution process for 7 consecutive days, forming a total of 50,000 training samples to constitute the training dataset.

[0010] The training steps of the temperature field sparse coding recognition model specifically involve initializing an overcomplete dictionary using the K-SVD dictionary learning algorithm, iteratively executing the sparse coding stage and the dictionary update stage. In the sparse coding stage, the overcomplete dictionary is fixed, and the orthogonal matching pursuit algorithm is used to solve for the sparse coefficients of all training samples in the training dataset. In the dictionary update stage, the sparse coefficients are fixed, and singular value decomposition is used to update the atomic vectors in the overcomplete dictionary one by one to minimize the reconstruction error. After 300 iterations, the overcomplete dictionary converges.

[0011] The sparsity constraint parameters in the sparse coding layer are determined based on the grid density ratio of the multi-scale adaptive grid system, the trace value of the grid-by-grid analysis error covariance matrix, and the energy concentration of the anomalous temperature field. The sparsity constraint parameters are positively correlated with the grid density ratio, negatively correlated with the trace value of the grid-by-grid analysis error covariance matrix, and positively correlated with the energy concentration of the anomalous temperature field.

[0012] The adaptive threshold is calculated based on the mean of historical reconstruction errors, the standard deviation of historical reconstruction errors, and the current density of observation stations in the sea area. The adaptive adjustment coefficient is equal to the ratio of the standard deviation of historical reconstruction errors to the standard deviation of standard reconstruction errors, multiplied by the square root of the ratio of the current density of observation stations in the sea area to the standard density of observation stations, and then multiplied by the ratio of the mean of historical reconstruction errors to the mean of standard reconstruction errors. The product of these three terms is then cubed.

[0013] The fractal dimension mutation identification algorithm is implemented after the temperature field sparse coding identification model identifies a suspected abnormal region. Specifically, it extracts the temperature field data of the suspected abnormal region, converts the temperature value in the temperature field data into a three-dimensional surface where the third dimension is the temperature value, sets a series of box sizes that increase from twice the minimum grid spacing to one-quarter of the region scale, calculates the minimum number of boxes required to cover the entire three-dimensional surface for each box size, fits a straight line in a logarithmic coordinate system to the logarithm of the box size and the logarithm of the minimum number of boxes, and uses the negative value of the slope of the straight line as the fractal dimension at the current moment.

[0014] Specifically, a 72-hour sliding time window was used to calculate the fractal dimension time series. The first difference of the fractal dimension time series was used to obtain the fractal dimension change rate. Based on the distribution of the fractal dimension change rate during historical normal periods, fractal control limits at a 99% confidence level were determined.

[0015] Specifically, the implementation of the partial differential equation inverse problem back derivation algorithm after determining the endogenous anomaly of the ocean involves using the surface temperature observation of the region corresponding to the endogenous anomaly of the ocean as known conditions, establishing a three-dimensional heat conduction forward model, initializing the internal heat source distribution to zero, setting the boundary condition as the historical average heat flux, running the three-dimensional heat conduction forward model to obtain the simulated temperature field, and calculating the weighted sum of squares of the residuals between the surface temperature observation and the simulated temperature field as the objective functional.

[0016] The process involves constructing an adjoint equation and solving for the adjoint variables by inverse integration. The gradient information of the objective functional with respect to the internal heat source distribution is calculated using the adjoint variables. The BFGS algorithm, which uses a quasi-Newton method, is used to update the internal heat source distribution and the boundary heat flux. A Tikhonov regularization term is added to constrain the sum of squares of the spatial second derivative of the internal heat source distribution.

[0017] In the multi-scale adaptive grid system, the grid spacing of the fine grid is one-quarter of the base grid spacing, and the gradient threshold is taken as the 90th percentile of the statistical distribution of the temperature gradient field under historical normal sea conditions as the fine grid trigger condition.

[0018] The random forest classifier consists of multiple decision trees, each trained on a random subset of training samples and a subset of features. The final classification result is determined by voting from all decision trees.

[0019] The early warning decision function calculates the warning level index by comprehensively considering the reconstruction error, the fractal dimension change rate, the duration of the anomaly, and the intensity of the anomaly. The warning response level is determined according to the interval to which the warning level index belongs. The warning response levels are divided into Level 1, Level 2, Level 3, and Level 4 warnings.

[0020] This invention establishes a complete technical chain from data fusion to anomaly identification and mechanism tracing by constructing an adaptive grid system for multi-source data fusion, integrating a multi-level discrimination mechanism of sparse coding anomaly identification and fractal dimension mutation detection, and combining a physical mechanism inversion method based on inverse deduction of partial differential equations. Sparse coding utilizes the inherent low-dimensional manifold structure of the normal sea surface temperature field to achieve efficient anomaly detection; fractal dimension mutation identification captures nonlinear changes in spatial structure; and inverse deduction of partial differential equations maps surface observations to three-dimensional anomaly sources. These three methods mutually verify each other at different physical levels, reducing the false alarm rate. The adaptive grid system dynamically adjusts the resolution according to the temperature gradient, and the integration of Kalman filtering data assimilation technology with multi-source observation data improves the quality of the initial field, providing high-precision input for subsequent anomaly detection. Simultaneously, multi-scale processing strategies and parallel algorithms ensure real-time performance. In summary, this invention solves the technical problem mentioned in the background art of simultaneously ensuring real-time performance and accuracy in sea surface temperature anomaly detection. Attached Figure Description

[0021] Figure 1 This is a flowchart of the method of the present invention.

[0022] Figure 2 This is a high-resolution spatial distribution map of the sea surface temperature field.

[0023] Figure 3 This is a time series plot of the instantaneous frequency of the intrinsic mode functions.

[0024] Figure 4 This is a graph showing the time evolution of the fractal dimension. Detailed Implementation

[0025] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0026] like Figure 1 The diagram shown is a flowchart of an early warning method for detecting sea surface temperature anomalies provided by this invention. This method includes the following steps:

[0027] S01. Construct a multi-source sea surface temperature observation data fusion system, collect satellite remote sensing data, buoy measured data and Argo profile data, perform quality control on the satellite remote sensing data, buoy measured data and Argo profile data to remove outliers, and form an initial observation dataset after unifying the spatiotemporal reference.

[0028] S02. Establish a multi-scale adaptive grid system, calculate the temperature gradient field based on the initial observation dataset, set a fine grid in the region where the temperature gradient is greater than the gradient threshold, and set the grid spacing of the fine grid to 1 / 4 of the reference grid spacing. Use the reference grid in the region where the temperature gradient is less than the gradient threshold.

[0029] S03. Using ensemble Kalman filter data assimilation technology, the initial observation dataset is integrated into the background field of the numerical ocean model. The initial field of the numerical ocean model background field is adjusted by variational method to generate a physically consistent high-resolution temperature field, and the error covariance matrix is ​​output on a grid-by-grid basis.

[0030] S04. Perform ensemble empirical mode decomposition on the physically consistent high-resolution temperature field, decompose the physically consistent high-resolution temperature field into a sequence of intrinsic mode functions, extract the instantaneous frequency of each intrinsic mode function in the sequence of intrinsic mode functions through Hilbert-Huang transform, identify non-stationary anomalous components and reconstruct the anomalous temperature field.

[0031] S05. Input the abnormal temperature field into the temperature field sparse coding recognition model. The temperature field sparse coding recognition model outputs the reconstruction error field and sparse coding coefficients. When the reconstruction error in the reconstruction error field exceeds the adaptive threshold, it is marked as a suspected abnormal region.

[0032] S06. Identify fractal dimension mutations in suspected abnormal regions, calculate the fractal dimension time series using the differential box counting method, and trigger a primary warning when the rate of change of the fractal dimension in the fractal dimension time series exceeds the fractal control limit.

[0033] S07. Real-time access to meteorological forecast information to establish an atmospheric forcing event feature database, use random forest classifier to identify the anomaly type in the area corresponding to the primary warning, and determine the marine endogenous anomaly when the duration of the anomaly exceeds the atmospheric disturbance characteristic time scale and the spatial distribution does not conform to the wind field structure.

[0034] S08. Initiate the inverse deduction of the partial differential equation for the region corresponding to the endogenous anomaly in the ocean, invert the distribution of internal heat sources based on the heat conduction equation, and use the adjoint method and quasi-Newton method for iterative optimization to reconstruct the three-dimensional anomaly structure and estimate the anomaly intensity.

[0035] S09. Establish an early warning decision function, calculate the early warning level index by comprehensively considering the reconstruction error in the reconstruction error field, the rate of change of fractal dimension in the fractal dimension time series, the duration of the anomaly, and the intensity of the anomaly, determine the early warning response level based on the interval to which the early warning level index belongs, and issue early warning information.

[0036] The specific steps of the ensemble empirical mode decomposition include: adding normally distributed white noise to a physically consistent high-resolution temperature field to form an ensemble sample; performing empirical mode decomposition on each ensemble sample to obtain single-sample intrinsic mode functions (IMFs); averaging the IMFs corresponding to all ensemble samples to eliminate the influence of white noise, resulting in an IMF sequence. The Hilbert-Huang transform is used to obtain an analytic signal by performing a Hilbert transform on the IMFs in the IMF sequence; the instantaneous frequency is obtained by calculating the derivative of the instantaneous phase of the analytic signal with respect to time; and abnormal fluctuations in the instantaneous frequency indicate the occurrence of non-stationary temperature anomalies.

[0037] The structure of the temperature field sparse coding recognition model is as follows: the input layer receives the spatial distribution data of the abnormal temperature field; the sparse coding layer uses an overcomplete dictionary to perform sparse linear representation of the spatial distribution data; the overcomplete dictionary contains 256 atomic vectors, and the dimension of each atomic vector is consistent with the number of grid points in the abnormal temperature field; the sparse coding layer solves for the sparse coefficient vector using the orthogonal matching pursuit algorithm; the reconstruction layer reconstructs the temperature field using the atomic vectors and sparse coefficient vectors in the overcomplete dictionary; and the output layer calculates the root mean square error between the abnormal temperature field and the temperature field reconstructed by the reconstruction layer as the reconstruction error in the reconstruction error field. The steps for establishing the training dataset for the temperature field sparse coding recognition model specifically include: collecting 10 years of historical normal sea surface temperature observation data; removing data from periods when known abnormal events occurred; normalizing the removed normal sea surface temperature observation data so that the temperature values ​​are distributed between 0 and 1; dividing the normalized temperature field into training samples according to spatial regions; each training sample contains the temperature field evolution process for 7 consecutive days, forming a total of 50,000 training samples constituting the training dataset. The specific steps for training the temperature field sparse coding recognition model include: initializing an overcomplete dictionary using the K-SVD dictionary learning algorithm, iteratively executing the sparse coding stage and the dictionary update stage; fixing the overcomplete dictionary in the sparse coding stage and using the orthogonal matching pursuit algorithm to solve for the sparse coefficients of all training samples in the training dataset in the orthogonal matching pursuit algorithm; fixing the sparse coefficients in the dictionary update stage and using singular value decomposition to update the atomic vectors in the overcomplete dictionary one by one to minimize the reconstruction error; setting a sparsity constraint that a maximum of 10 atomic vectors are used for each training sample; after 300 iterations, the overcomplete dictionary converges; and the final overcomplete dictionary is saved for anomaly detection.

[0038] The sparsity constraint parameters in the sparse coding layer are determined based on the grid density ratio of the multi-scale adaptive grid system, the trace value of the grid-by-grid analysis error covariance matrix, and the energy concentration of the anomalous temperature field. The grid density ratio is defined as the ratio of the area of ​​the fine grid region to the total area. The trace value of the grid-by-grid analysis error covariance matrix reflects the overall uncertainty level after data assimilation. The energy concentration of the anomalous temperature field is calculated by wavelet transform, determining the energy proportion of the top five largest coefficients. The sparsity constraint parameters are positively correlated with the grid density ratio, negatively correlated with the trace value of the grid-by-grid analysis error covariance matrix, and positively correlated with the energy concentration of the anomalous temperature field. A weighted combination of the grid density ratio, the trace value of the grid-by-grid analysis error covariance matrix, and the energy concentration of the anomalous temperature field is mapped to an integer range of 5 to 15 as the final sparsity constraint parameters.

[0039] The adaptive threshold adjustment function is used to adjust the adaptive threshold of the temperature field sparse coding recognition model. The adaptive adjustment coefficient is calculated based on the mean of historical reconstruction errors, the standard deviation of historical reconstruction errors, and the current observation station density in the sea area. The adaptive adjustment coefficient is expressed as follows: the adaptive adjustment coefficient equals the ratio of the standard deviation of historical reconstruction errors to the standard deviation of standard reconstruction errors, multiplied by the square root of the ratio of the current observation station density to the standard observation station density, and then multiplied by the ratio of the mean of historical reconstruction errors to the mean of standard reconstruction errors. The cube root of this product is then obtained. When the adaptive adjustment coefficient a∈[0.5, 0.8), the adaptive threshold is set to 0.7 times the baseline threshold; when the adaptive adjustment coefficient a∈[0.8, 1.2], the adaptive threshold is set to the baseline threshold; when the adaptive adjustment coefficient a∈(1.2, 2.0], the adaptive threshold is set to 1.5 times the baseline threshold; and when the adaptive adjustment coefficient a>2.0, the adaptive threshold is set to 2.0 times the baseline threshold.

[0040] The principle of the sparse representation dictionary learning anomaly coding algorithm is as follows: It is assumed that the normal sea surface temperature field is represented by a few atomic vectors from an overcomplete dictionary through a sparse linear combination. The number of atomic vectors in the overcomplete dictionary is greater than the dimension of the temperature field, giving the overcomplete dictionary representation flexibility. The spatial distribution pattern of the normal temperature field corresponds to sparse coding vectors in the dictionary space. Anomaly temperature fields, deviating from the normal pattern, cannot be effectively reconstructed using the few atomic vectors in the overcomplete dictionary, resulting in large residuals during reconstruction. The dictionary learning process adaptively extracts typical temperature field structures from historical normal data as atomic vectors in the overcomplete dictionary. Each atomic vector captures a spatial distribution pattern. The orthogonal matching pursuit algorithm greedily selects the atomic vectors most relevant to the residuals during encoding to gradually approximate the target temperature field. Sparse constraints limit the number of atomic vectors used to ensure the discriminative power of the encoding. The implementation of the sparse representation dictionary learning anomaly coding algorithm in the temperature field sparse coding recognition model is as follows: After inputting the anomaly temperature field, it is first normalized to a unit norm vector, the residual vector is initialized to be equal to the unit norm vector, the atomic vector with the largest inner product with the current residual vector in the overcomplete dictionary is iteratively selected, the coefficient of the atomic vector is calculated so that the projection of the residual vector in the direction of the atomic vector is zero, the residual vector is updated and the atomic vector is repeatedly selected until the norm of the residual vector is less than the residual threshold or the sparsity constraint parameter is reached, and finally the selected atomic vector index and corresponding coefficient constitute the sparse coding coefficient. The temperature field is reconstructed using the selected atomic vector and corresponding coefficient and the reconstruction error in the reconstruction error field is calculated. The sparse representation dictionary learning anomaly coding algorithm establishes a compact dictionary representation of the spatial pattern of the temperature field by learning the inherent sparse representation structure of the normal temperature field. During anomaly detection, only fast sparse coding and error calculation are required, avoiding statistical modeling and parameter estimation. The physical interpretability of the atomic vectors in the overcomplete dictionary enables the traceability of anomaly sources. The overcomplete dictionary is updated online based on new observation data to adapt to the long-term evolution of the climate system. Sparse coding also achieves data compression to reduce storage and transmission burden. The linear complexity and parallelization characteristics of the sparse representation dictionary learning anomaly coding algorithm meet the real-time processing requirements of high-dimensional sea surface temperature fields. Compared with traditional methods, it reduces computational resource consumption and response latency while maintaining detection sensitivity.

[0041] The principle of the fractal dimension mutation identification algorithm is as follows: the spatial distribution of the sea surface temperature field exhibits self-similar fractal characteristics. Under normal conditions, the temperature gradient structure at different observation scales follows a power-law relationship, and the fractal dimension is quantified by scale invariance. The occurrence of abnormal events disrupts the original fractal structure of the temperature field, introducing new spatial scale characteristics, leading to a jump in the fractal dimension. The differential box counting method covers the temperature field with boxes of different sizes and counts the minimum number of boxes required for complete coverage. In a log-log coordinate system, the logarithm of the box size and the logarithm of the number of boxes have a linear relationship. The negative value of the slope of the linear relationship line is the fractal dimension. The fractal dimension is calculated step by step in a sliding time window to form a fractal dimension time series. Abnormal fluctuations in the rate of change of the fractal dimension in the fractal dimension time series indicate a mutation in the spatial structure of the temperature field. The fractal dimension mutation identification algorithm is implemented after the temperature field sparse coding identification model identifies a suspected abnormal region as follows: extract the temperature field data of the suspected abnormal region, convert the temperature values ​​in the temperature field data into a three-dimensional surface where the third dimension is the temperature value, set a series of box sizes increasing from twice the minimum grid spacing to 1 / 4 of the region scale, calculate the minimum number of boxes required to cover the entire three-dimensional surface for each box size, fit a straight line in a logarithmic coordinate system to the logarithm of the box size and the minimum number of boxes, and use the negative value of the slope of the straight line as the fractal dimension at the current moment, calculate the fractal dimension time series using a sliding time window of 72 hours, obtain the fractal dimension change rate by taking the first difference of the fractal dimension time series, determine the fractal control limit at a 99% confidence level based on the distribution of the fractal dimension change rate during historical normal periods, and trigger a primary warning when the fractal dimension change rate exceeds the fractal control limit. The fractal dimension mutation identification algorithm utilizes the inherent fractal geometric properties of the temperature field to establish a geometric criterion for anomaly detection. It is sensitive to anomaly patterns that are invariant in spatial scale. The fractal dimension, as a single scalar index, simplifies feature extraction of multidimensional temperature fields. The difference box counting method requires only simple counting operations to avoid spatial correlation calculations. The computational complexity of the fractal dimension mutation identification algorithm is lower than that of traditional spatial autocorrelation analysis. The fractal dimension mutation detection does not depend on the absolute temperature value but focuses on changes in spatial structure. It is robust to instrument drift and system bias. The sliding time window mechanism enables continuous monitoring and early identification of anomalies. The adaptive setting of the fractal control limit balances detection sensitivity and false alarm rate based on the statistical characteristics of historical data.

[0042] The principle of the inverse derivation algorithm for the partial differential equation inverse problem is as follows: sea surface temperature distribution is controlled by the physical processes described by the heat conduction equation. The temperature anomalies observed at the surface are the result of the combined effects of deep heat sources, boundary heat flux, and three-dimensional heat transport. By inverting the unknown internal heat source distribution and boundary conditions through known surface temperature observations, the inverse problem of the partial differential equation is constructed. The forward model describes the evolution of the temperature field based on the three-dimensional heat conduction equation. The adjoint method constructs the adjoint equation, and its solution provides gradient information of the parameters to be inverted in the objective functional. The quasi-Newton method uses the gradient information to iteratively update the internal heat source distribution and boundary conditions, minimizing the weighted sum of squared residuals between the observed temperature and the simulated temperature. The regularization term constrains the spatial smoothness and temporal continuity of the solution to prevent ill-posedness of the inverse problem from leading to oscillating solutions. The implementation of the inverse derivation algorithm for the partial differential equation after determining the endogenous anomaly in the ocean is as follows: using the surface temperature observation of the region corresponding to the endogenous anomaly as known conditions, a three-dimensional forward heat conduction model is established, where the partial derivative of temperature with respect to time is equal to the thermal diffusivity multiplied by the Laplace operator of temperature plus the internal heat source term. The internal heat source distribution is initialized to zero, and the boundary condition is set as the historical average heat flux. The three-dimensional forward heat conduction model is run to obtain the simulated temperature field. The weighted sum of squares of the residuals between the surface temperature observation and the simulated temperature field is calculated as the objective functional. The adjoint equation is constructed and the adjoint variables are solved by inverse integration. The gradient information of the objective functional with respect to the internal heat source distribution is calculated using the adjoint variables. The BFGS algorithm with quasi-Newton method is used to update the internal heat source distribution and the boundary heat flux. A Tikhonov regularization term is added to constrain the sum of squares of the spatial second derivative of the internal heat source distribution. Iterative optimization is performed until the change of the objective functional is less than the convergence threshold, and the optimal heat source distribution and three-dimensional anomaly structure are output. The inverse partial differential equation algorithm fully integrates the physical laws of heat conduction and observational data, deepening anomaly detection from superficial analysis to the level of physical mechanisms. The internal heat source distribution obtained through inversion reveals the true origin and intensity of the anomaly. The three-dimensional anomaly structure reconstruction provides information on the vertical and horizontal range of the anomaly to support early warning decisions. Regularization constraints ensure the physical rationality and numerical stability of the inversion solution. The adjoint method efficiently calculates gradient information in the high-dimensional parameter space, avoiding the huge computational burden of finite difference. The multi-scale acceleration strategy quickly obtains a preliminary solution on a coarse grid and then performs local optimization on a fine grid. The parallel algorithm decomposes the three-dimensional space into multiple computing nodes for simultaneous solution, significantly shortening the computation time. The inverse partial differential equation algorithm is applicable to seabed topography and non-uniform water structures. The inversion results provide quantitative basis for the physical attribution and impact assessment of anomaly events.

[0043] The intrinsic mode functions (EMFs) are oscillating components with a single instantaneous frequency obtained from empirical mode decomposition (EMD). Each EMF represents the fluctuation characteristics of the original signal over time. The overcomplete dictionary refers to a dictionary matrix with a number of atomic vectors greater than the signal dimension, providing redundant representation bases to achieve sparse signal representation. The orthogonal matching pursuit algorithm is a greedy iterative algorithm. In each iteration, it selects the atomic vector most relevant to the current residual vector and updates the sparse coefficients and residual vector through orthogonal projection. The K-SVD dictionary learning algorithm alternately optimizes sparse coding and dictionary updates. In the dictionary update stage, singular value decomposition is used to update atomic vectors column-by-column to minimize reconstruction error. The differential box counting method covers fractal objects with grids of different scales, counts the number of non-empty grids, and obtains the fractal dimension by fitting a linear relationship between scale and number in logarithmic coordinates. The adjoint method constructs the adjoint equation corresponding to the forward equation. The solution of the adjoint variable provides gradient information of the objective functional with respect to the control variable, avoiding the huge computational cost of directly calculating high-dimensional gradient information. The quasi-Newton method is an iterative optimization algorithm that constructs the search direction by approximating the inverse of the Hessian matrix, avoiding the direct calculation and storage of the Hessian matrix. The Tikhonov regularization adds a penalty term to the norm or derivative norm of the solution in the objective functional, constraining the smoothness and stability of the solution. The random forest classifier consists of multiple decision trees, each trained on a random subset of training samples and a subset of features. The final classification result is determined by voting from all decision trees.

[0044] The atmospheric disturbance characteristic timescale refers to the duration of the influence of typhoon and cold wave atmospheric forcing events on sea surface temperature, determined by statistically analyzing the sea surface temperature response process of historical atmospheric events. The wind field structure refers to the spatial distribution pattern of sea surface temperature under wind stress, manifested as wind-direction-related temperature gradients and frontal characteristics. The grid-by-grid analysis error covariance matrix quantifies the uncertainty and spatial correlation of temperature estimates at each grid point after data assimilation; matrix elements represent the covariance of temperature errors at different locations. The baseline grid spacing is determined based on the computational power and target resolution of the numerical oceanography model, typically set to 10 to 50 kilometers. The gradient threshold... Based on the statistical distribution of the temperature gradient field under historical normal sea conditions, the 90th quantile is used as the fine-grid trigger condition. The variational method finds the initial model field that best matches the observations and is physically reasonable by minimizing the cost function that includes background field errors and observation errors. The fractal control limit is determined based on the mean of the rate of change of fractal dimension during historical normal periods plus three times the standard deviation, corresponding to a 99.7% confidence level. The standard reconstruction error standard deviation, standard observation station density, and standard reconstruction error mean are the long-term averages of the standard deviation, observation station density, and reconstruction error mean of historical global sea area reconstruction errors, respectively. The benchmark threshold is determined based on the statistical distribution of reconstruction errors in the reconstruction error field under historical normal sea conditions, using the 95th quantile. The residual threshold is set to 1% of the norm of the initial residual vector.

[0045] Optionally, the present invention also provides a computer-based early warning system for detecting sea surface temperature anomalies, wherein the computer is equipped with a readable storage medium storing program instructions, and the program instructions execute the aforementioned early warning method for detecting sea surface temperature anomalies when the computer is run.

[0046] The specific implementation methods of the above steps are described in detail below.

[0047] The specific implementation of step S01 is as follows: First, sea surface temperature remote sensing image data is acquired from a satellite remote sensing system. Atmospheric correction and cloud detection are performed using a radiative transfer model to remove cloud-covered pixels and sensor saturation values. Then, in-situ temperature measurement data is collected from an ocean observation buoy network. Data jumps, out-of-limit values, and physically unreasonable values ​​are detected and removed using a quality control algorithm. Next, temperature profile data is acquired from the global Argo buoy system to extract surface temperature information. Finally, the three types of data are uniformly projected onto the same geographic coordinate system. A time interpolation method is used to synchronize data from different sampling times to a unified time reference, forming an initial observation dataset with spatiotemporal registration. The quality control algorithm is based on a dual standard of climatological range detection and gradient detection. Climatological range detection removes data that exceeds the historical temperature extreme range, and gradient detection removes data whose temperature difference with neighboring observation points exceeds 5°C. The time interpolation method uses linear interpolation to maintain the temporal continuity of temperature.

[0048] The specific implementation of step S02 involves performing spatial difference operations on the initial observation dataset to calculate the temperature gradient field. The temperature gradient at each grid point is calculated using a central difference scheme, and a gradient threshold is set. Set to 0.02℃ / km, identify temperature gradients exceeding a gradient threshold in a temperature gradient field. Within the identified area, finer grid points are inserted. The grid spacing of the finer grid is set to 1 / 4 of the base grid spacing. The base grid spacing is typically set to 20 kilometers, and the finer grid spacing is 5 kilometers. This is done when the temperature gradient is less than the gradient threshold. The system maintains a baseline grid in the region, and this adaptive grid division provides a high-resolution description in frontal and vortex regions with drastic temperature changes, while reducing computation in stable regions. The central difference scheme calculates the temperature gradient by dividing the difference in temperature values ​​of adjacent grid points by the grid spacing. The adaptive grid system achieves a balance between computational accuracy and efficiency by dynamically adjusting the spatial resolution.

[0049] The specific implementation of step S03 involves using an ensemble Kalman filter algorithm to integrate the initial observation dataset into the background field of the numerical oceanography model. First, the numerical oceanography model is run to generate a background field with 50 ensemble members, each corresponding to different initial perturbations and parameter perturbations. Then, the mean and covariance of the ensemble background field are calculated. The background field is corrected using the Kalman gain matrix combined with the observation data to obtain the analysis field. The Kalman gain matrix is ​​calculated from the background field error covariance, the observation operator, and the observation error covariance. The analysis field is a physically consistent high-resolution temperature field. Simultaneously, a grid-by-grid analysis error covariance matrix is ​​output to quantify the uncertainty of the analysis field. Next, a variational method is used to construct a cost function, which includes background field error terms and observation error terms. By minimizing the cost function, the initial field of the model is adjusted to ensure that the analysis field is consistent with the observation data and physical laws. The ensemble Kalman filter algorithm avoids the direct calculation of the high-dimensional covariance matrix by estimating the error covariance through ensemble samples. The variational method ensures the dynamic rationality of the assimilation results by introducing physical constraints.

[0050] The specific implementation of step S04 involves adding white noise with an amplitude of 10% of the temperature standard deviation to a physically consistent high-resolution temperature field to generate 100 ensemble samples. Empirical mode decomposition (EMD) is performed on each ensemble sample. EMD iteratively extracts extreme points through a screening process to determine local means and separates single oscillation components from the original signal. This screening process is repeated until all intrinsic mode functions (EMFs) are extracted. The ensemble average of the EMFs of the same order corresponding to the 100 ensemble samples is taken to eliminate the influence of white noise, resulting in an EMF sequence. Then, each EMF in the EMF sequence is processed... The Hilbert transform constructs an analytic signal, where the real part of the analytic signal is the original intrinsic mode function and the imaginary part is its Hilbert transform. The arctangent of the amplitude of the analytic signal is calculated to obtain the instantaneous phase. The instantaneous phase is differentiated with respect to time to obtain the instantaneous frequency. The time evolution of the instantaneous frequency is analyzed to identify the time of anomalous frequency jumps and the corresponding intrinsic mode functions. The anomalous temperature field is reconstructed by linearly superimposing the intrinsic mode functions containing the frequency anomalies. The ensemble empirical mode decomposition improves the stability of the decomposition by adding white noise to solve the mode aliasing problem. The Hilbert-Huang transform captures the local features of non-stationary signals through instantaneous frequency analysis.

[0051] The specific implementation of step S05 involves inputting the spatial distribution data of the abnormal temperature field into the input layer of the temperature field sparse coding recognition model. The input layer normalizes the data to make its norm 1. The sparse coding layer loads a pre-trained overcomplete dictionary containing 256 atomic vectors. The orthogonal matching pursuit algorithm is used to sparsely encode the normalized data. The orthogonal matching pursuit algorithm initializes the residual vector to be equal to the input vector. Iteratively, the atomic vector with the largest inner product with the residual vector is selected and added to the support set. The coefficients of all atomic vectors in the support set are calculated using the least squares method to minimize the residual. The residual vector is then updated to be the input vector minus a linear combination of the support set. The iteration is repeated until the norm of the residual vector is less than 1% of the norm of the input vector or the number of atomic vectors in the support set reaches the sparsity constraint parameter. The sparsity constraint parameter is determined based on the grid density ratio, the trace value of the error covariance matrix of grid-by-grid analysis, and the energy concentration of the anomalous temperature field. The reconstruction layer reconstructs the temperature field using the atomic vectors and corresponding coefficients in the support set. The output layer calculates the root mean square error between the input anomalous temperature field and the reconstructed temperature field as the reconstruction error. When the reconstruction error exceeds the adaptive threshold, the corresponding region is marked as a suspected anomalous region. The adaptive threshold is dynamically adjusted based on the statistical characteristics of historical reconstruction errors and the current observation density. The baseline threshold is set to 0.5℃.

[0052] The specific implementation of step S06 involves extracting temperature field data from suspected anomaly areas to construct a three-dimensional surface. The horizontal coordinates of the three-dimensional surface represent geographical location, while the vertical coordinates represent temperature values. The fractal dimension is calculated using the box-difference counting method. The box size sequence starts at twice the minimum grid spacing and increases by a factor of two to one-quarter of the region scale. For each box size, the three-dimensional surface is covered with a regular grid. The minimum number of boxes containing the surface is counted. Linear fitting is performed in a double logarithmic coordinate system with the logarithm of the box size as the abscissa and the logarithm of the minimum number of boxes as the ordinate. The negative value of the slope of the fitted line is the fractal dimension at the current moment. A 72-hour sliding time window is used to repeatedly calculate the fractal dimension time series for consecutive moments. The change in fractal dimension between adjacent moments is calculated by first-order difference to obtain the fractal dimension change rate. The fractal control limit is determined based on the mean and standard deviation of the fractal dimension change rate during historical normal periods. The fractal control limit is set at a 99.7% confidence level corresponding to the mean plus three times the standard deviation. When the fractal dimension change rate exceeds the fractal control limit, a primary warning is triggered. The differential box counting method extracts the fractal features of the temperature field through multi-scale box coverage. The abrupt change in fractal dimension reflects the essential change in the spatial structure of the temperature field.

[0053] The specific implementation of step S07 involves establishing an atmospheric forcing event feature database by real-time access to typhoon path, intensity, wind field distribution, and cold wave outbreak forecast information from the meteorological forecasting system. This database stores the spatiotemporal characteristics of historical atmospheric events and corresponding sea surface temperature response patterns. A random forest classifier is used to classify temperature anomalies in the areas corresponding to the primary warning. The random forest classifier contains 100 decision trees, each trained based on random sampling of training samples and random feature selection. Classification features include anomaly duration, spatial distribution pattern, the angle between the temperature gradient direction and wind direction, and the temporal evolution of anomaly intensity. The rate and random forest classifier output anomaly type as either atmospheric-driven or ocean-endogenous. When the anomaly duration exceeds the atmospheric disturbance characteristic timescale and the correlation coefficient between the spatial distribution and wind field structure is less than 0.3, it is judged as an ocean-endogenous anomaly. The atmospheric disturbance characteristic timescale is determined based on the statistical analysis of the sea surface temperature response duration of historical typhoons and cold waves. The characteristic timescale of typhoon impact is 3 to 7 days, and the characteristic timescale of cold wave impact is 5 to 10 days. The random forest classifier improves classification accuracy and robustness through ensemble learning, and the anomaly type discrimination avoids false warnings caused by atmospheric disturbances.

[0054] The specific implementation of step S08 involves establishing a three-dimensional forward heat conduction model for the region corresponding to the endogenous anomaly in the ocean. The heat conduction equation describes the change in temperature over time as equal to the thermal diffusivity multiplied by the second spatial derivative of temperature plus an internal heat source term. The internal heat source distribution is initialized as a zero-field boundary condition, set to the historical average heat flux. The forward model is run to solve for the spatiotemporal evolution of the temperature field, resulting in a simulated temperature field. The sum of squared residuals between the observed surface temperature and the simulated surface temperature field is calculated as the target functional. An adjoint equation is constructed to transform the variation of the target functional with respect to the internal heat source distribution into a partial differential equation satisfied by the adjoint variables. The adjoint equation is integrated in reverse from the final time to the initial time to solve for the adjoint variables. The gradient information of the target functional with respect to the internal heat source distribution is calculated using the adjoint variables. The BFGS algorithm in the quasi-Newton method is used to update the internal heat sources based on the gradient information. The BFGS algorithm constructs the search direction by approximating the inverse of the Hessian matrix by storing a finite number of historical gradient vectors. A Tikhonov regularization term is added to the objective functional to constrain the spatial smoothness of the internal heat source distribution. The regularization term is the integral of the square of the second-order spatial derivative of the internal heat source distribution, with a regularization coefficient of 0.01. Iterative optimization is performed until convergence occurs when the relative change of the objective functional is less than 0.001. The algorithm outputs the optimal internal heat source distribution and the three-dimensional anomaly structure. The three-dimensional anomaly structure describes the spatial distribution of the anomalous temperature in the horizontal and vertical directions. The anomaly intensity is defined as the spatial average value of the temperature deviation in the anomalous region. The inverse derivation of the partial differential equation is achieved through physical model constraints to invert deep anomaly sources. The adjoint method efficiently calculates the gradient in the high-dimensional parameter space, avoiding the accumulation of errors in numerical differentiation.

[0055] The specific implementation of step S09 involves establishing an early warning decision function that integrates multiple anomaly indicators to calculate an early warning level index. Input parameters include the reconstruction error in the reconstruction error field, the rate of change of the fractal dimension in the fractal dimension time series, the duration of the anomaly, and the anomaly intensity. The reconstruction error is normalized to the 0-1 interval by dividing it by the historical maximum reconstruction error; the rate of change of the fractal dimension is normalized to the 0-1 interval by dividing it by the fractal control limit; the duration of the anomaly is normalized to the 0-1 interval by dividing it by the maximum monitoring duration of 72 hours; and the anomaly intensity is normalized to the 0-1 interval by dividing it by the historical maximum anomaly intensity. The early warning level index is calculated as the sum of the four normalized indicators. The weights are 0.3, 0.3, 0.2, and 0.2, respectively. The warning response level is determined based on the numerical range of the warning level index. When the warning level index is less than 0.4, it is a normal state and no warning is issued. When the warning level index is between 0.4 and 0.6, it is a blue warning. When the warning level index is between 0.6 and 0.8, it is a yellow warning. When the warning level index is greater than 0.8, it is a red warning. The warning information includes the geographical range, prediction intensity, duration, and three-dimensional structure of the abnormal area. The early warning decision function improves the reliability of the warning through multi-indicator fusion. The hierarchical warning mechanism supports the differentiated deployment of different response measures.

[0056] It should be noted that the key technical ideas of this invention include a fusion mechanism of multi-scale adaptive grid and data assimilation, the synergistic application of sparse representation dictionary learning anomaly encoding algorithm and fractal dimension mutation identification algorithm, and physical constraint anomaly tracing through inverse derivation of partial differential equation inverse problems. The fusion of multi-scale adaptive grid and data assimilation dynamically densifies the grid in regions with severe temperature gradients and uses a coarse grid in stable regions. Combined with ensemble Kalman filtering, sparse observation data is optimally fused with the numerical model background field. This ensures high-resolution characterization of mesoscale phenomena such as fronts and vortices while controlling computational resource consumption. The grid-by-grid analysis error covariance matrix output by data assimilation provides a quantitative uncertainty assessment for subsequent anomaly detection, resolving the contradiction between low interpolation accuracy in sparse observation regions and high computational resource requirements for high-resolution global coverage in traditional fixed grid systems. A multi-level anomaly detection mechanism is established through the collaborative application of a sparse representation dictionary learning anomaly coding algorithm and a fractal dimension mutation identification algorithm. Sparse coding achieves rapid initial screening of anomalies by learning the inherent sparse structure of the normal temperature field. The physical interpretability of dictionary atoms makes anomaly patterns traceable. Fractal dimension mutation identification captures changes in the spatial structure of the temperature field from a geometric topological perspective, making it sensitive to anomaly patterns that are invariant in spatial scale. The combination of the two algorithms cross-validates anomalies from both statistical and geometric dimensions, reducing false alarm rates and improving detection robustness. Inverse derivation of partial differential equations combines surface temperature observations with the physical laws of heat conduction to invert the distribution of deep heat sources and the three-dimensional anomaly structure. This overcomes the limitations of traditional methods that rely solely on surface observations, enabling physical mechanism-level tracing of anomalies. Adjoint methods and regularization constraints ensure computational efficiency and solution stability, providing quantitative physical evidence for the causal analysis and impact range prediction of anomalies. The synergistic effect of three key technological approaches constructs a complete technology chain from data fusion to anomaly detection and physical tracing. Multi-scale grids and data assimilation provide high-quality input data, while sparse coding and fractal dimension detection enable sensitive identification and accurate location of anomalies. Inverse problem deduction deepens the physical understanding of anomalies. The overall technical approach organically integrates data-driven methods with physical model constraints, significantly reducing false alarm rates and shortening early warning response time while maintaining high detection sensitivity. Compared with traditional single threshold detection or pure statistical methods, it demonstrates systematic advantages in dealing with key challenges such as multi-scale ocean processes, background trend separation, and extreme event interference.

[0057] It should be noted that this invention also solves the following technical problem: the difficulty in effectively distinguishing between endogenous marine anomalies and anomalies caused by atmospheric forcing. Traditional methods mainly rely on the spatiotemporal distribution characteristics of temperature fields for anomaly identification, but cannot determine the physical causes of anomalies, leading to misjudging wind-driven surface temperature changes as anomalies caused by internal marine heat sources. This invention establishes an atmospheric forcing event feature database by real-time access to meteorological forecast information, uses a random forest classifier to extract the temporal scale and spatial distribution characteristics of anomalies, compares the duration of anomalies with the temporal scale of atmospheric disturbance characteristics, and matches the spatial distribution pattern with the wind field structure. When the duration of anomalies exceeds the influence cycle of atmospheric disturbances and the spatial distribution does not conform to the wind stress-driven pattern, it is determined to be an endogenous marine anomaly, effectively distinguishing anomaly events with different physical causes, and providing accurate anomaly type identification for subsequent physical mechanism inversion and early warning decision-making.

[0058] Specifically, the principle of this invention is as follows: The core of this invention's solution, which addresses the difficulty of simultaneously ensuring real-time performance and accuracy, lies in the organic combination of a multi-level anomaly discrimination mechanism and physical constraint inversion. The sparse coding recognition model learns a sparse representation dictionary of the normal sea surface temperature field, transforming anomaly detection into low-dimensional manifold deviation judgment. Its computational complexity is linearly related only to sparsity, not the temperature field dimension, thus guaranteeing real-time performance. Fractal dimension mutation recognition utilizes the inherent scale invariance of the temperature field, exhibiting geometric sensitivity to spatial structural anomalies. A single scalar index simplifies feature extraction and avoids high-dimensional correlation calculations. The inverse derivation of partial differential equations uses the physical laws of heat conduction as strong constraints. The accompanying method efficiently calculates gradient information, and the regularization term ensures the physical rationality of the solution, verifying the authenticity of anomalies from a physical mechanism perspective. These three methods are cross-validated from the perspectives of statistical learning, geometric features, and physical laws, forming multiple discrimination barriers to reduce false alarm rates. Simultaneously, an adaptive grid system and data assimilation technology provide high-quality input data, achieving a unified technical architecture for real-time performance and accuracy.

[0059] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0060] The specific implementation of step S01 involves constructing a multi-source sea surface temperature (SST) observation data fusion system. First, three types of data sources are collected: satellite remote sensing data providing broad coverage, buoy-measured data providing high-precision point observations, and Argo profile data providing vertical structure information. The quality control process removes outliers from each data source, with the removal criterion being data points exceeding the historical mean plus or minus four standard deviations. The spatiotemporal reference is uniformly adopted using the World Geodetic Coordinate System (WGS84) and Coordinated Universal Time (UTC) standards. Spatial resolution interpolation is uniformly set to a 0.25-degree grid, and temporal resolution is uniformly set to daily data. The resulting initial observation dataset covers the complete spatiotemporal range of the study area.

[0061] The specific implementation of step S02 involves calculating the temperature gradient field based on the initial observation dataset when establishing a multi-scale adaptive grid system. The formula for calculating the temperature gradient field is as follows:

[0062] ;

[0063] In the formula, For position Temperature gradient at a location, in °C / km; This refers to the temperature field, with units of °C. and The coordinates are horizontal, and the unit is km. and The temperatures are respectively at direction and Partial derivative in direction, in °C / km. Gradient threshold. Based on the statistical distribution of the temperature gradient field under historical normal sea states, the 90th quantile is taken, typically ranging from 0.02℃ / km to 0.05℃ / km. When the temperature gradient is greater than... A fine grid is set for the region, with a grid spacing of 1 / 4 of the baseline grid spacing. The baseline grid spacing is determined based on the computational capabilities of the numerical oceanography model, typically ranging from 10 to 50 km, in areas with temperature gradients less than [value missing]. The region uses a reference grid.

[0064] The specific implementation of step S03 involves integrating the initial observation dataset into the numerical oceanographic model background field when using ensemble Kalman filter data assimilation technology. The calculation formula for the ensemble Kalman filter analysis field is as follows:

[0065] ;

[0066] In the formula, To analyze the temperature vector of the field, the dimension is... The unit is ℃. This represents the total number of grid points. The background field temperature vector has a dimension of The unit is ℃; The Kalman gain matrix has a dimension of . Dimensionless This represents the number of observation points; The observed temperature vector has a dimension of The unit is ℃; The observation operator matrix has a dimension of . Dimensionless, it maps the mode space to the observation space. Kalman gain matrix. The calculation formula is expressed as follows:

[0067] ;

[0068] In the formula, The background error covariance matrix has dimensions of . The unit is ℃ ; Let be the transpose of the observation operator, with dimension . Dimensionless; Let be the observation error covariance matrix, with dimension . The unit is ℃ The initial field of the mode is adjusted using a variational method, and the variational cost function is expressed as follows:

[0069] ;

[0070] In the formula, The cost function is dimensionless. Let be the temperature field vector to be optimized, with dimension . The unit is °C. Minimize the cost function. Generates a physically consistent high-resolution temperature field, and outputs the grid-by-grid analysis error covariance matrix. The calculation formula is expressed as follows:

[0071] ;

[0072] In the formula, It is an identity matrix with dimension 1. Dimensionless; To analyze the error covariance matrix, the dimension is... The unit is ℃ .

[0073] The specific implementation of step S04 involves performing ensemble empirical mode decomposition on a physically consistent high-resolution temperature field. First, white noise following a normal distribution is added to the high-resolution temperature field to form ensemble samples. The standard deviation of the white noise is taken as 0.1 to 0.3 times the standard deviation of the temperature field. Empirical mode decomposition is then performed on each sample in the ensemble samples to obtain the individual sample intrinsic mode functions (IMFs). The IMFs corresponding to all ensemble samples are averaged to eliminate the influence of white noise, resulting in an IMF sequence. ,in Indicates the first One eigenmode function This represents the total number of intrinsic mode functions, typically 5 to 12. Instantaneous frequencies are extracted using the Hilbert-Huang transform; the formula for calculating the instantaneous frequency is as follows:

[0074] ;

[0075] In the formula, For the first The instantaneous frequency of each intrinsic mode function, in Hz; The instantaneous phase is calculated by performing a Hilbert transform on the eigenmode functions to obtain the analytic signal, and the unit is rad. Time is expressed in hours (h). The criterion for identifying non-stationary anomalous components is that the instantaneous frequency standard deviation exceeds twice the instantaneous frequency standard deviation of historical normal periods. The identified non-stationary anomalous components are then used to reconstruct the anomalous temperature field.

[0076] The specific implementation of step S05 is that when the abnormal temperature field is input into the temperature field sparse coding recognition model, the model input layer receives the spatial distribution data of the abnormal temperature field, and the sparse coding layer uses an overcomplete dictionary. Spatial distributed data can be represented sparsely linearly, and the sparse representation formula is expressed as follows:

[0077] ;

[0078] In the formula, Let be the anomalous temperature field vector, with dimension . The unit is ℃. This represents the number of grid points. For an overcomplete dictionary matrix, the dimension is... Dimensionless, containing 256 atomic vectors, each with a dimension of . Consistent with the number of grid points; It is a sparse coefficient vector with dimension . The unit is ℃, which represents the weight coefficient of each atomic vector; To reconstruct the error vector, the dimension is The unit is ℃. The sparse coding layer solves for the sparse coefficient vector using the orthogonal matching pursuit algorithm. The reconstruction layer reconstructs the temperature field using an overcomplete dictionary and sparse coefficients. The output layer calculates the reconstruction error in the reconstruction error field, and the calculation formula is expressed as follows:

[0079] ;

[0080] In the formula, The reconstruction error is dimensionless. To reconstruct the temperature field vector, the dimension is The unit is ℃; Let be the Euclidean norm of the vector. When Exceeding the adaptive threshold The area was marked as a suspected abnormal region.

[0081] The specific implementation of step S06 involves using the difference box counting method to calculate the fractal dimension when identifying fractal dimension mutations in suspected abnormal regions. The formula for calculating the fractal dimension is as follows:

[0082] ;

[0083] In the formula, It is the fractal dimension, which is dimensionless; For size The minimum number of boxes required to cover a three-dimensional surface of a temperature field, dimensionless; Here are the box dimensions, in km. In actual calculations, a series of discrete box dimensions are used. ,in , The box size series is typically 8 to 15, and the slope of a straight line is fitted in a logarithmic coordinate system to obtain the desired value. The fractal dimension time series was calculated using a sliding time window of 72 hours. ,in Indicates the first At any given moment, the formula for calculating the rate of change of the fractal dimension is as follows:

[0084] ;

[0085] In the formula, The rate of change of the fractal dimension is dimensionless. For the previous moment. Fractal control limit Based on the mean of the rate of change of fractal dimension during historical normal periods Add 3 standard deviations Confirmed, the calculation formula is as follows: The empirical value is 0.05 to 0.15, when A primary warning is triggered at this time.

[0086] The specific implementation method of step S07 is the same as described above, and will not be repeated in detail here.

[0087] The specific implementation of step S08 involves establishing a forward model based on the three-dimensional heat conduction equation when initiating the inverse derivation of the partial differential equation for the region corresponding to the endogenous anomaly in the ocean. The heat conduction equation is expressed as follows:

[0088] ;

[0089] In the formula, This is the thermal diffusivity, in units of... Typical value to ; The coordinates are vertical, in meters, with the positive direction upwards. This represents the internal heat source term, with units of ℃ / s. Objective functional. The statement is as follows:

[0090] ;

[0091] In the formula, The objective functional is dimensionless; The number of observation points is dimensionless. For the first The observed temperatures at each observation point are in °C. For the first Simulated temperatures at each observation point, in °C; For reference temperature, 10℃ is used; For the first The standard deviation of the observation error for each observation point, in °C, is typically taken as 0.1 to 0.5 °C. The regularization parameter is dimensionless, and its empirical value is [value missing]. to ; The characteristic length scale is in meters, and it is taken as a typical scale of the study area. For reference heat source intensity, take ℃ / s; For volume infinitesimal elements, the unit is . The BFGS algorithm, based on a quasi-Newton method, iteratively updates the internal heat source distribution. The iterative formula is as follows:

[0092] ;

[0093] In the formula, For the first The internal heat source distribution field of the next iteration is expressed in °C / s. This represents the number of iterations. For the first The step size of each iteration is dimensionless and is determined by line search, typically ranging from 0.01 to 1.0. For the first The approximate Hessian matrix inverse of the nth iteration, in units of ; The gradient field of the objective functional with respect to the heat source distribution is dimensionless and is calculated using the adjoint method.

[0094] The specific implementation of step S09 is to establish the early warning decision function, and the early warning level index. The calculation formula is expressed as follows:

[0095] ;

[0096] In the formula, This is a dimensionless index representing the early warning level. The weighting coefficients are dimensionless and satisfy the following conditions: The values ​​are typically 0.25, 0.25, 0.25, and 0.25. The baseline reconstruction error is dimensionless and is taken as the 95th percentile of the historical normal sea state reconstruction error. The duration of the anomaly is expressed in hours (h). The time scale for atmospheric disturbance characteristics is in hours. The impact of a typhoon is usually 24 to 72 hours, and the impact of a cold wave is usually 48 to 96 hours. The maximum internal heat source intensity obtained by inversion is expressed in °C / s; For reference heat source intensity, the unit is ℃ / s, and it is usually taken as a value. ℃ / s. The warning response level is determined based on the range to which the warning level index belongs. There was no warning at the time, when The alert was in blue at the time. At the time, a yellow alert was issued. At the time, an orange alert was in effect. The alert level was red.

[0097] The specific implementation of the adaptive threshold adjustment function is based on the statistical characteristics of historical reconstruction errors and the current observation conditions, dynamically adjusting the threshold and the adaptive adjustment coefficient. The calculation formula is expressed as follows:

[0098] ;

[0099] In the formula, The adaptive adjustment coefficient is dimensionless. The standard deviation of the historical reconstruction error is dimensionless. The standard deviation of the reconstruction error is dimensionless and is the long-term average of the historical global sea area reconstruction error standard deviation. This represents the current density of observation stations in the sea area, expressed in units per... ; Standard observation station density, in units of [number] stations / The historical average density of global marine observation stations is taken. The mean of the historical reconstruction error is dimensionless. The standard reconstruction error mean is dimensionless and is the long-term average of historical global ocean reconstruction errors. Adaptive threshold. according to Value segmentation setting, when hour ,when hour ,when hour ,when hour ,in The baseline threshold is dimensionless and is determined based on the 95th percentile of the historical normal sea state reconstruction error, with an empirical value of 0.15 to 0.30.

[0100] To better understand and implement this invention, a specific application scenario of the invention is provided below as Example 2: To verify the effectiveness of the invention, technicians built a test environment and selected 15 consecutive days of sea surface temperature (SST) observation data from a certain sea area for anomaly detection and early warning testing. The test area covered a sea area of ​​approximately 500 km × 400 km. A suspected endogenous marine temperature anomaly event occurred in this area from day 9 to day 12 during the test period. Technicians first constructed a multi-source SST observation data fusion system to collect satellite remote sensing data, buoy measurement data, and Argo profile data for the sea area. The satellite remote sensing data came from MODIS and VIIRS sensors, with a temporal resolution of 4 times per day and a spatial resolution of 1 km. The buoy measurement data included real-time observations from 23 drifting buoys and 12 fixed buoys, with a sampling frequency of once per hour. The Argo profile data covered 8 profile buoys, with a vertical resolution ranging from 5 meters to 200 meters in depth. Quality control was performed on the three types of data to remove outliers exceeding the physical range (-2 to 35℃). After unifying the spatiotemporal reference, an initial observation dataset was formed, with a total data volume reaching [amount missing]. One effective observation point was established. Technicians built a multi-scale adaptive grid system to calculate the temperature gradient field of the initial observation dataset, setting the gradient threshold G_0 to 0.08℃ / km. A fine grid was set in the frontal region and vortex edge where the temperature gradient was greater than 0.08℃ / km, with a grid spacing of 5 km, accounting for 28% of the total area. A baseline grid was used in areas with smaller temperature gradients, with a grid spacing of 20 km. The entire region was divided into 625×500 grid points.

[0101] Technicians employed ensemble Kalman filtering data assimilation techniques, setting the ensemble size to 50 members, to integrate the initial observation dataset into the ROMS numerical oceanography model background field. The initial field was adjusted using a variational method, and after 15 iterations of optimization, a physically consistent high-resolution temperature field was generated, simultaneously outputting a grid-by-grid analysis error covariance matrix. The trace value of the grid-by-grid analysis error covariance matrix is... This reflects the overall level of uncertainty after data assimilation. For example... Figure 2 As shown, the high-resolution temperature field clearly reveals the structure of ocean fronts and the characteristics of mesoscale eddies. Technicians performed ensemble empirical mode decomposition on the high-resolution temperature field, adding white noise with an amplitude of 10% of the temperature standard deviation to generate 100 ensemble samples. Empirical mode decomposition was performed on each ensemble sample, yielding 8 intrinsic mode functions (EMFs). The EMFs of all ensemble samples were averaged to eliminate the white noise effect, resulting in an EMF sequence. The instantaneous frequencies of each EMF were extracted using the Hilbert-Huang transform, such as... Figure 3 As shown, the third intrinsic mode function exhibits abnormal instantaneous frequency fluctuations from day 9 to day 12, from the normal... leap to The abnormal temperature field was identified as a non-stationary anomaly. Technicians reconstructed the abnormal temperature field, finding that the maximum anomalous amplitude reached 2.8℃, with a spatial range of approximately 120 km × 95 km.

[0102] Technicians input the abnormal temperature field into a sparse coding model for temperature field identification. This model uses an overcomplete dictionary containing 256 atomic vectors, each with a dimension of 312,500, consistent with the number of grid points. The training dataset uses 10 years of historical normal sea surface temperature observation data, excluding known anomalous event periods, and after normalization, it is divided into 50,000 training samples, each containing the temperature field evolution over 7 consecutive days. The K-SVD dictionary learning algorithm is used to train the overcomplete dictionary, converging after 300 iterations. The sparsity constraint parameters are based on a grid density ratio of 0.28 and the trace value of the error covariance matrix analyzed on a grid-by-grid basis. The energy concentration of the anomalous temperature field was calculated to be 0.76, resulting in a value of 8. The model outputs a reconstruction error field and sparse coding coefficients, with the reconstruction error in the anomalous region reaching 0.58℃. Technicians calculated the adaptive adjustment coefficients, with a historical reconstruction error standard deviation of 0.18℃ and a standard reconstruction error standard deviation of 0.15℃. The current density of observation stations in the sea area is [missing value]. There are 2.3 stations, with a standard observation station density of [number missing] per [location missing]. There are 1.8 sites, with a historical reconstruction error mean of 0.12℃ and a standard reconstruction error mean of 0.10℃. The adaptive adjustment coefficient is calculated to be 1.15, falling within the interval [0.8, 1.2]. The adaptive threshold is set to the baseline threshold of 0.35℃. Since the reconstruction error of 0.58℃ exceeds the adaptive threshold of 0.35℃, it is marked as a suspected abnormal area.

[0103] Technicians performed fractal dimension abrupt change identification on suspected anomaly regions, extracted temperature field data for these regions, and converted it into a three-dimensional surface. Box sizes were set at eight scale levels, increasing from 10 km to 125 km. The minimum number of boxes at each scale was calculated using the differential box counting method, as shown in Table 1. A straight line was fitted in a logarithmic coordinate system, and the negative value of the slope was taken as the fractal dimension. A 72-hour sliding time window was used to calculate the fractal dimension time series, as shown in Table 1. Figure 4As shown, the fractal dimension was 2.38 before the anomaly and jumped to 2.67 during the anomaly. The first difference of the fractal dimension time series yielded the rate of change, with a maximum change of 0.082 / day. Based on historical normal period data, a fractal control limit of 0.045 / day with a 99% confidence level was determined. Since the rate of change of fractal dimension of 0.082 / day exceeded the fractal control limit of 0.045 / day, a primary warning was triggered. Technicians accessed real-time meteorological forecast information and established an atmospheric forcing event feature database, including the spatiotemporal characteristics of typhoons, cold waves, and strong winds. A random forest classifier was used to identify the anomaly type. The classifier consisted of 200 decision trees and was trained based on atmospheric wind field data, sea surface heat flux data, and a historical case database. The classification results showed that the anomaly lasted for 96 hours, exceeding the typical timescale of 48 hours for typhoon impacts, and its spatial distribution did not conform to the elliptical structure driven by the wind field, thus classifying it as an endogenous marine anomaly.

[0104] Table 1. Statistical results of different scale levels using the box-of-differences counting method.

[0105]

[0106] Technicians initiated a reverse derivation of the partial differential equation inverse problem for the region corresponding to endogenous marine anomalies, establishing a three-dimensional forward heat conduction model using surface temperature observations as known conditions. The thermal diffusivity was set as... The internal heat source distribution is initialized to zero, and the boundary condition is set to the historical average heat flux. The simulated temperature field is obtained by running the forward model, and the weighted sum of squared residuals is calculated as the initial value of the objective functional. The adjoint equation was constructed and the adjoint variables were solved by inverse integration. The spatial distribution of the adjoint variables revealed the sensitive region of the temperature field to heat source disturbances. The gradient of the objective functional with respect to the internal heat source distribution was calculated using the adjoint variables, and the heat source distribution and boundary heat flux were updated using the BFGS quasi-Newton method. A Tikhonov regularization term was added with a regularization parameter of 0.001 to constrain the sum of squares of the spatial second derivatives of the heat source distribution. After 42 iterations of optimization, the objective functional converged to 127, with a rate of change less than the convergence threshold of 0.1%. The optimal heat source distribution obtained by inversion showed that the anomalous source was located in the water layer at a depth of 85 to 150 meters, with a central intensity reaching [value missing]. The anomaly spans approximately 80 km × 60 km horizontally. The three-dimensional anomaly structure reconstruction reveals vertical profile characteristics of the temperature anomaly, with the core region showing a temperature 2.8℃ higher than normal, gradually decreasing towards the periphery.

[0107] Technicians established an early warning decision function, comprehensively considering factors such as a reconstruction error of 0.58℃, a fractal dimension change rate of 0.082 / day, anomaly duration of 96 hours, and anomaly intensity. The warning level index was calculated. The weights for reconstruction error, fractal dimension change rate, anomaly duration, and anomaly intensity were set to 0.25 and 0.20 respectively. The calculated warning level index was 7.8, falling within the range [7.0, 9.0), corresponding to a Level II warning response. Technical personnel issued a warning, advising relevant departments to monitor the development of the abnormal temperature in the sea area, increase the frequency of observations, and assess the potential impact on marine ecology and fishery resources. Subsequent monitoring showed that the abnormal event peaked on day 13 and gradually weakened, essentially dissipating on day 16, consistent with the warning information and verifying the effectiveness of the invention.

[0108] It should be noted that the variables involved in this invention are explained in detail in Tables 2 and 3.

[0109] Table 2. Variable Explanation Table (Part 1)

[0110]

[0111] Table 3. Variable Explanation Table (Part Two)

[0112]

[0113] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. An early warning method for detecting sea surface temperature anomalies, characterized in that, Includes the following steps: A multi-source sea surface temperature (SST) observation data fusion system was constructed, collecting satellite remote sensing data, buoy-measured data, and Argo profiling data. Outliers were removed through quality control, and an initial observation dataset was formed after unifying the spatiotemporal reference. A multi-scale adaptive grid system was established, calculating the temperature gradient field based on the initial observation dataset. A fine grid was set in regions where the temperature gradient was greater than a gradient threshold, while a baseline grid was used in regions where the temperature gradient was less than the gradient threshold. An ensemble Kalman filter data assimilation technique was employed to integrate the initial observation dataset into the background field of the numerical oceanography model. A variational method was used to adjust the model's initial field to generate a physically consistent high-resolution temperature field. Simultaneously, it outputs the grid-by-grid analysis error covariance matrix; it performs ensemble empirical mode decomposition on the physically consistent high-resolution temperature field, decomposing it into an intrinsic mode function sequence, extracting instantaneous frequencies through Hilbert-Huang transform, identifying non-stationary anomalous components, and reconstructing the anomalous temperature field; it inputs the anomalous temperature field into the temperature field sparse coding identification model, outputting the reconstruction error field and sparse coding coefficients, and marking it as a suspected anomalous region when the reconstruction error exceeds the adaptive threshold; it performs fractal dimension mutation identification on the suspected anomalous region, using the difference box counting method to calculate the fractal dimension time series, and triggering a primary warning when the rate of change of the fractal dimension exceeds the fractal control limit; Real-time access to meteorological forecast information is used to establish an atmospheric forcing event feature database. Random forest classifiers are used to identify the anomaly type in the area corresponding to the primary warning. When the duration of the anomaly exceeds the characteristic time scale of atmospheric disturbance and the spatial distribution does not conform to the wind field structure, it is determined to be an endogenous marine anomaly. For regions corresponding to endogenous marine anomalies, inverse partial differential equations are used to infer the distribution of internal heat sources based on the heat conduction equation. The adjoint method and quasi-Newton method are used for iterative optimization to reconstruct the three-dimensional anomaly structure and estimate the anomaly intensity. An early warning decision function is established, and the early warning level index is calculated by comprehensively considering the reconstruction error, the rate of change of fractal dimension, the duration of the anomaly, and the anomaly intensity. The early warning response level is determined according to the interval to which the early warning level index belongs, and the early warning information is issued.

2. The method according to claim 1, characterized in that, The steps of the ensemble empirical mode decomposition are as follows: a normally distributed white noise is added to a physically consistent high-resolution temperature field to form an ensemble sample; empirical mode decomposition is performed on each ensemble sample to obtain the single-sample intrinsic mode function; the single-sample intrinsic mode functions corresponding to all ensemble samples are averaged to eliminate the influence of white noise, resulting in an intrinsic mode function sequence.

3. The method according to claim 2, characterized in that, The Hilbert-Huang transform obtains an analytic signal by performing a Hilbert transform on the eigenmode functions in the eigenmode function sequence. The instantaneous frequency is obtained by calculating the derivative of the instantaneous phase of the analytic signal with respect to time. Abnormal fluctuations in the instantaneous frequency indicate the occurrence of non-stationary temperature anomalies.

4. The method according to claim 3, characterized in that, The structure of the temperature field sparse coding recognition model is as follows: the input layer receives the spatial distribution data of the abnormal temperature field, the sparse coding layer uses an overcomplete dictionary to perform sparse linear representation of the spatial distribution data, the sparse coding layer solves the sparse coefficient vector through the orthogonal matching pursuit algorithm, the reconstruction layer reconstructs the temperature field using the atomic vectors and sparse coefficient vectors in the overcomplete dictionary, and the output layer calculates the root mean square error between the abnormal temperature field and the temperature field reconstructed by the reconstruction layer as the reconstruction error.

5. The method according to claim 4, characterized in that, The overcomplete dictionary contains 256 atomic vectors, each with a dimension consistent with the number of grid points in the anomalous temperature field.

6. The method according to claim 5, characterized in that, The steps for establishing the training dataset for the sparse coding recognition model of the temperature field are as follows: collect normal sea surface temperature observation data from the past 10 years, remove data from the periods when known abnormal events occurred in the normal sea surface temperature observation data, normalize the normal sea surface temperature observation data after removal so that the temperature values ​​are distributed between 0 and 1, divide the normalized temperature field into training samples according to spatial regions, and each training sample contains the temperature field evolution process for 7 consecutive days, forming a total of 50,000 training samples to constitute the training dataset.

7. The method according to claim 6, characterized in that, The training steps of the temperature field sparse coding recognition model are as follows: an overcomplete dictionary is initialized using the K-SVD dictionary learning algorithm, and the sparse coding stage and dictionary update stage are executed iteratively. In the sparse coding stage, the overcomplete dictionary is fixed and the orthogonal matching pursuit algorithm is used to solve the sparse coefficients of all training samples in the training dataset. In the dictionary update stage, the sparse coefficients are fixed and the atomic vectors in the overcomplete dictionary are updated one by one using singular value decomposition to minimize the reconstruction error. After 300 iterations, the overcomplete dictionary converges.

8. The method according to claim 7, characterized in that, The sparsity constraint parameters in the sparse coding layer are determined based on the grid density ratio of the multi-scale adaptive grid system, the trace value of the grid-by-grid analysis error covariance matrix, and the energy concentration of the anomalous temperature field. The sparsity constraint parameters are positively correlated with the grid density ratio, negatively correlated with the trace value of the grid-by-grid analysis error covariance matrix, and positively correlated with the energy concentration of the anomalous temperature field.

9. The method according to claim 8, characterized in that, The adaptive threshold is calculated based on the mean of historical reconstruction errors, the standard deviation of historical reconstruction errors, and the current density of observation stations in the sea area. The adaptive adjustment coefficient is equal to the ratio of the standard deviation of historical reconstruction errors to the standard deviation of standard reconstruction errors, multiplied by the square root of the ratio of the current density of observation stations in the sea area to the standard density of observation stations, and then multiplied by the ratio of the mean of historical reconstruction errors to the mean of standard reconstruction errors. The product of these three terms is then cubed.

10. The method according to claim 9, characterized in that, The fractal dimension mutation identification algorithm is implemented after the temperature field sparse coding identification model identifies a suspected abnormal region. Specifically, it extracts the temperature field data of the suspected abnormal region, converts the temperature value in the temperature field data into a three-dimensional surface where the third dimension is the temperature value, sets a series of box sizes that increase from twice the minimum grid spacing to one-quarter of the region scale, calculates the minimum number of boxes required to cover the entire three-dimensional surface for each box size, fits a straight line in a logarithmic coordinate system to the logarithm of the box size and the logarithm of the minimum number of boxes, and uses the negative value of the slope of the straight line as the fractal dimension at the current moment.

Images