Submarine cable low-frequency vibration signal identification and monitoring method and system

Abstract

The application provides a submarine optical cable low-frequency vibration signal identification and monitoring method and system, relates to the technical field of underwater safety monitoring, and comprises the following steps: acquiring an acoustic wave vibration signal of a submarine optical cable; extracting a low-frequency component and converting the low-frequency component into strain time series data; analyzing a time difference and amplitude attenuation to calculate a medium sound velocity distribution; constructing a propagation path equation to determine a vibration source three-dimensional coordinate and a motion vector field; and predicting a trajectory envelope area and an approaching distance of a sensitive area to determine early warning. The application uses the optical cable as a distributed sensing network to realize accurate positioning and early warning of underwater threat targets and improve the safety supervision efficiency of sea areas.

Description

Technical Field

[0001] This invention relates to the field of underwater safety monitoring technology, and in particular to a method and system for identifying and monitoring low-frequency vibration signals of submarine optical cables. Background Technology

[0002] Submarine optical cables are a vital infrastructure for global communication networks, carrying a large number of international data transmission tasks. With the increasing development of marine resources and the growing number of seabed activities, the safety of submarine optical cables faces many threats, including human interference, natural disasters, and marine biological activities. To ensure the normal operation of submarine optical cables, real-time monitoring and identification of low-frequency vibration signals around the cables are of great significance for timely detection of potential threats and prevention of cable damage.

[0003] Traditional submarine fiber optic cable monitoring technologies primarily rely on periodic shipborne inspections or passive sonar systems at fixed monitoring stations. These methods typically only provide post-incident detection of damage and cannot offer real-time early warning of potential threats approaching the cable. With the development of distributed fiber optic sensing technology, it has become possible to utilize the cable itself as a vibration sensor, providing a new technological means for the protection of submarine fiber optic cables. Summary of the Invention

[0004] This invention provides a method and system for identifying and monitoring low-frequency vibration signals of submarine optical cables, which can solve the problems in the prior art.

[0005] A first aspect of the present invention provides a method for identifying and monitoring low-frequency vibration signals of submarine optical cables, comprising:

[0006] Acquire distributed acoustic vibration signals of submarine optical cables in the sea area to be monitored;

[0007] The acoustic vibration signal is decomposed in the frequency domain to extract the low-frequency components within the target frequency band, and the low-frequency components are converted into strain time series data based on the fiber strain response characteristics.

[0008] By analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, the medium sound velocity distribution between the vibration source and the optical cable is calculated by inversion, and the refraction path and number of reflections of the vibration wave at the seabed interface are derived based on the medium sound velocity distribution.

[0009] Based on the refraction path and the number of reflections, a propagation path equation system is constructed. By solving the least squares optimization solution of the propagation delay in the multi-path propagation equation system, the three-dimensional spatial coordinates of the vibration source are determined, and the motion vector field of the vibration source is fitted using the sequence of three-dimensional spatial coordinates at consecutive time moments.

[0010] Based on the velocity vector and heading rate of change of the motion vector field, the trajectory envelope region of the vibration source within the future time window is predicted by extrapolation of the kinematic equations, and the shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area is calculated to determine whether an early warning response is triggered.

[0011] The acoustic vibration signal is decomposed in the frequency domain to extract low-frequency components within the target frequency band, and the low-frequency components are converted into strain time-series data based on the fiber optic strain response characteristics, including:

[0012] The acoustic vibration signal is subjected to time-frequency transformation processing to obtain a three-dimensional distribution matrix of frequency time amplitude. The three-dimensional distribution matrix is ​​then filtered in the frequency domain according to a preset frequency cutoff boundary to extract the frequency domain energy components within the target frequency band as low-frequency components.

[0013] A strain transfer function model for optical fiber is established, which includes a frequency-dependent strain sensitivity coefficient and a phase delay coefficient. The strain sensitivity coefficient characterizes the unit strain caused by sound waves of different frequencies in the optical fiber, and the phase delay coefficient characterizes the phase lag of sound waves of different frequencies during propagation in the optical fiber.

[0014] Based on the frequency distribution of the low-frequency component, the strain sensitivity coefficient and the phase delay coefficient at the corresponding frequency in the fiber strain transfer function model are queried. The amplitude of the low-frequency component is compensated for using the strain sensitivity coefficient, and the phase of the low-frequency component is corrected for time delay using the phase delay coefficient, so as to obtain the compensated and corrected strain equivalent signal.

[0015] The strain equivalent signal is subjected to inverse time-frequency transformation to obtain the strain amplitude sequence of each spatial position of the optical fiber at continuous time, and arranged in chronological order to form strain time series data.

[0016] By analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, the medium sound velocity distribution between the vibration source and the optical cable is calculated by inversion. Based on the medium sound velocity distribution, the refraction path and number of reflections of the vibration wave at the seabed interface are derived, including:

[0017] The strain time series data is subjected to wave field separation processing. Based on the timing relationship of waveform arrival and the waveform energy concentration characteristics, the direct wave component, seabed reflected wave component and sea surface reflected wave component are separated and extracted from the aliased signal, and the arrival time difference of each wave component at different spatial locations of the optical cable is calculated.

[0018] Based on the arrival time difference of the direct wave component and the spatial coordinates of the optical cable, a direct wave time delay equation is constructed. Based on the arrival time difference of the seabed reflected wave component and the geometric relationship of the reflection path, a reflected wave time delay equation is constructed. Based on the amplitude attenuation gradient of each wave field component, a medium attenuation constraint equation is established.

[0019] By jointly solving the direct wave time delay equation, the reflected wave time delay equation, and the medium attenuation constraint equation, the sound velocity field of the seawater layer and the sound velocity field of the seabed layer are obtained by inversion. The position of the seawater-seabed interface is determined by using the numerical difference between the attenuation coefficient of the seawater layer and the attenuation coefficient of the seabed layer obtained by inversion.

[0020] Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, the vibration wave propagation trajectory is calculated using a ray tracing algorithm. The number of intersections between the vibration wave propagation trajectory and the seawater-seabed interface is counted as the number of reflections. The incident vector and refraction vector of the ray trajectory at the interface are calculated to determine the refraction path.

[0021] Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, a ray tracing algorithm is used to calculate the propagation trajectory of the vibration wave. The number of intersections between the vibration wave propagation trajectory and the seawater-seabed interface is counted as the number of reflections. The incident vector and refraction vector of the ray trajectory at the interface are calculated to determine the refraction path, including:

[0022] Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, the travel time field distribution from the vibration source to various points in space is calculated using the fast travel method. The travel time field represents the shortest time required for the vibration wave to propagate from the vibration source to any spatial location. Spatial points with equal values ​​in the travel time field are extracted to form travel time isosurfaces, which represent the wavefront position of the vibration wave at a specific moment.

[0023] Starting from each monitoring position of the optical cable, reverse tracing is performed along the negative gradient direction of the travel time field. The negative gradient direction of the travel time field points to the direction where the travel time decreases the fastest. The reverse tracing path eventually converges to the vibration source position. The reverse tracing path is the actual propagation path of the vibration wave.

[0024] During the reverse tracing process, it is determined whether the tracing path crosses the seawater-seabed interface. When the tracing path crosses the seawater-seabed interface, the spatial coordinates of the crossing point are recorded, and the number of interface crossing points on a single tracing path is counted to determine the number of reflections. At each interface crossing point, the tangent direction of the path is calculated based on the travel time field gradient as the incident vector, and the tangent direction of the refracted path is calculated based on the ratio of the sound speeds on both sides of the interface and the interface normal as the refraction vector, thus obtaining the complete refraction path.

[0025] Based on the refraction path and the number of reflections, a propagation path equation system is constructed. By solving the least-squares optimization solution of the propagation time delay in the multi-path propagation equation system, the three-dimensional spatial coordinates of the vibration source are determined. The motion vector field of the vibration source is then fitted using a sequence of three-dimensional spatial coordinates at consecutive time points, including:

[0026] For each monitoring location on the optical cable, the segmented path length between the seawater layer and the seabed layer is calculated based on the refraction path. The additional time delay introduced by the interface reflection is calculated based on the number of reflections. The theoretical propagation time delay is obtained by summing the quotient of each segmented path length divided by the corresponding sound velocity field value. The observation time delay of each monitoring location is extracted from the strain time series data. With the three-dimensional spatial coordinates of the vibration source as unknown variables, a multi-propagation path equation set is established. The multi-propagation path equation set contains multiple time delay equations.

[0027] The time delay deviation between the observation time delay and the theoretical propagation time delay is calculated, a residual vector is constructed, the multipath propagation equations are solved by the least squares method, and the optimal solution of the three-dimensional spatial coordinates of the vibration source is obtained by minimizing the square of the second norm of the residual vector. The three-dimensional spatial coordinate sequence of the vibration source is obtained by repeatedly locating the source at consecutive time points.

[0028] A piecewise cubic polynomial function is used to fit the three-dimensional spatial coordinate sequence. The first derivative of the piecewise cubic polynomial function is calculated. The velocity vector is obtained by sampling the first derivative function at each time step of the three-dimensional spatial coordinate sequence. The velocity vector is then organized with each three-dimensional spatial coordinate in chronological order to form the motion vector field of the vibration source.

[0029] Based on the velocity vector and heading rate of change of the motion vector field, the trajectory envelope region of the vibration source within a future time window is predicted by extrapolation of the kinematic equations. The shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area is calculated to determine whether a warning response should be triggered, including:

[0030] Based on the motion vector field, calculate the two-dimensional projection of the velocity vector in the horizontal plane at each moment and perform principal component analysis to obtain the principal axis direction and the secondary axis direction. Calculate the standard deviation of the velocity projection in the principal axis direction and the secondary axis direction, and construct a velocity uncertainty ellipse with the current velocity vector as the center.

[0031] Extract the velocity vector direction angles at consecutive moments and calculate the angular velocity sequence. Calculate the mean and standard deviation of the angular velocity sequence to determine the heading change trend and heading fluctuation amplitude, and calculate the heading change range. Construct a fan-shaped region with the current position as the vertex, the current heading as the central axis, and the heading change range as the subtended angle. Within the fan-shaped region, uniformly sample the angles of the velocity uncertainty ellipse to generate a set of discrete velocity vectors.

[0032] Set the prediction time length, calculate the product of each discrete velocity vector and the prediction time length as the displacement vector, and superimpose it onto the corresponding three-dimensional spatial coordinates to obtain the set of predicted endpoint positions. Connect the outer endpoints of the set of predicted endpoint positions to form the boundary of the trajectory envelope region.

[0033] Obtain the coordinates of the polygon vertices of the boundary of the sensitive sea area, calculate the vertical distance from each sampling point on the boundary of the trajectory envelope area to each line segment of the boundary of the sensitive sea area, take the minimum value as the shortest approximation distance, and trigger the early warning response when the shortest approximation distance is less than the intrusion warning distance threshold.

[0034] A second aspect of the present invention provides a system for identifying and monitoring low-frequency vibration signals of submarine optical cables, comprising:

[0035] The first unit is used to acquire distributed acoustic vibration signals of submarine optical cables in the sea area to be monitored.

[0036] The second unit is used to perform frequency domain decomposition on the acoustic vibration signal, extract the low-frequency components within the target frequency band, and convert the low-frequency components into strain time series data based on the fiber strain response characteristics.

[0037] The third unit is used to invert and calculate the medium sound velocity distribution between the vibration source and the optical cable by analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, and to deduce the refraction path and reflection number of the vibration wave at the seabed interface based on the medium sound velocity distribution.

[0038] The fourth unit is used to construct a propagation path equation set based on the refraction path and the number of reflections. By solving the least squares optimization solution of the propagation delay in the multi-path propagation equation set, the three-dimensional spatial coordinates of the vibration source are determined, and the motion vector field of the vibration source is fitted using the sequence of three-dimensional spatial coordinates at consecutive time moments.

[0039] The fifth unit is used to predict the trajectory envelope region of the vibration source within a future time window by extrapolating the kinematic equations based on the velocity vector and heading rate of change of the motion vector field, and to calculate the shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area to determine whether to trigger an early warning response.

[0040] A third aspect of the embodiments of the present invention,

[0041] An electronic device is provided, comprising:

[0042] processor;

[0043] Memory used to store processor-executable instructions;

[0044] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0045] Fourth aspect of the present invention,

[0046] A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0047] The beneficial effects of this application are as follows:

[0048] By performing frequency domain decomposition on the distributed acoustic vibration signal of submarine optical cables and extracting low-frequency components, the characteristic signals generated by specific vibration sources in the sea area can be effectively identified, improving the accuracy of vibration source identification and anti-interference capability. Based on the fiber strain response characteristics, the low-frequency components are converted into strain time-series data, realizing accurate conversion of physical quantities and providing a reliable data foundation for subsequent analysis.

[0049] By analyzing the arrival time difference and amplitude attenuation gradient of strain time series data at different spatial locations of the optical cable, the medium sound velocity distribution between the vibration source and the optical cable was calculated by inversion, thus solving the problem of accurate characterization of sound wave propagation characteristics in complex seawater-seabed environments.

[0050] Based on the obtained sound velocity distribution in the medium, the refraction path and reflection number of the vibration wave at the sea-seabed interface were derived, and a multipath propagation equation set that better conforms to the actual physical environment was constructed, improving the positioning accuracy. The least squares optimization method was used to solve the multipath propagation equation set to determine the three-dimensional spatial coordinates of the vibration source. Innovatively, the motion vector field of the vibration source was fitted using a sequence of three-dimensional spatial coordinates at consecutive time moments, achieving high-precision positioning and tracking of targets in the sea area. Attached Figure Description

[0051] Figure 1 This is a flowchart illustrating the method for identifying and monitoring low-frequency vibration signals of submarine optical cables according to an embodiment of the present invention.

[0052] Figure 2 This is a flowchart illustrating the method for deriving the refraction path and reflection number of vibration waves in an embodiment of the present invention. Detailed Implementation

[0053] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0054] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0055] Figure 1 This is a flowchart illustrating the method for identifying and monitoring low-frequency vibration signals of submarine optical cables according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0056] Acquire distributed acoustic vibration signals of submarine optical cables in the sea area to be monitored;

[0057] The acoustic vibration signal is decomposed in the frequency domain to extract the low-frequency components within the target frequency band, and the low-frequency components are converted into strain time series data based on the fiber strain response characteristics.

[0058] By analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, the medium sound velocity distribution between the vibration source and the optical cable is calculated by inversion, and the refraction path and number of reflections of the vibration wave at the seabed interface are derived based on the medium sound velocity distribution.

[0059] Based on the refraction path and the number of reflections, a propagation path equation system is constructed. By solving the least squares optimization solution of the propagation delay in the multi-path propagation equation system, the three-dimensional spatial coordinates of the vibration source are determined, and the motion vector field of the vibration source is fitted using the sequence of three-dimensional spatial coordinates at consecutive time moments.

[0060] Based on the velocity vector and heading rate of change of the motion vector field, the trajectory envelope region of the vibration source within the future time window is predicted by extrapolation of the kinematic equations, and the shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area is calculated to determine whether an early warning response is triggered.

[0061] In one optional implementation, the acoustic vibration signal is decomposed in the frequency domain to extract low-frequency components within the target frequency band, and the low-frequency components are converted into strain time-series data based on the fiber optic strain response characteristics, including:

[0062] The acoustic vibration signal is subjected to time-frequency transformation processing to obtain a three-dimensional distribution matrix of frequency time amplitude. The three-dimensional distribution matrix is ​​then filtered in the frequency domain according to a preset frequency cutoff boundary to extract the frequency domain energy components within the target frequency band as low-frequency components.

[0063] A strain transfer function model for optical fiber is established, which includes a frequency-dependent strain sensitivity coefficient and a phase delay coefficient. The strain sensitivity coefficient characterizes the unit strain caused by sound waves of different frequencies in the optical fiber, and the phase delay coefficient characterizes the phase lag of sound waves of different frequencies during propagation in the optical fiber.

[0064] Based on the frequency distribution of the low-frequency component, the strain sensitivity coefficient and the phase delay coefficient at the corresponding frequency in the fiber strain transfer function model are queried. The amplitude of the low-frequency component is compensated for using the strain sensitivity coefficient, and the phase of the low-frequency component is corrected for time delay using the phase delay coefficient, so as to obtain the compensated and corrected strain equivalent signal.

[0065] The strain equivalent signal is subjected to inverse time-frequency transformation to obtain the strain amplitude sequence of each spatial position of the optical fiber at continuous time, and arranged in chronological order to form strain time series data.

[0066] In this specific embodiment, it is necessary to obtain the distributed acoustic vibration signals of the submarine optical cable in the sea area to be monitored. Specifically, this can be achieved by deploying a coherent optical time domain reflectometer (COTDR) or a distributed acoustic sensing (DAS) system on the submarine optical cable to perform distributed fiber optic sensing monitoring. By transmitting probe light pulses to the optical cable and receiving the backscattered light signals generated by the slight deformation of the optical fiber caused by the sound waves, the received scattered light signals can be obtained after coherent demodulation to obtain the acoustic vibration signals at different spatial locations along the optical cable. These signals typically include acoustic characteristics generated by various vibration sources such as underwater vehicles, biological activities, and ocean current disturbances.

[0067] The acquired acoustic vibration signal is processed by time-frequency transformation. The time-domain vibration signal on the fiber segment with spatial location index i is denoted as x_i(t), where t represents the time variable. Short-time Fourier transform (STFT) is used to perform time-frequency analysis on the signal. A suitable window function w(t) is selected, such as Hanning window or Hamming window. The window length is selected to meet the minimum length required for frequency resolution, usually 128-512 sampling points. For each time segment, its Fourier transform is calculated to obtain the three-dimensional distribution matrix S_i(f,t) of frequency-time-amplitude, where f represents the frequency variable. This process can be described as performing a Fourier transform on the signal at each spatial location i within the time window [tT / 2, t+T / 2], where T is the window length.

[0068] The three-dimensional distribution matrix is ​​frequency-domain filtered according to the preset frequency cutoff boundary. Considering that low-frequency sound waves in the seabed environment have the characteristics of long propagation distance and small attenuation, and that underwater targets usually generate sound waves with characteristic frequencies in the range of several hertz to hundreds of hertz, the frequency cutoff boundary is set to [f_low, f_high], with typical values ​​of [5Hz, 200Hz]. The S_i(f,t) matrix is ​​bandpass filtered in the frequency domain to retain the energy components in this frequency band, thus obtaining the low-frequency component S_i_low(f,t). This processing can effectively eliminate the interference of marine environmental noise such as high-frequency thermal noise and water flow turbulence noise.

[0069] A strain transfer function model for optical fibers is established. When a submarine optical cable is subjected to sound waves, the sound wave energy needs to pass through the outer protective structure of the cable to reach the optical fiber. This process exhibits frequency-dependent transfer characteristics. The established transfer function model H(f) includes two key parameters: the strain sensitivity coefficient α(f) and the phase delay coefficient φ(f). The strain sensitivity coefficient characterizes the unit strain caused by sound waves of different frequencies in the optical fiber. Generally, the strain caused by low-frequency sound waves is smaller than that caused by mid-to-high-frequency sound waves. The phase delay coefficient characterizes the phase lag of sound waves of different frequencies during propagation in the optical fiber. It generally increases with increasing frequency. This model can be calibrated experimentally by measuring the relationship between the optical fiber feedback signal and the actual strain under the action of a standard vibration source with known frequency and amplitude, and fitting the coefficient values ​​at different frequencies.

[0070] Based on the frequency distribution of the low-frequency components, the fiber strain transfer function model is consulted. For each frequency component f in the low-frequency component Si_i_low(f,t), the strain sensitivity coefficient α(f) and phase delay coefficient φ(f) for the corresponding frequency are looked up in the transfer function model. During frequency compensation, the amplitude of Si_i_low(f,t) is divided by the strain sensitivity coefficient α(f) for the corresponding frequency to obtain the compensated amplitude. During phase correction, the phase delay coefficient φ(f) for the corresponding frequency is subtracted from the phase of Si_i_low(f,t) to obtain the corrected phase. This yields the compensated and corrected strain equivalent signal Si_i_compensated(f,t), which more accurately reflects the actual strain state generated by external acoustic waves acting on the fiber.

[0071] The strain equivalent signal is subjected to inverse time-frequency transformation. The compensated and corrected frequency domain signal S_i_compensated(f,t) is transformed back to the time domain through inverse short-time Fourier transform (ISTFT) to obtain the strain amplitude sequence ε_i(t) of each spatial position i in the fiber at continuous time t. In the inverse transformation process, the same window function type and length as the forward transformation are used, and appropriate window overlap is performed to reduce edge effects. The weighted averaging method is used to synthesize a continuous time domain signal in the overlapping area.

[0072] The strain amplitudes at each spatial location are arranged in chronological order to form strain time series data. For each spatial location i, the strain values ​​ε_i(t) at different times t are organized into time series data. To facilitate subsequent analysis, the time series data of all spatial locations are combined into a two-dimensional matrix E, where the row index represents the spatial location and the column index represents the time. This matrix form is convenient for observing the propagation characteristics of sound waves on optical cables and the spatial-temporal correlation.

[0073] In practical applications, low-frequency vibration signals can be feature extracted and identified based on strain time-series data. For example, by analyzing the amplitude, frequency characteristics, and propagation characteristics of strain time-series data, the presence and trajectory of underwater vehicles can be identified. Taking a submarine optical cable monitoring as an example, the acoustic vibration signal obtained by processing the above method revealed periodic low-frequency vibration characteristics in the 5-15Hz frequency band in the strain time-series data. These characteristics propagate along the optical cable at a speed of approximately 1.5 kilometers per second. By comparing with feature templates in the database, it was confirmed to be the acoustic wave characteristics of the propulsion system of a specific type of submersible, thus achieving effective detection of underwater targets.

[0074] Temperature variations and water flow disturbances in the seabed environment can also induce low-frequency strain in optical cables. To improve identification accuracy, an adaptive threshold detection algorithm can be applied to the strain time-series data to eliminate background disturbances. Specifically, by calculating the strain mean and standard deviation over a long time window, a dynamic detection threshold is set. Only when the strain amplitude exceeds the mean plus three times the standard deviation is it considered a valid vibration signal. This technical solution achieves accurate conversion from acoustic vibration signals to optical fiber strain time-series data, solves the problem of inconsistent coupling response characteristics between acoustic waves and optical fibers at different frequencies, and improves the accuracy of optical fiber sensing systems in the field of acoustic vibration monitoring.

[0075] Figure 2 This is a flowchart illustrating the method for deriving the refraction path and reflection number of a vibration wave according to an embodiment of the present invention. In an optional implementation, by analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, the sound velocity distribution of the medium between the vibration source and the optical cable is calculated by inversion, and the refraction path and reflection number of the vibration wave at the seabed interface are derived based on the sound velocity distribution of the medium, including:

[0076] The strain time series data is subjected to wave field separation processing. Based on the timing relationship of waveform arrival and the waveform energy concentration characteristics, the direct wave component, seabed reflected wave component and sea surface reflected wave component are separated and extracted from the aliased signal, and the arrival time difference of each wave component at different spatial locations of the optical cable is calculated.

[0077] Based on the arrival time difference of the direct wave component and the spatial coordinates of the optical cable, a direct wave time delay equation is constructed. Based on the arrival time difference of the seabed reflected wave component and the geometric relationship of the reflection path, a reflected wave time delay equation is constructed. Based on the amplitude attenuation gradient of each wave field component, a medium attenuation constraint equation is established.

[0078] By jointly solving the direct wave time delay equation, the reflected wave time delay equation, and the medium attenuation constraint equation, the sound velocity field of the seawater layer and the sound velocity field of the seabed layer are obtained by inversion. The position of the seawater-seabed interface is determined by using the numerical difference between the attenuation coefficient of the seawater layer and the attenuation coefficient of the seabed layer obtained by inversion.

[0079] Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, the vibration wave propagation trajectory is calculated using a ray tracing algorithm. The number of intersections between the vibration wave propagation trajectory and the seawater-seabed interface is counted as the number of reflections. The incident vector and refraction vector of the ray trajectory at the interface are calculated to determine the refraction path.

[0080] In this specific embodiment, wavefield separation processing is required for the strain time series data. By analyzing the vibration signal characteristics recorded on the submarine optical cable, wavefield components of different propagation paths are identified. Acoustic vibration signals in the seabed environment typically contain multiple components such as direct waves, seabed reflected waves, and sea surface reflected waves. These components superimpose on the optical cable to form a composite wavefield. Using time-frequency analysis methods, wavelet transform or short-time Fourier transform is performed on the strain time series data to obtain the time-frequency characteristic distribution of the signal. According to the physical characteristics of wavefield propagation, direct waves usually arrive first and have the most concentrated energy; seabed reflected waves arrive second, with more significant energy attenuation than direct waves; sea surface reflected waves arrive last, and due to the scattering effect of sea surface undulations, their energy is the most dispersed. Based on these differences in physical characteristics, adaptive filtering technology is used to design a pattern-matched filter, taking arrival time and energy distribution as key features to achieve the separation and extraction of direct wave components, seabed reflected wave components, and sea surface reflected wave components.

[0081] Calculate the time difference of arrival of each wave component at different spatial locations on the optical cable. Determine the precise arrival time of each separated wave component signal at each measurement point on the optical cable using peak detection or cross-correlation analysis. Assume there are N measurement points on the optical cable, labeled P1 to P2. N For each wave component (direct wave, seabed reflected wave, and sea surface reflected wave), an arrival time series containing N elements can be obtained. Taking the first measurement point P1 as the reference point, the arrival time differences of other measurement points relative to the reference point are calculated, forming the direct wave time difference matrix T_direct, the seabed reflected wave time difference matrix T_bottom, and the sea surface reflected wave time difference matrix T_surface. This time difference information directly reflects the propagation characteristics of sound waves in the seabed environment and is the key input data for subsequent sound velocity field inversion.

[0082] Based on the arrival time difference of the direct wave component and the spatial coordinates of the optical cable, a direct wave time delay equation is constructed. Submarine optical cables typically exhibit a curved distribution. The three-dimensional spatial coordinates of each measurement point on the cable can be obtained through cable laying data or measured by an inertial navigation system. For any two measurement points P... i and P jThe propagation delay of the direct wave between two points depends on the spatial distance between them and the average speed of sound along the propagation path. Considering the variation of seawater sound speed with factors such as depth, temperature, and salinity, a layered sound speed model is introduced, representing the seawater sound speed field as a piecewise function of depth z. Based on Fermat's principle, sound waves propagate along the minimum time path in a non-uniform medium. The propagation path of the direct wave between each measurement point is solved using optimization methods, establishing a functional relationship between the sound wave propagation delay and the sound speed field parameters, thus forming a set of equations for the direct wave delay.

[0083] Based on the geometric relationship between the arrival time difference of the seabed reflected wave components and the reflection path, a time delay equation for the reflected wave is constructed. The propagation path of the seabed reflected wave includes the incident path from the sound source to the seabed and the reflection path from the seabed to the receiving point. According to Snell's law, the incident angle equals the reflection angle. Based on this geometric constraint, combined with the spatial coordinates of the optical cable and the arrival time difference information, a time delay equation for the seabed reflected wave is established. This equation needs to consider the problem of the unknown location of the reflection point. An iterative optimization method is adopted to simultaneously solve for the location of the reflection point and the sound velocity parameters of the seabed by minimizing the difference between the observed time delay and the theoretical time delay. Similarly, a corresponding time delay equation is also established for the sea surface reflected wave, forming a complete set of time delay equations for the reflected wave.

[0084] Based on the amplitude attenuation gradient of each wavefield component, a medium attenuation constraint equation is established. During sound wave propagation, energy attenuation follows an exponential decay law, and the degree of attenuation is related to the medium properties and propagation distance. For each separated wavefield component, the amplitude variation trend at various measurement points on the optical cable is analyzed, and the normalized amplitude attenuation curve is calculated. There is a significant difference in the attenuation coefficients of seawater and seabed media, with seabed media typically exhibiting a higher attenuation rate. Based on this difference, a model relating amplitude attenuation to medium parameters is established, forming a medium attenuation constraint equation, providing additional constraints for sound velocity field inversion.

[0085] By jointly solving the direct wave time delay equation, the reflected wave time delay equation, and the medium attenuation constraint equation, the sound velocity fields of the seawater layer and the seabed layer are obtained through inversion. Using a Bayesian inversion framework, the sound velocity inversion problem is transformed into a posterior probability maximization problem. A prior model is constructed, and the sound velocity field parameters are initialized using historical oceanographic data or empirical models. A likelihood function is defined to quantify the degree of fit between the observed data and the model predictions. The sound velocity field parameters are iteratively solved using optimization algorithms such as the conjugate gradient method until the optimal solution is converged. The inversion results include the sound velocity field C_water(z) of the seawater layer and the sound velocity field C_bottom(z) of the seabed layer, where z represents the depth variable. At the same time, the inversion process also yields the seawater layer attenuation coefficient α_water and the seabed layer attenuation coefficient α_bottom.

[0086] By utilizing the numerical difference between the attenuation coefficients of the seawater layer and the seabed layer obtained through inversion, the location of the seawater-seabed interface can be determined. At the interface, the attenuation coefficients exhibit significant changes. By analyzing the depth gradient of the attenuation coefficients, the interface location can be accurately identified. Specifically, the first derivative of the attenuation coefficients is calculated along the depth direction; the point where the derivative reaches its maximum value is the interface location. This method is more robust than relying solely on changes in sound velocity, and is particularly suitable for soft seabed environments where sound velocity changes are relatively gradual but attenuation characteristics differ significantly. The seawater-seabed interface location determined in this way provides important geological boundary conditions for subsequent sound wave propagation analysis.

[0087] Based on the sound velocity field obtained through inversion, a ray tracing algorithm is used to calculate the propagation trajectory of the vibration wave. The ray tracing algorithm is based on Snell's law; in a medium with a changing sound velocity gradient, the sound ray trajectory curves in an arc. The ray equation is solved using numerical integration to obtain the precise path of the sound wave from the vibration source to each point on the optical cable. The number of intersections between the vibration wave propagation trajectory and the seabed interface is counted to determine the number of reflections. At each intersection, the reflection coefficient and transmission coefficient are calculated based on the incident angle and the ratio of the sound velocities of the two media, thus determining the energy distribution between reflection and transmission. For the refraction path, the incident direction vector and the refracted direction vector when the sound ray intersects the interface are recorded; these vectors define the refraction path of the sound wave at the interface.

[0088] In practical applications, submarine optical cable monitoring systems collect acoustic signals emitted by distant ships. By processing strain time-series data using the methods described above, the direct wave and reflected wave components were successfully separated, and the sound velocity of seawater and seabed was obtained through inversion. Based on the inversion results, the seabed reflection and sea surface reflection experienced by the sound wave during propagation were calculated, and the refraction path of the sound wave at the seabed interface was determined.

[0089] In one optional implementation, based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, a ray tracing algorithm is used to calculate the propagation trajectory of the vibration wave. The number of intersections between the vibration wave propagation trajectory and the seawater-seabed interface is counted as the number of reflections. The incident vector and refraction vector of the ray trajectory at the interface are calculated to determine the refraction path, including:

[0090] Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, the travel time field distribution from the vibration source to various points in space is calculated using the fast travel method. The travel time field represents the shortest time required for the vibration wave to propagate from the vibration source to any spatial location. Spatial points with equal values ​​in the travel time field are extracted to form travel time isosurfaces, which represent the wavefront position of the vibration wave at a specific moment.

[0091] Starting from each monitoring position of the optical cable, reverse tracing is performed along the negative gradient direction of the travel time field. The negative gradient direction of the travel time field points to the direction where the travel time decreases the fastest. The reverse tracing path eventually converges to the vibration source position. The reverse tracing path is the actual propagation path of the vibration wave.

[0092] During the reverse tracing process, it is determined whether the tracing path crosses the seawater-seabed interface. When the tracing path crosses the seawater-seabed interface, the spatial coordinates of the crossing point are recorded, and the number of interface crossing points on a single tracing path is counted to determine the number of reflections. At each interface crossing point, the tangent direction of the path is calculated based on the travel time field gradient as the incident vector, and the tangent direction of the refracted path is calculated based on the ratio of the sound speeds on both sides of the interface and the interface normal as the refraction vector, thus obtaining the complete refraction path.

[0093] In this specific embodiment, it is necessary to obtain the sound velocity field of the seawater layer and the sound velocity field of the seabed layer. The travel time field distribution from the vibration source to each point in space is calculated using the fast travel method. The monitored sea area is divided into a three-dimensional grid structure. The grid spacing is determined according to the required accuracy, usually 10-20 meters in the horizontal direction and 5-10 meters in the vertical direction. The sound velocity value at the corresponding position is assigned to each grid node to form a discretized sound velocity field model.

[0094] Assuming the vibration source is located at a certain grid node, its travel time is initialized to 0, while the travel times of all other nodes are initialized to infinity. The core idea of ​​the fast travel method is to start from a node with a known travel time and gradually update the travel times of its neighboring nodes, similar to the process of water waves spreading outward from the source. In practice, a narrow band region is maintained, containing the nodes to be processed. Each time, the node with the smallest travel time is selected from the narrow band for processing, the travel times of its neighboring nodes are updated, and these updated nodes are added to the narrow band. The update of node travel times follows the Eckner equation, taking into account the influence of spatial variations in sound speed on wavefront propagation. The iterative process continues until the travel times of all nodes are determined, finally obtaining the travel time field distribution of the entire computational domain.

[0095] Spatial points with equal values ​​in the travel time field are extracted to form travel time isosurfaces. After the travel time field calculation is completed, a series of time points T1, T2, ..., T are selected. n This is typically chosen at uniform time intervals, such as every 0.1 seconds or every 0.2 seconds. For each time point T... i Find the travel time value equal to T in three-dimensional space. iAll grid nodes that form the travel-time isosurface. In actual operation, due to the limitations of the discretized grid, there are few nodes whose travel time is exactly equal to the specified value. Therefore, interpolation methods are used to determine the intersection points of the isosurface and the grid edges. Commonly used interpolation methods include linear interpolation or cubic spline interpolation. By connecting these intersection points, a triangular grid representation of the travel-time isosurface is formed. These isosurfaces have clear physical meanings and represent the wavefront positions of vibration waves during propagation. The i-th isosurface represents the instantaneous position of the vibration wave after traveling for time T i After time, the shape of the isosurface is affected by the inhomogeneity of the sound speed field. In regions with large sound speed gradients, the isosurface will bend or deform.

[0096] Starting from each monitoring position of the optical cable, reverse tracking is carried out along the negative gradient direction of the travel-time field. After each monitoring point on the optical cable receives a vibration signal, it is necessary to determine the propagation path and the source location of the signal. The gradient vector of the travel-time field points in the direction where the travel time increases fastest, and its negative gradient points in the direction where the travel time decreases fastest, that is, the actual propagation direction of the vibration wave. At each monitoring point, the discrete gradient of the travel-time field is calculated, and the central difference format is used to improve the calculation accuracy. Along the calculated negative gradient direction, move forward with an appropriate step size. After each step forward, recalculate the gradient of the travel-time field at the current position and adjust the tracking direction. This iterative process is essentially to solve the characteristic line equation. The characteristic line is orthogonal to the travel-time isosurface and represents the energy propagation path. During the reverse tracking process, as approaching the vibration source, the travel time value continuously decreases and finally approaches 0. When the travel time of the tracking point is less than the preset threshold, it is considered that the vibration source location has been found and the tracking is terminated. If multiple tracking paths converge to a similar position, then this position is the actual vibration source location.

[0097] During the reverse tracking process, it is necessary to judge whether the tracking path crosses the seawater-seabed interface. The seawater-seabed interface is an interface where the acoustic properties change significantly. When the vibration wave crosses this interface, reflection and refraction will occur. During the tracking process, it is necessary to monitor the relationship between the position of the tracking point and the seawater-seabed interface in real time. The specific method is that after each step of tracking, check whether the two positions before and after the tracking point are on different sides of the interface. The seawater-seabed interface function z = f(x, y) determined in the early stage can be used, where z is the depth and x and y are the horizontal coordinates. For the tracking point (x, y, z), if z < f(x, y), the point is located in the water layer; if z > f(x, y), the point is located in the seabed layer. When the tracking point enters the seabed layer from the water layer or enters the water layer from the seabed layer, it means that the tracking path has crossed the interface. At this time, record the spatial coordinates of the crossing point, and these coordinates can be determined by linear interpolation: Assume the previous point is (x1, y1, z1) and the current point is (x2, y2, z2), then the coordinates of the interface crossing point are (x1, y1, z1) + t·((x2, y2, z2) - (x1, y1,z1)), where t is the interpolation parameter that makes the crossing point exactly lie on the interface.

[0098] The number of reflections is determined by counting the number of interface crossing points on a single tracking path. Each reverse tracking path starting from a monitoring point will cross the sea-seabed interface multiple times, and each crossing corresponds to a reflection or transmission phenomenon. By counting all the crossing points recorded on the tracking path, the number of interface crossings for that path, or the number of reflections, can be obtained. The number of reflections is an important basis for determining the type and location of the vibration source. For example, if the number of reflections for most tracking paths is 0, it indicates that the vibration source is located in the water and there is no obstruction between it and the optical cable; if the number of reflections is 1, it is a vibration source located on the seabed; if the number of reflections is greater than 1, it is a vibration signal propagating over a long distance. By analyzing the distribution of the number of reflections obtained from different monitoring points, the location range of the vibration source can be further narrowed down.

[0099] At each interface crossing point, the tangent direction of the path is calculated based on the travel time field gradient as the incident vector. When the tracking path crosses the interface, the vibration wave is refracted, and the propagation direction changes. At the interface crossing point, the tangent direction of the incident path is calculated, which is the direction vector of the tracking path before the crossing point. This direction vector can be calculated by the coordinate difference between the two tracking points before and after the crossing point, or by directly using the negative gradient direction of the travel time field at the crossing point. At the same time, the normal vector of the interface is calculated, which is the perpendicular direction of the seawater-seabed interface at the crossing point. The interface normal vector can be obtained by the gradient of the interface function. The normal vector of the interface function z=f(x,y) is (-∂f / ∂x,-∂f / ∂y,1). After normalization, the unit normal vector is obtained. The angle between the incident vector and the normal vector is the incident angle, which is the basis for calculating the refraction angle.

[0100] The tangent direction of the refracted path is calculated based on the ratio of the sound velocities on both sides of the interface and the interface normal. According to Snell's law, the sine ratio of the angle of refraction to the angle of incidence is equal to the ratio of the sound velocities on both sides: sin(θ2) / sin(θ1)=c2 / c1, where θ1 is the angle of incidence, θ2 is the angle of refraction, and c1 and c2 are the sound velocities on the incident and refracted sides, respectively. At the crossing point, the sound velocity values ​​on both sides are extracted from the sound velocity field model based on its position. Combined with the known angle of incidence, the angle of refraction is calculated. Using the angle of refraction and the interface normal, the direction of the refracted vector is calculated. Specifically, the incident vector is decomposed into two components: parallel to the interface and perpendicular to the interface. The perpendicular component is transformed according to the angle of refraction and recombined with the unchanged parallel component to obtain the refracted vector.

[0101] After obtaining the complete refraction path, path verification and optimization are performed. The calculated refraction paths are compared with the actual observed arrival times of the vibration signals to verify the accuracy of the tracking results. If there is a significant difference between the theoretical propagation time corresponding to the calculated path and the observation time, the sound velocity field model needs to be adjusted or the travel time field calculation needs to be re-performed. At the same time, the convergence of the tracking paths obtained from different monitoring points is analyzed. Multiple paths should converge to similar source point positions. If the convergence positions are scattered, it indicates that there are multiple vibration sources or that the sound velocity field model is not accurate enough. Through an iterative optimization process, the most reliable vibration source position and propagation path are finally determined.

[0102] In a submarine optical cable monitoring embodiment, a low-frequency vibration signal was detected. Based on the previously retrieved sound velocity field, the sound velocity in the water layer was 1520 m / s, the sound velocity at the seabed surface was 1650 m / s, and it gradually increased to 1780 m / s at depth. After calculating the travel time field using the rapid travel method, a reverse tracing was performed starting from 10 monitoring points. It was found that 8 paths converged to the same area about 30 meters below the seabed, and most paths crossed an interface once. Analysis of the refraction path showed that the incident angle of the vibration wave at the interface was about 45 degrees, and the refraction angle was about 40 degrees, which is consistent with the theoretical calculation. Combining the number of reflections and refraction characteristics, the vibration source was determined to be located inside the seabed, and it was a low-frequency signal generated by local geological activity. This method successfully achieved the precise location of the submarine vibration source and provided an effective tool for submarine environmental monitoring.

[0103] In one optional implementation, a propagation path equation set is constructed based on the refraction path and the number of reflections. The three-dimensional spatial coordinates of the vibration source are determined by solving the least-squares optimization solution of the propagation time delay in the multi-path propagation equation set. The motion vector field of the vibration source is then fitted using a sequence of three-dimensional spatial coordinates at consecutive time points, including:

[0104] For each monitoring location on the optical cable, the segmented path length between the seawater layer and the seabed layer is calculated based on the refraction path. The additional time delay introduced by the interface reflection is calculated based on the number of reflections. The theoretical propagation time delay is obtained by summing the quotient of each segmented path length divided by the corresponding sound velocity field value. The observation time delay of each monitoring location is extracted from the strain time series data. With the three-dimensional spatial coordinates of the vibration source as unknown variables, a multi-propagation path equation set is established. The multi-propagation path equation set contains multiple time delay equations.

[0105] The time delay deviation between the observation time delay and the theoretical propagation time delay is calculated, a residual vector is constructed, the multipath propagation equations are solved by the least squares method, and the optimal solution of the three-dimensional spatial coordinates of the vibration source is obtained by minimizing the square of the second norm of the residual vector. The three-dimensional spatial coordinate sequence of the vibration source is obtained by repeatedly locating the source at consecutive time points.

[0106] A piecewise cubic polynomial function is used to fit the three-dimensional spatial coordinate sequence. The first derivative of the piecewise cubic polynomial function is calculated. The velocity vector is obtained by sampling the first derivative function at each time step of the three-dimensional spatial coordinate sequence. The velocity vector is then organized with each three-dimensional spatial coordinate in chronological order to form the motion vector field of the vibration source.

[0107] In this specific embodiment, for each monitoring location on the optical cable, the segmented path lengths of the seawater layer and the seabed layer are calculated based on the refraction path. The refraction path is obtained through inverse tracing of the travel time field. Each path contains multiple segments located within the seawater layer and the seabed layer, respectively. For each refraction path, the length of each segment needs to be calculated. The entire refraction path is discretized into multiple small line segments. The starting and ending coordinates of each line segment are known, and the length of the line segment can be calculated using Euclidean distance. The medium layer in which the line segment is located is determined based on the positions of the starting and ending points. If both the starting and ending points are located within the seawater layer, then the line segment belongs to the seawater layer path; if both are located within the seabed layer, then it belongs to the seabed layer path; if one point is in the seawater layer and the other point is in the seabed layer, then the line segment needs to be divided into two parts, and the corresponding medium layer path lengths are counted separately. All line segments are categorized and summarized to obtain the total length of the seawater layer path and the total length of the seabed layer path.

[0108] The additional time delay introduced by interface reflection is calculated based on the number of reflections. When a vibration wave reflects at the seawater-seabed interface, a phase change and energy loss occur, introducing an additional time delay. The number of reflections is obtained by statistically analyzing the number of times the wave crosses the interface in the reverse tracing path. For each reflection, the calculation of the additional time delay considers the acoustic impedance difference between the media on both sides of the interface and the incident angle. Generally, the additional time delay is related to the reflection coefficient, which can be calculated from the density and sound velocity ratio of the media on both sides of the interface. In practical applications, an empirical relationship between the number of reflections and the additional time delay can be established based on experimental data. For example, each reflection introduces an additional time delay of 0.01-0.05 seconds, the specific value depending on the incident angle and interface characteristics. For multiple reflections, the additional time delays introduced by each reflection are accumulated to obtain the total additional time delay.

[0109] The theoretical propagation delay is obtained by summing the quotients of each segment's path length divided by the corresponding sound velocity field value. For a refraction path from the vibration source to the monitoring point, the theoretical propagation delay consists of the path propagation delay and the additional reflection delay. The formula for calculating the path propagation delay is: Path propagation delay = Seawater layer path length / Seawater layer sound velocity + Seabed layer path length / Seabed layer sound velocity. The seawater layer sound velocity and seabed layer sound velocity are obtained from the sound velocity field model. Since the sound velocity field in the actual marine environment is spatially variable, this variation must be considered during calculation. Each path segment can be further subdivided to make the sound velocity within each small segment approximately constant, and the propagation delays of each small segment are accumulated to obtain the theoretical propagation delay = path propagation delay + additional reflection delay. A theoretical propagation delay value is calculated for each monitoring location on the optical cable.

[0110] The observation delay at each monitoring location is extracted from strain time-series data. The distributed fiber optic sensing system of the submarine optical cable can collect strain data at various points on the cable in real time, forming a strain time-series. After the vibration signal propagates to the optical cable, it will form characteristic waveforms in the strain time-series. By analyzing these waveforms, the arrival time of the vibration signal at each monitoring location can be extracted. Specifically, methods such as threshold detection or cross-correlation analysis are used to identify abrupt changes or characteristic peaks in the strain signal. Let the time of the vibration event be t0, and the time of the vibration signal arriving at a certain monitoring location be t1. Then the observation delay at that location is t1-t0. If the time of the vibration event is unknown, the point on the optical cable that first receives the signal can be taken as a reference point, and the relative time delay of other points with respect to this reference point can be calculated.

[0111] Using the three-dimensional spatial coordinates of the vibration source as unknown variables, a system of multiple propagation path equations is established. Let the three-dimensional spatial coordinates of the vibration source be (x_s, y_s, z_s), the coordinates of the i-th monitoring position on the optical cable be (x_i, y_i, z_i), the observation delay be T_i^obs, and the theoretical propagation delay be T_i^theo(x_s, y_s, z_s). Then the time delay equation is: T_i^obs = T_i^theo(x_s, y_s, z_s). For n monitoring positions on the optical cable, n time delay equations can be established, forming a system of multiple propagation path equations. The theoretical propagation delay is a nonlinear function of the vibration source coordinates. Because it involves travel time field calculation and path tracing, it is difficult to obtain an analytical expression. In actual solutions, a numerical method can be used to construct the mapping relationship between the vibration source coordinates and the theoretical propagation delay.

[0112] Calculate the time delay deviation between the observed time delay and the theoretical propagation time delay, construct a residual vector, and for a given predicted value of the vibration source coordinates (x_s, y_s, z_s), calculate the theoretical propagation time delay T_i^theo from that point to each monitoring location, and compare it with the observed time delay T_i^obs to obtain the time delay deviation δT_i = T_i^obs - T_i^theo. Construct a residual vector δT = [δT1, δT2, ..., δT] for all monitoring locations. n The magnitude of the residual vector reflects the accuracy of the current guessed coordinates.

[0113] The least squares method is used to solve the multipath propagation equations. The goal of the least squares method is to find the vibration source coordinate solution that minimizes the squared L2 norm of the residual vector. The objective function is defined as J(x_s,y_s,z_s)=||δT|| 2 =Σ(δT_i) 2 The goal is to find the equation (x_s, y_s, z_s) that minimizes J. Due to the nonlinear nature of the equation system, an iterative optimization algorithm is typically used. An initial guess is selected, and the gradient of the objective function is calculated based on prior information or a rough location result to determine the descent direction. The guess is then updated, and this process is repeated until convergence. During the solution process, a regularization term can be introduced to improve the stability of the solution, especially when the monitoring points are unevenly distributed or the measurement noise is high.

[0114] The optimal solution for the three-dimensional spatial coordinates of the vibration source is obtained by minimizing the squared L2 norm of the residual vector. When the iterative process converges, the obtained vibration source coordinates are the optimal solution. The quality of the optimal solution can be evaluated by the residual size, the covariance matrix of the parameter estimation, or the confidence interval. If the residual is still large, it is necessary to re-examine whether the sound velocity field model, reflection time delay estimation, and other aspects are accurate. For complex seabed environments, there are multiple local optima. A multi-starting point initialization strategy is required to start optimization from different initial points, compare the residual size of each solution, and select the global optimal solution. In practical applications, computational efficiency also needs to be considered. Parallel computing or optimization algorithms can be used to accelerate the solution process.

[0115] To obtain the three-dimensional spatial coordinate sequence of the vibration source, repeated localization is performed at consecutive time points. Since the vibration source is in motion, continuous tracking is required. The acquired strain time-series data is divided into fixed time windows. Within each time window, vibration source localization is performed once to obtain the vibration source coordinates at that moment. The length of the time window depends on the sampling rate of the strain signal and the desired tracking accuracy, typically 0.5-2 seconds. Adjacent time windows can overlap appropriately to improve the temporal resolution of the coordinate sequence. For each time window, the aforementioned localization process is repeated to obtain a series of three-dimensional coordinates of the vibration source at various moments, forming a coordinate sequence {(x_s(t_j),y_s(t_j),z_s(t_j))}, j=1,2,...,m.

[0116] Piecewise cubic polynomial functions were used to fit the three-dimensional spatial coordinate sequence. Since noise and errors are unavoidable in the positioning process, directly using the original coordinate sequence for subsequent analysis leads to unstable results. Piecewise cubic polynomial fitting smooths the coordinate sequence while preserving the main characteristics of the vibration source's motion. The x, y, and z coordinate components were fitted separately, and the time axis was divided into several intervals. Each interval was represented by a cubic polynomial function indicating the change of coordinates over time. The fitted functions of adjacent intervals remained continuous at the connection points, and the first and second derivatives were also continuous, forming a smooth transition. The least squares method was used in the fitting process to minimize the deviation between the fitted curve and the original data points.

[0117] Calculate the first derivative of the piecewise cubic polynomial function. The cubic polynomial is f(t) = a·t. 3 +b·t 2 The first derivative of +c·t+d is f'(t)=3a·t 2 +2b·t+c, for the fitted function in each time interval, calculate its first derivative expression. Since piecewise cubic spline fitting is used, the first derivative is guaranteed to be continuous at the interval connection points. Therefore, the obtained derivative function is globally smooth. This step does not require actual data and is based solely on the fitted function obtained in the previous step for analytical derivation.

[0118] The velocity vector is obtained by sampling the first derivative function at each time point in the three-dimensional spatial coordinate sequence. For each time point t_j in the time series, it is substituted into the first derivative function of the corresponding interval to calculate the instantaneous velocity components in the x, y, and z directions, forming the velocity vector v(t_j) = [v_x(t_j), v_y(t_j), v_z(t_j)]. The velocity vector represents the direction and speed of the vibration source at that moment; its magnitude reflects the velocity, and its direction represents the trend of motion.

[0119] The motion vector field of the vibration source is formed by organizing the three-dimensional spatial coordinates in chronological order. The position coordinates and velocity vectors at each moment are combined to form the motion vector field {(x_s(t_j),y_s(t_j),z_s(t_j),v_x(t_j),v_y(t_j),v_z(t_j))},j=1,2,...,m. The motion vector field comprehensively describes the motion state of the vibration source, including the changes in position and velocity over time. By analyzing the motion vector field, the motion mode of the vibration source can be identified, such as uniform motion, accelerated motion, oscillating motion, etc., providing a basis for vibration source type identification and behavior prediction.

[0120] In one optional implementation, based on the velocity vector and heading rate of change of the motion vector field, the trajectory envelope region of the vibration source within a future time window is predicted by extrapolation of kinematic equations, and the shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area is calculated to determine whether to trigger an early warning response, including:

[0121] Based on the motion vector field, calculate the two-dimensional projection of the velocity vector in the horizontal plane at each moment and perform principal component analysis to obtain the principal axis direction and the secondary axis direction. Calculate the standard deviation of the velocity projection in the principal axis direction and the secondary axis direction, and construct a velocity uncertainty ellipse with the current velocity vector as the center.

[0122] Extract the velocity vector direction angles at consecutive moments and calculate the angular velocity sequence. Calculate the mean and standard deviation of the angular velocity sequence to determine the heading change trend and heading fluctuation amplitude, and calculate the heading change range. Construct a fan-shaped region with the current position as the vertex, the current heading as the central axis, and the heading change range as the subtended angle. Within the fan-shaped region, uniformly sample the angles of the velocity uncertainty ellipse to generate a set of discrete velocity vectors.

[0123] Set the prediction time length, calculate the product of each discrete velocity vector and the prediction time length as the displacement vector, and superimpose it onto the corresponding three-dimensional spatial coordinates to obtain the set of predicted endpoint positions. Connect the outer endpoints of the set of predicted endpoint positions to form the boundary of the trajectory envelope region.

[0124] Obtain the coordinates of the polygon vertices of the boundary of the sensitive sea area, calculate the vertical distance from each sampling point on the boundary of the trajectory envelope area to each line segment of the boundary of the sensitive sea area, take the minimum value as the shortest approximation distance, and trigger the early warning response when the shortest approximation distance is less than the intrusion warning distance threshold.

[0125] In this specific embodiment, after obtaining the motion vector field data of the vibration source, velocity vector analysis needs to be performed. The three-dimensional velocity vector is projected onto the horizontal plane to obtain two-dimensional velocity components. Principal component analysis is used to process these two-dimensional velocity components, calculate the covariance matrix, and obtain its eigenvalues ​​and eigenvectors. The two directions corresponding to the eigenvectors are the principal axis direction and the secondary axis direction. The square root of the eigenvalue corresponds to the standard deviation in these two directions. With the current velocity vector as the center and the principal axis and secondary axis directions as the coordinate axes, a velocity uncertainty ellipse is constructed. For example, if the standard deviation in the principal axis direction is 2.5 m / s and the standard deviation in the secondary axis direction is 1.2 m / s, then a velocity uncertainty ellipse of the corresponding size can be constructed.

[0126] Extract velocity vector direction angle data at continuous time points, calculate the angle difference between adjacent time points, divide by the time interval to obtain the angular velocity sequence, and calculate the mean and standard deviation of this angular velocity sequence. The mean reflects the overall trend of heading change, and the standard deviation reflects the degree of heading fluctuation. According to the 3σ principle, the heading change range is set to mean ± 3 times the standard deviation. Using the current position of the vibration source as the vertex, the current heading as the central axis, and the heading change range as the angle, construct a fan-shaped prediction region. For example, if the mean angular velocity is 0.5 degrees / second, the standard deviation is 0.2 degrees / second, and the prediction time window is 10 minutes, then the heading change range is approximately ±36 degrees.

[0127] Within the constructed sector-shaped region, the velocity uncertainty ellipse is sampled at uniform angles. Typically, 8 to 16 uniformly distributed angles are selected, and a velocity vector is extracted from the ellipse boundary in each angular direction. This generates a discrete set of velocity vectors, which represent the motion direction and velocity combination of the vibration source.

[0128] Set the prediction time length, such as 10 minutes or 30 minutes, depending on the application scenario and early warning requirements. For each discrete velocity vector, multiply it by the prediction time length to obtain the corresponding displacement vector. Superimpose these displacement vectors onto the current three-dimensional spatial coordinates of the vibration source to obtain a series of predicted endpoint positions. These endpoint positions form a set, and connect its outer points to form the boundary of the trajectory envelope region. Usually, a convex hull algorithm (such as Graham scan method) is used to determine the outer points and construct a closed polygon.

[0129] After obtaining the vertex coordinates of the boundary polygon of the sensitive area in the sea area, calculate the shortest distance between the trajectory envelope and the sensitive area. For the sampling points on the boundary of the trajectory envelope, calculate the perpendicular distance from them to each line segment of the boundary of the sensitive area. If the projection of the sampling point on the line segment is not on the line segment, calculate the distance from the point to the two endpoints of the line segment and take the minimum value. Traverse all combinations of sampling points and line segments, and the minimum distance value obtained is the shortest approximation distance.

[0130] When the calculated shortest approach distance is less than the preset intrusion warning distance threshold, a warning response is triggered. The warning threshold can be dynamically adjusted according to the type of vibration source, the importance of the sensitive area, and the response time requirements. For example, a larger warning distance threshold (e.g., 5000 meters) can be set for high-speed moving vibration sources; a smaller warning distance threshold (e.g., 2000 meters) can be set for low-speed moving vibration sources.

[0131] To improve prediction accuracy, the prediction model can be optimized by combining historical trajectory data. By comparing the deviation between historical prediction results and actual trajectories, the parameters of the velocity uncertainty ellipse and the calculation method of the heading change range can be adjusted. At the same time, the influence of environmental factors such as ocean currents and wind on the motion of the vibration source can be considered, and corresponding correction terms can be introduced into the calculation of velocity vector and heading change.

[0132] Furthermore, different kinematic models can be used for trajectory prediction for different types of vibration sources. For vibration sources that are mainly in uniform linear motion, a constant velocity model can be used. For vibration sources that frequently change direction, it is necessary to increase the range of heading changes and use curve fitting methods to predict the trajectory.

[0133] In practical applications, warning levels can be divided into multiple grades, determined by the ratio of the shortest approach distance to a threshold. When predictions indicate that a vibration source is entering a sensitive area, the system will issue a warning message of the corresponding level, including the estimated approach time, approach location, and recommended measures.

[0134] The present invention provides a system for identifying and monitoring low-frequency vibration signals of submarine optical cables, comprising:

[0135] The first unit is used to acquire distributed acoustic vibration signals of submarine optical cables in the sea area to be monitored.

[0136] The second unit is used to perform frequency domain decomposition on the acoustic vibration signal, extract the low-frequency components within the target frequency band, and convert the low-frequency components into strain time series data based on the fiber strain response characteristics.

[0137] The third unit is used to invert and calculate the medium sound velocity distribution between the vibration source and the optical cable by analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, and to deduce the refraction path and reflection number of the vibration wave at the seabed interface based on the medium sound velocity distribution.

[0138] The fourth unit is used to construct a propagation path equation set based on the refraction path and the number of reflections. By solving the least squares optimization solution of the propagation delay in the multi-path propagation equation set, the three-dimensional spatial coordinates of the vibration source are determined, and the motion vector field of the vibration source is fitted using the sequence of three-dimensional spatial coordinates at consecutive time moments.

[0139] The fifth unit is used to predict the trajectory envelope region of the vibration source within a future time window by extrapolating the kinematic equations based on the velocity vector and heading rate of change of the motion vector field, and to calculate the shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area to determine whether to trigger an early warning response.

[0140] A third aspect of the present invention provides an electronic device, comprising:

[0141] processor;

[0142] Memory used to store processor-executable instructions;

[0143] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0144] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0145] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0146] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for identifying and monitoring low-frequency vibration signals of submarine optical cables, characterized in that, include: Acquire distributed acoustic vibration signals of submarine optical cables in the sea area to be monitored; The acoustic vibration signal is decomposed in the frequency domain to extract the low-frequency components within the target frequency band, and the low-frequency components are converted into strain time series data based on the fiber strain response characteristics. By analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, the medium sound velocity distribution between the vibration source and the optical cable is calculated by inversion, and the refraction path and number of reflections of the vibration wave at the seabed interface are derived based on the medium sound velocity distribution. Based on the refraction path and the number of reflections, a propagation path equation system is constructed. By solving the least squares optimization solution of the propagation delay in the multi-path propagation equation system, the three-dimensional spatial coordinates of the vibration source are determined, and the motion vector field of the vibration source is fitted using the sequence of three-dimensional spatial coordinates at consecutive time moments. Based on the velocity vector and heading rate of change of the motion vector field, the trajectory envelope region of the vibration source within the future time window is predicted by extrapolation of the kinematic equations, and the shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area is calculated to determine whether an early warning response is triggered.

2. The method according to claim 1, characterized in that, The acoustic vibration signal is decomposed in the frequency domain to extract low-frequency components within the target frequency band, and the low-frequency components are converted into strain time-series data based on the fiber optic strain response characteristics, including: The acoustic vibration signal is subjected to time-frequency transformation processing to obtain a three-dimensional distribution matrix of frequency time amplitude. The three-dimensional distribution matrix is ​​then filtered in the frequency domain according to a preset frequency cutoff boundary to extract the frequency domain energy components within the target frequency band as low-frequency components. A strain transfer function model for optical fiber is established, which includes a frequency-dependent strain sensitivity coefficient and a phase delay coefficient. The strain sensitivity coefficient characterizes the unit strain caused by sound waves of different frequencies in the optical fiber, and the phase delay coefficient characterizes the phase lag of sound waves of different frequencies during propagation in the optical fiber. Based on the frequency distribution of the low-frequency component, the strain sensitivity coefficient and the phase delay coefficient at the corresponding frequency in the fiber strain transfer function model are queried. The amplitude of the low-frequency component is compensated for using the strain sensitivity coefficient, and the phase of the low-frequency component is corrected for time delay using the phase delay coefficient, so as to obtain the compensated and corrected strain equivalent signal. The strain equivalent signal is subjected to inverse time-frequency transformation to obtain the strain amplitude sequence of each spatial position of the optical fiber at continuous time, and arranged in chronological order to form strain time series data.

3. The method according to claim 1, characterized in that, By analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, the medium sound velocity distribution between the vibration source and the optical cable is calculated by inversion. Based on the medium sound velocity distribution, the refraction path and number of reflections of the vibration wave at the seabed interface are derived, including: The strain time series data is subjected to wave field separation processing. Based on the timing relationship of waveform arrival and the waveform energy concentration characteristics, the direct wave component, seabed reflected wave component and sea surface reflected wave component are separated and extracted from the aliased signal, and the arrival time difference of each wave component at different spatial locations of the optical cable is calculated. Based on the arrival time difference of the direct wave component and the spatial coordinates of the optical cable, a direct wave time delay equation is constructed. Based on the arrival time difference of the seabed reflected wave component and the geometric relationship of the reflection path, a reflected wave time delay equation is constructed. Based on the amplitude attenuation gradient of each wave field component, a medium attenuation constraint equation is established. By jointly solving the direct wave time delay equation, the reflected wave time delay equation, and the medium attenuation constraint equation, the sound velocity field of the seawater layer and the sound velocity field of the seabed layer are obtained by inversion. The position of the seawater-seabed interface is determined by using the numerical difference between the attenuation coefficient of the seawater layer and the attenuation coefficient of the seabed layer obtained by inversion. Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, the vibration wave propagation trajectory is calculated using a ray tracing algorithm. The number of intersections between the vibration wave propagation trajectory and the seawater-seabed interface is counted as the number of reflections. The incident vector and refraction vector of the ray trajectory at the interface are calculated to determine the refraction path.

4. The method according to claim 3, characterized in that, Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, a ray tracing algorithm is used to calculate the propagation trajectory of the vibration wave. The number of intersections between the vibration wave propagation trajectory and the seawater-seabed interface is counted as the number of reflections. The incident vector and refraction vector of the ray trajectory at the interface are calculated to determine the refraction path, including: Based on the sound velocity field of the seawater layer and the sound velocity field of the seabed layer, the travel time field distribution from the vibration source to various points in space is calculated using the fast travel method. The travel time field represents the shortest time required for the vibration wave to propagate from the vibration source to any spatial location. Spatial points with equal values ​​in the travel time field are extracted to form travel time isosurfaces, which represent the wavefront position of the vibration wave at a specific moment. Starting from each monitoring position of the optical cable, reverse tracing is performed along the negative gradient direction of the travel time field. The negative gradient direction of the travel time field points to the direction where the travel time decreases the fastest. The reverse tracing path eventually converges to the vibration source position. The reverse tracing path is the actual propagation path of the vibration wave. During the reverse tracing process, it is determined whether the tracing path crosses the seawater-seabed interface. When the tracing path crosses the seawater-seabed interface, the spatial coordinates of the crossing point are recorded, and the number of interface crossing points on a single tracing path is counted to determine the number of reflections. At each interface crossing point, the tangent direction of the path is calculated based on the travel time field gradient as the incident vector, and the tangent direction of the refracted path is calculated based on the ratio of the sound speeds on both sides of the interface and the interface normal as the refraction vector, thus obtaining the complete refraction path.

5. The method according to claim 1, characterized in that, Based on the refraction path and the number of reflections, a propagation path equation system is constructed. By solving the least-squares optimization solution of the propagation time delay in the multi-path propagation equation system, the three-dimensional spatial coordinates of the vibration source are determined. The motion vector field of the vibration source is then fitted using a sequence of three-dimensional spatial coordinates at consecutive time points, including: For each monitoring location on the optical cable, the segmented path length between the seawater layer and the seabed layer is calculated based on the refraction path. The additional time delay introduced by the interface reflection is calculated based on the number of reflections. The theoretical propagation time delay is obtained by summing the quotient of each segmented path length divided by the corresponding sound velocity field value. The observation time delay of each monitoring location is extracted from the strain time series data. With the three-dimensional spatial coordinates of the vibration source as unknown variables, a multi-propagation path equation set is established. The multi-propagation path equation set contains multiple time delay equations. The time delay deviation between the observation time delay and the theoretical propagation time delay is calculated, a residual vector is constructed, the multipath propagation equations are solved by the least squares method, and the optimal solution of the three-dimensional spatial coordinates of the vibration source is obtained by minimizing the square of the second norm of the residual vector. The three-dimensional spatial coordinate sequence of the vibration source is obtained by repeatedly locating the source at consecutive time points. A piecewise cubic polynomial function is used to fit the three-dimensional spatial coordinate sequence. The first derivative of the piecewise cubic polynomial function is calculated. The velocity vector is obtained by sampling the first derivative function at each time step of the three-dimensional spatial coordinate sequence. The velocity vector is then organized with each three-dimensional spatial coordinate in chronological order to form the motion vector field of the vibration source.

6. The method according to claim 1, characterized in that, Based on the velocity vector and heading rate of change of the motion vector field, the trajectory envelope region of the vibration source within a future time window is predicted by extrapolation of the kinematic equations. The shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area is calculated to determine whether a warning response should be triggered, including: Based on the motion vector field, calculate the two-dimensional projection of the velocity vector in the horizontal plane at each moment and perform principal component analysis to obtain the principal axis direction and the secondary axis direction. Calculate the standard deviation of the velocity projection in the principal axis direction and the secondary axis direction, and construct a velocity uncertainty ellipse with the current velocity vector as the center. Extract the velocity vector direction angles at consecutive moments and calculate the angular velocity sequence. Calculate the mean and standard deviation of the angular velocity sequence to determine the heading change trend and heading fluctuation amplitude, and calculate the heading change range. Construct a fan-shaped region with the current position as the vertex, the current heading as the central axis, and the heading change range as the subtended angle. Within the fan-shaped region, uniformly sample the angles of the velocity uncertainty ellipse to generate a set of discrete velocity vectors. Set the prediction time length, calculate the product of each discrete velocity vector and the prediction time length as the displacement vector, and superimpose it onto the corresponding three-dimensional spatial coordinates to obtain the set of predicted endpoint positions. Connect the outer endpoints of the set of predicted endpoint positions to form the boundary of the trajectory envelope region. Obtain the coordinates of the polygon vertices of the boundary of the sensitive sea area, calculate the vertical distance from each sampling point on the boundary of the trajectory envelope area to each line segment of the boundary of the sensitive sea area, take the minimum value as the shortest approximation distance, and trigger the early warning response when the shortest approximation distance is less than the intrusion warning distance threshold.

7. A system for identifying and monitoring low-frequency vibration signals of submarine optical cables, used to implement the method as described in any one of claims 1-6, characterized in that, include: The first unit is used to acquire distributed acoustic vibration signals of submarine optical cables in the sea area to be monitored. The second unit is used to perform frequency domain decomposition on the acoustic vibration signal, extract the low-frequency components within the target frequency band, and convert the low-frequency components into strain time series data based on the fiber strain response characteristics. The third unit is used to invert and calculate the medium sound velocity distribution between the vibration source and the optical cable by analyzing the arrival time difference and amplitude attenuation gradient of the strain time series data at different spatial locations of the optical cable, and to deduce the refraction path and reflection number of the vibration wave at the seabed interface based on the medium sound velocity distribution. The fourth unit is used to construct a propagation path equation set based on the refraction path and the number of reflections. By solving the least squares optimization solution of the propagation delay in the multi-path propagation equation set, the three-dimensional spatial coordinates of the vibration source are determined, and the motion vector field of the vibration source is fitted using the sequence of three-dimensional spatial coordinates at consecutive time moments. The fifth unit is used to predict the trajectory envelope region of the vibration source within a future time window by extrapolating the kinematic equations based on the velocity vector and heading rate of change of the motion vector field, and to calculate the shortest approximation distance between the trajectory envelope region and the boundary of the sensitive sea area to determine whether to trigger an early warning response.

8. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images