A sound speed profile assimilation inversion method and system based on physical data double driving

Abstract

The application discloses a sound velocity profile assimilation inversion method and system based on physical data double driving, and the method comprises the following steps: adaptive observation and sound propagation path sparse representation are performed to obtain a sound propagation path sparse representation result; long-time noise cross-correlation calculation and passive time delay self-calibration are performed to obtain a virtual sound source arrival time sequence corrected by an absolute time delay reference; constraint verification is performed through a physical-data double driving double flow coupling generation type inversion network; a wide-band electromagnetic response time sequence is collected, cross-medium electromagnetic-acoustic joint inversion constraint is performed, and a full-water-depth sound velocity profile after strong constraint calibration is obtained; and set Kalman filtering assimilation and uncertainty propagation based on POD reduction are performed. The application can realize uncertainty quantification of an inversion result. The application, as a sound velocity profile assimilation inversion method and system based on physical data double driving, can be widely applied to the technical field of underwater acoustic engineering.

Description

Technical Field

[0001] This invention relates to the field of underwater acoustic engineering technology, and in particular to a method and system for assimilation and inversion of sound velocity profiles based on dual physical data driving. Background Technology

[0002] The speed of sound in seawater is the most fundamental physical parameter determining the characteristics of underwater sound propagation. It directly affects the propagation path, propagation loss, effective range, and accuracy of target detection and localization. Whether for underwater target detection and localization, underwater acoustic communication, underwater navigation, marine environmental monitoring, or seabed topography mapping, accurate understanding of the seawater sound speed field is indispensable. The distribution of sound speed in seawater is not uniform but varies with temperature, salinity, and hydrostatic pressure, exhibiting a complex vertically layered structure known as the Sound Speed ​​Profile (SSP). The morphology of the SSP—especially characteristic structures such as the surface mixing layer, the principal transition layer, the deep-sea acoustic axis, and the deep-sea isotherm—directly determines the curvature of sound rays underwater, the location of convergence zones, and the propagation efficiency of sound energy. Therefore, rapidly and accurately obtaining a full-depth sound speed profile is of crucial strategic importance for improving the effectiveness of underwater information systems.

[0003] In recent years, with the development of artificial intelligence technology, deep learning methods have been introduced into the field of sound velocity profile inversion, such as Long Short-Term Memory (LSTM) networks, Convolutional Neural Networks (CNN), and Graph Attention Networks (GAT), showing potential in capturing the nonlinear mapping relationship between sea surface data and underwater sound velocity. However, existing deep learning inversion methods still follow the basic paradigm of "point-by-point inversion" or "sea surface information driven," and have significant shortcomings in terms of physical interpretability, uncertainty quantification, and multiphysics fusion.

[0004] Current mainstream sound velocity profile inversion methods have the following limitations: EOF statistical regression extracts orthogonal modes from historical data and establishes a mapping between sea surface parameters and coefficients, but mode truncation leads to information loss, point-by-point inversion efficiency is low, and relying solely on sea surface data cannot capture complex disturbances such as internal waves and vortices. Acoustic tomography matching field method utilizes sound propagation time / angle inversion, with a solid physical foundation, but it is computationally intensive, sensitive to model accuracy, and prone to getting trapped in local optima. While recent integration with deep learning has improved efficiency, it still relies on active sound sources and has limited ray coverage. Deep learning methods (LSTM, 3D-CNN, GAT) can establish end-to-end mappings of observed profiles, improving accuracy and efficiency, but their black-box nature leads to poor physical interpretability and a lack of uncertainty quantification.

[0005] In summary, the gaps in existing technologies are mainly reflected in the following aspects: the observation system and inversion algorithm are open-loop, making it impossible to dynamically adjust the path; it relies on active sound sources and lacks passive inversion and time delay benchmark schemes; the integration of physical laws and data-driven approaches is insufficient, and the modes and physical space are separated; there is a lack of cross-medium (seawater-seabed) constraints; the inversion results lack system uncertainty quantification; and the physical characteristics such as the deep-sea acoustic channel axis are not used as constraints, resulting in distortion of deep-sea inversion. Summary of the Invention

[0006] To address the aforementioned technical problems, the present invention aims to provide a sound velocity profile assimilation and inversion method and system based on dual physical data driving, which can quantify the uncertainty of the inversion results.

[0007] The first technical solution adopted in this invention is: a sound velocity profile assimilation and inversion method based on dual-drive physical data, comprising the following steps:

[0008] The planned path of the underwater unmanned swarm is generated and adaptive observation and sparse representation of the sound propagation path are performed to obtain the sparse representation result of the sound propagation path.

[0009] Marine environmental noise data was collected and long-term noise cross-correlation calculations and passive time delay self-calibration were performed to obtain the virtual sound source arrival time series after absolute time delay reference correction.

[0010] By combining the sparse characterization results of the sound propagation path with the arrival time series of the virtual sound source, and performing constraint verification through a physical-data dual-drive dual-stream coupled generative inversion network, the verified full-depth sound velocity is obtained.

[0011] Broadband electromagnetic response time series were collected, and combined with the verified full-depth sound velocity, cross-medium electromagnetic-acoustic joint inversion constraints were performed to obtain a strongly constrained calibrated full-depth sound velocity profile.

[0012] The full-depth sound velocity profile after strong constraint calibration is assimilated by ensemble Kalman filtering based on POD order reduction and uncertainty propagation to obtain the full-depth sound velocity mean profile and the full-depth sound velocity standard deviation profile.

[0013] Furthermore, the step of generating the planned path of the underwater unmanned swarm and performing adaptive observation and sparse representation of the sound propagation path to obtain the sparse representation result of the sound propagation path specifically includes:

[0014] Based on the historical background field or the inversion results of the previous moment, the sensitivity distribution of the full-depth sound velocity field to each mode of the empirical orthogonal function is calculated.

[0015] Based on the sensitivity distribution, the curve with the largest information entropy gradient is selected by locating high sensitivity distributions above a preset threshold, and used as the planning path for the underwater unmanned swarm.

[0016] Based on the planned path of the underwater unmanned swarm, while the swarm is sailing along the planned path, each node acts as a transmitter and receiver for each other, forming a sparse acoustic link network.

[0017] Based on a sparse acoustic link network, the multipath arrival times of several links are extracted and sparsely represented using compressed sensing methods to obtain sparse representation results of the sound propagation path.

[0018] Furthermore, the step of collecting marine environmental noise data and performing long-term noise cross-correlation calculations and passive time delay self-calibration to obtain the virtual sound source arrival time series after absolute time delay reference correction specifically includes:

[0019] Based on a deep-sea fixed vertical array, marine environmental noise data is collected and superimposed with long-term noise for cross-correlation calculation to obtain the arrival time of virtual sound sources.

[0020] The theoretical arrival time of the sound source is determined based on the depth difference of seawater and the compressibility coefficient of seawater;

[0021] Based on the deviation between the arrival time of the virtual sound source and the arrival time of the theoretical sound source, the absolute time drift and the correction value for the micro-perturbation of the array element position are determined.

[0022] Absolute time calibration is performed based on the absolute time drift and the micro-perturbation correction value of the array element position. The propagation time delay structure of the virtual sound source from the sea surface to the seabed is derived, and the arrival time sequence of the virtual sound source after absolute time delay reference correction is obtained.

[0023] Furthermore, the step of combining the sparse characterization results of the sound propagation path with the arrival time series of the virtual sound source, and performing constraint verification through a physical-data dual-drive dual-stream coupled generative inversion network to obtain the verified full-depth sound velocity specifically includes:

[0024] We acquire sea surface satellite remote sensing data, and combine the sparse characterization results of sound propagation paths with the arrival time series of virtual sound sources after absolute time delay benchmark correction to perform data preprocessing and construct a feature dataset.

[0025] Based on a multi-layer fully connected sensing mechanism, the modal flow is modeled, and the physical flow is constructed by combining a deep residual network and a physical information neural network. The modal flow and physical flow are integrated and a cross-attention mechanism is introduced to construct a dual-stream parallel network.

[0026] The feature dataset is input into a dual-stream parallel network. Modal orthogonality constraints are performed based on modal flow, and empirical orthogonal function coefficients are output. Sound velocity values ​​are calculated based on physical flow, and a full-depth sound velocity profile is output.

[0027] A generative adversarial network (GAN) discriminant is introduced and pre-trained to construct a trained GAN discriminant.

[0028] The trained generative adversarial network discriminator performs a physical consistency check on the empirical orthogonal function coefficients and the full-depth sound velocity profile to obtain the verified full-depth sound velocity.

[0029] Furthermore, the step of integrating modal streams and physical streams and introducing a cross-attention mechanism to construct a dual-stream parallel network specifically includes:

[0030] The last layer of the modal flow is mapped to a query vector through a linear transformation, and the feature vectors corresponding to the intermediate feature layers of the physical flow are extracted as key vectors and value vectors.

[0031] By calculating the dot product similarity between the query vector and the key vector, and then normalizing it with Softmax to obtain the attention weight, the value vector is then weighted and summed using the attention weight to obtain the fused feature vector, thus realizing the fusion of modal flow to physical flow.

[0032] The feature vectors corresponding to the intermediate feature layers of the physical flow are mapped to key vectors through a linear transformation. The feature map of the current layer of the modal flow is used as the query vector, and the attention weights between the query vector and the key vector are calculated.

[0033] The current feature map of the modal flow is weighted and corrected by attention weights to obtain a correction vector, which is then added to the original features of the modal flow and continued to propagate forward to achieve the fusion of physical flow and modal flow.

[0034] Furthermore, the step of acquiring a broadband electromagnetic response time series, combining it with the verified full-depth sound velocity, and performing cross-medium electromagnetic-acoustic joint inversion constraints to obtain a strongly constrained calibrated full-depth sound velocity profile specifically includes:

[0035] Broadband electromagnetic response time series and electric and magnetic field component data are collected and Fourier transformed to obtain the frequency domain electromagnetic response;

[0036] The resistivity structure of the frequency domain electromagnetic response is inverted using a three-dimensional electromagnetic inversion algorithm to obtain three-dimensional resistivity grid data.

[0037] Based on three-dimensional resistivity grid data and combined with sediment grain size parameters, sound velocity conversion was performed through a rock physics model to obtain a sound velocity dispersion curve table and a reflection coefficient matrix.

[0038] Based on the verified full-depth sound velocity, combined with the seawater density and bottom seawater sound velocity gradient obtained from the temperature-salt-pressure state equation, the acoustic impedance of the seawater side of the seabed interface is calculated.

[0039] Obtain the sound velocity and density of the seabed sediment layer to determine the acoustic impedance on the sediment layer side;

[0040] The ratio of the acoustic impedance on the seawater side to the acoustic impedance on the sedimentary layer side of the seabed interface is obtained to calculate the reflection coefficient at vertical incidence. The result is compared with the reflection coefficient matrix and the relative deviation signal is calculated.

[0041] By recalculating the acoustic impedance and reflection coefficient using the perturbation gradient value, the sensitivity of the bottom seawater sound velocity gradient to the reflection coefficient deviation is obtained, and the Jacobian matrix of the gradient change on the reflection coefficient is obtained.

[0042] The relative deviation signal is converted into a corrected gradient of the hidden layer weights related to the bottom seawater in the core inversion network. Combined with the Jacobian matrix, it is backpropagated to the inversion network through the backpropagation algorithm for iterative correction, resulting in a fully constrained calibrated sound velocity profile.

[0043] Furthermore, the step of converting sound velocity based on three-dimensional resistivity grid data, combined with sediment grain size parameters, and using a rock physics model to obtain the sound velocity dispersion curve table and reflection coefficient matrix specifically includes:

[0044] Based on three-dimensional resistivity grid data and combined with sediment grain size auxiliary parameters obtained through historical geological sampling, an improved Biot-Stoll two-phase medium model framework was constructed. Combined with pore water resistivity, sediment porosity was obtained by back-calculating grid cell by grid cell using the Archie formula.

[0045] Based on sediment porosity, combined with the solid skeleton bulk modulus and shear modulus obtained from the table and the pore fluid bulk modulus calculated from temperature, salinity and pressure, the data are input into the improved Biot-Stoll two-phase medium model framework. Considering the Biot dissipation mechanism, the longitudinal wave velocity dispersion curve, transverse wave velocity and attenuation coefficient are calculated.

[0046] Based on the longitudinal wave velocity dispersion curve, transverse wave velocity and attenuation coefficient, the modulus and phase of the plane wave reflection coefficient at different grazing angles are calculated using a multilayer medium reflection coefficient algorithm, resulting in a velocity dispersion curve table and a reflection coefficient matrix.

[0047] Furthermore, the step of performing ensemble Kalman filtering assimilation and uncertainty propagation based on POD-based order reduction on the full-depth sound velocity profile after strong constraint calibration to obtain the full-depth sound velocity mean profile and the full-depth sound velocity standard deviation profile specifically includes:

[0048] The verified full-depth sound velocity is arranged in the same depth order, a matrix is ​​constructed and singular value decomposition is performed to obtain the left singular vector, singular value and right singular vector. The left singular vector is regarded as the POD mode basis.

[0049] Based on the cumulative contribution rate of singular values, the pre-preset order POD mode basis is selected and retained, and a mode basis matrix is ​​constructed.

[0050] New observation data are acquired and combined with the modal basis matrix. Kalman filtering is then performed in a low-dimensional space to assimilate the data and obtain a set of sound velocity profiles for the analysis field.

[0051] By combining the set of sound velocity profiles in the analysis field with the full-depth sound velocity profile after strong constraint calibration, the uncertainty of water depth adaptive quantification is performed to obtain the mean sound velocity profile and the standard deviation sound velocity profile of the full depth.

[0052] Furthermore, the step of acquiring new observation data, combining it with the modal basis matrix, and performing Kalman filtering assimilation in a low-dimensional space to obtain a set of sound velocity profiles for the analysis field specifically includes:

[0053] Acquire new observation data, which includes sparse field profile measurement data and multiple-path arrival times of several links;

[0054] Based on the low-dimensional POD coefficient space, the mean coefficient vector is calculated according to the modal basis matrix, and the deviation of each mean coefficient vector from the mean is calculated to construct the background field covariance matrix.

[0055] Based on the new observation data, the corresponding observation operator is constructed, and the full-depth sound velocity profile is reconstructed by combining it with the modal basis matrix. Then, the sound velocity value at the observation depth is extracted, and the theoretical arrival time of the sound wave is calculated.

[0056] Based on the background field covariance matrix, the theoretical arrival time of sound waves, and the error covariance of the observation data, the Kalman gain matrix is ​​calculated to determine the correction weight of the observation data on the background field.

[0057] Based on the Kalman gain matrix and the corrected weights, the POD coefficients are updated and calculated, multiplied by the POD mode basis, and projected back to the original high-dimensional physical space to obtain the set of sound velocity profiles of the analysis field.

[0058] The second technical solution adopted in this invention is: a sound velocity profile assimilation and inversion system based on dual physical data driving, comprising:

[0059] The first module is used to generate the planned path of the underwater unmanned swarm and perform adaptive observation and sparse representation of the sound propagation path to obtain the sparse representation result of the sound propagation path.

[0060] The second module is used to collect marine environmental noise data and perform long-term noise cross-correlation calculation and passive time delay self-calibration to obtain the virtual sound source arrival time sequence after absolute time delay reference correction.

[0061] The third module is used to combine the sparse characterization results of the sound propagation path with the arrival time series of the virtual sound source, and perform constraint verification through a physical-data dual-drive dual-flow coupled generative inversion network to obtain the verified full-depth sound velocity.

[0062] The fourth module is used to acquire broadband electromagnetic response time series, and combined with the verified full-depth sound velocity, to perform cross-medium electromagnetic-acoustic joint inversion constraints, and obtain a strongly constrained calibrated full-depth sound velocity profile.

[0063] The fifth module is used to perform ensemble Kalman filtering assimilation and uncertainty propagation on the full-depth sound velocity profile after strong constraint calibration, based on POD order reduction, to obtain the full-depth sound velocity mean profile and the full-depth sound velocity standard deviation profile.

[0064] The beneficial effects of the method and system of this invention are as follows: This invention generates a planned path for an underwater unmanned swarm and performs adaptive observation and sparse characterization of the sound propagation path to obtain the sparse characterization result of the sound propagation path. It directly links the motion planning of the observation platform with the "information requirement" for the sound velocity mode within the inversion algorithm, achieving a leap from "data-driven inversion" to "inversion-driven observation." Furthermore, it collects marine environmental noise data and performs long-term noise cross-correlation calculations and passive time delay self-calibration to obtain a virtual sound source arrival time series corrected for an absolute time delay reference. This transforms the deep-sea background noise field into a self-calibrated, continuously usable "passive sound source." Moreover, it utilizes the immutable physical laws of the deep-sea isotherm to solve the long-standing problem of the absolute time delay reference in passive acoustic tomography, enabling the system to possess all-weather, full-depth, zero-energy passive inversion capabilities. Further, it combines the sparse characterization result of the sound propagation path with the virtual sound source arrival time series and performs constraint verification through a physical-data dual-drive, dual-flow coupled generative inversion network to obtain the verified full-depth sound velocity. By combining dual-stream coupling and generative adversarial verification, the global statistical laws of modal space and the local differential constraints of physical space are organically unified, significantly improving the physical realism and reliability of the inversion results. Broadband electromagnetic response time series are acquired, and combined with verified full-depth sound velocity, cross-medium electromagnetic-acoustic joint inversion constraints are applied to obtain a strongly constrained calibrated full-depth sound velocity profile. The cross-medium electromagnetic detection data is transformed into quantifiable acoustic boundary conditions through a rock physics model, and a closed-loop feedback from boundary conditions to inversion network weights is established, enabling physical consistency coupling between seawater sound velocity inversion and seabed sediment characteristics. Finally, the strongly constrained calibrated full-depth sound velocity profile is assimilated and uncertainty propagated using ensemble Kalman filtering based on POD order reduction, yielding a full-depth sound velocity mean profile and a full-depth sound velocity standard deviation profile. Real-time observation data is used to update the profile in low-dimensional space using ensemble Kalman filtering, and the updated results are projected back into physical space. Simultaneously, the probability density function of the sound velocity value at each depth layer is output, realizing the quantification of uncertainty in the inversion results. Attached Figure Description

[0065] Figure 1 This is a flowchart of the steps of a sound velocity profile assimilation and inversion method based on dual physical data driving according to the present invention.

[0066] Figure 2 This is a structural block diagram of a sound velocity profile assimilation and inversion system based on dual physical data driving according to the present invention. Detailed Implementation

[0067] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. The step numbers in the following embodiments are only for ease of explanation and do not limit the order of the steps. The execution order of each step in the embodiments can be adapted according to the understanding of those skilled in the art.

[0068] First, it's important to note that the ocean is a highly dynamic, nonlinear system. Ocean dynamics processes such as mesoscale eddies, ocean fronts, and internal waves constantly alter the water's temperature, salinity, and structure, leading to significant spatiotemporal variations in the sound velocity field. For example, mesoscale eddies can cause localized shifts in the sound velocity profile of a given area of ​​tens of meters per second, while internal wave activity can induce high-frequency oscillations at the depth of the thermocline, posing significant challenges to real-time monitoring of the sound velocity field. In the complex deep-sea environment, seawater temperature fluctuates irregularly over time, resulting in strong nonlinearity and uncertainty in the sound velocity profile across time and space. Simultaneously, the increasing depth of human ocean activities—from deep-sea resource exploration and subsea pipeline laying to the widespread application of underwater unmanned platforms—places higher demands on sound velocity profile acquisition: not only high precision but also wide coverage, high temporal resolution, and real-time performance are required. However, directly obtaining full-depth sound velocity profiles through on-site observation is costly and limited by the spatiotemporal coverage capabilities of observation platforms, making it difficult to meet these requirements.

[0069] Domestic and international scholars have conducted long-term research on the acquisition and inversion of sound velocity profiles, resulting in three main technical approaches. The first is the direct measurement method, which includes direct on-site measurements using a sound velocity profiler (CTD) or a discardable temperature-depth meter (XBT), and indirect calculation of sound velocity using temperature, salinity, and depth data combined with empirical formulas. This method offers the highest accuracy, but is limited by ship time costs and the sparsity of observation stations, making it difficult to construct a large-scale, high-frequency full-depth sound velocity field. The second is the statistical regression reconstruction method, which uses empirical orthogonal function (EOF) regression as its core framework. It establishes statistical relationships between sea surface satellite remote sensing data (such as sea surface temperature SST and sea surface height anomaly SLA) and historical sound velocity profile data, thereby achieving the goal of extrapolating the full-depth profile from sea surface information. This method improves convenience, but its accuracy is highly dependent on the representativeness of the historical database, and EOF mode truncation can lead to information loss. The third method is acoustic inversion, which utilizes the path integral information carried by the sound propagation signal (such as acoustic tomography or matched field processing) to invert the equivalent sound velocity profile through a sound field propagation model. This type of method has strong physical interpretability, but it is computationally intensive, sensitive to the accuracy of sound field modeling, and difficult to cover areas that sonar systems cannot reach.

[0070] Therefore, the existing technology has the following shortcomings, for example:

[0071] 1) The observation system and inversion algorithm are in an open-loop state, and the observation efficiency has not been fully utilized. Currently, whether based on Argo buoys, moorings, or underwater unmanned vehicles, the station deployment and route planning are independent of the inversion algorithm itself. Observation platforms typically perform observation tasks according to preset regular grids or fixed depth layers, failing to adaptively adjust to the dynamic changes in the information requirements of each mode of the sound velocity field according to the inversion algorithm. This leads to a severe mismatch between the spatial distribution of observation data and the sensitivity distribution of the inversion algorithm to information entropy. There is insufficient observation in key depth layers requiring dense observation (such as the thermocline and near the sound channel axis), while the same amount of observation resources are invested in weakly changing areas with low information gain. The lack of a closed-loop mechanism for real-time interaction between observation and inversion prevents the limited resources of the observation platform from maximizing information capture.

[0072] 2) Strong dependence on active sound sources and a severe lack of passive inversion capabilities. Existing acoustic inversion methods (such as acoustic tomography and matched field inversion) generally rely on actively emitted sound sources or specifically designed sound propagation experiments, requiring continuous energy consumption and exposure of their own location, making it difficult to meet the needs of long-term, continuous, and covert ocean monitoring. Although ocean background noise contains rich propagation information, current technologies have not yet been able to effectively convert it into usable inversion observations. The key challenges are: how to extract coherent propagation structures from random noise fields, how to solve the problem of absolute time delay reference under passive observation, and how to utilize inherent physical characteristics such as the deep-sea isotherm for self-calibration of the observation system. These are all technical blind spots that current technologies have not yet addressed.

[0073] 3) Insufficient integration of physical laws and data-driven approaches leads to poor physical interpretability of inversion results. While emerging deep learning methods (such as LSTM, CNN, and graph neural networks) have improved inversion accuracy to some extent, their "black box" nature makes it difficult to interpret the output results from a physical perspective and to guarantee compliance with the inherent laws of ocean dynamics. Existing methods either rely entirely on data-driven approaches, ignoring the physical constraints that sound velocity profiles should satisfy (such as the continuity of the thermocline gradient and the temperature-salinity relationship implicit in the seawater equation of state); or they adhere rigidly to physical models (such as the matched field method), making it difficult to fully utilize the statistical regularities contained in massive historical data. In particular, there is a general information gap between the modal space (EOF coefficients) and the physical space (depth-sound velocity values), failing to achieve collaborative optimization and mutual verification between the two spaces at the algorithm architecture level. This results in physically unreasonable oscillations in the inversion results at depths lacking observational data.

[0074] 4) The unique physical structures of the deep sea have not been fully explored and utilized, making it difficult to guarantee the accuracy of deep-sea inversion. The deep-sea sonic duct axis (SOFAR) is the most physically anchoring feature in the sound velocity profile, exhibiting a "capture" effect on sound rays and directly determining the energy distribution of long-range sound propagation. However, existing inversion methods typically treat the entire water depth as a continuous layer for uniform processing, failing to incorporate characteristic parameters such as the depth, intensity, and thickness of the sonic duct axis as physical constraints into the inversion process. In the absence of measured deep-sea data, existing methods often lack physical anchoring for inversion results below the deep-sea isotherm and sonic duct axis, leading to an accumulation and amplification of deep-sea sound velocity errors with increasing depth, severely impacting the reliability of deep-sea sound propagation predictions.

[0075] 5) The inversion results lack systematic uncertainty quantification, making it difficult to support high-reliability applications. Existing technologies generally output a single "optimal" sound velocity profile, ignoring the non-uniform uncertainty distribution characteristics of the inversion results as a function of water depth, observation density distribution, and physical constraint strength. This "deterministic output" cannot provide reliable error envelope information for downstream applications (such as underwater target location, weapon guidance, and acoustic propagation prediction), limiting the practical value of the inversion results in mission planning and decision support. Even in the few studies that mention uncertainty assessment, simplified assumptions are often used (such as uncertainty monotonically increasing with depth), failing to establish a dynamic uncertainty propagation mechanism linked to physical constraint strength.

[0076] 6) Lack of cross-medium physical field constraints leads to significant ambiguity in seabed boundary layer inversion. Existing technologies treat seawater sound velocity inversion and seabed sedimentary layer acoustic parameter inversion as two independent problems, failing to consider the physical continuity constraint of the solid-liquid interface. The seabed, as a crucial boundary for sound propagation, directly influences the sound velocity gradient inversion of the bottom seawater due to its acoustic properties (reflection coefficient, sound velocity dispersion characteristics). However, existing inversion methods rely solely on acoustic observations in seawater, leading to severe inversion ambiguity in near-seabed regions due to limited sound ray penetration—multiple combinations of bottom seawater sound velocity gradients can fit the same acoustic propagation observation data, resulting in multiple solutions for seabed boundary layer inversion. Currently, there is no effective solution for introducing seabed geophysical information (such as resistivity and porosity) to create a reverse constraint for seawater sound velocity inversion.

[0077] 7) Data assimilation frameworks struggle to balance high-dimensional state spaces with real-time computation requirements. The full-depth sound velocity field is a typical high-dimensional nonlinear state space. While traditional data assimilation methods such as Ensemble Kalman Filtering (EnKF) are theoretically mature, they face the "curse of dimensionality" in practical applications: if assimilation and updates are performed directly in a high-dimensional grid space, the computational cost increases exponentially with the number of water depth layers, making it difficult to meet the needs of real-time or near-real-time inversion. Existing methods either employ coarse-grid dimensionality reduction, sacrificing inversion accuracy, or abandon the assimilation framework altogether, relying solely on single inversion results, failing to achieve continuous fusion of observational data and the background field.

[0078] 8) The system architecture is fragmented and lacks integrated design. The existing technology system consists of independent components such as the observation platform, data preprocessing, inversion algorithm, and uncertainty assessment, lacking a unified data flow and control flow architecture. The observation platform cannot perceive the information requirements of the inversion algorithm, and the inversion algorithm cannot provide optimization instructions to the observation system; multiple physical field data (acoustic, electromagnetic, satellite remote sensing) have not been integrated on the same platform; and the uncertainty information of the inversion results is difficult to transmit to the subsequent application layer. This fragmented architecture not only reduces overall efficiency but also limits the system's adaptability in complex marine environments.

[0079] Reference Figure 1 This invention provides a sound velocity profile assimilation and inversion method based on dual physical data driving, the method comprising the following steps:

[0080] S100. Generate the planned path of the underwater unmanned swarm and perform adaptive observation and sparse representation of the sound propagation path to obtain the sparse representation result of the sound propagation path.

[0081] S110. Based on the historical background field or the inversion results of the previous moment, calculate the sensitivity distribution of the full-depth sound velocity field to each mode of the empirical orthogonal function.

[0082] S120. Based on the sensitivity distribution, select the high sensitivity distribution above the preset threshold and connect them to form the curve with the largest information entropy gradient, which is used as the planning path for the underwater unmanned swarm.

[0083] In this embodiment, the sensitivity distribution of the full-depth sound velocity field to each mode of the empirical orthogonal function (EOF) is first calculated based on the historical background field or the inversion results of the previous moment. Specifically, instead of calculating the conventional spatial variance, the system calculates the point where adding an observation sample to a certain spatial location would maximize the reduction in uncertainty of the first three EOF mode coefficients under the existing observation configuration. The system connects these highly sensitive points into a curve with the "maximum information entropy gradient," which serves as the planning path for the underwater unmanned swarm. This path is characterized by the fact that it may not be a straight line or a regular grid, but rather concentrated in key physical transition zones such as the top and bottom interfaces of the thermocline and the vicinity of the acoustic duct axis. Furthermore, the path density dynamically changes with water depth—sparse in the weak gradient region of the deep sea and dense in the main transition layer.

[0084] S130. Based on the underwater unmanned swarm planning path, while the swarm sails along the planned path, each node acts as a transmitter and receiver for each other, forming a sparse acoustic link network.

[0085] S140. Based on the sparse acoustic link network, the multipath arrival times of several links are extracted and sparsely represented by the compressed sensing method to obtain the sparse representation results of the sound propagation path.

[0086] In this embodiment, while the cluster navigates along the planned path, each node acts as a transmitter and receiver for the others, forming a sparse acoustic link network. The system extracts the multipath arrival time of each link and uses compressed sensing technology to perform sparse characterization of the sound propagation path, that is, to identify the sound rays most sensitive to the vertical distribution of the sound velocity field (such as the convergence zone sound rays captured by the sound duct axis and the seabed reflection sound rays). Based on the sparse characterization results, the system evaluates the "modal observability" of the current link set for the full-depth sound velocity field in real time. If the observability of a certain EOF mode (such as the mode corresponding to the deep-sea isotherm) is found to be below the threshold, the system immediately instructs a node in the cluster to change its course or depth and actively add a new type of sound propagation link sensitive to that mode (for example, diving below the sound duct axis to add seabed reflection paths), forming a real-time dynamic closed loop between observation and inversion algorithms. The innovation of this process lies in that it directly links the motion planning of the observation platform with the "information requirements" of the sound velocity mode within the inversion algorithm, realizing a leap from "data-driven inversion" to "inversion-driven observation".

[0087] S200: Collect marine environmental noise data and perform long-term noise cross-correlation calculation and passive time delay self-calibration to obtain the virtual sound source arrival time sequence after absolute time delay reference correction;

[0088] S210. Based on a deep-sea fixed vertical array, marine environmental noise data is collected and superimposed with long-term noise for cross-correlation calculation to obtain the arrival time of the virtual sound source.

[0089] In this embodiment, marine environmental noise data is continuously collected for more than 24 hours on a fixed vertical array in the deep sea (e.g., a 32-element hydrophone array deployed on the seabed). Cross-correlation calculations are performed on the noise signals received by any two hydrophones in the array, and these signals are superimposed for a sufficiently long time (e.g., 72 hours). Symmetrical arrival structures with clear physical meaning will appear on the positive and negative time axes of the cross-correlation function. These structures correspond to the propagation paths excited by the virtual sound source between the two hydrophone locations, effectively constructing a "virtual tomography" system without active transmission.

[0090] S220. Based on the seawater depth difference and seawater compressibility coefficient, the theoretical arrival time of the sound source is determined;

[0091] S230. Based on the deviation between the arrival time of the virtual sound source and the arrival time of the theoretical sound source, determine the absolute time drift and the correction value for the micro-perturbation of the array element position.

[0092] In this embodiment, the depth below 3500 meters in the deep sea is generally an isothermal layer, where the speed of sound is mainly controlled by hydrostatic pressure, exhibiting a stable trend of slowly increasing with depth. This invention utilizes this physical characteristic: in the deep layers of the vertical array (such as the lower 4000-5000 meter section), the theoretical sound propagation time between any two hydrophones can be accurately calculated using the known depth difference and seawater compressibility, unaffected by temperature and salinity variations. Therefore, by comparing the virtual arrival time obtained from actual cross-correlation with this theoretically calculated value, the absolute time drift of the entire observation system and the correction value for minor perturbations in array element positions can be obtained. This correction value is used as an "anchor point" to eliminate systematic errors introduced by clock asynchrony and array offset in all subsequent sound propagation measurements.

[0093] S240. Based on the absolute time drift and the correction value of the micro-perturbation of the array element position, the absolute time is calibrated, and the propagation time delay structure of the virtual sound source from the sea surface to the seabed is derived, thus obtaining the arrival time sequence of the virtual sound source after correction by the absolute time delay reference.

[0094] In this embodiment, after absolute time calibration, the propagation time delay structure of the virtual sound source from the sea surface to the seabed is derived by using the cross-correlation results of different depth combinations on the vertical array. This structure includes not only surface sound rays, but more importantly, deep-sea sound rays propagating along the acoustic duct axis and bottom-sea sound rays reflected from the seabed. All these passively extracted arrival delays are uniformly fed into the core inversion network as observation data independent of active sound sources. The innovation of this embodiment lies in the fact that it transforms the deep-sea background noise field into a self-calibrated and continuously usable "passive sound source" for the first time, and solves the long-standing problem of absolute time delay reference in passive acoustic tomography by utilizing the immutable physical laws of the deep-sea isotherm, enabling the system to have all-weather, all-ocean-depth, and zero-energy passive inversion capabilities.

[0095] S300, combining the sparse characterization results of the sound propagation path with the arrival time series of the virtual sound source, and performing constraint verification through a physical-data dual-drive dual-flow coupled generative inversion network, the verified full-depth sound velocity is obtained;

[0096] S310. Acquire sea surface satellite remote sensing data, and combine the sparse characterization results of sound propagation paths with the arrival time series of virtual sound sources after absolute time delay benchmark correction, perform data preprocessing, and construct a feature dataset.

[0097] In this embodiment, the input data includes three categories: The first category is from a sparse sound propagation arrival time set, specifically the multipath arrival time values ​​of several acoustic links (e.g., between nodes in an unmanned cluster or between a node and a fixed array). Each link includes depth and horizontal distance information for the transmitting and receiving nodes, and the data format is a time-series vector with a sampling interval typically in milliseconds. The second category is from passive delay self-calibration data, i.e., virtual sound source arrival time series corrected for absolute delay reference. This series covers multiple depth layers from the sea surface to the seabed, with each depth layer corresponding to a set of virtual arrival delays. The third category is sea surface satellite remote sensing data (sea surface temperature SST, sea surface height anomaly SSH) and historical Argo climatological profiles (depth-indexed temperature-salinity grid data with a spatial resolution of 1°×1° and vertically divided into standard layers of 0~2000 meters). Before being input into the network, all data needs to be preprocessed: the arrival time data is normalized to the 0~1 interval, the satellite remote sensing data is interpolated to the same geographical location as the target profile, and the Argo profile is matched to the background field of the target point through bilinear interpolation.

[0098] S320: Based on a multi-layer fully connected sensing mechanism, the modal flow is modeled, and the physical flow is constructed by combining a deep residual network and a physical information neural network. The modal flow and physical flow are integrated and a cross-attention mechanism is introduced to construct a dual-stream parallel network.

[0099] Specifically, the last layer of the modal flow is mapped to a query vector through a linear transformation. The feature vectors corresponding to the intermediate feature layers of the physical flow are extracted and used as key and value vectors. The dot product similarity between the query vector and the key vector is calculated and normalized using Softmax to obtain attention weights. The value vectors are then weighted and summed using the attention weights to obtain a fused feature vector, thus achieving the fusion of modal flow and physical flow. The feature vectors corresponding to the intermediate feature layers of the physical flow are mapped to key vectors through a linear transformation. The feature map of the current layer of the modal flow is used as the query vector, and the attention weights between the query vector and the key vector are calculated. The current feature map of the modal flow is weighted and corrected using the attention weights to obtain a corrected vector. This corrected vector is then added to the original features of the modal flow and propagated forward to achieve the fusion of physical flow and modal flow.

[0100] In this embodiment, the network adopts a dual-stream parallel structure. The modal stream consists of six fully connected multilayer perceptrons (MLPs), with each layer containing 256, 128, 64, 32, 16, and M neurons (M is the order of the empirical orthogonal function; in this invention, M=6). The input vector is a one-dimensional vector (dimension approximately 200-500, depending on the number of acoustic links) concatenated from the three types of preprocessed observation data mentioned above. The output consists of M EOF coefficients (typically ranging from -3 to 3). The loss function of the modal stream includes two terms: the first term is the root mean square error, which calculates the difference between the EOF coefficients output by the network and the true coefficients obtained from historical profile statistics; the second term is the modal orthogonality constraint term, which specifically calculates the inner product between the sound velocity profile reconstructed by the output coefficients and each order of EOF modes in the background field, forcing the off-diagonal elements to approach zero, thereby ensuring that the modes corresponding to the coefficients maintain orthogonality.

[0101] The physical flow employs a deep residual network (ResNet) combined with a physically-informed neural network (PINN). The input is a normalized depth value (mapping 0-6000 meters to 0-1), and the output is the speed of sound at that depth (in m / s). The main body of the physical flow consists of eight residual blocks, each containing two convolutional layers (kernel size 3) and skip connections. In the PINN part, two additional physical constraints are added to the loss function: the first is a smoothing constraint on the second derivative of the speed of sound with respect to depth. Specifically, this is achieved by automatically differentiating the network output to calculate the second derivative with respect to the input depth and applying a sparse regularization term to this derivative. Large abrupt changes in the second derivative are allowed at the upper and lower bounds of the thermocline (depth ranges pre-defined by historical climatological data), while at other depths, the second derivative is forced to approach zero, thus reflecting the gradient abrupt change characteristics of the thermocline. The second term is the implicit temperature-salinity constraint in the seawater state equation: while outputting the speed of sound, the network simultaneously outputs temperature T and salinity S (again, with depth as input) through an auxiliary branch. Then, it calculates the theoretical speed of sound using empirical formulas for seawater speed of sound (such as the Chen-Millero formula), forcing the network to minimize the difference between the output speed of sound and the theoretical speed of sound, thereby ensuring that the output speed of sound, temperature, and salinity meet physical consistency. The final output of the physical flow is a complete full-depth sound speed profile, that is, a sequence of sound speed values ​​at each discrete depth point (usually at 1-meter intervals or standard layers) from the sea surface to the seabed.

[0102] Furthermore, a cross-attention mechanism is introduced between specific layers of the modal flow and physical flow to achieve bidirectional information fusion. Specifically, the last layer of the modal flow (i.e., the feature layer before the output EOF coefficient, with a dimension of 32) and the intermediate feature layer of the physical flow (e.g., the feature map output after the 4th residual block, which contains depth and channel dimensions, with a dimension of depth points × 64) are selected. First, the fusion from the modal flow to the physical flow is performed. The 32-dimensional feature vector of the last layer of the modal flow is mapped to a query vector through a linear transformation. At the same time, the feature vectors of the intermediate feature layer of the physical flow at three key depth layers (thermocline upper boundary, acoustic duct axis depth, and seabed boundary layer) are extracted. Each key depth layer corresponds to a 64-dimensional feature vector, which serves as the key and value. By calculating the dot product similarity between the query and each key, the attention weight is obtained after Softmax normalization. Then, the value vector is weighted and summed using this weight to obtain a fused feature vector. The fused feature vector is added to the input of subsequent layers of the physical flow, thereby injecting global information of the modal space (i.e., the overall profile shape represented by the EOF coefficients) into the local depth features of the physical flow.

[0103] Conversely, the fusion from the physical flow to the modal flow involves concatenating the sound velocity gradient features (i.e., the derivative of sound velocity with respect to depth, calculated by automatic differentiation) output by the physical flow at the three key depth layers into a 192-dimensional feature vector, which is then mapped to a key vector through a linear transformation. Simultaneously, the features of the current layer of the modal flow are used as queries. Attention weights between the queries and keys are calculated, and then the current feature map of the modal flow is weighted and corrected to obtain a correction vector. This correction vector is added to the original features of the modal flow and then propagated forward. The purpose of this reverse correction mechanism is that if the EOF coefficient combination generated by the modal flow leads to physically unreasonable gradients at the key depth layers (e.g., incorrect gradient sign near the vocal tract axis), the gradient features extracted by the physical flow amplify this inconsistency through attention weights, thereby forcing the modal flow to adjust its coefficients. This bidirectional feature exchange is dynamically executed in each training iteration until the difference between the outputs of the two flows (e.g., the root mean square error between the profile reconstructed by the modal flow coefficients and the profile directly output by the physical flow) is less than a preset threshold.

[0104] S330. Input the feature dataset into the dual-stream parallel network, perform modal orthogonality constraints based on modal flow, output empirical orthogonal function coefficients, calculate sound velocity values ​​based on physical flow, and output full-depth sound velocity profile.

[0105] S340. Introduce a generative adversarial network discriminant and pre-train it to construct the trained generative adversarial network discriminant.

[0106] S350: Based on the trained generative adversarial network discriminator, the physical consistency between the empirical orthogonal function coefficients and the full-depth sound velocity profile is verified to obtain the verified full-depth sound velocity.

[0107] In this embodiment, after forward computation by the dual-flow coupled network, two sets of outputs are obtained: the EOF coefficients of the modal flow output and the full-depth sound velocity profile of the physical flow output. To further ensure the physical realism of the output, this embodiment introduces a lightweight generative adversarial network (GAN) discriminator. The discriminator is a five-layer convolutional neural network. The input is a complete sound velocity profile (length equal to the number of discrete points in the water depth, e.g., 1000 points), and the output is a scalar score representing the probability (between 0 and 1) that the profile belongs to real physical reality. The discriminator needs to be independently trained in advance: the training dataset contains two classes of samples. Positive samples are real sound velocity profiles obtained from a large number of measurements or high-precision numerical simulations (covering typical profiles of different sea areas and seasons), and negative samples are artificially constructed physically unreasonable profiles (e.g., artificially adding negative gradient oscillations in the thermocline, deleting the acoustic duct axis to make the sound velocity monotonically decrease with depth, introducing periodic fluctuations in the deep-sea isotherm, etc.). The discriminator is trained to convergence using binary cross-entropy loss, enabling it to distinguish between positive and negative samples with an accuracy of over 95%.

[0108] In the online inversion phase, the sound velocity profile output by the two-stream network is input into a pre-trained discriminator, which outputs a confidence score. If the score is higher than 0.8, the profile is considered physically reasonable and is directly output as the final result. If the score is lower than 0.8, a backpropagation is triggered: the discriminator's loss gradient is backpropagated to the two-stream network while keeping the discriminator's own weights frozen (only the parameters of the two-stream network are updated). This gradient signal forces the two-stream network to adjust its internal weights, making the profile generated in the next iteration closer to the discriminator's perception of the true profile. This process can be iterated multiple times (usually 3-5 times) until the discriminator score reaches a threshold or the number of iterations is exhausted. The final output sound velocity profile not only fits the observation data but also passes the physical consistency check, maintaining a physically reasonable gradient even in depth ranges lacking observation data (e.g., below the deep-sea acoustic duct axis to the seabed). The core innovation of this embodiment lies in the fact that by combining dual-stream coupling and generative adversarial verification, the global statistical laws of modal space and the local differential constraints of physical space are organically unified, which significantly improves the physical realism and reliability of the inversion results.

[0109] S400: Acquire broadband electromagnetic response time series, combine with verified full-depth sound velocity, perform cross-medium electromagnetic-acoustic joint inversion constraint, and obtain a strongly constrained calibrated full-depth sound velocity profile.

[0110] S410: Acquire broadband electromagnetic response time series and electric and magnetic field component data and perform Fourier transform to obtain frequency domain electromagnetic response;

[0111] S420. Using a three-dimensional electromagnetic inversion algorithm, resistivity structure inversion is performed on the frequency domain electromagnetic response to obtain three-dimensional resistivity grid data.

[0112] In this embodiment, the input data consists of broadband electromagnetic response time series collected by seabed electromagnetic receiving stations deployed in the target sea area (node ​​spacing is typically 2-5 km, covering the inversion area), combined with electric and magnetic field component data recorded at multiple transmit-receive offsets by a controlled-source electromagnetic transmission system (transmission frequency range of 0.01 Hz to 10 Hz, covering seabed exploration depths on a decades-scale). The acquisition duration is at least 30 minutes per transmission frequency to suppress environmental noise. The raw data is converted into frequency domain electromagnetic response using Fourier transform, including apparent resistivity and phase curves. The processing employs a three-dimensional electromagnetic inversion algorithm, specifically a Gauss-Newton inversion framework based on finite element or finite difference forward modeling. First, an initial resistivity model is constructed, treating the seawater layer as a known highly conductive homogeneous layer (resistivity approximately 0.3 Ω·m), and the area below the seabed as a horizontally layered or homogeneous half-space initial model. The inversion objective function includes a data fit difference term and a model smoothing constraint term. During the iteration process, the theoretical electromagnetic response at each measurement point and frequency is calculated in each forward modeling iteration, and the root mean square error is calculated between this calculation and the actual observed response. The gradient of the objective function with respect to the model parameters is calculated using the adjoint method to update the resistivity distribution. The iteration convergence condition is set to the change in the fitting difference between two adjacent iterations being less than 1% or reaching the preset maximum number of iterations (e.g., 20 iterations). The final output is three-dimensional resistivity grid data within a depth range of 500 meters below the seabed, with a horizontal grid spacing of 100-200 meters and vertical stratification of 10-20 layers from the seabed surface downwards (the top layer is 2 meters thick, gradually increasing in thickness downwards). This resistivity distribution visually reflects the spatial variation of sediment porosity, saturation, and fluid properties.

[0113] S430. Based on three-dimensional resistivity grid data and combined with sediment grain size parameters, sound velocity conversion was performed through a rock physics model to obtain a sound velocity dispersion curve table and a reflection coefficient matrix.

[0114] Specifically, based on three-dimensional resistivity grid data and sediment grain size auxiliary parameters obtained through historical geological sampling, an improved Biot-Stoll two-phase medium model framework is constructed. Using pore water resistivity, sediment porosity is calculated cell-by-cell using the Archie formula. Based on sediment porosity, the solid skeleton bulk modulus and shear modulus obtained from tables, as well as the pore fluid bulk modulus calculated from temperature, salinity, and pressure, are input into the improved Biot-Stoll two-phase medium model framework. Considering the Biot dissipation mechanism, the P-wave velocity dispersion curve, S-wave velocity, and attenuation coefficient are calculated. Based on the P-wave velocity dispersion curve, S-wave velocity, and attenuation coefficient, the modulus and phase of the plane wave reflection coefficient at different grazing angles are calculated using a multilayer medium reflection coefficient algorithm, resulting in a sound velocity dispersion curve table and a reflection coefficient matrix.

[0115] In this embodiment, the input is three-dimensional resistivity grid data and auxiliary sediment grain size parameters (obtainable from historical geological sampling or regional geological maps, including average grain size, sorting factor, clay content, etc.). The processing adopts an improved Biot-Stoll equivalent medium model suitable for deep-sea sediments. This model treats sediments as a two-phase medium composed of a solid skeleton and pore fluid, with seawater as the pore fluid. The core relationship of the model is: resistivity and porosity are linked through the Archie formula, that is, the ratio of resistivity to pore water resistivity is equal to the negative power of porosity (the exponent m is called the cementation index, usually taken as 1.8~2.2); porosity and P-wave velocity are linked through the Wood formula or Gassmann equation, while considering the frequency-dependent Biot dissipation mechanism. In specific implementation, porosity is first calculated cell by cell based on the resistivity grid: porosity is inversely calculated using known pore water resistivity (which can be estimated through in-situ sampling or seawater conductivity) and the Archie formula. Then, the porosity, bulk modulus of the solid skeleton, shear modulus (obtained from tables based on mineral composition and grain size parameters), and bulk modulus of the pore fluid (calculated from temperature, salinity, and pressure) are input into the Biot-Stoll model. Within a preset frequency range (10 Hz ~ 10 kHz, covering the underwater acoustic application frequency band), the P-wave velocity dispersion curve and S-wave velocity (considering only shear waves; these may be zero for unconsolidated sediments) are calculated. Simultaneously, a multilayer media reflection coefficient algorithm (such as Kennett's reflectivity method) is used to calculate the plane wave reflection coefficient modulus and phase at different grazing angles (0° ~ 90°, step size 1°). The final output consists of two datasets at the seabed interface: one is a table of sound velocity dispersion curves, formatted as frequency-P-wave velocity-S-wave velocity-attenuation coefficient; the other is a reflection coefficient matrix, formatted as frequency-grazing angle-reflection coefficient modulus. These data characterize the intrinsic properties of the seabed as an acoustic boundary.

[0116] S440. Based on the verified full-depth sound velocity, combined with the seawater density and bottom seawater sound velocity gradient obtained from the temperature-salt-pressure state equation, the acoustic impedance of the seawater side of the seabed interface is calculated.

[0117] S450. Obtain the sound velocity and density of the seabed sediment layer to determine the acoustic impedance on the sediment layer side.

[0118] S460. Obtain the ratio of the acoustic impedance on the seawater side to the acoustic impedance on the sedimentary layer side of the seabed interface, calculate the reflection coefficient at vertical incidence, compare it with the reflection coefficient matrix, and calculate the relative deviation signal.

[0119] S470. Recalculate the acoustic impedance and reflection coefficient by perturbation gradient value, obtain the sensitivity of the bottom seawater sound velocity gradient to the reflection coefficient deviation, and obtain the Jacobian matrix of the gradient change on the reflection coefficient.

[0120] S480. The relative deviation signal is converted into the corrected gradient of the hidden layer weights related to the bottom seawater in the core inversion network. Combined with the Jacobian matrix, it is backpropagated to the inversion network through the backpropagation algorithm for iterative correction, and the full-depth sonic velocity profile after strong constraint calibration is obtained.

[0121] In this embodiment, the input consists of two parts: first, the sound velocity dispersion curve and reflection coefficient matrix of the seabed interface, serving as strong physical anchors; second, the full-depth sound velocity profile output by the core inversion network, with particular attention to the sound velocity gradient distribution of the bottom seawater (within 10 meters of depth above the seabed interface). The processing first performs a consistency check: based on the bottom seawater sound velocity gradient output by the core inversion network, combined with the seawater density (calculated from temperature, salinity, and pressure using the equation of state), the acoustic impedance (sound velocity multiplied by density) on the seawater side of the seabed interface is calculated. Simultaneously, the acoustic impedance on the sedimentary layer side is calculated from the sound velocity and sedimentary layer density (which can be calculated from porosity and skeleton density). The ratio of the impedances on both sides determines the reflection coefficient (the reflection coefficient at vertical incidence equals the impedance difference divided by the sum of impedances). The calculated vertical reflection coefficient is compared with the reflection coefficient modulus at the corresponding frequency to calculate the relative deviation. If the deviation exceeds a preset threshold, a feedback correction mechanism is triggered. The specific process of feedback correction is as follows: The system calculates the sensitivity of the bottom seawater sound velocity gradient to the reflection coefficient deviation (by perturbing the gradient value, recalculating the impedance and reflection coefficient, and obtaining the Jacobian matrix of the gradient change on the reflection coefficient), and then converts the deviation signal into a correction gradient for the hidden layer weights related to the bottom seawater in the core inversion network. This gradient is backpropagated to the dual-flow network architecture through the backpropagation algorithm, forcing the network to adjust the sound velocity output of the bottom seawater in the next iteration, so that the reflection coefficient calculated by the new gradient gradually approaches the reflection coefficient determined by the electromagnetic inversion method. The correction process adopts a gradual update, adjusting only 0.1~1% of the network weights in each iteration to avoid disrupting the continuity of the entire profile. Iteration continues until the reflection coefficient deviation is less than the threshold, or the maximum number of corrections is reached (e.g., 10 times). The final output is a full-depth sound velocity profile after strong constraint calibration at the seabed interface. Its sound velocity gradient at the seabed boundary layer and the seabed acoustic characteristics derived from the electromagnetic-rock physics model satisfy the acoustic impedance matching constraint, thereby eliminating the inversion ambiguity of relying solely on acoustic observations of the bottom seawater. The core innovation of this embodiment lies in transforming cross-medium electromagnetic detection data into quantifiable acoustic boundary conditions through a rock physics model, and establishing a closed-loop feedback from boundary conditions to inversion network weights, thereby achieving physical consistency coupling between seawater sound velocity inversion and seabed sediment characteristics.

[0122] S500 performs ensemble Kalman filtering assimilation and uncertainty propagation on the full-depth sound velocity profile after strong constraint calibration, based on POD order reduction, to obtain the full-depth sound velocity mean profile and the full-depth sound velocity standard deviation profile.

[0123] S510. Arrange the verified full-depth sound velocity in the same depth order, construct a matrix and perform singular value decomposition to obtain the left singular vector, singular value and right singular vector. The left singular vector is regarded as the POD mode basis.

[0124] S520. Based on the cumulative contribution rate of singular values, select and retain the pre-preset order POD mode basis, and construct the mode basis matrix;

[0125] In this embodiment, the input data is a set of full-depth sonic velocity profiles accumulated historically or output by a two-stream coupled generative inversion network. This set contains multiple profile samples, typically between one hundred and five hundred, each covering the entire depth from the sea surface to the seabed in the vertical direction. The number of discrete depth points depends on the vertical resolution (e.g., at one-meter intervals, there would be six thousand depth points at a depth of six thousand meters). These profiles are arranged in the same depth order to form a large matrix, where rows correspond to depth layers and columns correspond to different profile samples.

[0126] The processing begins with singular value decomposition (SVD) of the matrix, a matrix factorization technique that decomposes the original data into three parts: left singular vectors, singular values, and right singular vectors. The left singular vectors represent the main spatial patterns of sound velocity profile variation. Each left singular vector is a sequence of lengths equal to the number of depth points, called a POD mode basis, which describes the typical shape of sound velocity variation with depth. For example, the first mode might represent the overall shift, the second mode the thermocline intensity variation, and the third mode the vertical movement of the deep-sea acoustic duct axis. The singular values ​​reflect the contribution of each mode to the overall data variance; the larger the singular value, the more important the mode. Subsequently, the system determines the number of modes to retain based on the cumulative contribution rate of the singular values, typically retaining the first five to ten modes, ensuring these modes can explain more than 99.5% of the overall variance. The contribution of truncated modes is very small and can be ignored. After this reduction process, the original high-dimensional sound velocity profile (e.g., 6,000-dimensional) can be approximated by a low-dimensional coefficient vector (only five to ten dimensions). Specifically, the preserved POD mode basis is arranged in columns to form a mode basis matrix. Any sound velocity profile is equal to the product of this matrix and a coefficient vector. Each element of the coefficient vector represents the weight of the corresponding mode in the profile.

[0127] The output of this embodiment includes a POD modal basis matrix, the size of which is the number of depth points multiplied by the order of the preserved modes, with each column being a modal basis vector; and a coefficient matrix of the historical profile set in the POD space, the size of which is the order of the preserved modes multiplied by the number of profile samples, with each element corresponding to a specific modal coefficient of a specific profile. Singular values ​​corresponding to each mode are also preserved for subsequent uncertainty scaling. This order reduction process drastically reduces the computational complexity of subsequent Kalman filtering operations from being proportional to the cube of the number of depth points to being only proportional to the cube of the order of the preserved modes, enabling real-time computation.

[0128] S530. Acquire new observation data, combine it with the modal basis matrix, and perform Kalman filtering assimilation in low-dimensional space to obtain a set of sound velocity profiles of the analysis field.

[0129] Specifically, new observation data is acquired, including sparse field profile measurement data and multipath arrival times of several links; based on the low-dimensional POD coefficient space, mean coefficient vectors are calculated according to the modal basis matrix, and the deviation of each mean coefficient vector relative to the mean is calculated to construct the background field covariance matrix; corresponding observation operators are constructed based on the new observation data, and the full-depth sound velocity profile is reconstructed in combination with the modal basis matrix, and the sound velocity values ​​at the observation depth are extracted to calculate the theoretical arrival time of sound waves; based on the background field covariance matrix, the theoretical arrival time of sound waves, and the error covariance attached to the observation data, the Kalman gain matrix is ​​calculated to determine the correction weight of the observation data to the background field; based on the Kalman gain matrix and the correction weight, the POD coefficient update calculation is performed, multiplied by the POD modal basis, and projected back to the original high-dimensional physical space to obtain the set of sound velocity profiles of the analysis field.

[0130] In this embodiment, the input includes two parts: first, a set of background field sound velocity profiles obtained at the previous moment or the previous assimilation period, which has been stored in the form of POD coefficients in step S100, i.e., each set member corresponds to a set of low-dimensional coefficient vectors; second, newly provided observation data. These observations can be sparse field profile measurements, such as temperature and salinity data covering only a portion of the depth layer obtained by a CTD instrument at a certain point, or arrival time data of the sound propagation link, including multiple arrival times. The observation data includes measurement error information, usually expressed as the variance of each observation, reflecting the instrument accuracy and noise level.

[0131] The processing is carried out in a low-dimensional POD coefficient space, and the specific steps are as follows: First, all POD coefficient vectors of the background field set are aggregated, their average value is calculated, and a mean coefficient vector is obtained; simultaneously, the deviation of each coefficient relative to the mean is calculated, forming a perturbation matrix. Then, based on these deviations, the background field covariance matrix is ​​calculated. This matrix is ​​only the product of the preserved mode order, which is very small. Next, a mapping relationship from POD coefficients to observations needs to be established, called the observation operator. For CTD profile observations, the observation operator's role is to first reconstruct the full-depth sound velocity profile by multiplying the POD coefficients by the modal basis matrix, and then extract the sound velocity values ​​at the observation depth. For sound propagation time observations, the observation operator needs to call the ray tracing model to calculate the theoretical arrival time of the sound wave from the emission point to the receiving point based on the reconstructed sound velocity profile. Since the preserved mode order is very small, the computational efficiency of this mapping is very high. Then, the system calculates the Kalman gain matrix based on the background field covariance, the linearized approximation of the observation operator, and the observation error covariance. This matrix determines the correction weights of the observation data to the background field during the update process. To maintain the statistical consistency of the set, the system adds a random perturbation conforming to the measurement error distribution to each observation, generating a set of perturbed observation vectors, with each set member corresponding to a perturbed observation. Finally, for each set member, its POD coefficients are updated using the Kalman gain matrix and the corresponding perturbed observation: the updated coefficients are equal to the background coefficients plus the gain matrix multiplied by the observation residual (i.e., the perturbed observation minus the theoretical observation calculated from the background coefficients using the observation operator). After all members are updated, the set of POD coefficients for the analysis field is obtained. Multiplying these coefficients by the POD modal basis matrix allows projection back into the original high-dimensional physical space, yielding a new set of full-depth sound velocity profiles.

[0132] The output of this embodiment is a set of sound velocity profiles of the analysis field, containing the same number of profile samples as the background field, with each profile covering the entire depth layer; at the same time, the corresponding set of POD coefficients of the analysis field is also output for use in the next assimilation cycle.

[0133] S540. By combining the set of sound velocity profiles in the analysis field with the full-depth sound velocity profile after strong constraint calibration, the uncertainty of water depth adaptive quantification is performed to obtain the mean profile and the standard deviation profile of the full-depth sound velocity.

[0134] In this embodiment, the input is the set of sound velocity profiles of the analysis field output in step S200, and the pre-calculated physical constraint strength indices for each depth layer. The physical constraint strength indices include three categories: the first is the sound ray penetration density, obtained from the statistics of the sound propagation links, indicating how many sound rays pass through or are reflected at a certain depth layer. After normalization, the value is between zero and one; the higher the density, the stronger the acoustic constraint on that layer. The second is the electromagnetic constraint confidence, derived from the inversion results, reflecting the ability of electromagnetic data near the seabed interface to constrain sound velocity, also normalized to the zero-to-one range. The third is the indicator function near the sound channel axis depth, located by the background field, and a Gaussian window centered on the sound channel axis is set. The closer to the sound channel axis, the closer the index is to one; the farther away, the closer it is to zero.

[0135] The processing first extracts the sound velocity values ​​from all profiles at each depth layer, forming a numerical set. The average of this set is calculated as the most probable sound velocity value for that depth layer. The standard deviation of the set is then calculated to reflect the dispersion of sound velocity at that layer. Assuming that the sound velocity follows a normal distribution, a 95% confidence interval is obtained. However, the directly calculated standard deviation does not reflect the additional dispersion modulation caused by the varying strengths of physical constraints at different depth layers. Therefore, the system introduces a depth-adaptive set expansion mechanism. Specifically, for each depth layer, the three types of physical constraint strength indices are weighted and summed to obtain a comprehensive physical constraint strength factor. The weights can be adjusted according to the actual application scenario; for example, the weight for ray penetration density is 0.5, the weight for electromagnetic constraint confidence is 0.3, and the weight for the acoustic tract axis indicator function is 0.2. Then, this strength factor is used to modulate the original standard deviation: the modulated standard deviation is equal to the original standard deviation multiplied by one, minus a modulation strength coefficient, and then multiplied by the comprehensive strength factor. The modulation intensity factor is an adjustable parameter ranging from 0.5 to 0.8. Its physical meaning is that the stronger the physical constraints, the more the ensemble dispersion is compressed. When the overall intensity factor is close to 1, the modulated standard deviation can be compressed to 20% to 50% of the original standard deviation. When the overall intensity factor is close to zero (e.g., in the deep-sea mesosphere with no observations and no physical constraints), the modulated standard deviation retains its original value, i.e., it preserves the natural variability of the background field. It is important to note that the modulated standard deviation should not fall below a minimum threshold to avoid unrealistic accuracy due to overconfidence.

[0136] After modulation, the system recalculates the confidence interval for each depth layer and outputs a non-uniform uncertainty field that varies with water depth. Furthermore, the system can also output a probability density function histogram or fitted normal distribution parameters for the sound velocity values ​​at each depth layer for downstream applications. The final output of this embodiment is a data structure containing: a full-depth sound velocity mean profile, i.e., the most probable sound velocity value for each depth layer; a full-depth sound velocity standard deviation profile, i.e., the standard deviation sequence after depth-adaptive modulation; and the upper and lower bounds of the 95% confidence interval for each depth layer. This uncertainty output reflects both the sparsity of the observation data (fewer observations result in greater dispersion of the original set) and the physical anchoring capability at different depth layers (stronger physical constraints compress more uncertainty), thus providing reliable error envelope information for underwater equipment applications such as path planning, target detection, and positioning.

[0137] Reference Figure 2 A sound velocity profile assimilation and inversion system based on dual physical data driving, comprising:

[0138] The first module 201 is used to generate the planned path of the underwater unmanned swarm and perform adaptive observation and sparse representation of the sound propagation path to obtain the sparse representation result of the sound propagation path.

[0139] The second module 202 is used to collect marine environmental noise data and perform long-term noise cross-correlation calculation and passive time delay self-calibration to obtain the virtual sound source arrival time sequence after absolute time delay reference correction.

[0140] The third module 203 is used to combine the sparse characterization results of the sound propagation path with the arrival time series of the virtual sound source, and perform constraint verification through a physical-data dual-drive dual-stream coupled generative inversion network to obtain the verified full-depth sound velocity.

[0141] The fourth module 204 is used to acquire broadband electromagnetic response time series, and combined with the verified full-depth sound velocity, to perform cross-medium electromagnetic-acoustic joint inversion constraint to obtain a strongly constrained and calibrated full-depth sound velocity profile.

[0142] Module 5, 205, is used to perform ensemble Kalman filtering assimilation and uncertainty propagation on the full-depth sound velocity profile after strong constraint calibration, based on POD order reduction, to obtain the full-depth sound velocity mean profile and the full-depth sound velocity standard deviation profile.

[0143] The content of the above method embodiments is applicable to this system embodiment. The specific functions implemented in this system embodiment are the same as those in the above method embodiments, and the beneficial effects achieved are also the same as those achieved in the above method embodiments.

[0144] The above is a detailed description of the preferred embodiments of the present invention. However, the present invention is not limited to the embodiments described. Those skilled in the art can make various equivalent modifications or substitutions without departing from the spirit of the present invention. All such equivalent modifications or substitutions are included within the scope defined by the claims of this application.

Claims

1. A sound velocity profile assimilation and inversion method based on dual-driven physical data, characterized in that, Includes the following steps: Based on the historical background field or the inversion results of the previous moment, the sensitivity distribution of the full-depth sound velocity field to each mode of the empirical orthogonal function is calculated. Based on the sensitivity distribution, the curve with the largest information entropy gradient is selected by locating high sensitivity distributions above a preset threshold, and used as the planning path for the underwater unmanned swarm. Based on the planned path of the underwater unmanned swarm, while the swarm is sailing along the planned path, each node acts as a transmitter and receiver for each other, forming a sparse acoustic link network. Based on a sparse acoustic link network, the multipath arrival times of several links are extracted and sparsely represented by compressed sensing method to obtain the sparse representation results of the sound propagation path. Based on a deep-sea fixed vertical array, marine environmental noise data is collected and superimposed with long-term noise for cross-correlation calculation to obtain the arrival time of virtual sound sources. The theoretical arrival time of the sound source is determined based on the depth difference of seawater and the compressibility coefficient of seawater; Based on the deviation between the arrival time of the virtual sound source and the arrival time of the theoretical sound source, the absolute time drift and the correction value for the micro-perturbation of the array element position are determined. Absolute time calibration is performed based on the absolute time drift and the correction value of the micro-perturbation of the array element position. The propagation time delay structure of the virtual sound source from the sea surface to the seabed is derived, and the arrival time sequence of the virtual sound source after absolute time delay reference correction is obtained. We acquire sea surface satellite remote sensing data, and combine the sparse characterization results of sound propagation paths with the arrival time series of virtual sound sources after absolute time delay benchmark correction to perform data preprocessing and construct a feature dataset. Based on a multi-layer fully connected sensing mechanism, the modal flow is modeled, and the physical flow is constructed by combining a deep residual network and a physical information neural network. The modal flow and physical flow are integrated and a cross-attention mechanism is introduced to construct a dual-stream parallel network. The method of integrating modal streams and physical streams and introducing a cross-attention mechanism to construct a dual-stream parallel network also includes: The last layer of the modal flow is mapped to a query vector through a linear transformation, and the feature vectors corresponding to the intermediate feature layers of the physical flow are extracted as key vectors and value vectors. By calculating the dot product similarity between the query vector and the key vector, and then normalizing it with Softmax to obtain the attention weight, the value vector is then weighted and summed using the attention weight to obtain the fused feature vector, thus realizing the fusion of modal flow to physical flow. The feature vectors corresponding to the intermediate feature layers of the physical flow are mapped to key vectors through a linear transformation. The feature map of the current layer of the modal flow is used as the query vector, and the attention weights between the query vector and the key vector are calculated. The current feature map of the modal flow is weighted and corrected by attention weights to obtain a correction vector, which is then added to the original features of the modal flow and propagated forward to achieve the fusion of physical flow to modal flow. The feature dataset is input into a dual-stream parallel network. Modal orthogonality constraints are performed based on modal flow, and empirical orthogonal function coefficients are output. Sound velocity values ​​are calculated based on physical flow, and a full-depth sound velocity profile is output. A generative adversarial network (GAN) discriminant is introduced and pre-trained to construct a trained GAN discriminant. The physical consistency between the empirical orthogonal function coefficients and the full-depth sound velocity profile is verified based on the trained generative adversarial network discriminator, and the verified full-depth sound velocity is obtained. Broadband electromagnetic response time series and electric and magnetic field component data are collected and Fourier transformed to obtain the frequency domain electromagnetic response; The resistivity structure of the frequency domain electromagnetic response is inverted using a three-dimensional electromagnetic inversion algorithm to obtain three-dimensional resistivity grid data. Based on three-dimensional resistivity grid data and combined with sediment grain size parameters, sound velocity conversion was performed through a rock physics model to obtain a sound velocity dispersion curve table and a reflection coefficient matrix. Based on the verified full-depth sound velocity, combined with the seawater density and bottom seawater sound velocity gradient obtained from the temperature-salt-pressure state equation, the acoustic impedance of the seawater side of the seabed interface is calculated. Obtain the sound velocity and density of the seabed sediment layer to determine the acoustic impedance on the sediment layer side; The ratio of the acoustic impedance on the seawater side to the acoustic impedance on the sedimentary layer side of the seabed interface is obtained to calculate the reflection coefficient at vertical incidence. The result is compared with the reflection coefficient matrix and the relative deviation signal is calculated. By recalculating the acoustic impedance and reflection coefficient using the perturbation gradient value, the sensitivity of the bottom seawater sound velocity gradient to the reflection coefficient deviation is obtained, and the Jacobian matrix of the gradient change on the reflection coefficient is obtained. The relative deviation signal is converted into the corrected gradient of the hidden layer weights related to the bottom seawater in the core inversion network. Combined with the Jacobian matrix, it is backed to the inversion network through the backpropagation algorithm for iterative correction, and the full-depth sound velocity profile after strong constraint calibration is obtained. The verified full-depth sound velocity is arranged in the same depth order, a matrix is ​​constructed and singular value decomposition is performed to obtain the left singular vector, singular value and right singular vector. The left singular vector is regarded as the POD mode basis. Based on the cumulative contribution rate of singular values, the pre-preset order POD mode basis is selected and retained, and a mode basis matrix is ​​constructed. New observation data are acquired and combined with the modal basis matrix. Kalman filtering is then performed in a low-dimensional space to assimilate the data and obtain a set of sound velocity profiles for the analysis field. By combining the set of sound velocity profiles in the analysis field with the full-depth sound velocity profile after strong constraint calibration, the uncertainty of water depth adaptive quantification is performed to obtain the mean sound velocity profile and the standard deviation sound velocity profile of the full depth.

2. The sound velocity profile assimilation and inversion method based on dual physical data driving according to claim 1, characterized in that, The step of converting sound velocity based on three-dimensional resistivity grid data, combined with sediment grain size parameters, and using a rock physics model to obtain a sound velocity dispersion curve table and a reflection coefficient matrix specifically includes: Based on three-dimensional resistivity grid data and combined with sediment grain size auxiliary parameters obtained through historical geological sampling, an improved Biot-Stoll two-phase medium model framework was constructed. Combined with pore water resistivity, sediment porosity was obtained by back-calculating grid cell by grid cell using the Archie formula. Based on sediment porosity, combined with the solid skeleton bulk modulus and shear modulus obtained from the table and the pore fluid bulk modulus calculated from temperature, salinity and pressure, the data are input into the improved Biot-Stoll two-phase medium model framework. Considering the Biot dissipation mechanism, the longitudinal wave velocity dispersion curve, transverse wave velocity and attenuation coefficient are calculated. Based on the longitudinal wave velocity dispersion curve, transverse wave velocity and attenuation coefficient, the modulus and phase of the plane wave reflection coefficient at different grazing angles are calculated using a multilayer medium reflection coefficient algorithm, resulting in a velocity dispersion curve table and a reflection coefficient matrix.

3. The sound velocity profile assimilation and inversion method based on dual physical data driving according to claim 2, characterized in that, The step of acquiring new observation data, combining it with the modal basis matrix, and performing Kalman filtering assimilation in low-dimensional space to obtain a set of sound velocity profiles for the analysis field specifically includes: Acquire new observation data, which includes sparse field profile measurement data and multiple-path arrival times of several links; Based on the low-dimensional POD coefficient space, the mean coefficient vector is calculated according to the modal basis matrix, and the deviation of each mean coefficient vector from the mean is calculated to construct the background field covariance matrix. Based on the new observation data, the corresponding observation operator is constructed, and the full-depth sound velocity profile is reconstructed by combining it with the modal basis matrix. Then, the sound velocity value at the observation depth is extracted, and the theoretical arrival time of the sound wave is calculated. Based on the background field covariance matrix, the theoretical arrival time of sound waves, and the error covariance of the observation data, the Kalman gain matrix is ​​calculated to determine the correction weight of the observation data on the background field. Based on the Kalman gain matrix and the corrected weights, the POD coefficients are updated and calculated, multiplied by the POD mode basis, and projected back to the original high-dimensional physical space to obtain the set of sound velocity profiles of the analysis field.

4. A system for a sound velocity profile assimilation and inversion method based on dual physical data as described in claim 1, characterized in that, Includes the following modules: The first module is used to generate the planned path of the underwater unmanned swarm and perform adaptive observation and sparse representation of the sound propagation path to obtain the sparse representation result of the sound propagation path. The second module is used to collect marine environmental noise data and perform long-term noise cross-correlation calculation and passive time delay self-calibration to obtain the virtual sound source arrival time sequence after absolute time delay reference correction. The third module is used to combine the sparse characterization results of the sound propagation path with the arrival time series of the virtual sound source, and perform constraint verification through a physical-data dual-drive dual-flow coupled generative inversion network to obtain the verified full-depth sound velocity. The fourth module is used to acquire broadband electromagnetic response time series, and combined with the verified full-depth sound velocity, to perform cross-medium electromagnetic-acoustic joint inversion constraints, and obtain a strongly constrained calibrated full-depth sound velocity profile. The fifth module is used to perform ensemble Kalman filtering assimilation and uncertainty propagation on the full-depth sound velocity profile after strong constraint calibration, based on POD order reduction, to obtain the full-depth sound velocity mean profile and the full-depth sound velocity standard deviation profile.

Images