Wave parameter inversion method under strong tidal environment based on pressure sensor array

Abstract

The application belongs to the technical field of ocean wave observation, and discloses a wave parameter inversion method under strong tidal environment based on a pressure sensor array. The method comprises the following steps: synchronously collecting pressure time sequences and average water depths of the pressure sensor array; for each sensor, a dynamic localized tidal model is constructed and a tidal prediction sequence is generated; a common mode of tidal prediction sequences of all sensors of the array is calculated and deducted from original data, and a singular value decomposition is performed on the deducted data matrix to obtain pure wave pressure signals; based on an array cross-spectrum matrix and by fusing a tidal current main flow direction extracted from the local tidal model as a priori constraint, an iterative maximum likelihood method with a penalty function is used to invert a high-resolution wave direction spectrum. The application can effectively separate wave signals under the conditions of short observation sequences, strong tidal interference and data discontinuity, completely retain long-period wave information, eliminate tidal residual inconsistency between sensors, and significantly improve wave direction inversion accuracy and stability.

Description

Technical Field

[0001] This invention belongs to the field of ocean wave observation technology, and specifically relates to a wave parameter inversion method based on a pressure sensor array under strong tidal conditions. Background Technology

[0002] Waves are a core element of the marine dynamic environment, and high-precision wave parameters (especially wave direction spectra) are crucial for marine engineering design, shipping safety, coastal zone management, and military activities. Currently, wave observation mainly relies on equipment such as wave radar, acoustic Doppler current profilers (ADCP), or wave buoys. While these devices offer high accuracy, they generally suffer from high costs, complex deployment and maintenance, and vulnerability to damage.

[0003] Underwater pressure sensors have become a promising alternative or supplementary observation method due to their advantages such as low cost, high reliability, and ease of concealed deployment. However, the signals recorded by pressure sensors are the superposition of dynamic wave pressure and hydrostatic pressure changes caused by low-frequency processes such as tides. In nearshore or shallow water areas, the tidal amplitude is much greater than the pressure changes caused by waves, resulting in significant background noise.

[0004] To eliminate tidal interference, existing technologies typically employ high-pass filtering, which involves high-pass filtering the pressure time series to remove low-frequency tidal components. However, this method simultaneously filters out long-period waves containing important oceanographic information, leading to a severe underestimation of low-frequency wave energy. Furthermore, it cannot handle the spectral aliasing problem caused by the nonlinear interaction between tides and waves. Additionally, simple time-domain filtering disrupts the original phase relationship between sensors, introducing fundamental errors for subsequent wave direction inversion based on phase differences.

[0005] Another commonly used method is the standard harmonic analysis method, which performs tidal harmonic analysis based on long-term series (usually requiring at least one month) and subtracts tidal components from the forecast model. However, in practical applications, it often faces problems such as short observation sequences, data discontinuities, or contamination by severe sea conditions, leading to unstable harmonic analysis results and decreased forecast accuracy. More importantly, this method is essentially a single-point processing method. When applied to pressure sensor arrays, it cannot handle the inconsistencies in the residual phase and amplitude of tidal signals caused by slight differences in installation height, local topography, or instrument errors among the sensors. These residual inconsistencies severely contaminate the cross-spectral phase calculation between sensors, thus introducing significant errors into wave direction inversion.

[0006] In summary, existing wave spectrum inversion techniques largely rely on ADCP or array radar. Methods for high-precision wave spectrum inversion from low-cost pressure sensor arrays, especially under conditions of strong tidal interference, remain lacking. Therefore, there is an urgent need to develop a method that can fuse prior information to stably construct a local tidal model and eliminate residual tidal inconsistencies between sensors through array collaborative processing. This would enable high-fidelity extraction of wave information from pressure data with strong tidal pollution, achieving high-resolution two-dimensional wave spectrum inversion. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention provides a wave parameter inversion method based on a pressure sensor array under strong tidal conditions. This method aims to extract wave information from pressure sensor array data with high fidelity under conditions of strong tidal interference, short observation sequences, and imperfect arrays, thereby achieving high-precision inversion of wave height, wave period, and high-resolution wave spectrum.

[0008] To achieve the above objectives, the technical solution of the present invention is as follows:

[0009] A wave parameter inversion method based on a pressure sensor array under strong tidal conditions includes the following steps:

[0010] S1: Synchronously collect pressure data from each sensor in the pressure sensor array deployed in the target sea area, and obtain the average water depth;

[0011] S2: For the pressure data of each sensor, combined with the regional tidal prior information, a dynamic localized tidal model for each sensor location is constructed through harmonic analysis with regularization constraints and robust estimation, and tidal forecast sequences are generated respectively.

[0012] S3: Calculate the common pattern of the tide forecast sequences of all sensors and subtract it from the raw pressure data of each sensor; perform singular value decomposition on the pressure data array matrix after subtracting the common pattern, extract and eliminate the first principal component, and obtain the pure wave pressure signal of each sensor.

[0013] S4: Based on the pure wave pressure signal, the wave height and wave period are inverted using the pressure-wave surface transfer function; based on the pure wave pressure signal, the cross spectrum matrix is ​​calculated and the wave direction spectrum is inverted using the direction spectrum estimation algorithm to obtain the wave direction information.

[0014] In the above scheme, the pressure sensor array consists of a triangular array of three underwater pressure sensors.

[0015] In the above scheme, S2 specifically includes:

[0016] S2.1: Construct observation equations to represent the tidal components observed by each sensor as a linear superposition of multiple tidal constituents, and establish the relationship between the harmonic constant to be determined and the observed values;

[0017] S2.2: Introduce prior constraints, obtain prior harmonic constants from the regional tidal model, and construct an objective function containing regularization terms to constrain the parameters to be determined from deviating too much from the prior values ​​under the condition of short observation sequences;

[0018] S2.3: Robust iterative solution: Based on the objective function, iterative reweighted least squares method is used to solve the problem. During the iteration process, the weights of the observed data are dynamically adjusted according to the residual size to suppress the influence of outliers. After convergence, a stable local harmonic constant is obtained.

[0019] S2.4: Generate local tide forecasts by using the stable local harmonic constant to generate a tide forecast sequence for each sensor at any time.

[0020] In a further technical solution, in S2.1, for the pressure data of the i-th sensor... Its tidal components are represented as the superposition of multiple tidal constituents:

[0021] ;

[0022] Where i is the sensor number. Let j be the tide observed by the i-th sensor, and j be the tide fraction number. M is the selected number of tidal constituents. For the nodal factor of the j-th tidal constituent, Let be the astronomical initial phase angle of the j-th tidal constituent. The correction angle for the j-th tidal constituent node. , Let be the amplitude and lag angle of the i-th sensor at the j-th tidal fraction, respectively, i.e., the harmonic constant. This is the average sea level.

[0023] In a further technical solution, S2.2, the objective function containing the regularization term is constructed as follows:

[0024] ;

[0025] in, The design matrix is ​​K×(2M+1), and its k-th row is:

[0026] ;

[0027] in, For the Mth tidal component The node factor at time, For the Mth tidal component Phase angle at time, , Let ω be the angular velocity of the Mth tidal constituent. For the Mth tidal component The node factor at time, For the Mth tidal component The nodal correction angle at time , It is the mean sea level.

[0028] There are 2M+1 parameter vectors to be determined. ,in, Let be the cosine component coefficient of the i-th sensor at the j-th tidal fraction. Let be the sinusoidal component coefficient of the i-th sensor at the j-th tidal fraction. , , Let be the amplitude and lag angle of the i-th sensor at the j-th tidal fraction, respectively. It is the mean sea level.

[0029] Let K×1 be the observation vector, consisting of tidal observations from K time sampling points. Composition, where K is the time sampling point;

[0030] This is a regularization parameter used to balance the strength of data fitting and prior constraints;

[0031] The weight matrix is ​​(2M+1)×(2M+1), and the diagonal matrix is ​​taken, with the diagonal elements being the reciprocal of the confidence level of each prior parameter;

[0032] These are prior parameters;

[0033] ;

[0034] in, Let be the prior cosine component coefficients of the i-th sensor at the j-th tidal fraction. This represents the prior sinusoidal component coefficient of the i-th sensor at the j-th tidal constituent. , It is the amplitude of the i-th sensor at the j-th tidal constituent, obtained from the regional tidal model. It is the lag angle of the i-th sensor at the j-th tidal constituent, obtained from the regional tidal model.

[0035] In a further technical solution, the specific steps of S2.3 are as follows:

[0036] (1) Initialize the parameter vector Initialize the observation weight matrix , Let K be a K×K identity matrix, and let the iteration counter p=0;

[0037] (2) Iterative calculation: in the first... In the next iteration, the current weight matrix is ​​used. Solve for parameters :

[0038] ;

[0039] in, For the i-th sensor at the th The parameter vector after the nth iteration. For the first The observation weight matrix after the next iteration;

[0040] (3) Calculate the observation residuals under the current iteration:

[0041] ;

[0042] in, Indicates the first Observation residuals under the next iteration;

[0043] And the robust standard deviation is estimated using the median absolute deviation:

[0044] ;

[0045] in, This represents the median;

[0046] (4) Update the observation weights using the Huber weight function:

[0047] ;

[0048] in, Represents a diagonal matrix; express The observation weight matrix after the next iteration has diagonal elements representing the weights of each observation in the current iteration. express After the nth iteration Weights of each observation:

[0049] ;

[0050] in, This represents the residual of the k-th observation. The tuning constant is... For robust standard deviation;

[0051] (5) Convergence judgment:

[0052] ;

[0053] when If the threshold value is less than ε, the iteration terminates; otherwise, the iteration continues. ε is the threshold value.

[0054] (6) After the iteration converges, stable linear parameter estimates are obtained. Then, the estimated harmonic constants are obtained through inverse transformation:

[0055] ;

[0056] ;

[0057] in, This represents the amplitude estimate of the j-th tidal fraction from the i-th sensor; This represents the estimated lag angle of the j-th tidal constituent from the i-th sensor; This represents the estimated cosine component of the j-th tidal constituent from the i-th sensor. This represents the estimated sinusoidal component of the j-th tidal constituent from the i-th sensor. This is an overestimate of the mean sea level.

[0058] In a further technical solution, in S2.4, the harmonic constant obtained in S2.3 is used. and mean sea level Generate the tidal forecast value for the i-th sensor at time t:

[0059] ;

[0060] in, Let be the tidal forecast value of the i-th sensor at time t. For the j-th tidal component The node factor at time, Let be the astronomical initial phase angle of the j-th tidal constituent. The correction angle for the j-th tidal constituent.

[0061] In the above scheme, S3 specifically includes:

[0062] S3.1: Common tidal pattern extraction and preliminary subtraction: that is, the arithmetic mean of the tidal forecast sequences of all sensors is calculated as the common tidal pattern, and the common pattern is subtracted from the original pressure data of each sensor to obtain the signal after preliminary de-tidal processing.

[0063] S3.2: Cooperative extraction and elimination of residual common low-frequency noise: The signals from each sensor after initial de-tidal processing are organized into a data matrix. Singular value decomposition is performed on the matrix to extract the first principal component with the highest energy and remove it from the matrix to obtain the pure wave pressure signal from each sensor, which is then used for subsequent wave parameter inversion.

[0064] In a further technical solution, S3.2, an observation matrix is ​​first constructed to initially remove the tidal signals. A data matrix organized as N×K Where N is the number of sensors and K is the time sampling point; matrix elements This represents the initial detidal signal of the i-th sensor at the k-th sampling time;

[0065] Again Perform singular value decomposition:

[0066] ;

[0067] in, It is an N×N left singular vector matrix, where each column represents a spatial mode; It is an N×K diagonal matrix; It is a K×K right singular vector matrix, where each column represents the time coefficient of the corresponding spatial mode;

[0068] Subtracting the reconstructed component of the first principal component from the original data matrix yields the data matrix after eliminating common noise:

[0069] ;

[0070] in, The maximum singular value obtained from singular value decomposition. The first column of the left singular vector matrix U, This is the transpose of the first row of the right singular vector matrix V;

[0071] The i-th row vector is the pure wave pressure signal of the i-th sensor after array collaborative noise reduction.

[0072] In the above scheme, the directional spectrum estimation algorithm described in S4 is an iterative maximum likelihood method, and the local power flow information obtained from S2 is incorporated as a priori constraint during the iteration process, as detailed below:

[0073] Extracting the main current direction from the local tidal model Construct a penalty function :

[0074] ;

[0075] in, Indicates wave direction;

[0076] The penalty function is introduced into the objective function of the iterative maximum likelihood method to form a constrained directional spectrum estimation criterion:

[0077] ;

[0078] in, The final frequency-direction spectrum; This is a mathematical operator that represents finding the minimum value of the objective function within the parentheses. ML(E) is the original maximum likelihood objective function, and α is the regularization parameter; The frequency-direction spectrum is continuously updated during the iteration process;

[0079] During each iteration of updating the directional spectrum, a penalty is imposed on wave direction solutions that deviate significantly from the mainstream trend of the current, ultimately yielding a directional spectrum that incorporates prior information about the current. .

[0080] Through the above technical solution, the wave parameter inversion method based on a pressure sensor array under strong tidal conditions provided by the present invention has the following beneficial effects:

[0081] 1. High-precision tidal modeling with short-sequence data:

[0082] This invention introduces prior harmonic constants from regional tidal models or historical data to construct an objective function with regularization constraints. Combined with robust iterative solution, it can stably estimate the local harmonic constants of each sensor location under conditions of short observation sequences, data discontinuity, or pollution from severe sea conditions, thereby generating high-precision tidal forecast sequences and overcoming the dependence of traditional harmonic analysis on long-sequence data.

[0083] 2. Fully preserve long-cycle wave information:

[0084] Unlike high-pass filtering, which simply truncates low-frequency signals, this invention employs a two-step array collaborative processing method: "common tidal mode subtraction and principal component analysis / singular value decomposition principal component elimination." This method eliminates tidal and common low-frequency noise while fully preserving long-period wave components such as swells, avoiding the loss of low-frequency wave energy and providing more realistic wave signals for subsequent wave height and direction inversion.

[0085] 3. Precisely eliminate residual tidal inconsistencies between sensors:

[0086] To address the inconsistency in tidal residual phase and amplitude between array sensors caused by differences in installation height, local terrain changes, or instrument errors, which cannot be handled by standard harmonic analysis, this invention effectively eliminates common low-frequency noise between channels by extracting and subtracting common tidal modes and the first principal component. This provides high-fidelity phase information for wave direction inversion and significantly improves the inversion accuracy of wave direction spectrum.

[0087] 4. Significantly improves the accuracy of wave direction inversion in strong current sea areas:

[0088] This invention extracts the tidal current elliptical elements (such as the main current direction) of the major tidal constituents from the local tidal model and introduces them as physical prior constraints into the objective function of directional spectrum inversion. By penalizing wave direction solutions that are seriously deviated from the known tidal current direction, it effectively suppresses false wave direction interference under strong current background and improves wave direction resolution and inversion stability.

[0089] 5. Low cost and easy deployment:

[0090] This invention is based on an underwater pressure sensor array. Compared with traditional equipment such as wave radar, ADCP or wave buoys, it has the advantages of low cost, high reliability, easy concealment and deployment and simple maintenance, providing an economical and effective technical solution for wave observation in nearshore and shallow water areas. Attached Figure Description

[0091] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below.

[0092] Figure 1 This is an overall flowchart of the method of the present invention;

[0093] Figure 2 This is a schematic diagram of the layout of a pressure sensor array (triangle). Detailed Implementation

[0094] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.

[0095] This invention provides a wave parameter inversion method based on a pressure sensor array under strong tidal conditions, such as... Figure 1 As shown, it includes the following steps:

[0096] Step S1: Synchronous acquisition and preprocessing of array data

[0097] In this embodiment, a triangular spatial array consisting of three high-precision underwater pressure sensors is deployed on the seabed of the target sea area (e.g., a strong tidal sea area with a near-shore water depth of approximately 20m). Figure 2As shown, to ensure uniformity of wave resolution in all directions, the three sensors (S1, S2, S3) are arranged in an equilateral triangle with a baseline length of 10 meters. All sensors are strictly synchronized via a synchronization cable or a built-in high-precision clock, and the sampling frequency is uniformly set to 2Hz.

[0098] Pressure time series P1(t), P2(t), and P3(t) from three sensors were continuously collected over 20 days (T≥15 days). Simultaneously, the average water depth d (e.g., 19.8 meters) of the array coverage area was obtained through high-precision depth measurement or from regional nautical charts for subsequent calculation of the pressure-wave surface transfer function.

[0099] Step S2: Dynamic localized tidal modeling based on prior constraints and robust estimation

[0100] For each sensor's pressure data, the following sub-steps are performed to build a high-precision local tidal model.

[0101] S2.1: Constructing the observation equation

[0102] For the pressure data of the i-th sensor The tidal components it contains are represented as the superposition of M tidal constituents, and the observation equation is as follows:

[0103] ;

[0104] Where i is the sensor number. Let j be the tide observed by the i-th sensor, and j be the tide fraction number. M is the selected number of tidal constituents. For the nodal factor of the j-th tidal constituent, Let be the astronomical initial phase angle of the j-th tidal constituent. The correction angle for the j-th tidal constituent node. , Let be the amplitude and lag angle of the i-th sensor at the j-th tidal fraction, respectively, i.e., the harmonic constant. This is the average sea level.

[0105] Taking sensor 1 as an example, its pressure time series Tidal components It can be represented as The linear superposition of the major tidal constituents. This embodiment selects... =4 main tidal constituents, including the semi-diurnal tidal constituents , and diurnal tides 1. O1.

[0106] S2.2: Introducing Prior Constraints

[0107] Extract the data from the high-precision regional tidal model (such as TPXO8 or NAO.99b) covering the target sea area, focusing on the same batch of data at sensor location i. The a prior harmonic constant of each tidal component , And convert it to a linear parameter:

[0108] ;

[0109] ;

[0110] Construct the objective function that includes the regularization term:

[0111] ;

[0112] in, The design matrix is ​​K×(2M+1), and its k-th row is:

[0113] ;

[0114] in, For the Mth tidal component The node factor at time, For the Mth tidal component Phase angle at time, , Let ω be the angular velocity of the Mth tidal constituent. For the Mth tidal component The node factor at time, For the Mth tidal component The nodal correction angle at time , It is the mean sea level.

[0115] There are 2M+1 parameter vectors to be determined. ,in, Let be the cosine component coefficient of the i-th sensor at the j-th tidal fraction. Let be the sinusoidal component coefficient of the i-th sensor at the j-th tidal fraction. , , Let be the amplitude and lag angle of the i-th sensor at the j-th tidal fraction, respectively. It is the mean sea level.

[0116] Let K×1 be the observation vector, consisting of tidal observations from K time sampling points. Composition, where K is the time sampling point; + , For observation error, the first term is to make the prediction... Closer observation The sum of squared residuals;

[0117] This is a regularization parameter used to balance the strength of data fitting and prior constraints, and can be set empirically; in this embodiment, it is set to 0.1.

[0118] The weight matrix is ​​(2M+1)×(2M+1), and the diagonal matrix is ​​taken, with the diagonal elements being the reciprocal of the confidence level of each prior parameter;

[0119] These are prior parameters.

[0120] ;

[0121] in, Let be the prior cosine component coefficients of the i-th sensor at the j-th tidal fraction. This represents the prior sinusoidal component coefficient of the i-th sensor at the j-th tidal constituent. , It is the amplitude of the i-th sensor at the j-th tidal constituent, obtained from the regional tidal model. It is the lag angle of the i-th sensor at the j-th tidal constituent, obtained from the regional tidal model.

[0122] S2.3: Robust Iterative Solution

[0123] The objective function is solved using an iterative reweighted least squares method to suppress outliers in the observed data. During the iteration process, the weights of the observed data are dynamically adjusted according to the residuals to suppress the influence of outliers, and a stable local harmonic constant is obtained after convergence.

[0124] The specific steps are as follows:

[0125] (1) Initialize the parameter vector Initialize the observation weight matrix , Let K be a K×K identity matrix, and let the iteration counter p=0;

[0126] (2) Iterative calculation: in the first... In the next iteration, the current weight matrix is ​​used. Solve for parameters :

[0127] ;

[0128] in, For the i-th sensor at the th The parameter vector after the nth iteration. For the first The observation weight matrix after the next iteration;

[0129] (3) Calculate the observation residuals under the current iteration:

[0130] ;

[0131] in, Indicates the first Observation residuals under the next iteration;

[0132] And the robust standard deviation is estimated using the median absolute deviation:

[0133] ;

[0134] in, This represents the median;

[0135] (4) Update the observation weights using the Huber weight function:

[0136] ;

[0137] in, Represents a diagonal matrix; express The observation weight matrix after the next iteration has diagonal elements representing the weights of each observation in the current iteration. express After the nth iteration Weights of each observation:

[0138] ;

[0139] in, This represents the residual of the k-th observation. The tuning constant is... For robust standard deviation;

[0140] (5) Convergence judgment:

[0141] ;

[0142] when If the threshold value is less than ε, the iteration terminates; otherwise, the iteration continues. ε is the threshold value.

[0143] (6) After the iteration converges, stable linear parameter estimates are obtained. Then, the estimated harmonic constants are obtained through inverse transformation:

[0144] ;

[0145] ;

[0146] in, This represents the amplitude estimate of the j-th tidal fraction from the i-th sensor; This represents the estimated lag angle of the j-th tidal constituent from the i-th sensor; This represents the estimated cosine component of the j-th tidal constituent from the i-th sensor. This represents the estimated sinusoidal component of the j-th tidal constituent from the i-th sensor. This is an overestimate of the mean sea level.

[0147] Repeat steps S2.1 to S2.3 above for the three sensors to obtain their respective local harmonic constants.

[0148] S2.4: Generate local tide forecast

[0149] The harmonic constant obtained using S2.3 and mean sea level Generate the tidal forecast value for the i-th sensor at time t:

[0150] ;

[0151] in, Let be the tidal forecast value of the i-th sensor at time t. For the j-th tidal component The node factor at time, Let be the astronomical initial phase angle of the j-th tidal constituent. The correction angle for the j-th tidal constituent.

[0152] The same method was used to obtain local tidal forecasts from the three sensors.

[0153] Step S3: Array-coordinated tidal noise reduction for wave direction inversion

[0154] S3.1: Extraction and Preliminary Deduction of Common Tidal Patterns

[0155] Calculate the arithmetic mean of the tidal forecast sequences from the three sensors as a common tidal pattern:

[0156] ;

[0157] Then from the raw pressure of each sensor Subtracting the common tidal pattern from the signal yields the signal after preliminary tidal removal:

[0158] ;

[0159] This operation eliminates the main tidal component shared by all sensors in the array, but it still has drawbacks due to local water depth differences, installation height deviations, different bottom reflection characteristics, and harmonic constant estimation errors among the sensors. There may still be residual low-frequency tidal margins or common low-frequency noise from different sensors.

[0160] S3.2: Cooperative Extraction and Elimination of Residual Common Low-Frequency Noise

[0161] To eliminate the aforementioned residual common low-frequency noise, this step employs singular value decomposition (SVD) and adaptive filtering utilizing the spatial correlation of the array's multi-channel array. The preliminary signals from the three sensors are then processed. , , (Each signal has K sampling points) organized into a 3×K data matrix Matrix elements This represents the initial detidal signal of the i-th sensor at the k-th sampling time;

[0162] Perform singular value decomposition (SVD) on it:

[0163] ;

[0164] in, It is an N×N left singular vector matrix, where each column represents a spatial mode; It is an N×K diagonal matrix, with diagonal elements being singular values, arranged in descending order of energy; It is a K×K right singular vector matrix, where each column represents the time coefficient of the corresponding spatial mode.

[0165] Because common noise such as tides and low-frequency residual water levels is highly correlated across all sensors, and its energy is typically much higher than that of wave signals, its corresponding modal behavior exhibits the largest singularity. and its corresponding spatial vector (First column) and time vector (First column). This mode reflects the prevalent in-phase low-frequency variation in the array. Removing this principal component from the data matrix eliminates residual common noise, yielding the noise-free data matrix:

[0166] ;

[0167] in, The maximum singular value obtained from singular value decomposition. The first column of the left singular vector matrix U, This is the transpose of the first row of the right singular vector matrix V; The i-th row vector is the pure wave pressure signal from the i-th sensor after array-based collaborative noise reduction (i.e., the signal required by subsequent S4). This signal has largely eliminated tidal interference and residual inconsistencies between sensors, and can be used for subsequent high-precision wave inversion.

[0168] This operation is equivalent to performing an adaptive high-pass filter in the spatial dimension, which removes the low-frequency components with the strongest energy that are common to all sensors, thereby preserving the unique signals of each sensor (mainly dynamic pressure changes caused by waves and random noise from the sensors themselves).

[0169] Step S4: High-precision inversion of wave parameters

[0170] S4.1: Wave Height and Period Inversion

[0171] Dynamic pressure time series acquired by any pressure sensor A pressure response function is established using linear wave theory, and then transformed into a time-wave equation for water surface ripples. Then, wave parameters such as effective wave height and average wave period can be obtained through spectrum analysis.

[0172] (1) Frequency domain pressure-wavefront transfer function

[0173] With the still water surface as the origin and vertical upward as the positive direction, the seabed position z = -d. According to linear theory, dynamic pressure and free surface displacement... The following transfer relationship is satisfied in the frequency domain:

[0174] ;

[0175] in, For the Fourier transform of dynamic pressure For the dynamic pressure response transfer function:

[0176] ;

[0177] In the formula, k is the wavenumber, which is determined by the linear dispersion relation:

[0178] ;

[0179] For each frequency Solving the above equation yields the corresponding wave number. , thus calculating .

[0180] (2) Reconstruction of the time equation for water surface fluctuation

[0181] Dynamic pressure The water surface displacement spectrum was obtained after processing. Then, the time-series water surface fluctuation sequence is obtained through inverse Fourier transform. If the water is shallow and the wave frequency is low, then the following conditions are met. ,but Then the transfer function is approximately 1, which can be simplified to:

[0182] ;

[0183] This condition only applies to shallow water conditions; for general water depths, a complete transfer function must be used.

[0184] (3) Calculation of energy spectral density

[0185] For the reconstructed water surface wave sequence Its power spectral density was estimated using Welch's method. .Will The data is divided into N segments, each with a length L (a power of 2, such as 1024), with 50% overlap between adjacent segments. A Hanning window is applied to each segment to reduce spectral leakage. A Fast Fourier Transform is performed on each segment to calculate the one-sided power spectrum. Finally, the spectra of all segments are averaged to obtain a smoothed power spectral density estimate.

[0186] ;

[0187] final Let be the power spectral density of the wave surface displacement, whose integral is equal to the variance of the water surface ripple.

[0188] (4) Wave parameter extraction

[0189] Calculate the spectral moment from the power spectral density:

[0190] ;

[0191] The significant wave height is calculated using the following formula:

[0192] ;

[0193] In the formula, the coefficient 4.004 is accurately obtained based on the relationship between wave height and variance under the narrow spectrum assumption, and 4 can also be used in practical applications. The average wave period adopts the zero-period transition. :

[0194] ;

[0195] in, The zeroth moment of the spectrum is obtained by weighted integration of the wave spectrum; is the second moment of the spectrum.

[0196] By following the steps above, a high-precision wave spectrum and characteristic parameters can be derived from the time series of a single pressure sensor.

[0197] S4.2: Wave Spectrum Inversion

[0198] This step is based on the pure wave pressure signals extracted from each sensor by S3. By combining the spatial geometry of the sensor array, a high-resolution frequency-directional spectrum can be retrieved. Then, wave direction parameters such as the main wave direction and the directional distribution width are extracted.

[0199] (1) Construction of cross-spectral matrix

[0200] For N pressure sensors (N=3 in this scheme), calculate the wave pressure signals of any two sensors m and n. and cross power spectral density .

[0201] right and Perform Fourier transform to obtain and Thus, the cross-power spectral density can be calculated.

[0202] ;

[0203] in, for The complex conjugate;

[0204] The cross-spectral matrix of all sensors is formed in the following form. :

[0205] ;

[0206] (2) Definition of array manifold vector

[0207] Define frequency and wave direction The array manifold vector :

[0208] ;

[0209] in, For the first The horizontal coordinates of each sensor (with the array center as the origin), where , It can be obtained from step one of S4.1.

[0210] (3) Directional spectrum estimation algorithm

[0211] The iterative maximum likelihood method is used as the core inversion algorithm. For each frequency... directional spectrum The solution is obtained through the following iterative process. Initial direction spectrum is set. For uniform distribution, or as a preliminary estimate obtained from traditional beamforming methods:

[0212] ;

[0213] in, This is the conjugate transpose of the array manifold vector. It is a cross-spectral matrix. For array manifold vectors;

[0214] Iterative update: For the p-th iteration, the covariance matrix is ​​first reconstructed based on the directional spectrum of the previous round.

[0215] ;

[0216] in, Let be the covariance matrix obtained from the p-th iteration reconstruction. This is the directional spectrum at the current frequency obtained in the (p-1)th iteration;

[0217] Then update the directional spectrum:

[0218] ;

[0219] The iteration terminates when the difference in spectral estimation between two adjacent iterations is less than a preset value or when the maximum number of iterations is reached.

[0220] (4) Fusion of trend prior constraints

[0221] This method utilizes the tidal information extracted from the local tidal model in step S2 as prior conditions and improves the inversion accuracy through a Bayesian framework. According to tidal harmonic analysis theory, tidal currents can be represented as the superposition of several tidal constituents. For a single tidal constituent, the endpoint of the tidal velocity vector on the horizontal plane traces an elliptical trajectory over time, i.e., the tidal ellipse. The harmonic constants of each tidal constituent obtained from S2... The main current direction of the tidal current ellipse is calculated based on tidal dynamics. The average mainstream direction within the array coverage area As a priori wave direction reference, a penalty term is introduced into the objective function of the iterative maximum likelihood algorithm to construct a new directional spectrum estimation criterion: ;

[0222] in, The final frequency-direction spectrum; This is a mathematical operator that represents finding the minimum value of the objective function within the parentheses. ML(E) is the original maximum likelihood objective function. This is the regularization parameter (set through cross-validation or empirically, such as 0.1-0.3). The penalty function is defined as follows:

[0223] ;

[0224] in, The function represents the wave direction; it takes a maximum of 1 when the wave direction is perpendicular to the tidal current direction and 0 when it is parallel, thus penalizing wave direction solutions that deviate significantly from the tidal current direction. This constrained objective function is embedded into the iterative process of the maximum likelihood objective function. In each iteration of calculating the directional spectrum, soft constraints are applied to the estimation results, ultimately fusing the directional spectrum with fused tidal current prior information to obtain the final directional spectrum. .

[0225] (5) Wave direction parameter extraction

[0226] Directional spectrum obtained from inversion Calculate the following wave direction parameters.

[0227] Integrating the directional spectrum with respect to frequency yields the directional distribution. .

[0228] Integrating the directional spectrum with respect to the direction yields the frequency spectrum. .

[0229] For each frequency Main wave direction Pick The angle corresponding to the maximum value. Overall main wave direction. Pick The angle corresponding to the maximum value.

[0230] Step S5: Output and Visualization

[0231] Finally, the effective wave height, average wave period, main wave direction, directional spectrum, and other information obtained from the inversion, as well as intermediate results in the processing (such as local harmonic constants), are formatted and graphically displayed for user analysis.

[0232] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A wave parameter inversion method based on a pressure sensor array under strong tidal conditions, characterized in that, Includes the following steps: S1: Synchronously collect pressure data from each sensor in the pressure sensor array deployed in the target sea area, and obtain the average water depth; S2: For the pressure data of each sensor, combined with the regional tidal prior information, a dynamic localized tidal model for each sensor location is constructed through harmonic analysis with regularization constraints and robust estimation, and tidal forecast sequences are generated respectively. S3: Calculate the common pattern of the tide forecast sequences of all sensors and subtract it from the raw pressure data of each sensor; perform singular value decomposition on the pressure data array matrix after subtracting the common pattern, extract and eliminate the first principal component, and obtain the pure wave pressure signal of each sensor. S4: Based on the pure wave pressure signal, the wave height and wave period are inverted using the pressure-wave surface transfer function; based on the pure wave pressure signal, the cross spectrum matrix is ​​calculated and the wave direction spectrum is inverted using the direction spectrum estimation algorithm to obtain the wave direction information; S2 specifically includes: S2.1: Construct observation equations to represent the tidal components observed by each sensor as a linear superposition of multiple tidal constituents, and establish the relationship between the harmonic constant to be determined and the observed values; S2.2: Introduce prior constraints, obtain prior harmonic constants from the regional tidal model, and construct an objective function containing regularization terms to constrain the parameters to be determined from deviating too much from the prior values ​​under the condition of short observation sequences; S2.3: Robust iterative solution: Based on the objective function, iterative reweighted least squares method is used to solve the problem. During the iteration process, the weights of the observed data are dynamically adjusted according to the residual size to suppress the influence of outliers. After convergence, a stable local harmonic constant is obtained. S2.4: Generate local tide forecasts by using the stable local harmonic constant to generate a tide forecast sequence for each sensor at any time. S3 specifically includes: S3.1: Common tidal pattern extraction and preliminary subtraction: that is, the arithmetic mean of the tidal forecast sequences of all sensors is calculated as the common tidal pattern, and the common pattern is subtracted from the original pressure data of each sensor to obtain the signal after preliminary de-tidal processing. S3.2: Cooperative extraction and elimination of residual common low-frequency noise: The signals from each sensor after initial de-tidal processing are organized into a data matrix. Singular value decomposition is performed on the matrix to extract the first principal component with the highest energy and remove it from the matrix to obtain the pure wave pressure signal from each sensor, which is then used for subsequent wave parameter inversion.

2. The wave parameter inversion method based on a pressure sensor array under strong tidal conditions according to claim 1, characterized in that, The pressure sensor array consists of a triangular array of three underwater pressure sensors.

3. The wave parameter inversion method based on a pressure sensor array under strong tidal conditions according to claim 1, characterized in that, In S2.1, for the pressure data of the i-th sensor Its tidal components are represented as the superposition of multiple tidal constituents: ； Where i is the sensor number. Let j be the tide observed by the i-th sensor, and j be the tide fraction number. M is the selected number of tidal constituents. For the nodal factor of the j-th tidal constituent, Let be the astronomical initial phase angle of the j-th tidal constituent. The correction angle for the j-th tidal constituent node. , Let be the amplitude and lag angle of the i-th sensor at the j-th tidal fraction, respectively, i.e., the harmonic constant. This is the average sea level.

4. The wave parameter inversion method based on a pressure sensor array under strong tidal conditions according to claim 3, characterized in that, In S2.2, the objective function containing the regularization term is constructed as follows: ； in, The design matrix is ​​K×(2M+1), and its k-th row is: ； in, For the Mth tidal component The node factor at time, For the Mth tidal component Phase angle at time, , Let ω be the angular velocity of the Mth tidal constituent. For the Mth tidal component The node factor at time, For the Mth tidal component The nodal correction angle at time , It is the mean sea level. There are 2M+1 parameter vectors to be determined. ,in, Let be the cosine component coefficient of the i-th sensor at the j-th tidal fraction. Let be the sinusoidal component coefficient of the i-th sensor at the j-th tidal fraction. , , Let be the amplitude and lag angle of the i-th sensor at the j-th tidal fraction, respectively. It is the mean sea level. Let K×1 be the observation vector, consisting of tidal observations from K time sampling points. Composition, where K is the time sampling point; This is a regularization parameter used to balance the strength of data fitting and prior constraints; The weight matrix is ​​(2M+1)×(2M+1), and the diagonal matrix is ​​taken, with the diagonal elements being the reciprocal of the confidence level of each prior parameter; For prior parameters, ； in, Let be the prior cosine component coefficients of the i-th sensor at the j-th tidal fraction. This represents the prior sinusoidal component coefficient of the i-th sensor at the j-th tidal constituent. , It is the amplitude of the i-th sensor at the j-th tidal constituent, obtained from the regional tidal model. It is the lag angle of the i-th sensor at the j-th tidal constituent, obtained from the regional tidal model.

5. The wave parameter inversion method based on a pressure sensor array under strong tidal conditions according to claim 4, characterized in that, The specific steps of S2.3 are as follows: (1) Initialize the parameter vector Initialize the observation weight matrix , Let K be a K×K identity matrix, and let the iteration counter p=0; (2) Iterative calculation: in the first... In the next iteration, the current weight matrix is ​​used. Solve for parameters : ； in, For the i-th sensor at the th The parameter vector after the nth iteration. For the first The observation weight matrix after the next iteration; (3) Calculate the observation residuals under the current iteration: ； in, Indicates the first Observation residuals under the next iteration; And the robust standard deviation is estimated using the median absolute deviation: ； in, This represents the median; (4) Update the observation weights using the Huber weight function: ； in, Represents a diagonal matrix; express The observation weight matrix after the next iteration has diagonal elements representing the weights of each observation in the current iteration. express After the nth iteration Weights of each observation: ； in, This represents the residual of the k-th observation. The tuning constant is... For robust standard deviation; (5) Convergence judgment: ； when If the threshold value is less than ε, the iteration terminates; otherwise, the iteration continues. ε is the threshold value. (6) After the iteration converges, stable linear parameter estimates are obtained. Then, the estimated harmonic constants are obtained through inverse transformation: ； ； in, This represents the amplitude estimate of the j-th tidal fraction from the i-th sensor; This represents the estimated lag angle of the j-th tidal constituent from the i-th sensor; This represents the estimated cosine component of the j-th tidal constituent from the i-th sensor. This represents the estimated sinusoidal component of the j-th tidal constituent from the i-th sensor. This is an overestimate of the mean sea level.

6. The wave parameter inversion method based on a pressure sensor array under strong tidal conditions according to claim 5, characterized in that, In S2.4, the harmonic constant is obtained using S2.

3. and mean sea level Generate the tidal forecast value for the i-th sensor at time t: ； in, Let be the tidal forecast value of the i-th sensor at time t. For the j-th tidal component The node factor at time, Let be the astronomical initial phase angle of the j-th tidal constituent. The correction angle for the j-th tidal constituent.

7. The wave parameter inversion method based on a pressure sensor array under strong tidal conditions according to claim 1, characterized in that, In S3.2, the observation matrix is ​​first constructed, and the signal after preliminary tidal removal is... A data matrix organized as N×K Where N is the number of sensors and K is the time sampling point; matrix elements This represents the initial detidal signal of the i-th sensor at the k-th sampling time; Again Perform singular value decomposition: ； in, It is an N×N left singular vector matrix, where each column represents a spatial mode; It is an N×K diagonal matrix; It is a K×K right singular vector matrix, where each column represents the time coefficient of the corresponding spatial mode; Subtracting the reconstructed component of the first principal component from the original data matrix yields the data matrix after eliminating common noise: ； in, The maximum singular value obtained from singular value decomposition. The first column of the left singular vector matrix U, This is the transpose of the first row of the right singular vector matrix V; The i-th row vector is the pure wave pressure signal of the i-th sensor after array collaborative noise reduction.

8. The wave parameter inversion method based on a pressure sensor array under strong tidal conditions according to claim 1, characterized in that, The directional spectrum estimation algorithm described in S4 is an iterative maximum likelihood method, and it incorporates local power flow information obtained from S2 as a priori constraint during the iteration process, as detailed below: Extracting the main current direction from the local tidal model Construct a penalty function : ； in, Indicates wave direction; The penalty function is introduced into the objective function of the iterative maximum likelihood method to form a constrained directional spectrum estimation criterion: ； in, The final frequency-direction spectrum; This is a mathematical operator that represents finding the minimum value of the objective function within the parentheses. ML(E) is the original maximum likelihood objective function, and α is the regularization parameter; The frequency-direction spectrum is continuously updated during the iteration process; During each iteration of updating the directional spectrum, a penalty is imposed on wave direction solutions that deviate significantly from the mainstream trend of the current, ultimately yielding a directional spectrum that incorporates prior information about the current. .

Images