Underwater target cooperative communication guidance method and system based on multi-modal perception

Abstract

The application provides a kind of underwater target cooperative communication guidance method and system based on multi-modal perception, it is related to underwater environment detection and communication technical field, first acquire the multi-modal perception flow such as sound wave, optical and magnetic anomaly etc. that underwater multi-modal perception node array gathers, generate spatial topological relation descriptor by cross-modal correlation feature analysis, then according to spatial topological relation descriptor identification target motion trend and shelter distribution boundary parameter, construct dynamic channel occupancy state evolution model, call dynamic channel occupancy state evolution model to deduce candidate communication link quality, obtain quality fluctuation characteristic sequence and interference overlapping area identification, finally based on these data generates cooperative communication guidance instruction set, realizes the efficient cooperative communication guidance of underwater target.

Description

Technical Field

[0001] This invention relates to the field of underwater environment detection and communication technology, and more specifically, to an underwater target cooperative communication guidance method and system based on multimodal perception. Background Technology

[0002] In the field of underwater environment detection and communication, accurate perception and efficient communication guidance of underwater targets are crucial issues. Existing underwater target perception and communication methods have many limitations. On the one hand, traditional underwater perception methods often rely solely on single-modal perception data, such as using only acoustic or optical sensing. However, the underwater environment is complex and variable, and single-modal perception data is easily affected by environmental factors. For example, acoustic waves are reflected, refracted, and scattered in complex seabed topography, leading to a decrease in perception accuracy; optical imaging underwater is affected by factors such as water turbidity and light intensity, resulting in unstable image quality. On the other hand, in terms of communication guidance, existing methods do not fully consider the dynamic channel occupancy status between underwater targets and communication nodes. Obstacles and target movement in the underwater environment cause the channel status to change continuously. Existing communication guidance strategies cannot adjust communication parameters in real time according to these dynamic changes, thus affecting communication quality and guidance effectiveness, and failing to meet the high requirements for target cooperative communication guidance in complex underwater environments. Summary of the Invention

[0003] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide an underwater target cooperative communication guidance method based on multimodal perception, the method comprising: The multimodal sensing stream acquired by the underwater multimodal sensing node array includes a set of acoustic wave propagation waveform samples, a set of optical imaging image frames, and a set of magnetic anomaly change time series records. The multimodal sensing flow is subjected to cross-modal correlation feature parsing processing to generate a spatial topological relationship descriptor that reflects the environmental elements of the target water area. The spatial topological relationship descriptor contains a feature identification mapping chain that points to the same underwater target body in different modal sensing data. The motion trend parameters and cover distribution boundary parameters of the underwater target are identified based on the spatial topology descriptor, and a dynamic channel occupancy state evolution model between the underwater target and each communication node is constructed based on the motion trend parameters and the cover distribution boundary parameters. The dynamic channel occupancy state evolution model is invoked to perform link quality evolution deduction on candidate communication links between each communication node, thereby obtaining the quality fluctuation feature sequence of each candidate communication link in the time dimension and the interference overlap area identifier in the spatial dimension. Based on the quality fluctuation characteristic sequence and the interference overlap region identifier, a set of cooperative communication guidance instructions for the underwater target is generated. The set of cooperative communication guidance instructions includes the transmit power adjustment amount, beam pointing angle adjustment amount, and communication timing scheduling identifier allocated to each communication node.

[0004] Furthermore, embodiments of the present invention also provide an underwater target cooperative communication guidance system based on multimodal perception, comprising: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to execute the above-described underwater target cooperative communication guidance method based on multimodal perception by executing the machine-executable instructions.

[0005] In another aspect, embodiments of the present invention also provide a computer program product, the computer program product including machine-executable instructions stored in a computer-readable storage medium, a processor of a multimodal perception-based underwater target cooperative communication guidance system reading the machine-executable instructions from the computer-readable storage medium, the processor executing the machine-executable instructions, causing the multimodal perception-based underwater target cooperative communication guidance system to execute the above-described multimodal perception-based underwater target cooperative communication guidance method.

[0006] Based on the above, by acquiring sensing flows containing multiple modes such as acoustic, optical, and magnetic anomalies collected by an underwater multimodal sensing node array, and performing cross-modal correlation feature analysis, a spatial topological relationship descriptor reflecting the environmental elements of the target water area is generated. This fully utilizes the complementarity of different modal sensing data, improving the accuracy and comprehensiveness of underwater target sensing. Based on the spatial topological relationship descriptor, the motion trend of underwater targets and the boundary parameters of obstruction distribution are identified, and a dynamic channel occupancy state evolution model is constructed. This model can reflect the changes in channel state between underwater targets and communication nodes in real time and accurately. The dynamic channel occupancy state evolution model is then used to perform link quality evolution deduction processing on candidate communication links, obtaining quality fluctuation feature sequences and interference overlap area identifiers. Based on these data, a set of cooperative communication guidance instructions is generated. This allows for real-time adjustment of the transmission power, beam pointing angle, and communication timing of each communication node according to dynamic changes in the underwater environment, effectively improving the quality and stability of underwater communication and achieving efficient cooperative communication guidance for underwater targets. Attached Figure Description

[0007] Figure 1 This is a schematic diagram of the execution flow of the underwater target cooperative communication guidance method based on multimodal perception provided in an embodiment of the present invention.

[0008] Figure 2This is a schematic diagram of exemplary hardware and software components of an underwater target cooperative communication guidance system based on multimodal perception provided in an embodiment of the present invention. Detailed Implementation

[0009] Figure 1 This is a flowchart illustrating an embodiment of the underwater target cooperative communication guidance method based on multimodal perception provided by the present invention, which will be described in detail below.

[0010] Step S110: Obtain the multimodal sensing stream collected by the underwater multimodal sensing node array. The multimodal sensing stream includes a set of acoustic wave propagation waveform samples, a set of optical imaging image frames, and a set of magnetic anomaly change time series records.

[0011] Specifically, in a specific underwater monitoring area, an underwater multimodal sensing node array consisting of multiple sensing nodes is pre-deployed. Each sensing node integrates an acoustic wave acquisition module, an optical imaging module, and a magnetic anomaly detection module. Each sensing node synchronously initiates data acquisition according to a preset sampling frequency. The acoustic wave acquisition module continuously acquires underwater sound field signals through a hydrophone, performs analog-to-digital conversion on the received analog sound wave signals, and generates a sound wave propagation waveform sampling set. The data format in this sound wave propagation waveform sampling set is a discrete-time series containing timestamps and amplitude values. Each sampling point records the corresponding sampling time and sound pressure amplitude quantization value. The optical imaging module captures underwater environmental images at fixed time intervals, generating a set of optical imaging image frames. Each image frame consists of a pixel matrix, and each pixel contains the brightness values ​​of the red, green, and blue channels, as well as the pixel's coordinates in the image plane coordinate system. The magnetic anomaly detection module uses a magnetometer to continuously record the total field strength of the geomagnetic field or the three-component fluxgate data, generating a set of time-series records of magnetic anomaly changes. This set consists of a series of magnetic induction intensity vectors arranged in chronological order, with each record containing a timestamp and the magnetic induction intensity component values ​​in three-dimensional space.

[0012] Step S120: Perform cross-modal correlation feature parsing processing on the multimodal sensing flow to generate a spatial topological relationship descriptor that reflects the environmental elements of the target water area. The spatial topological relationship descriptor contains a feature identification mapping chain that points to the same underwater target body in different modal sensing data.

[0013] Step S121: Perform pulse feature extraction processing on the sound wave propagation waveform sampling set to obtain a sound wave pulse arrival time distribution set and a waveform distortion coefficient set. The sound wave pulse arrival time distribution set includes the time offset of the sound wave pulse received by each sensing node, and the waveform distortion coefficient set includes the pulse broadening parameter and the pulse energy attenuation gradient parameter.

[0014] Step S1211: Perform energy envelope detection processing on each waveform sample in the sound wave propagation waveform sampling set to obtain the waveform energy envelope curve. The horizontal axis of the waveform energy envelope curve is the time sampling point, and the vertical axis is the normalized energy amplitude.

[0015] Specifically, for each discrete-time waveform sequence in the sound wave propagation waveform sampling set, it is represented as a function of amplitude changing with time. A Hilbert transform is applied to this sequence to construct an analytic signal. The magnitude of the analytic signal is the instantaneous amplitude envelope of the original signal. The calculated instantaneous amplitude values ​​are arranged in chronological order to form an initial energy envelope curve. To eliminate the influence of amplitude differences between different waveforms, the initial energy envelope curve is normalized; that is, the maximum amplitude value on the curve is found, and the amplitude value at each point on the curve is divided by this maximum amplitude value, resulting in a normalized energy envelope curve with a ordinate ranging from zero to one. The abscissa of this normalized energy envelope curve retains the same time sampling point index as the original waveform, thus obtaining a waveform energy envelope curve that can be used for subsequent peak detection.

[0016] Step S1212: Perform peak detection processing on the waveform energy envelope curve, identify local maxima points that exceed a preset energy threshold, and use the time sampling points corresponding to the local maxima points as candidate pulse arrival time points to obtain a set of candidate pulse arrival time points.

[0017] Specifically, the waveform energy envelope curve generated in step S1211 is scanned point by point. For each point on the curve, its energy amplitude is compared with the energy amplitudes of its preceding and following points. If the amplitude of this point is greater than the amplitudes of both its preceding and following points, then this point is marked as a local maximum. Simultaneously, an energy threshold is set, which can be dynamically adjusted according to the background noise level, for example, three times the average energy of the entire envelope curve. Only those local maxima whose energy amplitude exceeds this energy threshold are retained. The time sampling points corresponding to the retained local maxima points are recorded as candidate pulse arrival times. All candidate points are arranged in chronological order to form a set of candidate pulse arrival times.

[0018] Step S1213: Perform cluster analysis on the candidate pulse arrival time point set, group adjacent candidate pulse arrival time points with time intervals less than a preset time interval threshold into the same pulse cluster, each pulse cluster corresponds to one sound wave pulse, extract the point with the highest energy in each pulse cluster as the arrival time point of the pulse, and generate a set of sound wave pulse arrival time points.

[0019] Specifically, the points in the candidate pulse arrival time set obtained in step S1212 are sorted chronologically. A time interval threshold is set, which is usually slightly smaller than the expected minimum pulse interval. Starting from the first candidate point, the time difference between the current candidate point and the next candidate point is calculated sequentially. If the time difference is less than the preset time interval threshold, the two points are considered to belong to the same pulse cluster, and subsequent points that meet the condition are included in the current cluster. When the time difference between the next candidate point and the current point is greater than or equal to the threshold, the current cluster ends and a new cluster begins. After traversing all candidate points, each cluster represents an independent acoustic pulse. For each pulse cluster, the energy amplitude corresponding to all candidate points within the cluster is compared, and the candidate point with the largest energy amplitude is selected, and its time sampling point is determined as the precise arrival time of the pulse. The arrival time points determined for all pulse clusters are arranged in chronological order to generate the final set of acoustic pulse arrival time points.

[0020] Step S1214: Based on the waveform sampling values ​​within a preset time window before and after the arrival time of each sound wave pulse, extract pulse waveform segments, perform Fourier transform processing on the pulse waveform segments to obtain the pulse spectrum distribution, and calculate the proportion of the main frequency energy in the pulse spectrum distribution to the total energy as the pulse energy concentration parameter.

[0021] Specifically, for each arrival time point in the set of acoustic pulse arrival time points generated in step S1213, a preset number of sampling points are extracted from the original acoustic wave propagation waveform sampling sequence, centered on that arrival time point, both forward and backward, forming a continuous pulse waveform segment. This preset time window size should be large enough to completely cover the duration of a pulse. Subsequently, a Fast Fourier Transform is performed on this pulse waveform segment to convert it from a time-domain signal to a frequency-domain signal, obtaining the pulse spectrum distribution. This pulse spectrum distribution shows the distribution of signal energy at different frequencies, specifically the amplitude or power spectral density of each frequency component. The frequency corresponding to the maximum amplitude in the spectrum distribution is identified as the dominant frequency. A bandwidth range centered on the dominant frequency is determined, and the sum of the energy of all frequency components within this bandwidth range is calculated to obtain the dominant frequency energy. Simultaneously, the total energy of all frequency components within the entire spectrum range is calculated to obtain the total energy. The ratio of the dominant frequency energy to the total energy is calculated, and this ratio is used as the pulse energy concentration parameter. The larger this parameter value, the more concentrated the pulse energy is near the dominant frequency, and the smaller the waveform distortion.

[0022] Step S1215: Perform pulse width measurement processing on the pulse waveform segment, measure the time span from the start of the rising edge to the end of the falling edge of the pulse, obtain the pulse duration parameter, compare the pulse duration parameter with the preset standard pulse duration, calculate the duration deviation, and use it as the pulse width parameter.

[0023] Specifically, for each pulse waveform segment captured in step S1214, the start and end points of the pulse are first determined. A threshold method can be used: searching backward from the pulse arrival time to find the point where the energy amplitude first falls below a preset starting threshold as the rising edge start point; searching backward from the pulse arrival time to find the point where the energy amplitude last falls below a preset ending threshold as the falling edge end point. The time span from the rising edge start point to the falling edge end point is calculated, which is the product of the difference between the two time sampling points and the sampling interval, to obtain the pulse duration parameter. This measured pulse duration parameter is compared with a pre-stored standard pulse duration (e.g., the nominal width of the transmitted pulse), and the difference between the two is calculated. If the pulse duration parameter is greater than the standard pulse duration, the difference is positive, indicating that the pulse is broadened; if it is less, the difference is negative, indicating that the pulse is compressed. This difference, or the ratio of the difference to the standard pulse duration, is used as the pulse broadening degree parameter.

[0024] Step S1216: Perform energy integration processing on the pulse waveform segment to obtain the total pulse energy value. Compare the total pulse energy value with the reference energy value of the transmitting end and calculate the energy attenuation factor as the pulse energy attenuation gradient parameter. The pulse energy attenuation gradient parameter reflects the energy loss on the sound wave propagation path.

[0025] Specifically, for each pulse waveform segment captured in step S1214, the squares of the amplitudes of all sampling points within the segment are summed. For more precise energy calculation, the squares of the amplitudes can be integrated over time. The result is the total pulse energy value. The reference energy value of the transmitting end emitting the sound pulse is obtained from the system configuration parameters. This reference energy value is measured at close range under standard conditions. The ratio of the transmitting end reference energy value to the current total pulse energy value is calculated to obtain the energy attenuation factor. For example, if the reference energy value is E_ref and the total pulse energy value is E_pulse, then the energy attenuation factor is equal to E_ref divided by E_pulse. This energy attenuation factor directly reflects the degree of total energy loss caused by factors such as seawater absorption, scattering, and geometric spread during the propagation of the sound wave from the transmitting source to the receiving node. This energy attenuation factor is used as the pulse energy attenuation gradient parameter.

[0026] Step S1217: Perform hyperbolic positioning calculation on the arrival time points of acoustic pulses received by multiple sensing nodes corresponding to the same underwater target to obtain the acoustic positioning coordinates of the underwater target. Combine the coordinates of each sensing node to generate a set of acoustic pulse arrival time distributions. The set of acoustic pulse arrival time distributions includes the acoustic propagation delay corresponding to the distance between each sensing node and the underwater target.

[0027] For example, in step S12171, select any three arrival time points of sensing nodes from the set of arrival time points of the sound wave pulse, calculate the arrival time difference between each pair, and obtain a first time difference, a second time difference, and a third time difference. The first time difference is the arrival time difference between the first sensing node and the second sensing node, the second time difference is the arrival time difference between the first sensing node and the third sensing node, and the third time difference is the arrival time difference between the second sensing node and the third sensing node.

[0028] Specifically, assuming that a certain pulse sequence has been confirmed to originate from the same underwater target from the set of acoustic pulse arrival times generated in step S1213, three nodes are selected from the numerous sensing nodes that received the pulse, denoted as node A, node B, and node C. Their precise times of receiving the pulse are recorded as t_A, t_B, and t_C. Using the arrival time of node A as a reference, the arrival time difference between node B and node A is calculated, i.e., Δt_AB = t_B - t_A, and this is defined as the first time difference. The arrival time difference between node C and node A is calculated, i.e., Δt_AC = t_C - t_A, and this is defined as the second time difference. The arrival time difference between node C and node B is calculated, i.e., Δt_BC = t_C - t_B, and this is defined as the third time difference. These three time differences constitute the basic observation data for subsequent positioning calculations.

[0029] Step S12172: Based on the first time difference and the propagation speed of sound waves in water, calculate the distance difference between the first sensing node and the second sensing node to the underwater target to obtain the first distance difference. Similarly, based on the second time difference, calculate the distance difference between the first sensing node and the third sensing node to the underwater target to obtain the second distance difference.

[0030] Specifically, the speed of sound in water is known to be a constant c. Multiplying the first time difference Δt_AB obtained in step S12171 by the speed of sound c yields the distance difference between node A and node B to the underwater target, denoted as Δd_AB = c × Δt_AB. This distance difference is an algebraic value, and its sign depends on which node is closer to the target. Similarly, multiplying the second time difference Δt_AC by the speed of sound c yields the distance difference between node A and node C to the target, denoted as Δd_AC = c × Δt_AC. These two distance differences will serve as constant terms in constructing the hyperbolic (surface) equation.

[0031] Step S12173: Using the spatial coordinates of the first sensing node and the second sensing node as the focus and the first distance difference as a constant, construct the first hyperbolic surface equation; using the spatial coordinates of the first sensing node and the third sensing node as the focus and the second distance difference as a constant, construct the second hyperbolic surface equation.

[0032] Specifically, let the spatial coordinates of node A be (x_A, y_A, z_A), and the spatial coordinates of node B be (x_B, y_B, z_B). For any point P(x, y, z) in space, its distance to node A is d_A = √[(x-x_A)² + (y-y_A)² + (z-z_A)²], and its distance to node B is d_B = √[(x-x_B)² + (y-y_B)² + (z-z_B)²]. The first hyperbolic surface equation is defined as the set of all points P that satisfy d_A - d_B = Δd_AB. This is a hyperboloid of revolution. Similarly, with nodes A and C as foci, the set of all points P that satisfy d_A - d_C = Δd_AC constitutes the second hyperbolic surface equation.

[0033] Step S12174: Solve for the intersection of the first hyperbolic surface equation and the second hyperbolic surface equation to obtain two candidate positioning points, which are located on both sides of the plane where the focus is located.

[0034] Specifically, the two hyperboloid equations constructed in step S12173 are solved simultaneously. This is a nonlinear system of equations. Numerical methods, such as the Newton-Raphson iterative method, can be used to iterate from an initial guess point. Since the intersection of two hyperboloids in space usually produces a curve, but given sufficient constraints, their intersection is usually two points. These two points are symmetric about a plane formed by three focal nodes, and are located on opposite sides of this plane. The coordinates of these two points are calculated and used as candidate location points for the underwater target.

[0035] Step S12175: Based on the arrival time of the fourth sensing node, calculate the arrival time difference between the first sensing node and the fourth sensing node to obtain the third time difference, and then obtain the distance difference between the first sensing node and the fourth sensing node to the underwater target as the third distance difference.

[0036] Specifically, a fourth sensing node, denoted as node D, is introduced, with spatial coordinates (x_D, y_D, z_D). The time at which it receives the same pulse is t_D. The arrival time difference between node D and node A, Δt_AD = t_D - t_A, is calculated. Multiplying this by the speed of sound c, we obtain the distance difference between node A and node D from the target, Δd_AD = c × Δt_AD, which is defined as the third distance difference.

[0037] Step S12176: Using the spatial coordinates of the first sensing node and the fourth sensing node as the focus and the third distance difference as a constant, construct a third hyperbolic surface equation. Substitute the two candidate positioning points into the third hyperbolic surface equation for verification. The candidate positioning points that meet the accuracy requirements of the equation are determined as the preliminary acoustic positioning coordinates of the underwater target.

[0038] Specifically, with nodes A and D as foci and Δd_AD as a constant, a third hyperbolic surface equation is constructed, which is the set of points P that satisfy d_A - d_D = Δd_AD. The coordinates of the two candidate positioning points obtained in step S12174 are substituted into the equation d_A - d_D to calculate the left-hand side value. The deviation of this value from Δd_AD is calculated. The candidate point with the smaller absolute value of the deviation is selected as the possible location of the true target. If the deviations of both points are large, it may be necessary to reselect nodes or consider positioning errors. The verified candidate point is determined as the preliminary acoustic positioning coordinates of the underwater target.

[0039] Step S12177: Iteratively optimize the preliminary acoustic positioning coordinates by using the arrival time of all available sensing nodes and employing the least squares method to iteratively adjust the preliminary acoustic positioning coordinates, so as to minimize the sum of squared residuals between the theoretical arrival time and the actual arrival time of all sensing nodes, thereby obtaining the optimized acoustic positioning coordinates.

[0040] Specifically, the preliminary acoustic positioning coordinates obtained in step S12176 are used as initial values. Now, all sensing nodes that received the pulse are utilized (assuming there are N nodes, N>3). For each node i, its theoretical arrival time t_theory_i can be calculated by dividing the distance between the initial positioning coordinates and the node coordinates by the speed of sound c. The actual arrival time is t_actual_i. The residual r_i is defined as t_theory_i - t_actual_i. The objective function is constructed as the sum of squares of the residuals of all nodes, i.e., F=Σ(r_i)². Then, a nonlinear least squares optimization algorithm, such as the Gauss-Newton method or the Levenberg-Marquardt method, is used to iteratively adjust the target coordinates. In each iteration, the correction amount of the coordinates is calculated based on the residual vector and the Jacobian matrix, the target coordinates are updated, and the residuals are recalculated. This process is repeated until the decrease in the objective function is less than a preset threshold or the number of iterations reaches the upper limit. The target coordinates obtained at this time are the optimized acoustic positioning coordinates, which have higher accuracy than the preliminary positioning results using only four nodes.

[0041] Step S12178: Spatial registration is performed between the optimized acoustic positioning coordinates and the optical target contour position detected in the optical imaging image frame set at the same moment. If the registration error is less than a preset error threshold, the optimized acoustic positioning coordinates are output as the final acoustic positioning coordinates of the underwater target.

[0042] Specifically, after obtaining the optimized acoustic positioning coordinates, optical information is needed for final confirmation. From the optical target contour sequence obtained in step S122, the optical image frame closest to the arrival time of the current acoustic pulse is identified. In this frame, the optical target contour has been detected, and the coordinates of the optical target in three-dimensional space (or a line-of-sight ray emanating from the camera) can be calculated using stereo vision or known target depth information. The optimized acoustic positioning coordinates are projected onto this optical image, or conversely, the spatial coordinates of the optical target are compared with the acoustic positioning coordinates. The Euclidean distance between the two is calculated as the spatial registration error. A registration error threshold is set. If the calculated registration error is less than the threshold, it indicates that the acoustic positioning is highly consistent with the optical observation, and the positioning result is highly reliable and can be used as the final acoustic positioning coordinates output for the underwater target. If the error is too large, it is necessary to check each step of the acoustic positioning or optical detection, and it may trigger a repositioning or data fusion process. This step utilizes the high angular resolution characteristics of optical data to verify and confirm the acoustic positioning results.

[0043] Step S1218: The pulse energy concentration parameter, the pulse broadening parameter, and the pulse energy attenuation gradient parameter are associated and stored according to each pulse to form the waveform distortion coefficient set. Each coefficient vector in the waveform distortion coefficient set corresponds one-to-one with the arrival time point of the corresponding sound wave pulse.

[0044] Specifically, for each successfully detected and located acoustic pulse, the pulse energy concentration parameter calculated in step S1214, the pulse broadening parameter calculated in step S1215, and the pulse energy attenuation gradient parameter calculated in step S1216 need to be integrated. These three parameters are combined into a three-dimensional coefficient vector. The first component of this three-dimensional coefficient vector is the pulse energy concentration parameter, the second component is the pulse broadening parameter, and the third component is the pulse energy attenuation gradient parameter. A unique identifier is assigned to each of these coefficient vectors, and the arrival time of the corresponding acoustic pulse and the identifier of the sensing node from which the pulse originates are recorded. Then, the coefficient vectors of all pulses, their corresponding time points, and node identifiers are organized into a data structure in chronological order. This data structure can be in list form, where each element is a mapping entry, with the pulse arrival time as the key and a structure containing the node identifier and the three-dimensional coefficient vector as the value. This forms a set of waveform distortion coefficients. This set of waveform distortion coefficients ensures that the quantization characteristics of the waveform distortion of each received acoustic pulse can precisely correspond to the specific receiving time and receiving location.

[0045] Step S122: Perform target edge detection processing on the optical imaging image frame set to obtain an optical target contour sequence and an optical target contour motion offset set. The optical target contour sequence contains a pixel coordinate chain of the same optical target contour in consecutive image frames, and the optical target contour motion offset set contains the contour centroid displacement vector between adjacent image frames.

[0046] Specifically, each frame of the optical imaging image is read sequentially from the set of frames. For each frame, preprocessing is performed first, including noise reduction and contrast enhancement, to improve the sharpness of the target edges. Then, an edge detection algorithm, such as the Canny or Sobel operator, is applied to calculate the gradient magnitude and direction of each pixel in the image. Through non-maximum suppression and double thresholding, pixels with drastic brightness changes are extracted, forming a binarized edge image. Connectivity analysis is performed on this edge image to connect adjacent edge pixels into a contour. Based on the target's prior shape or motion features, the optical target contour belonging to the underwater target is identified from all detected contours. The image coordinates of all pixels on the identified contour are recorded in traversal order, forming the pixel coordinate chain of the optical target contour for that frame. For consecutive frames, a target matching algorithm, such as based on contour shape similarity or centroid position prediction, is used to determine contours belonging to the same physical target in different frames, and their pixel coordinate chains are arranged in chronological order to form an optical target contour sequence. Simultaneously, the centroid displacement of the target contour between adjacent image frames is calculated. First, calculate the centroid coordinates of the current frame's contour, which is the average of the coordinates of all pixels within the contour. Then, calculate the centroid coordinates of the same target contour in the next frame. The vector difference between the two centroid coordinates, including direction and magnitude, is the contour centroid displacement vector. Repeat this calculation for all adjacent frames to obtain a set of displacement vectors arranged in chronological order, forming a set of optical target contour motion offsets.

[0047] Step S123: Perform anomaly peak localization processing on the magnetic anomaly change time series record set to obtain a set of magnetic anomaly source location coordinates and a set of magnetic anomaly intensity change gradients. The set of magnetic anomaly source location coordinates contains the coordinate estimates of the magnetic anomaly sources detected by each sensing node in three-dimensional space, and the set of magnetic anomaly intensity change gradients contains the rate of change parameter of magnetic anomaly intensity over time.

[0048] Specifically, the magnetic anomaly change time series record set consists of magnetic field intensity time series recorded by multiple sensing nodes. For the time series record of each sensing node, preprocessing is first performed, including removing the geomagnetic background field and eliminating diurnal variation interference, to obtain a pure magnetic anomaly signal time series. Peak detection is performed on this series to identify the local maxima of the magnetic anomaly signal and their corresponding times. When multiple nodes detect magnetic anomaly peaks at similar times, these peaks are considered to be caused by the same magnetic anomaly source. Using the peak times and spatial coordinates of these nodes, hyperbolic or intersection localization methods similar to acoustic localization are used to estimate the spatial location of the magnetic anomaly source. Since the magnetic anomaly signal intensity decays with distance, the anomaly amplitude measured by multiple nodes can be combined to obtain a more accurate estimate of the magnetic anomaly source location coordinates by fitting a field source model. All coordinates located by the same source are summarized to form a set of magnetic anomaly source location coordinates. Simultaneously, a first-order difference calculation is performed on the magnetic anomaly signal time series of each node, that is, the change in anomaly intensity at adjacent time points is calculated and divided by the time interval to obtain the instantaneous rate of change. After smoothing the rate of change sequence, the rate of change values ​​corresponding to the time periods when the magnetic anomaly source appears are extracted, especially the rate values ​​at the rising and falling edges of the anomaly peak, forming a set of magnetic anomaly intensity change gradients. Each parameter in this set of magnetic anomaly intensity change gradients describes the rate of change of the magnetic anomaly intensity over time.

[0049] Step S124: Based on the acoustic pulse arrival time distribution set and the waveform distortion coefficient set, combined with the optical target contour sequence and the optical target contour motion offset set, perform acoustic-optical target association probability calculation processing to generate an acoustic-optical target association probability distribution map. The acoustic-optical target association probability distribution map contains the probability values ​​that the acoustic features and optical features at different spatial locations belong to the same underwater target.

[0050] Specifically, the spatial region is divided into a fine three-dimensional grid. For each grid point, firstly, based on the set of acoustic pulse arrival time distributions, the theoretical acoustic propagation delay from a target at that grid point to each sensing node is calculated. This theoretical delay is compared with the actual detected acoustic pulse arrival times to calculate the degree of matching. A higher degree of matching indicates a greater likelihood of an acoustic target at that grid point. Simultaneously, by incorporating parameters from the waveform distortion coefficient set, such as pulse broadening and energy attenuation, this probability can be further corrected, as different acoustic propagation paths at different locations lead to different distortion characteristics. Secondly, based on the optical target contour sequence, the target contour on the image plane is back-projected into three-dimensional space using a camera imaging model, forming a cone-shaped region traversing the contour from the camera. If a grid point falls within this cone-shaped region, it is likely to correspond to a target in the optical image. Simultaneously, by combining the set of optical target contour motion offsets, the target's direction and velocity in space can be inferred, thus constraining the motion state of the grid points. The degree to which each grid point simultaneously satisfies both acoustic and optical constraints is quantified. Probabilistic methods, such as Bayes' theorem, can be used. The acoustic matching degree can be used as the likelihood probability, and the existence probability of the optical back-projection region can be used as the prior probability. The posterior probability of each grid point simultaneously generating observed acoustic and optical features can then be calculated. This posterior probability is used as the acoustic-optical association probability value for that grid point. After traversing all grid points, a three-dimensional acoustic-optical target association probability distribution map is obtained. Each point on this map is marked with a probability value between zero and one, representing the probability that the location is simultaneously detected by acoustic and optical means and belongs to the same target.

[0051] Step S125: Based on the set of magnetic anomaly source location coordinates and the set of magnetic anomaly intensity change gradients, perform spatial location matching processing with the acoustic-optical target association probability distribution map to determine the same underwater target object that is pointed to by multiple modes, and obtain the multimodal identifier of the underwater target object. The multimodal identifier is used to uniquely index the feature records of the underwater target object in different modal data.

[0052] Specifically, each coordinate point in the set of magnetic anomaly source location coordinates generated in step S123 is mapped to the acoustic-optical target association probability distribution map generated in step S124. The acoustic-optical association probability value of the grid where the coordinate point is located is checked. If the probability value is higher than a preset high confidence threshold, it strongly indicates that the magnetic anomaly source and the high-probability target detected by acoustic-optical means are located in the same spatial location. At the same time, the gradient characteristics of the magnetic anomaly intensity change of the magnetic anomaly source are analyzed, such as whether its rate of change is consistent with the magnetic field disturbance characteristics caused by the target's motion velocity inferred by acoustic or optical means. If the two are highly overlapping in space and synchronized in time (i.e., the time of the magnetic anomaly occurrence is consistent with the time of the target detection by acoustic-optical means), it is determined that the acoustic, optical, and magnetic anomaly modal data all point to the same underwater target. A globally unique identifier, i.e., a multimodal identifier, is assigned to the jointly confirmed underwater target. This multimodal identifier can be a string composed of a timestamp and a random number. Afterwards, in all the original data records of all modalities, all data belonging to this target is tagged with this multimodal identifier. For example, a field can be added to the record of the set of arrival times of acoustic pulses to store this identifier; a similar identifier field can be added to the record of the optical target contour sequence. In this way, the multimodal identifier becomes a bridge for indexing all feature records of the same target in different modalities.

[0053] Step S126: Extract the acoustic subset corresponding to the underwater target in the acoustic wave propagation waveform sampling set, the optical image subset corresponding to the optical imaging image frame set, and the magnetic anomaly time series subset corresponding to the magnetic anomaly change time series record set, and construct the feature identification mapping chain. The feature identification mapping chain includes the correspondence between each modal feature and the underwater target and the mutual index pointers between each modal feature.

[0054] Specifically, based on the multimodal identifier determined in step S125, all data associated with the identifier are extracted back from the original multimodal sensing stream. From the acoustic wave propagation waveform sampling set, all arrival times of acoustic pulses marked with the identifier, and the corresponding original waveform segments, are identified and organized into an acoustic subset. From the optical imaging image frame set, all image frames marked with the identifier, and the corresponding target contour pixel coordinate chain in each frame, are identified and organized into an optical image subset. From the magnetic anomaly change time-series record set, magnetic anomaly signal segments corresponding to the target's occurrence time window and matching its spatial location, along with the corresponding magnetic anomaly source location coordinates and intensity change gradient, are identified and organized into a magnetic anomaly time-series subset. A data structure called the feature identifier mapping chain is constructed. The core of this feature identifier mapping chain is an association matrix or graph structure. The central node of the graph is the multimodal identifier of the underwater target. Starting from the central node, three edges point to the storage addresses of the acoustic subset, the optical image subset, and the magnetic anomaly time-series subset, respectively. Simultaneously, cross-indexing is established between data entries in these three subsets based on timestamps and spatial coordinates. For example, a record in the acoustic wave subset may contain a pointer to the identifier of the optical image frame closest to the moment of the acoustic pulse. Similarly, the record of the optical image frame may contain a pointer to the identifier of a magnetic anomaly signal detected before or after the moment the image was acquired. In this way, the feature identifier mapping chain not only establishes the attribution relationship between modes and targets but also establishes direct navigational links between feature data from different modes.

[0055] Step S127: Perform multipath time delay analysis on the set of arrival time distributions of the sound wave pulses, identify the reflected echoes caused by the water stratification interface in the sound wave propagation path, and generate the set of water stratification interface distribution coordinates as the boundary parameters of the first type of shielding object.

[0056] Specifically, for an acoustic pulse radiated from the same underwater target, its arrival time distribution set may include not only the arrival time of the direct wave but also the arrival time of reflected waves after reflection from the seabed, sea surface, or thermocline within the water body. The pulse arrival times recorded by each sensing node are sorted chronologically. The earliest arriving pulse is usually the direct wave. Subsequent pulses are potential reflected waves, i.e., multipath signals. The arrival time difference between each subsequent pulse and the direct wave is calculated. Using ray acoustics theory, a possible layered interface location is assumed, and the theoretical propagation time of the acoustic wave from the target, reflected from this interface, to the receiving node is calculated. Through iterative search, one or more interface locations are found such that, for multiple receiving nodes, the theoretically calculated arrival time difference of the reflected wave best matches the actual observed sequence of arrival time differences of multiple subsequent pulses. This process can be achieved by constructing a cost function, such as the mean square error between the theoretical and actual time delays, and using optimization algorithms (such as genetic algorithms or particle swarm optimization) to solve for the optimal interface location. Finally, the reflection interfaces that can explain the observation multipath delay are identified, and their three-dimensional spatial coordinates (such as the coordinates of several points on the interface or the fitted plane equation parameters) are recorded to generate a set of coordinates of the water body layer interface distribution. These interfaces will become shields that block or reflect sound wave signals in subsequent communications, and therefore serve as the boundary parameters of the first type of shield.

[0057] Step S128: Perform dark area range detection processing on the pixel region behind the underwater target in the optical imaging image frame set to generate a set of optical occlusion region boundary coordinates as the boundary parameters of the second type of occlusion.

[0058] Specifically, in optical image frames containing underwater targets, the region behind the target's outline is analyzed. Since illumination typically comes from above or in front, the target itself blocks light, creating a dark area behind it that is significantly less bright than the surrounding environment—a shadow. For each frame, image segmentation techniques are first used to precisely segment the target outline. Then, using the target outline as the boundary, region growing is performed in the direction away from the light source (estimated based on the image brightness gradient) to identify all consecutive pixel regions with brightness below the background brightness threshold; these regions are potential optical occlusion areas. Boundary extraction is then performed on these regions to obtain a series of boundary pixel coordinates describing the shape of the dark area. Since a single image only provides two-dimensional information, to obtain the boundary of the occlusion in three-dimensional space, it is necessary to combine multiple images taken from different angles, or to perform triangulation using the geometric relationship between the target, the light source, and the camera. By using a multi-view geometry method, the two-dimensional boundary coordinates of the same shadow area detected in multiple images are matched and three-dimensionally reconstructed to obtain the set of coordinate points of the shadow area in three-dimensional space. The set of coordinate points constitutes the boundary of the optical occlusion area, reflecting the spatial range of the target being occluded. Although objects in the above area (such as reefs or shipwrecks behind the target) may not be directly illuminated, they will block communication signals. Therefore, the above three-dimensional coordinate set is used as the boundary parameter of the second type of occlusion.

[0059] Step S129: Merge the set of coordinates of the water body layer interface distribution and the set of coordinates of the optical occlusion region boundary, perform boundary overlap region elimination and boundary continuity enhancement processing, and generate the occlusion distribution boundary parameters. The occlusion distribution boundary parameters include the closed boundary surface descriptor of the occlusion in three-dimensional space.

[0060] Specifically, the coordinate set of water stratification interface distribution (first type of shading boundary) obtained in step S127 and the coordinate set of optical shading region boundary (second type of shading boundary) obtained in step S128 are spatially registered and fused. First, the two types of boundary point cloud data are placed in the same three-dimensional coordinate system. Spatial overlap detection is performed. If a region is found to be covered by both types of boundaries simultaneously, for example, an optical shadow region happens to fall on an acoustic reflection interface, then it is necessary to discard or use a weighted average based on the confidence level of the data source to eliminate overlapping regions and avoid duplicate counting. Then, the fused discrete boundary point cloud is processed for surface reconstruction. The Poisson surface reconstruction algorithm or the moving cube algorithm can be used to fit the discrete point cloud into a continuous, closed triangular mesh surface. During the fitting process, for regions with sparse boundaries or voids, interpolation or deduction based on a physical model (e.g., assuming that the surface of the shading object has a certain smoothness) is used to enhance the boundary continuity, fill the voids, and make the final generated surface descriptor continuous and closed. The closed boundary surface descriptor exists in the form of a three-dimensional mesh model. Each mesh vertex is represented by three-dimensional coordinates, and the mesh patch is composed of three vertex indices. It completely outlines the spatial shape and location of all physical obstacles in the underwater environment that can hinder or significantly affect signal propagation. The above closed boundary surface descriptor constitutes the boundary parameters of the shielding object distribution.

[0061] Step S1210: Based on the changes in the optical target contour motion offset set and the magnetic anomaly source position coordinate set over time, calculate the instantaneous motion velocity vector and motion direction angle change rate of the underwater target, and generate the motion trend parameters; associate and encapsulate the motion trend parameters with the cover distribution boundary parameters to generate the spatial topology descriptor, which contains mutually independent but related motion trend parameter components and cover distribution boundary parameter components.

[0062] Specifically, the pixel displacement of the target on the image plane can be extracted from the set of optical target contour motion offsets obtained in step S122. Combining the camera's intrinsic and extrinsic parameters (i.e., position and attitude), the pixel displacement is converted into the actual displacement of the target in three-dimensional space through coordinate transformation. Simultaneously, the change in the magnetic anomaly source position coordinates obtained in step S123 over time can also directly yield the target's position change in three-dimensional space. By fusing the displacement information from these two sources and employing data fusion algorithms such as Kalman filtering or particle filtering, the optimal estimate of the underwater target's motion state is obtained. The difference between the fused position coordinates of adjacent time points is calculated and divided by the time interval to obtain the instantaneous velocity vector, which includes the velocity magnitude (rate) and direction. The velocity direction at multiple consecutive moments is differentially calculated to obtain the rate of change of the motion direction angle over time, i.e., the turning angular velocity. The instantaneous velocity vector and the rate of change of the motion direction angle are encapsulated together to form a motion trend parameter describing the target's dynamic characteristics. Subsequently, this motion trend parameter is associated and encapsulated with the occlusion distribution boundary parameter (i.e., the closed boundary surface descriptor) generated in step S129. The aforementioned association is not a simple data splicing, but rather the establishment of a composite data structure. In this structure, motion trend parameters serve as one independent data component, and occlusion distribution boundary parameters serve as another independent data component. These two components are linked through a common scene identifier or timestamp. In this way, the spatial topology descriptor contains information on two key environmental elements: "moving objects" and "stationary obstacles," and clearly defines the relative spatial relationship between them.

[0063] Step S130: Identify the motion trend parameters and cover distribution boundary parameters of the underwater target body according to the spatial topology descriptor, and construct a dynamic channel occupancy state evolution model between the underwater target body and each communication node based on the motion trend parameters and the cover distribution boundary parameters.

[0064] Step S131: Parse the spatial topology descriptor and extract the current three-dimensional spatial coordinate sequence of the underwater target and the closed boundary surface descriptor in the distribution boundary parameters of the cover. The current three-dimensional spatial coordinate sequence contains the spatial position coordinates of the underwater target at continuous time points.

[0065] Specifically, from the spatial topology descriptor generated in step S1210, the motion trend parameter component is first parsed. From this component, a series of three-dimensional spatial coordinates of the underwater target, recorded at fixed sampling intervals over a past period, are extracted. These coordinates are arranged in chronological order to form the current three-dimensional spatial coordinate sequence. Each element in the sequence contains a timestamp and its corresponding coordinate value. Simultaneously, from another component of the spatial topology descriptor, the boundary parameters of the obstruction distribution are extracted.

[0066] Step S132: Perform trajectory fitting processing on the current three-dimensional spatial coordinate sequence to generate the predicted motion trajectory curve of the underwater target. The predicted motion trajectory curve includes the set of spatial points that the underwater target may pass through in the future time period and the corresponding time labels.

[0067] Specifically, the current three-dimensional spatial coordinate sequence extracted in step S131 is used as the observation data points. Curve fitting or state estimation methods are employed, such as polynomial fitting, spline interpolation, or motion model prediction based on Kalman filtering. If polynomial fitting is used, it is assumed that the target's trajectory can be described by a polynomial function of time, for example, the three-dimensional coordinates are represented as cubic polynomials of time. The coefficients of the polynomial are solved using the least squares method to make the fitted curve as close as possible to the historical observation points. After obtaining the curve equation, this equation is extrapolated to future times to calculate the coordinate values ​​at a series of future time points. If Kalman filtering is used, after constructing the target's motion state equation (such as a uniform or uniformly accelerated model) and observation equation, the current state is updated using the historical observation sequence, and the state is iteratively predicted forward for multiple future steps to obtain the predicted position. The final generated predicted motion trajectory curve is a data structure containing a list of future time points, each time point corresponding to a spatial coordinate obtained by the prediction algorithm, and may include a covariance matrix representing the prediction uncertainty.

[0068] Step S133: Based on the closed boundary surface descriptor, the water environment around the underwater target is spatially rasterized to generate an environmental raster model containing multiple three-dimensional raster units. Each three-dimensional raster unit in the environmental raster model is marked with a shading attribute identifier.

[0069] Specifically, a three-dimensional bounding box is defined centered on the current position of the underwater target or based on the entire area of ​​interest. This bounding box is divided into a series of identical cubic units, or 3D grid units, along its length, width, and height according to a preset grid size. Each 3D grid unit is uniquely identified by its center point coordinates or index coordinates. Each 3D grid unit is traversed, and its spatial intersection with the closed boundary surface descriptor extracted in step S131 is determined. It is then determined whether the spatial cube represented by the grid unit overlaps or intersects with the triangular mesh model describing the occlusion. If the grid unit is completely inside any occlusion surface, it is marked as "completely occluded." If the grid unit partially intersects with an occlusion, the proportion of the intersection volume to the total volume of the grid unit is calculated. If the proportion exceeds a certain threshold, it is marked as "partially occluded," and an optional transmission coefficient is recorded, which is inversely proportional to the intersection volume. If the grid unit has no intersection with any occlusion, it is marked as "unoccluded." After traversing all grid cells, a three-dimensional environment grid model is obtained. This environment grid model is actually a three-dimensional array, where each element corresponds to a grid cell and stores the occlusion attribute identifier of that cell.

[0070] Step S134: Perform acoustic ray tracing processing on the straight-line propagation path between each communication node and the underwater target. Combine the occlusion attribute identifiers of each three-dimensional grid unit in the environmental grid model to determine the length of the obstruction segment and the position coordinates of the obstruction segment of each straight-line propagation path, and generate a path obstruction parameter set.

[0071] Specifically, for each communication node, its fixed spatial coordinates are known, along with the predicted coordinates of the underwater target at that moment obtained from step S132. A virtual straight line is constructed between these two points, representing the shortest possible propagation path of the sound wave. In the environmental grid model generated in step S133, three-dimensional digital differential analysis or the Bressenham line algorithm is used to calculate the sequence of all three-dimensional grid cells traversed by this straight line path. The obstruction attribute identifier of each grid cell on the path is checked sequentially. If a grid cell marked as "completely obstructed" is encountered, the coordinate range of the cell is recorded, and the path from entering the cell to leaving the cell is marked as a completely blocked segment. If a grid cell marked as "partially obstructed" is encountered, the path segment is marked as a partially blocked segment based on its transmission coefficient, and the attenuation coefficient is recorded. The sum of the lengths of all completely blocked segments on the entire path is calculated, along with the start and end coordinates of each blocked segment (which can be calculated from the boundaries of the grid cells it traverses). All this information, including the total number of blocking segments, the length of each segment, its location coordinates, and the attenuation coefficient of some blocking segments, is compiled into a set of path blocking parameters for this communication node.

[0072] Step S135: Based on the predicted motion trajectory curve and the path obstruction parameter set, calculate the straight-line distance change sequence and effective line-of-sight ratio sequence between each communication node and the underwater target at different time points. The effective line-of-sight ratio sequence is the proportion of the unobstructed path length to the total path length.

[0073] Specifically, a series of future discrete time points are selected along the predicted motion trajectory curve generated in step S132. For each time point, the predicted coordinates of the underwater target are obtained. For each communication node, the Euclidean distance between the node coordinates and the predicted target coordinates at that time is calculated to obtain the straight-line distance. The straight-line distances calculated for all these time points are arranged in chronological order to obtain a sequence of straight-line distance changes. Simultaneously, for each time point, the acoustic ray tracing method from step S134 is reused to calculate the total length L_total of the straight-line path from the node to the target at that time, and all completely blocked grid cells on the path are identified, and the sum of the path lengths within these grid cells, L_blocked, is calculated. The effective line-of-sight length L_los is equal to L_total minus L_blocked. The effective line-of-sight ratio R_los is equal to L_los divided by L_total. If there are partially blocked cells on the path, their impact can be first converted to an equivalent blocking length through the transmission coefficient, and then subtracted from the total length. The effective line-of-sight ratios calculated for each time point are arranged in chronological order to obtain a sequence of effective line-of-sight ratios. This effective line-of-sight ratio sequence visually reflects the degree to which the direct line-of-sight path between the target and the fixed node is obstructed by obstacles at various future moments as the target moves.

[0074] Step S136: Perform bidirectional channel detection simulation processing on the candidate communication links between each communication node. Based on the distribution of obstructions in the environmental grid model and the motion of the underwater target, generate an initial channel impulse response estimation sequence for each candidate communication link. The initial channel impulse response estimation sequence includes multipath delay distribution and amplitude attenuation factor.

[0075] Step S1361: Extract the occlusion attribute identifier of each 3D grid cell from the environmental grid model. The occlusion attribute identifier includes three states: complete occlusion, partial occlusion, and no occlusion. Record the transmission coefficient of the partially occluded cell.

[0076] Specifically, before performing link simulation, the environmental grid model constructed in step S133 needs to be loaded first. This environmental grid model is stored in the form of a three-dimensional array. The array is traversed, and for each three-dimensional grid cell, its pre-calculated and stored occlusion attribute identifier is read. This occlusion attribute identifier is an enumerated value that explicitly indicates whether the cell is "completely occluded," "partially occluded," or "unoccluded." For cells identified as "partially occluded," a transmission coefficient value also needs to be read from its associated data structure. This transmission coefficient is a real number between 0 and 1, representing the proportion of energy that can pass through the cell without attenuation. For example, if the transmission coefficient of a partially occluded cell is tau, then the energy of the signal after passing through the cell is multiplied by tau. The aforementioned identifier and coefficient constitute the basic environmental data for sound tracing and channel estimation.

[0077] Step S1362: For each candidate communication link, draw a three-dimensional straight path in the environmental grid model according to the spatial coordinates of the communication nodes at both ends, and determine all the three-dimensional grid cell sequences that the three-dimensional straight path passes through and the traversal length in each cell.

[0078] Specifically, considering all possible pairs of nodes participating in communication, a candidate communication link is formed between each pair. For one of these links, the 3D coordinates of nodes A and B are known. Starting from A and ending at B, a virtual 3D straight line is drawn in the environment grid model loaded in step S1361. A grid traversal algorithm is used to calculate each 3D grid cell traversed sequentially from the starting point to the ending point. The algorithm outputs an ordered list, where each element corresponds to a traversed grid cell, and records the geometric length of the line segment traversing within that grid cell. This length can be obtained by calculating the distance between the two intersection points where the line enters and exits the cell. This yields the "grid cell sequence" corresponding to the candidate communication link and the "traversal length" within each cell.

[0079] Step S1363: For the units with complete occlusion properties in the three-dimensional grid unit sequence, mark their corresponding path segments as signal completely blocked segments; for the units with partial occlusion properties, calculate the signal attenuation factor of the unit based on its transmission coefficient; and for the units with no occlusion properties, mark them as free propagation segments.

[0080] Specifically, the grid cell sequence obtained in step S1362 is traversed. For each grid cell in the sequence, its occlusion attribute identifier is checked. If it is identified as "completely occluded," the path segment corresponding to that cell (whose length has been determined in step S1362) is marked as a "completely blocked signal segment." This means that signal energy is considered to be completely attenuated and cannot penetrate in this path segment. If it is identified as "partially occluded," the transmission coefficient is read from the data of that cell. The attenuation factor of the cell for the signal is defined as the transmission coefficient itself, or more commonly, it is converted into an attenuation value in decibels, for example, the attenuation in decibels is equal to negative ten times the common logarithm of the transmission coefficient. The path segment corresponding to that cell is marked as a "partially attenuated segment," and the length of the segment and the calculated attenuation factor are recorded. If it is identified as "unoccluded," the corresponding path segment is marked as a "free propagation segment," meaning that the signal can propagate freely in this segment without additional attenuation, only following basic propagation losses such as spherical spread.

[0081] Step S1364: Divide the three-dimensional straight path into multiple continuous small segments, each segment corresponding to a three-dimensional grid unit. Add up the lengths of all free propagation segments to obtain the total free propagation distance. Add up the lengths of all partially blocked segments and multiply by the corresponding transmission coefficient to obtain the equivalent free propagation distance. The length of the completely blocked segment is the invalid propagation distance.

[0082] Specifically, following the processing result of step S1363, the entire link path has now been divided into multiple continuous segments, each corresponding to a grid cell, and each segment carries a label (complete block, partial attenuation, or free propagation). For all segments labeled "free propagation segment," their lengths are summed to obtain the total free propagation distance L_free. For all segments labeled "partial attenuation segment," their length is multiplied by the corresponding transmission coefficient to obtain the equivalent free propagation distance contributed by that segment. These products are summed to obtain the total equivalent free propagation distance L_eq. For example, a partially attenuated segment with length d and transmission coefficient tau has an equivalent free propagation distance of d multiplied by tau. For all segments labeled "complete block," their length is the ineffective propagation distance L_blocked, which contributes nothing to signal propagation. In subsequent calculations of the direct path amplitude, L_free plus L_eq will be used as the equivalent propagation distance.

[0083] Step S1365: Based on the total free propagation distance and the equivalent free propagation distance, and using the acoustic wave propagation attenuation model, calculate the amplitude attenuation factor of the direct path, and based on the spatial coordinates of the underwater target, determine the possible reflection path of the underwater target as a reflector, and calculate the length of the reflection path and the arrival time delay.

[0084] Specifically, based on the total free propagation distance L_free and the total equivalent free propagation distance L_eq calculated in step S1364, the total equivalent propagation distance L_total_eq of the direct path can be calculated, which is equal to L_free plus L_eq. Then, the propagation attenuation model of sound waves in water is applied. This propagation attenuation model typically includes two parts: spread loss and absorption loss. The spread loss is generally inversely proportional to the square of the distance, while the absorption loss is exponentially related to the product of the distance and the absorption coefficient. Therefore, the amplitude attenuation factor A_direct of the direct path can be expressed as a function of L_total_eq. For example, A_direct is equal to the amplitude at a certain reference distance divided by L_total_eq and multiplied by an exponential attenuation term related to the frequency and the seawater absorption coefficient. Meanwhile, considering that the underwater target itself may be a good acoustic reflector, possible reflection paths need to be calculated. Using the current position of the target as the reflection point, a path is constructed from communication node A to the target and then to communication node B. The total length L_ref of the path A-target-B is calculated. Dividing the propagation time corresponding to this length by the speed of sound c yields the arrival time delay t_ref of the reflection path, which is equal to L_ref divided by c minus the direct path length L_total_eq divided by c. If multiple reflection paths exist, the geometry and reflection coefficient of the target object must be considered, but in this simplified model, only the case of the target as a point reflection source is considered.

[0085] Step S1366: Track the reflection path and determine whether the reflection path is blocked by other obstructions. If it is not blocked, calculate the amplitude attenuation factor of the reflection path. The amplitude attenuation factor depends on the reflection coefficient of the reflecting surface and the total length of the reflection path.

[0086] Specifically, for the reflection path (A-target-B) constructed in step S1365, the ray tracing method from steps S1362 to S1364 needs to be reapplied. This time, however, the path starts at node A and ends at node B, but passes through the target point, thus the path consists of two straight lines: A to the target and the target to B. Grid cell traversal analysis is performed on these two straight lines to identify whether they are blocked by obstructions. It is particularly important to note that the target point's position may change with movement, therefore the obstruction status of the reflection path is dynamic. If any segment of the reflection path is marked as "completely blocked," the reflection path is considered blocked and ignored in the channel impulse response. If it is not blocked, the total equivalent propagation distance L_ref_eq of the reflection path is calculated (again, by accumulating the lengths of the free propagation segment and the partially attenuated segment corrected for the transmission coefficient). Then, the amplitude attenuation factor A_ref of the reflection path is calculated. A_ref depends not only on the propagation distance L_ref_eq but also on the reflectivity R_target of the target surface. Therefore, A_ref can be represented as the A_direct function applied to L_ref_eq and then multiplied by R_target. The reflection coefficient R_target is a value between zero and one, depending on the target material and the incident angle of the sound wave. In this model, it can be preset as a constant or estimated based on the incident angle.

[0087] Step S1367: Sort the direct path and all unblocked reflection paths in ascending order of arrival time delay. Each path corresponds to an impulse response tap. The amplitude of the tap is the amplitude attenuation factor of the corresponding path, and the delay of the tap is the arrival time difference of the path relative to the direct path. Construct the initial channel impulse response estimate.

[0088] Specifically, the direct path calculated in step S1365 and all unobstructed reflection paths calculated in step S1366 are combined. For each path, its total equivalent propagation distance L_eq_path has been calculated. Its absolute propagation time t_abs is equal to L_eq_path divided by the speed of sound c. The absolute propagation time t_abs_min of the path with the shortest propagation time among all paths (usually the direct path) is used as the benchmark. For each path, its arrival time delay delta_t relative to the benchmark path is equal to its absolute propagation time t_abs minus t_abs_min. The amplitude attenuation factor of each path is used as its gain coefficient. Then, the above paths are arranged into a sequence in ascending order of delta_t. This sequence is a multipath channel impulse response model. Each item in the sequence (called a tap) contains two values: time delay delta_t and gain (amplitude attenuation factor). This yields the initial channel impulse response estimate for the candidate communication link at the current time point.

[0089] Step S1368: At each discrete time point, based on the predicted motion trajectory curve of the underwater target, update the spatial position coordinates of the underwater target, recalculate the time delay and amplitude of the reflection path, and whether the direct path changes due to the target's motion, and generate a time-varying channel impulse response estimation sequence.

[0090] Specifically, a series of discrete future time points are taken along the predicted motion trajectory curve generated in step S132. For each time point, steps S1362 to S1367 are repeated. In each iteration, the spatial coordinates of the underwater target are first updated to the predicted coordinates for that time point. Due to the target's movement, its position as a reflection point changes, thus altering the length and obstruction of the reflection path A-target-B. Simultaneously, the target's movement may slightly alter the distance from the node to the target (for some specific links), but this usually has a minor impact. Crucially, for direct links between communication nodes, the target's movement does not directly affect the straight path from node A to B, but it significantly affects the path reflected by the target. Therefore, by recalculating at each future time point, a series of time-evolving channel impulse response estimates can be obtained. Arranging these estimates in chronological order constitutes the initial channel impulse response estimation sequence for the candidate communication link. This initial channel impulse response estimation sequence reflects the dynamic changes in the channel multipath structure caused by the target's movement over a future period.

[0091] Step S1369: Combine the channel impulse response estimates of each candidate communication link at different time points to form the initial channel impulse response estimation sequence. Each element in the initial channel impulse response estimation sequence contains a set of multipath tap coefficients and their corresponding time delay values.

[0092] Specifically, in step S1368, a channel impulse response estimate has been generated for each discrete future time point, such as time points t1, t2, t3, ..., tn. Each of these estimates is a data structure containing an ordered list, where each item corresponds to a distinguishable propagation path (a multipath tap). Each tap item contains two values: one is the arrival time delay of the path relative to the reference path, and the other is the complex gain coefficient or amplitude attenuation factor of the path. Now, it is necessary to integrate the above time-discrete estimates into a sequence. Specifically, a new data structure is created, which is essentially a time-indexed array. The first element of the array corresponds to the channel impulse response estimate for time point t1, the second element corresponds to the estimate for time point t2, and so on, up to time point tn. Each array element stores the complete list of multipath taps for that time. This yields the initial channel impulse response estimate sequence. This initial channel impulse response estimation sequence fully describes how the channel multipath structure of each candidate communication link dynamically evolves over time due to the motion of the underwater target and changes in the environment.

[0093] Step S137: Based on the predicted motion trajectory curve of the underwater target, dynamically update the initial channel impulse response estimation sequence of the candidate communication link to obtain a dynamic evolution sequence of the channel impulse response over time. Each time point in the dynamic evolution sequence of the channel impulse response corresponds to a set of multipath component parameters.

[0094] Specifically, the initial channel impulse response estimation sequence generated in step S1368 already includes dynamic changes in the time dimension, but it is a single prediction based on a grid model and static obstructions. Here, "dynamic update" can be understood as interpolating or smoothing the sequence obtained in step S1368 at a finer time scale, or by incorporating the continuous equations of target motion, to ensure its physical plausibility. For example, assuming step S1368 calculates the channel impulse response at future seconds T1 and T2, then for any time t between T1 and T2, the estimated value at time t can be obtained by linearly interpolating the delay and amplitude of the multipath components at times T1 and T2. This results in a dynamically evolving sequence of the channel impulse response that changes continuously in time. This dynamically evolving sequence of the channel impulse response is essentially a time function; at any given time point tau, it can output a set of multipath component parameters, each including its delay and amplitude.

[0095] Step S138: Couple the channel impulse response dynamic evolution sequence with the transmit power configuration parameters and receive sensitivity parameters of each communication node for analysis and processing to generate the signal-to-noise ratio dynamic change curve of the link between each communication node and the underwater target, as well as the interference power dynamic change curve of the link between each communication node.

[0096] Specifically, for the link between a communication node and an underwater target, the communication quality depends on the strength of the signal transmitted by the target and received by the node after passing through the channel. Let the transmit power of a communication node be P_tx, and its antenna or transducer gain be G_tx. The underwater target, as the receiver, has a receiving sensitivity of S_rx. From the dynamic evolution sequence of the channel impulse response from the node to the target, the amplitude attenuation factor A_direct(t) of the direct path can be extracted. Then, the signal power received by the target, P_rx_target(t), is equal to P_tx multiplied by G_tx and then by the square of A_direct(t). The noise power N is the ambient noise power at the target node. Therefore, the signal-to-noise ratio (SNR) of this uplink, SNR_up(t), is equal to P_rx_target(t) divided by N. Plotting this function over time yields the dynamic change curve of the SNR. Similarly, the downlink (target to node) can be calculated in a similar way. For links between communication nodes, a signal on one link may interfere with the receiver on another link. For example, consider node i sending a signal to node j, while node k is also sending a signal. The useful signal power received by node j from node i is P_useful_ij(t), which is equal to P_tx_i multiplied by G_tx_i multiplied by the square of the direct path amplitude attenuation factor A_ij_direct(t) in the channel impulse response from i to j. The interference signal power received by node j from node k is P_interf_kj(t), which is equal to P_tx_k multiplied by G_tx_k multiplied by the total energy of all paths (possibly including reflections) in the channel impulse response from k to j. This is usually obtained by summing the squares of the amplitudes of all taps in the channel impulse response. The total interference power P_interf_total_j(t) is obtained by summing the interference power generated at node j by all interference sources (other transmitting nodes besides i). Plotting the total interference power against time yields the dynamic variation curve of the interference power at node j. This curve reflects the change in mutual interference between communication nodes over time.

[0097] Step S139: Using the dynamic change curve of signal-to-noise ratio and the dynamic change curve of interference power as independent components, construct the dynamic channel occupancy state evolution model. The dynamic channel occupancy state evolution model includes a channel quality function with time as the independent variable and an inter-node interference function with time as the independent variable.

[0098] Specifically, the two dynamic curves calculated in step S138 are integrated to form a unified model that can be called upon in subsequent steps. This dynamic channel occupancy state evolution model is essentially a data structure or a set of function interfaces containing multiple independent components. The first component is the "channel quality function," which takes time and the specific communication node-target pair (or node-node pair) as input and outputs the predicted signal-to-noise ratio (SNR) of the link. The second component is the "inter-node interference function," which takes time and the specific interfered node as input and outputs the predicted total interference power from all other nodes at that node. Although these two components originate from the same physical process, they provide information in different dimensions. This dynamic channel occupancy state evolution model can be viewed as a digital twin, encapsulating the complex relationships between target movement, environmental obstruction, and channel propagation, and providing quantified, time-varying channel state predictions for the generation of communication guidance commands.

[0099] Step S140: Call the dynamic channel occupancy state evolution model to perform link quality evolution deduction processing on the candidate communication links between each communication node, and obtain the quality fluctuation feature sequence of each candidate communication link in the time dimension and the interference overlap area identifier in the spatial dimension.

[0100] Step S141: Extract the dynamic signal-to-noise ratio (SNR) change curve of the link between each communication node and the underwater target from the dynamic channel occupancy state evolution model, perform peak and valley detection processing on the dynamic SNR change curve, and generate a set of time points with maximum SNR and a set of time points with minimum SNR, which are used as the first fluctuation feature component in the quality fluctuation feature sequence.

[0101] Specifically, from the dynamic channel occupancy state evolution model constructed in step S139, its "channel quality function" interface is called to extract the dynamic signal-to-noise ratio (SNR) change curves within the predicted time range for all uplink or downlink connections between communication nodes and underwater targets. An extreme point detection algorithm is executed for each of these curves. The algorithm calculates the first derivative of the curve and identifies points where the derivative is zero and the second derivative is less than zero as local maxima. The times corresponding to these maxima are recorded, forming a set of SNR maxima time points. Similarly, points where the derivative is zero and the second derivative is greater than zero are identified as local minima, and their times are recorded, forming a set of SNR minima time points. These two sets together describe when the communication quality of the link is at its best (SNR peak) and when it is at its worst (SNR trough), constituting the first component of the link's quality fluctuation characteristic sequence.

[0102] Step S142: Extract the dynamic change curve of interference power of the link between each communication node from the dynamic channel occupancy state evolution model, perform sliding window energy accumulation calculation on the dynamic change curve of interference power, and generate the average interference power value in each time window as the second fluctuation feature component in the quality fluctuation feature sequence.

[0103] Specifically, from the dynamic channel occupancy state evolution model, its "inter-node interference function" interface is invoked to extract the dynamic change curve of the total interference power received by each communication node that may act as a receiver from all other nodes. A fixed-length sliding time window is set, for example, a window length of T_win. This window starts from the prediction start time and slides to the end time with a fixed step size. At each window position, the interference power curve segment falling within the window is integrated or averaged to obtain the average interference power value within that window. The above average values ​​are arranged in order of the window start time to form a sequence. This sequence reflects the magnitude of interference power at different time periods and is another important dimension for measuring link quality, because interference directly reduces the signal-to-noise ratio. This sequence serves as the second component of the quality fluctuation characteristic sequence.

[0104] Step S143: Input the dynamic change curve of signal-to-noise ratio and the dynamic change curve of interference power into the time series alignment module to perform time axis synchronization alignment processing to obtain the aligned signal-to-noise ratio-interference power joint curve. Each time point on the joint curve is marked with both the signal-to-noise ratio value and the interference power value.

[0105] Specifically, the signal-to-noise ratio (SNR) and interference power (EP) curves extracted in steps S141 and S142 may be based on different time sampling rates or different initial phases. For joint analysis, their time axes need to be precisely aligned. A new, unified time axis is created, with the higher of the two time resolutions. Using interpolation methods, such as linear interpolation or spline interpolation, both the SNR and EEP curves are resampled onto this unified time axis. After resampling, for each time point t_i, a corresponding SNR value SNR(t_i) and a corresponding EEP value P_intf(t_i) are obtained. These two values ​​are stored together to form a two-dimensional joint curve. Each point on this joint curve simultaneously carries information about the signal quality and interference level at the current moment.

[0106] Step S144: Perform abrupt change point detection processing on the joint curve, identify the time intervals in which the signal-to-noise ratio and interference power change drastically at the same time, and generate a set of unstable link quality period identifiers, which includes the start and end time points of the unstable period.

[0107] Step S1441: Perform first-order difference calculation on the signal-to-noise ratio curve in the joint curve to obtain the signal-to-noise ratio change rate curve, and perform first-order difference calculation on the interference power curve in the joint curve to obtain the interference power change rate curve.

[0108] Specifically, for the aligned joint curve generated in step S143, the signal-to-noise ratio (SNR) sequence and the interference power sequence are extracted respectively. For the SNR sequence, the difference between the SNR values ​​at adjacent time points is calculated, and then divided by the time interval to obtain the instantaneous rate of change at each time point, forming the SNR rate of change curve. Similarly, the same difference calculation is performed on the interference power sequence to obtain the interference power rate of change curve. These two rate of change curves reflect the severity of the fluctuations in signal quality and interference level over time.

[0109] Step S1442: Align the signal-to-noise ratio change rate curve and the interference power change rate curve by time point. For each time point, determine whether the absolute value of the signal-to-noise ratio change rate is greater than a preset signal-to-noise ratio change rate threshold, and at the same time determine whether the absolute value of the interference power change rate is greater than a preset interference power change rate threshold.

[0110] Specifically, two thresholds are pre-defined: one for the signal-to-noise ratio (SNR) change rate, denoted as Th_SNR_rate; and the other for the interference power change rate, denoted as Th_intf_rate. These two thresholds can be determined based on the statistical characteristics of historical data or the system's stability requirements. At each time point, both change rates are checked simultaneously. If the conditions are met—that is, the absolute value of the SNR change rate is greater than Th_SNR_rate, and the absolute value of the interference power change rate is also greater than Th_intf_rate—then it is considered that both indicators have changed drastically at that time point.

[0111] Step S1443: Mark the time points that simultaneously satisfy the absolute value of the rate of change of signal-to-noise ratio greater than the threshold and the absolute value of the rate of change of interference power greater than the threshold as candidate mutation points, and obtain the time series of candidate mutation points.

[0112] Specifically, all time points are iterated through, and those time points that are judged to satisfy both conditions in step S1442 are selected and recorded in chronological order. The selected time points constitute a preliminary candidate mutation point time series, which indicates that the channel state may have entered an unstable stage.

[0113] Step S1444: Perform continuity analysis on the time series of the candidate mutation points, and merge adjacent candidate mutation points with a time interval less than a preset merging time interval into a mutation time period. Each mutation time period includes a start time and an end time.

[0114] Specifically, a merging time interval threshold, denoted as T_merge, is set. After sorting the candidate mutation point time series obtained in step S1443 by time, the time difference between adjacent points is checked sequentially. If the time difference is less than or equal to T_merge, the two points are considered to belong to the same continuous unstable event, and they, along with all time points in between, are grouped into the same time period. This process continues until a point's time difference with the next point is greater than T_merge, at which point the current time period ends and a new time period begins. Ultimately, each time period is defined by the time of its first point as the start time and the time of its last point as the end time. This results in several mutation time periods.

[0115] Step S1445: Calculate the average signal-to-noise ratio and average interference power within each mutation time period, and compare them with the average signal-to-noise ratio and average interference power of the stable time periods before and after the mutation time period. If the difference exceeds the preset difference threshold, the mutation time period is confirmed as a period of unstable link quality.

[0116] Specifically, to eliminate potential false detections, each abrupt change time period needs to be verified. First, calculate the average signal-to-noise ratio (SNR) SNR_mean_seg and the average interference power P_intf_mean_seg across all time points within the abrupt change time period. Then, take a stable time period before and after the abrupt change time period (e.g., time periods of equal length that do not contain other abrupt change points), and calculate the average SNR SNR_mean_before and SNR_mean_after, and the average interference power P_intf_mean_before and P_intf_mean_after for these two stable time periods respectively. Next, set two difference thresholds: the SNR difference threshold D_SNR and the interference power difference threshold D_intf. If the absolute value of the difference between the mean signal-to-noise ratio during the abrupt change period and the mean of any preceding or following stable period is greater than D_SNR, and the absolute value of the difference between the mean interference power during the abrupt change period and the mean of any preceding or following stable period is greater than D_intf, then the abrupt change period is considered to be caused by a fundamental change in the channel state, and it is confirmed to be a period of unstable link quality.

[0117] Step S1446: Arrange all confirmed unstable link quality periods in chronological order to generate a set of unstable link quality period identifiers. Each entry in the set of unstable link quality period identifiers includes the start time, end time, signal-to-noise ratio fluctuation amplitude, and interference power fluctuation amplitude within the unstable period.

[0118] Specifically, all unstable periods identified in step S1445 are compiled into a list. Each entry in the list corresponds to an unstable period. Each entry contains at least three fields: the start time and end time of the unstable period. Furthermore, to provide richer information, two calculated values ​​can be appended to the entry: the difference between the maximum and minimum signal-to-noise ratio (SNR) values ​​within the period, representing the SNR fluctuation amplitude; and the difference between the maximum and minimum interference power values ​​within the period, representing the interference power fluctuation amplitude. This generates a complete set of link quality unstable period identifiers.

[0119] Step S1447: For each period of unstable link quality, backtrack the predicted motion trajectory of the environmental grid model and the underwater target to identify the physical cause of the instability during that period of unstable link quality, and record the physical cause as an additional attribute in the unstable link quality period identifier.

[0120] Specifically, this is an optional deepening step. For each confirmed unstable period, its time interval can be correlated with the predicted motion trajectory in step S132 and the environmental grid model in step S133. For example, it can be checked whether, during this time period, the underwater target happened to move to the edge of a certain obstruction, causing a drastic change in its path as a reflector; or whether a direct path between nodes was blocked by an obstruction in the newly entering path. Through the above causal backtracking, the physical causes of instability can be identified, such as "target entering the shadow area" or "target leaving the multipath area," and these causes can be recorded as additional attributes in the link quality unstable period identifier in the form of text or enumeration codes. This helps to make more targeted decisions when generating subsequent guidance instructions.

[0121] Step S145: In the spatial dimension, perform three-dimensional spatial mapping processing on the spatial position coordinates of each communication node to obtain the set of node grid coordinates of each communication node in the three-dimensional grid. Combined with the occlusion attribute identifiers of each three-dimensional grid unit in the environmental grid model, generate a node spatial distribution map.

[0122] Specifically, the actual three-dimensional spatial coordinates of each communication node are transformed using the origin and grid size of the environmental grid model defined in step S133. The grid index corresponding to the node coordinates is calculated by subtracting the origin coordinates from the node coordinates, dividing by the grid size, and rounding down to obtain the row, column, and layer index values ​​in the grid model. The grid index values ​​of all nodes constitute the node grid coordinate set. Then, a three-dimensional array with the same dimensions as the environmental grid model in step S133 is created as the basis for the node spatial distribution map. All nodes are traversed, and at their corresponding grid index positions, a communication node is marked within the grid cell, and the node's identifier can be recorded. For grid cells without nodes, the occlusion attribute identifier of that cell in step S133 is copied. In this way, a comprehensive spatial distribution map integrating node positions and environmental occlusion information is generated.

[0123] Step S146: Calculate the sequence of three-dimensional grid cells that the straight path of each candidate communication link passes through on the node spatial distribution map, and extract the number of other candidate communication links that exist simultaneously in each three-dimensional grid cell as the link density parameter of the straight path.

[0124] Specifically, for each candidate communication link (i.e., a pair of nodes), a straight path drawing similar to step S1362 is performed again. On the node spatial distribution map generated in step S145, the sequence of all three-dimensional grid cells traversed by this link from the starting node to the ending node is calculated. For each traversed grid cell in the sequence, the number of other candidate communication link paths that also traverse that cell, besides the endpoint node of the current link itself, is counted. This count can be achieved by pre-calculating the traversing cells of all links. Finally, the "link density parameter" of this link path can be defined as a statistical measure of the link density values ​​of all grid cells traversed by the path, such as the maximum value, average value, or cumulative sum. This parameter intuitively reflects the degree of spatial overlap between this communication link and other links; the greater the overlap, the greater the possibility of future interference in that area.

[0125] Step S147: Mark the three-dimensional grid cells whose link density parameters exceed the preset density threshold as potential interference overlap areas, and count the candidate communication link identifiers contained in each potential interference overlap area to generate the interference overlap area identifier. The interference overlap area identifier includes the area spatial coordinate range and the set of link identifiers within the area.

[0126] Specifically, a spatial link density threshold D_th is set. All three-dimensional grid cells in the environmental grid model are traversed. For each cell, the number of candidate communication links with straight-line paths passing through it (i.e., the link density of the cell calculated in step S146) is calculated. If this number exceeds D_th, the grid cell is marked as a potential interference overlap region. Then, for each marked grid cell, the identifiers of all candidate communication links passing through it are recorded, forming a link identifier set. Simultaneously, the spatial coordinate range of the grid cell is recorded, i.e., its minimum and maximum x, y, and z coordinate values. Thus, each potential interference overlap region corresponds to an information block containing the spatial range and a list of links. All of the above information blocks are summarized to form an interference overlap region identifier. This interference overlap region identifier spatially indicates which areas are hotspots with dense communication links and prone to mutual interference.

[0127] Step S148: Perform connectivity analysis on the potential interference overlapping regions, merge adjacent three-dimensional grid cells into continuous interference overlapping regions, and generate an interference overlapping region connectivity graph. The interference overlapping region connectivity graph includes the boundary coordinate sequence of each interference region and the internal link density distribution heatmap.

[0128] Specifically, the potential interference overlap regions marked in step S147 are discrete grid cells. To obtain continuous interference regions with practical physical meaning, connectivity analysis is required. A three-dimensional connected component marking algorithm is used to traverse all marked grid cells, merging those cells that are coplanar or adjacent in three-dimensional space into the same connected component. Each connected component represents a continuous interference overlap region. For each of the above continuous regions, its boundary is extracted. This can be done by finding all cells belonging to the region but not belonging to its 26-neighborhood that are completely filled, as boundary cells, and connecting the center point coordinates of these boundary cells to form a boundary coordinate sequence. Meanwhile, for each grid cell within the region, its link density value is known. The ratio of the link density within each cell to the average link density of the region can be calculated, or the original link density value can be directly used to generate a pseudo-color map within the region, forming a link density distribution heatmap. Finally, for all continuous interference overlap regions, they are organized into a graph structure, where each region is a node, and if two regions are spatially close or have overlapping links, there is an edge between them. This graph contains the boundary coordinate sequence of each region and the internal heatmap data, forming a connected graph of interfering and overlapping regions.

[0129] Step S149: Associate and map each fluctuation feature component in the quality fluctuation feature sequence with the interference overlap region connectivity graph to generate a correspondence table between the quality fluctuation feature sequence of each candidate communication link in the time dimension and the interference overlap region identifier in the spatial dimension.

[0130] Specifically, the quality fluctuation feature sequences (including temporal extreme points, average interference, etc.) generated in steps S141 and S142 are associated with the interference overlap region connectivity map (spatial interference hotspots) generated in step S148. For each candidate communication link, its quality fluctuation feature sequence already contains its signal-to-noise ratio and interference information at future time points. Simultaneously, from the identifiers in step S147 or S148, it can be determined which spatial interference overlap regions the link's path traverses. A correspondence table is established. The rows of the table are the identifiers of each candidate communication link. The columns of the table can contain two parts: one part is the temporal dimension characteristics of the link, i.e., the storage address or specific value list of its quality fluctuation feature sequence; the other part is the spatial dimension characteristics of the link, i.e., the list of identifiers of all interference overlap regions it traverses, and its link density value within each region.

[0131] Step S150: Based on the quality fluctuation feature sequence and the interference overlap region identifier, generate a set of cooperative communication guidance instructions for the underwater target. The set of cooperative communication guidance instructions includes the transmit power adjustment amount, beam pointing angle adjustment amount, and communication timing scheduling identifier allocated to each communication node.

[0132] Step S151: Analyze the quality fluctuation feature sequence, extract the set of time points with maximum signal-to-noise ratio (SNR), the set of time points with minimum SNR, and the benchmark quality score for each candidate communication link, and construct a link quality priority list. The link quality priority list sorts each candidate communication link from high to low according to the benchmark quality score.

[0133] Specifically, from the data generated in steps S141 to S144, key quality-related information is compiled for each candidate communication link. The sets of time points with maximum and minimum signal-to-noise ratios are directly taken from step S141. The baseline quality score can be a comprehensive indicator, such as by integrating or averaging the dynamic change curve of the link's signal-to-noise ratio over the entire prediction time domain, and then subtracting a penalty term related to the interference level. The magnitude of the penalty term can be proportional to the average interference power value obtained in step S142. After calculating the baseline quality score for each link, all candidate communication links are sorted from high to low according to the score, forming a link quality priority list. Links with higher scores should theoretically receive higher priority in subsequent resource allocation, such as being allocated better time slots or greater power.

[0134] Step S152: Based on the interference overlap region connectivity graph in the interference overlap region identifier, identify the set of candidate communication links located in the same interference overlap region, and generate a link conflict group list. The link conflict group list includes the communication node pair identifier and the corresponding spatial coordinate range of the interference region in each link conflict group.

[0135] Specifically, the interference overlap region connectivity graph generated in step S148 is read. Each consecutive interference overlap region node in the graph is traversed. For each of these regions, the correspondence table generated in step S149 is queried to find the identifiers of all candidate communication links that traverse that region. These identifiers are grouped together to form a link conflict group. All links in this link conflict group share the same physical area in space, therefore they have the potential for strong mutual interference. For each link conflict group, all link identifiers (i.e., communication node pair identifiers) that make up the group are recorded, along with the spatial coordinate range of the interference overlap region (which can be obtained from the connectivity graph nodes).

[0136] Step S153: For each link conflict group in the link conflict group list, perform time-domain peak-shifting scheduling analysis and processing. Based on the set of time points with maximum and minimum signal-to-noise ratio of each candidate communication link, determine the communication time slot allocation scheme that can be staggered on the time axis for different links, and generate a communication timing scheduling identifier. The communication timing scheduling identifier includes the start time and duration of the transmission time slot of each communication node in a continuous time frame.

[0137] Specifically, time-domain scheduling is performed for each conflict group in the link conflict group list. The set of member links of the conflict group is L_set. For each link l in L_set, the set of time points with maximum signal-to-noise ratio (SNR) Peaks_l and the set of time points with minimum SNR Troughs_l are obtained from its quality fluctuation characteristic sequence. The goal of scheduling is to stagger the active communication times of the above links as much as possible to avoid interference caused by them transmitting at the same time. One approach is to construct a time axis that divides the entire predicted communication period into continuous, fixed-length time frames. For each frame, one or more time slots are allocated to each link in the group. The allocation principle is: for link l, its time slots are preferentially allocated near the time points with maximum SNR, while avoiding the time points with minimum SNR. At the same time, it is ensured that at the same time point, only one link in the group is allocated a transmission time slot. This can be achieved by solving a constraint satisfaction problem or using a greedy algorithm. Finally, for each communication node (which may belong to multiple links), a communication timing scheduling identifier can be generated. The communication timing schedule identifier is a timetable that indicates which time slots (start time and duration) in which future time frames the node is authorized to send signals.

[0138] Step S154: According to the link quality priority list, allocate an initial transmit power adjustment amount to each communication node. The initial transmit power adjustment amount is negatively correlated with the baseline quality score of the corresponding link. The higher the baseline quality score, the smaller the transmit power adjustment amount is allocated to the link.

[0139] Specifically, based on the link quality priority list generated in step S151, an initial transmit power adjustment is assigned to the two communication nodes involved in each link. The adjustment refers to an increment or decrement relative to a certain reference power. The allocation logic is as follows: for links with high reference quality scores, it indicates that their channel conditions are already good, so a smaller transmit power adjustment can be assigned, or even the power can be reduced to save energy and reduce interference; for links with low reference quality scores, it indicates that their channel conditions are poor, and a larger transmit power adjustment may be assigned, i.e., the power needs to be increased to ensure basic communication quality. For example, a mapping function can be set to linearly map the reference quality score to a value between zero and the maximum adjustment; the higher the score, the smaller the mapped adjustment.

[0140] Step S155: Input the initial transmit power adjustment amount into the power optimization iterative model, and perform power redistribution processing in conjunction with the link conflict group list. Adjust the transmit power adjustment amount of each communication node in the same link conflict group. The transmit power adjustment amount is the absolute adjustment value of the power. The optimization objective is to minimize the sum of the transmit power adjustment amounts of all links in the same interference overlap area, and obtain the optimized transmit power adjustment amount.

[0141] Specifically, a power optimization iterative model is constructed. The input to this model is the initial transmit power adjustment allocated to each node in step S154, and the list of link conflict groups generated in step S152. The optimization objective of the model is to minimize the sum of transmit power adjustments for all links within the same interference overlap area, while meeting the minimum communication requirements of each link (e.g., ensuring its receiver signal-to-noise ratio is not lower than a certain threshold). Since the power adjustment directly reflects the power magnitude, minimizing their sum is equivalent to minimizing the total transmit power, thereby reducing interference and energy consumption. This is a typical convex optimization problem or linear programming problem. It can be solved using optimization techniques such as iterative water-filling algorithms or dual ascent methods. In each iteration, the power optimization iterative model reallocates the power adjustment based on the interference coupling relationship between links within the conflict group. For example, for links sharing the same area, if one link increases its power, it may cause greater interference to other links, requiring other links to also increase their power to counteract the interference, ultimately potentially increasing the total power instead of decreasing it. The optimization model aims to find a balanced set of power adjustments that minimizes the total power. After iterative convergence, the optimized transmit power adjustment for each node is obtained.

[0142] Step S156: Based on the predicted motion trajectory curve of the underwater target and the spatial coordinates of each communication node, calculate the instantaneous beam pointing angle of each communication node pointing to the underwater target, and generate an initial beam pointing angle set.

[0143] Specifically, for each communication node, its fixed spatial coordinates are known. From the predicted motion trajectory curve in step S132, the predicted coordinates of the underwater target at each future moment can be obtained. For a specific future moment t, the vector pointing from the node's coordinates to the target's coordinates is calculated. This vector is transformed from the Cartesian coordinate system to the spherical coordinate system, yielding two angles: an azimuth angle relative to true north and a pitch angle relative to the horizontal plane. These two angles are the instantaneous beam pointing angles of the node towards the target at moment t. The above calculation is performed for each node, and for a series of future time points, a set of pointing angles varying over time is obtained, forming an initial set of beam pointing angles.

[0144] Step S157: Perform beamwidth adaptive adjustment processing on the initial beam pointing angle set. Based on the distance change sequence between each communication node and the underwater target, determine the beamwidth adjustment factor for each communication node. The beamwidth adjustment factor is positively correlated with the distance, and the adjusted beam pointing angle adjustment amount is obtained.

[0145] For example, in step S1571, the spatial coordinates of each communication node are obtained from the positioning system of each communication node, and combined with the predicted motion trajectory curve of the underwater target, the Euclidean distance between each communication node and the underwater target at each future time point is calculated to generate a distance change sequence.

[0146] Specifically, similar to step S156, but the goal here is not angle, but distance. For each communication node, its own fixed coordinates and the target coordinates predicted in step S132 at various future time points are substituted into the Euclidean distance formula to calculate the straight-line distance at each time point. Arranging these distance values ​​in chronological order yields the distance change sequence for that node.

[0147] Step S1572: Perform smoothing filtering on the distance change sequence to remove measurement noise and obtain a smoothed distance sequence, wherein each distance value in the smoothed distance sequence corresponds to a time point.

[0148] Specifically, due to potential noise in the predicted trajectory, the directly calculated distance sequence may not be smooth. A low-pass filter, such as a moving average filter or a Kalman filter, is applied to the distance change sequence obtained in step S1571. The purpose of filtering is to remove high-frequency jitter caused by inaccurate predictions, resulting in a more stable and smoother distance sequence. Each point in the smoothed sequence still corresponds to a future time point.

[0149] Step S1573: Extract the maximum and minimum distance values ​​from the smoothed distance sequence, calculate the distance dynamic range, and divide the distance into three distance intervals: near distance interval, medium distance interval, and far distance interval based on the distance dynamic range.

[0150] Specifically, identify the maximum value D_max and the minimum value D_min in the smoothed distance sequence. Calculate the distance dynamic range D_range, which is D_max minus D_min. Then, divide the entire distance interval from D_min to D_max into three consecutive sub-intervals, which may be of equal or unequal length. For example, the near-distance interval can be defined as from D_min to D_min plus one-third of D_range, the medium-distance interval as from D_min plus one-third of D_range to D_min plus two-thirds of D_range, and the far-distance interval as from D_min plus two-thirds of D_range to D_max.

[0151] Step S1574: Preset a reference beamwidth value for each distance range. The near-distance range corresponds to a wide beamwidth, the medium-distance range corresponds to a medium beamwidth, and the far-distance range corresponds to a narrow beamwidth. The reference beamwidth value decreases as the distance increases.

[0152] Specifically, according to communication principles, when the distance is short, signal propagation loss is small, and a wider beam can be used to simplify tracking or cover a larger area; when the distance is long, the signal is weak, and a narrow beam is needed to concentrate energy on the target to improve the signal-to-noise ratio. Therefore, different reference beamwidth values ​​are preset for the three intervals divided in step S1573. A larger beamwidth value (wide beam) is preset for the short-range interval, a medium-sized beamwidth value is preset for the medium-range interval, and a smaller beamwidth value (narrow beam) is preset for the long-range interval.

[0153] Step S1575: For each communication node, select the corresponding reference beamwidth as the initial beamwidth according to the distance range it is in at the current moment, and generate the beam pointing angle, which is the azimuth and pitch angle of the communication node pointing to the current position of the underwater target.

[0154] Specifically, for a future time t, the smoothed distance value D_t of the node at that time is first determined. Based on which interval D_t falls within as defined in step S1573, the corresponding reference beamwidth value BW_base is retrieved from step S1574 as the initial beamwidth for that time. Simultaneously, the calculation results from step S156 are reused to obtain the azimuth angle Az_t and elevation angle El_t pointing towards the target at that time. Thus, the initial beam pointing parameters for that time are obtained: center pointing (Az_t, El_t) and beamwidth BW_base.

[0155] Step S1576: Calculate the Doppler frequency shift based on the velocity vector of the relative motion between the communication node and the underwater target. If the Doppler frequency shift exceeds a preset frequency shift threshold, further narrow the beamwidth to improve directivity and compensate for the frequency shift caused by the Doppler effect.

[0156] Specifically, in addition to distance, the radial motion of the target also causes a Doppler frequency shift. From the predicted trajectory in step S132, the radial velocity V_radial of the target relative to the node can be calculated. The Doppler frequency shift f_d is equal to V_radial divided by the wavelength of the sound wave. A Doppler frequency shift threshold Th_doppler is set. If the calculated |f_d| is greater than Th_doppler, it indicates that the relative motion is intense, which may cause a significant frequency shift, leading to demodulation difficulties at the receiver. To compensate for the above effects, the beamwidth can be further narrowed, allowing the energy to be more precisely directed towards the moving target, thereby stabilizing the received signal. Therefore, based on the initial beamwidth BW_base, a coefficient of one is multiplied to obtain the narrowed beamwidth BW_narrowed. The specific value of this coefficient can be related to the degree to which |f_d| exceeds the threshold.

[0157] Step S1577: Compare the narrowed beamwidth with the reference beamwidth and calculate the beamwidth adjustment amount. The beamwidth adjustment amount is expressed as the ratio of the narrowed beamwidth to the reference beamwidth, and is used as the beamwidth adjustment factor.

[0158] Specifically, if a narrowing operation is performed in step S1576, the beamwidth adjustment factor Factor_bw is equal to BW_narrowed divided by BW_base, and this value is less than 1. If no narrowing operation is performed, Factor_bw is equal to 1. This beamwidth adjustment factor quantifies the degree to which additional beam narrowing is required due to the Doppler effect.

[0159] Step S1578: If the communication node is also in a link collision group, adjust the beam pointing angle of this node according to the beam pointing of other nodes in the collision group, so that its main lobe direction avoids the receiving direction of other nodes and reduces mutual interference. This adjustment does not change the beamwidth adjustment factor, but only changes the beam pointing angle value.

[0160] Specifically, check whether the communication node belongs to a link collision group defined in step S152. If so, spatial interference avoidance needs to be considered. Assume the node is transmitting a signal to the target, and its main lobe beam direction is originally precisely pointing to the target. However, at the same time, the receiving antennas (or receiving beams) of other nodes in the collision group may also be located in a similar direction. If the main lobe of the node's transmitting beam happens to be aligned with the receiving direction of another node, it will cause strong interference. Therefore, it is necessary to fine-tune the pointing of the node's transmitting beam within an allowable range, so that its main lobe direction deviates slightly from the receiving direction of other nodes. The premise of the above fine-tuning is that it should not deviate too much from the target, causing the target's received signal to weaken. The amount of fine-tuning can be determined by solving a constrained optimization problem: while ensuring that the target point is still within a certain gain drop tolerance range of the main lobe beam (e.g., the gain drop does not exceed 3dB), maximize the angle between the transmitting beam direction and the receiving direction of other nodes. Finally, a new pointing angle (Az_t_adjusted, El_t_adjusted) is obtained, while the beamwidth adjustment factor Factor_bw remains unchanged.

[0161] Step S1579: Associate and store the beamwidth adjustment factor of each communication node with the corresponding beam pointing angle value to form the adjusted beam pointing angle adjustment amount, which includes the azimuth adjustment value, the elevation adjustment value and the beamwidth scaling factor.

[0162] Specifically, for each communication node, its final determined beam pointing angle adjustment at time t is encapsulated into a data structure containing three parts: the final azimuth angle (possibly adjusted for interference avoidance, Az_t_adjusted), the final elevation angle (El_t_adjusted), and the beamwidth scaling factor (Factor_bw). The azimuth and elevation angles define the direction the beam should point, and the scaling factor defines the beamwidth relative to the reference width. Storing this information for all nodes at all relevant times yields the complete adjusted beam pointing angle adjustment.

[0163] Step S158: Integrate the optimized transmit power adjustment, the adjusted beam pointing angle adjustment, and the communication timing scheduling identifier to generate a preliminary collaborative communication guidance instruction set. The preliminary collaborative communication guidance instruction set includes the target communication node identifier and the corresponding parameter adjustment value for each instruction.

[0164] Specifically, the optimized transmit power adjustment amount for each node obtained in step S155, the adjusted beam pointing angle adjustment amount for each node obtained in step S157, and the communication timing scheduling identifier for each node obtained in step S153 are integrated. For each communication node, these three types of information are combined. For example, node N is instructed at time t to use the transmit power adjustment amount P_adj_t, the beam pointing angle adjustment amount (Az_t, El_t, Factor_bw_t), and to transmit within the time slot identified as S_t. The above information is broken down according to the action of a node at a specific time corresponding to each instruction, forming a series of instruction entries. Each instruction entry explicitly contains the identifier of the target communication node and the parameter values ​​(power, angle, time slot) that need to be adjusted. The set of all these entries is the preliminary cooperative communication guidance instruction set.

[0165] Step S159: Perform instruction conflict detection processing on the preliminary cooperative communication guidance instruction set, check whether there are multiple conflicting transmit power adjustment amounts or beam pointing angle adjustment amounts assigned to the same communication node at the same time. If so, select the instruction with higher priority according to the link quality priority list and retain it. After resolving the conflict, generate the final cooperative communication guidance instruction set.

[0166] Specifically, the initial instruction set is checked for logical contradictions. The most common scenario is that a communication node may participate in multiple links simultaneously, and at the same time, guidance instructions from different links may require the node to perform different power or angle settings. For example, a instruction from link A requires node N to increase power at time t, while a instruction from link B requires node N to decrease power at the same time. This is an instruction conflict. To resolve conflicts, all instructions are iterated through to identify all conflicting instruction groups targeting the same node at the same time. For each conflict group, the link quality priority list from step S151 is consulted to identify the links associated with the conflicting instruction, and the priorities of these links are compared. Then, the instruction corresponding to the link with the highest priority is retained, while other conflicting instructions are discarded or shelved. This ensures that each node receives only one set of explicit, non-contradictory parameter adjustment instructions at any given time. After resolving all conflicts, the remaining instruction set constitutes the final cooperative communication guidance instruction set.

[0167] Step S1510: The final collaborative communication guidance instruction set is encapsulated according to the receiving protocol of each communication node to generate a binary instruction frame sequence containing instruction type field, parameter value field and timestamp field, which is output as a collaborative communication guidance instruction set for the underwater target.

[0168] Specifically, the final generated instructions are high-level logic instructions. To be distributed to each communication node via the underwater communication network, these instructions need to be encapsulated according to the communication protocols supported by each node. A fixed binary format is defined for each instruction type. For example, for the "transmit power adjustment" instruction, the instruction type field can be a fixed byte, indicating that this is a power adjustment command; the parameter value field can be two bytes, indicating the specific power adjustment amount (e.g., an integer in decibels); and the timestamp field can be four bytes, indicating the absolute time when the instruction took effect or the time offset relative to a certain reference. All instructions are converted into the aforementioned fixed-length binary frames one by one. Finally, these binary frames are arranged in chronological order of execution to form a set of cooperative communication guidance instructions.

[0169] In one exemplary embodiment, a multimodal perception-based underwater target cooperative communication guidance system is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 2As shown, the underwater target cooperative communication guidance system based on multimodal perception includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements an underwater target cooperative communication guidance method based on multimodal perception. The display unit is used to generate a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the shell of the underwater target cooperative communication guidance system based on multimodal perception, or an external keyboard, touchpad, or mouse, etc.

[0170] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A method for cooperative communication guidance of underwater targets based on multimodal perception, characterized in that, The method includes: The multimodal sensing stream acquired by the underwater multimodal sensing node array includes a set of acoustic wave propagation waveform samples, a set of optical imaging image frames, and a set of magnetic anomaly change time series records. The multimodal sensing flow is subjected to cross-modal correlation feature parsing processing to generate a spatial topological relationship descriptor that reflects the environmental elements of the target water area. The spatial topological relationship descriptor contains a feature identification mapping chain that points to the same underwater target body in different modal sensing data. The motion trend parameters and cover distribution boundary parameters of the underwater target are identified based on the spatial topology descriptor, and a dynamic channel occupancy state evolution model between the underwater target and each communication node is constructed based on the motion trend parameters and the cover distribution boundary parameters. The dynamic channel occupancy state evolution model is invoked to perform link quality evolution deduction on candidate communication links between each communication node, thereby obtaining the quality fluctuation feature sequence of each candidate communication link in the time dimension and the interference overlap area identifier in the spatial dimension. Based on the quality fluctuation characteristic sequence and the interference overlap region identifier, a set of cooperative communication guidance instructions for the underwater target is generated. The set of cooperative communication guidance instructions includes the transmit power adjustment amount, beam pointing angle adjustment amount, and communication timing scheduling identifier allocated to each communication node.

2. The underwater target cooperative communication guidance method based on multimodal perception according to claim 1, characterized in that, The process of performing cross-modal correlation feature parsing on the multimodal sensing flow to generate a spatial topological relationship descriptor reflecting the environmental elements of the target water area includes: The sound wave propagation waveform sampling set is subjected to pulse feature extraction processing to obtain a sound wave pulse arrival time distribution set and a waveform distortion coefficient set. The sound wave pulse arrival time distribution set includes the time offset of the sound wave pulse received by each sensing node, and the waveform distortion coefficient set includes the pulse broadening parameter and the pulse energy attenuation gradient parameter. The optical imaging image frame set is subjected to target edge detection processing to obtain an optical target contour sequence and an optical target contour motion offset set. The optical target contour sequence contains a pixel coordinate chain of the same optical target contour in consecutive image frames, and the optical target contour motion offset set contains the contour centroid displacement vector between adjacent image frames. The magnetic anomaly change time series record set is subjected to anomaly peak localization processing to obtain a set of magnetic anomaly source location coordinates and a set of magnetic anomaly intensity change gradients. The set of magnetic anomaly source location coordinates contains the coordinate estimates of the magnetic anomaly source detected by each sensing node in three-dimensional space. The set of magnetic anomaly intensity change gradients contains the change rate parameter of magnetic anomaly intensity over time. Based on the set of acoustic pulse arrival time distributions and the set of waveform distortion coefficients, combined with the optical target contour sequence and the set of optical target contour motion offsets, an acoustic-optical target association probability calculation is performed to generate an acoustic-optical target association probability distribution map. The acoustic-optical target association probability distribution map contains the probability values ​​that acoustic features and optical features at different spatial locations belong to the same underwater target. Based on the set of magnetic anomaly source location coordinates and the set of magnetic anomaly intensity change gradients, spatial location matching processing is performed with the acoustic-optical target association probability distribution map to determine the same underwater target object that is pointed to by multiple modes, and the multimodal identifier of the underwater target object is obtained. The multimodal identifier is used to uniquely index the feature records of the underwater target object in different modal data. Extract the acoustic subset corresponding to the underwater target in the acoustic wave propagation waveform sampling set, the optical image subset corresponding to the optical imaging image frame set, and the magnetic anomaly time series subset corresponding to the magnetic anomaly change time series record set, and construct the feature identification mapping chain. The feature identification mapping chain contains the correspondence between each modal feature and the underwater target and the mutual index pointers between each modal feature. Multipath delay analysis is performed on the set of arrival time distributions of the sound pulses to identify reflected echoes caused by water stratification interfaces in the sound propagation path, and a set of water stratification interface distribution coordinates is generated as the boundary parameters of the first type of shielding object. Dark area detection processing is performed on the pixel region behind the underwater target in the optical imaging image frame set to generate a set of optical occlusion region boundary coordinates, which serves as the boundary parameters of the second type of occlusion. By fusing the set of coordinates of the water body layer interface distribution and the set of coordinates of the optical occlusion region boundary, boundary overlap region elimination and boundary continuity enhancement processing are performed to generate the distribution boundary parameters of the occlusion object. The distribution boundary parameters of the occlusion object include the closed boundary surface descriptor of the occlusion object in three-dimensional space. Based on the changes in the optical target contour motion offset set and the magnetic anomaly source position coordinate set over time, the instantaneous motion velocity vector and motion direction angle change rate of the underwater target are calculated to generate the motion trend parameters. The motion trend parameters are then associated and encapsulated with the cover distribution boundary parameters to generate the spatial topology descriptor, which contains mutually independent but related motion trend parameter components and cover distribution boundary parameter components.

3. The underwater target cooperative communication guidance method based on multimodal perception according to claim 1, characterized in that, The step of identifying the motion trend parameters and cover distribution boundary parameters of the underwater target based on the spatial topology descriptor, and constructing a dynamic channel occupancy state evolution model between the underwater target and each communication node based on the motion trend parameters and cover distribution boundary parameters, includes: The spatial topology descriptor is parsed to extract the current three-dimensional spatial coordinate sequence of the underwater target and the closed boundary surface descriptor in the distribution boundary parameters of the cover. The current three-dimensional spatial coordinate sequence contains the spatial position coordinates of the underwater target at continuous time points. The current three-dimensional spatial coordinate sequence is subjected to trajectory fitting processing to generate a predicted motion trajectory curve of the underwater target. The predicted motion trajectory curve includes a set of spatial points that the underwater target may pass through in a future time period and the corresponding time labels. Based on the closed boundary surface descriptor, the water environment around the underwater target is spatially rasterized to generate an environmental raster model containing multiple three-dimensional raster units. Each three-dimensional raster unit in the environmental raster model is marked with a shading attribute identifier. Sound ray tracing processing is performed on the straight-line propagation path between each communication node and the underwater target. Combined with the occlusion attribute identifier of each three-dimensional grid unit in the environmental grid model, the length of the obstruction segment and the position coordinate of the obstruction segment of each straight-line propagation path are determined, and a path obstruction parameter set is generated. Based on the predicted motion trajectory curve and the path obstruction parameter set, calculate the straight-line distance change sequence and effective line-of-sight ratio sequence between each communication node and the underwater target at different time points. The effective line-of-sight ratio sequence is the proportion of the unobstructed path length to the total path length. A bidirectional channel sounding simulation is performed on the candidate communication links between each communication node. Based on the distribution of obstructions in the environmental grid model and the motion of the underwater target, an initial channel impulse response estimation sequence is generated for each candidate communication link. The initial channel impulse response estimation sequence includes multipath delay distribution and amplitude attenuation factor. Based on the predicted motion trajectory curve of the underwater target, the initial channel impulse response estimation sequence of the candidate communication link is dynamically updated to obtain the dynamic evolution sequence of the channel impulse response over time. Each time point in the dynamic evolution sequence of the channel impulse response corresponds to a set of multipath component parameters. The dynamic evolution sequence of the channel impulse response is coupled with the transmit power configuration parameters and receive sensitivity parameters of each communication node for analysis and processing to generate the dynamic change curve of the signal-to-noise ratio of the link between each communication node and the underwater target, as well as the dynamic change curve of the interference power of the link between each communication node. The dynamic channel occupancy state evolution model is constructed by treating the signal-to-noise ratio dynamic change curve and the interference power dynamic change curve as independent components. The dynamic channel occupancy state evolution model includes a channel quality function with time as the independent variable and an inter-node interference function with time as the independent variable.

4. The underwater target cooperative communication guidance method based on multimodal perception according to claim 1, characterized in that, The dynamic channel occupancy state evolution model is invoked to perform link quality evolution deduction processing on candidate communication links between communication nodes, obtaining the quality fluctuation feature sequence of each candidate communication link in the time dimension and the interference overlap region identifier in the spatial dimension, including: The dynamic signal-to-noise ratio (SNR) change curves of the links between each communication node and the underwater target are extracted from the dynamic channel occupancy state evolution model. Peak and valley detection processing is performed on the dynamic SNR change curves to generate a set of time points with maximum SNR and a set of time points with minimum SNR, which are used as the first fluctuation feature component in the quality fluctuation feature sequence. The dynamic interference power change curves of the links between communication nodes are extracted from the dynamic channel occupancy state evolution model. The dynamic interference power change curves are processed by sliding window energy accumulation calculation to generate the average interference power value in each time window, which is used as the second fluctuation feature component in the quality fluctuation feature sequence. The signal-to-noise ratio dynamic change curve and the interference power dynamic change curve are input into the time series alignment module for time axis synchronization alignment processing to obtain the aligned signal-to-noise ratio-interference power joint curve. Each time point on the joint curve is marked with both the signal-to-noise ratio value and the interference power value. The joint curve is subjected to abrupt change point detection processing to identify the time interval in which the signal-to-noise ratio and interference power change drastically at the same time, and a set of unstable link quality period identifiers is generated. The set of unstable link quality period identifiers includes the start time point and end time point of the unstable period. Based on the set of unstable link quality period identifiers, the average signal-to-noise ratio and average interference power during the stable period are extracted from the joint curve to generate a benchmark quality score for each candidate communication link, which is used as the third fluctuation feature component in the quality fluctuation feature sequence. In terms of spatial dimension, the spatial coordinates of each communication node are processed by three-dimensional spatial mapping to obtain the set of node grid coordinates of each communication node in the three-dimensional grid. Combined with the occlusion attribute identifiers of each three-dimensional grid unit in the environmental grid model, a node spatial distribution map is generated. Calculate the sequence of three-dimensional grid cells that the straight path of each candidate communication link passes through on the node spatial distribution map, and extract the number of other candidate communication links that exist simultaneously in each three-dimensional grid cell as the link density parameter of the straight path; The three-dimensional grid cells whose link density parameters exceed a preset density threshold are marked as potential interference overlap regions. A list of candidate communication link identifiers contained in each potential interference overlap region is compiled to generate the interference overlap region identifier. The interference overlap region identifier includes the region's spatial coordinate range and the set of link identifiers within the region. Connectivity analysis is performed on the potential interference overlapping regions to merge adjacent 3D grid cells into continuous interference overlapping regions, generating an interference overlapping region connectivity graph. The interference overlapping region connectivity graph includes the boundary coordinate sequence of each interference region and the internal link density distribution heatmap. The fluctuation feature components in the quality fluctuation feature sequence are associated with the interference overlap region connectivity graph to generate a correspondence table between the quality fluctuation feature sequence in the time dimension and the interference overlap region identifier in the spatial dimension for each candidate communication link.

5. The underwater target cooperative communication guidance method based on multimodal perception according to claim 1, characterized in that, The step of generating a set of cooperative communication guidance instructions for the underwater target based on the quality fluctuation feature sequence and the interference overlap region identifier includes: The quality fluctuation feature sequence is analyzed, and the set of time points with maximum signal-to-noise ratio (SNR), the set of time points with minimum SNR, and the benchmark quality score for each candidate communication link are extracted. A link quality priority list is constructed, and the link quality priority list sorts each candidate communication link from high to low according to the benchmark quality score. Based on the interference overlap region connectivity graph in the interference overlap region identifier, a set of candidate communication links located in the same interference overlap region is identified, and a list of link conflict groups is generated. The list of link conflict groups includes the communication node pair identifiers and the corresponding spatial coordinate range of the interference region in each link conflict group. For each link conflict group in the link conflict group list, time-domain peak-shifting scheduling analysis is performed. Based on the set of time points with maximum and minimum signal-to-noise ratio of each candidate communication link, a communication time slot allocation scheme that can be staggered on the time axis for different links is determined, and a communication timing scheduling identifier is generated. The communication timing scheduling identifier includes the start time and duration of the transmission time slot of each communication node in a continuous time frame. According to the link quality priority list, an initial transmit power adjustment is assigned to each communication node. The initial transmit power adjustment is negatively correlated with the baseline quality score of the corresponding link. Links with higher baseline quality scores are assigned smaller transmit power adjustments. The initial transmit power adjustment is input into the power optimization iterative model, and power redistribution is performed in conjunction with the link conflict group list. The transmit power adjustment of each communication node in the same link conflict group is adjusted. The transmit power adjustment is the absolute adjustment value of the power. The optimization objective is to minimize the sum of the transmit power adjustments of all links in the same interference overlap area, and thus obtain the optimized transmit power adjustment. Based on the predicted motion trajectory curve of the underwater target and the spatial coordinates of each communication node, the instantaneous beam pointing angle of each communication node pointing to the underwater target is calculated, and an initial beam pointing angle set is generated. The initial beam pointing angle set is subjected to beamwidth adaptive adjustment processing. Based on the distance change sequence between each communication node and the underwater target, the beamwidth adjustment factor of each communication node is determined. The beamwidth adjustment factor is positively correlated with the distance, and the adjusted beam pointing angle adjustment amount is obtained. By integrating the optimized transmit power adjustment, the adjusted beam pointing angle adjustment, and the communication timing scheduling identifier, a preliminary set of cooperative communication guidance instructions is generated. The preliminary set of cooperative communication guidance instructions includes the target communication node identifier and the corresponding parameter adjustment value for each instruction. The preliminary set of cooperative communication guidance instructions is subjected to instruction conflict detection processing to check whether there are multiple conflicting transmit power adjustment amounts or beam pointing angle adjustment amounts assigned to the same communication node at the same time. If so, the instruction with higher priority is selected and retained according to the link quality priority list. After resolving the conflict, the final set of cooperative communication guidance instructions is generated. The final set of cooperative communication guidance instructions is encapsulated according to the receiving protocol of each communication node to generate a sequence of binary instruction frames containing instruction type field, parameter value field and timestamp field, which is then output as a set of cooperative communication guidance instructions for the underwater target.

6. The underwater target cooperative communication guidance method based on multimodal perception according to claim 2, characterized in that, The step of performing pulse feature extraction processing on the sampled set of sound wave propagation waveforms to obtain a set of sound wave pulse arrival time distributions and a set of waveform distortion coefficients includes: Energy envelope detection processing is performed on each waveform sample in the sound wave propagation waveform sampling set to obtain a waveform energy envelope curve. The horizontal axis of the waveform energy envelope curve is the time sampling point, and the vertical axis is the normalized energy amplitude. Peak detection processing is performed on the waveform energy envelope curve to identify local maxima exceeding a preset energy threshold. The time sampling points corresponding to the local maxima are used as candidate pulse arrival times to obtain a set of candidate pulse arrival times. Cluster analysis is performed on the candidate pulse arrival time points. Adjacent candidate pulse arrival time points with a time interval less than a preset time interval threshold are grouped into the same pulse cluster. Each pulse cluster corresponds to a sound wave pulse. The point with the highest energy in each pulse cluster is extracted as the arrival time point of the pulse, thus generating a set of sound wave pulse arrival time points. Based on the waveform sampling values ​​within a preset time window before and after the arrival time of each sound pulse, a pulse waveform segment is extracted, and the pulse waveform segment is subjected to Fourier transform processing to obtain the pulse spectrum distribution. The proportion of the main frequency energy in the pulse spectrum distribution to the total energy is calculated as the pulse energy concentration parameter. The pulse waveform segment is processed by pulse width measurement. The time span from the start of the rising edge to the end of the falling edge is measured to obtain the pulse duration parameter. The pulse duration parameter is compared with the preset standard pulse duration to calculate the duration deviation, which is used as the pulse width parameter. The pulse waveform segment is subjected to energy integration processing to obtain the total pulse energy value. The total pulse energy value is compared with the reference energy value of the transmitting end, and the energy attenuation factor is calculated as the pulse energy attenuation gradient parameter. The pulse energy attenuation gradient parameter reflects the energy loss on the sound wave propagation path. Hyperbolic positioning calculation is performed on the arrival time points of acoustic pulses received by multiple sensing nodes corresponding to the same underwater target to obtain the acoustic positioning coordinates of the underwater target. Combined with the coordinates of each sensing node, an acoustic pulse arrival time distribution set is generated. The acoustic pulse arrival time distribution set includes the acoustic propagation delay corresponding to the distance between each sensing node and the underwater target. The pulse energy concentration parameter, the pulse broadening parameter, and the pulse energy attenuation gradient parameter are associated and stored for each pulse to form the waveform distortion coefficient set. Each coefficient vector in the waveform distortion coefficient set corresponds one-to-one with the arrival time of the corresponding acoustic pulse.

7. The underwater target cooperative communication guidance method based on multimodal perception according to claim 3, characterized in that, The process of performing bidirectional channel sounding simulation on candidate communication links between communication nodes, based on the distribution of obstructions in the environmental grid model and the motion of the underwater target, generates an initial channel impulse response estimation sequence for each candidate communication link, including: The occlusion attribute identifier of each 3D grid cell is extracted from the environmental grid model. The occlusion attribute identifier includes three states: complete occlusion, partial occlusion, and no occlusion. The transmission coefficient of the partially occluded cell is recorded. For each candidate communication link, a three-dimensional straight path is drawn in the environmental grid model based on the spatial coordinates of the communication nodes at both ends, and the sequence of all three-dimensional grid cells that the three-dimensional straight path passes through and the traversal length within each cell are determined. For units in the three-dimensional grid unit sequence that have complete occlusion properties, their corresponding path segments are marked as signal completely blocked segments. For units with partial occlusion properties, the signal attenuation factor of the unit is calculated based on its transmission coefficient. For units with no occlusion properties, they are marked as free propagation segments. The three-dimensional straight path is divided into multiple continuous small segments, each segment corresponding to a three-dimensional grid unit. The total free propagation distance is obtained by summing the lengths of all free propagation segments. The equivalent free propagation distance is obtained by summing the lengths of all partially blocked segments and multiplying them by the corresponding transmission coefficient. The length of the completely blocked segment is the invalid propagation distance. Based on the total free propagation distance and the equivalent free propagation distance, and based on the acoustic wave propagation attenuation model, the amplitude attenuation factor of the direct path is calculated. Based on the spatial coordinates of the underwater target, the possible reflection path of the underwater target as a reflector is determined, and the length of the reflection path and the arrival time delay are calculated. The reflection path is tracked to determine whether it is blocked by other obstructions. If it is not blocked, the amplitude attenuation factor of the reflection path is calculated. This amplitude attenuation factor depends on the reflection coefficient of the reflecting surface and the total length of the reflection path. The direct path and all unblocked reflection paths are sorted in ascending order of arrival time delay. Each path corresponds to an impulse response tap. The amplitude of the impulse response tap is the amplitude attenuation factor of the corresponding path, and the delay of the impulse response tap is the arrival time difference of the path relative to the direct path. The initial channel impulse response estimate is constructed. At each discrete time point, based on the predicted motion trajectory curve of the underwater target, the spatial coordinates of the underwater target are updated, the time delay and amplitude of the reflection path are recalculated, and whether the direct path changes due to the target's motion is determined, thereby generating a time-varying channel impulse response estimation sequence. The channel impulse response estimates of each candidate communication link at different time points are combined to form the initial channel impulse response estimation sequence. Each element in the initial channel impulse response estimation sequence contains a set of multipath tap coefficients and their corresponding time delay values.

8. The underwater target cooperative communication guidance method based on multimodal perception according to claim 4, characterized in that, The process of detecting abrupt changes in the joint curve identifies time intervals where both the signal-to-noise ratio and interference power change drastically, generating a set of identifiers for unstable link quality periods, including: The signal-to-noise ratio curve in the joint curve is calculated by first-order difference to obtain the signal-to-noise ratio change rate curve, and the interference power curve in the joint curve is calculated by first-order difference to obtain the interference power change rate curve. Align the signal-to-noise ratio change rate curve and the interference power change rate curve by time point. For each time point, determine whether the absolute value of the signal-to-noise ratio change rate is greater than a preset signal-to-noise ratio change rate threshold, and at the same time determine whether the absolute value of the interference power change rate is greater than a preset interference power change rate threshold. The time points that simultaneously satisfy the absolute value of the rate of change of signal-to-noise ratio being greater than a threshold and the absolute value of the rate of change of interference power being greater than a threshold are marked as candidate mutation points, thus obtaining the time series of candidate mutation points; A continuity analysis is performed on the time series of the candidate mutation points, and adjacent candidate mutation points with a time interval less than a preset merging time interval are merged into a mutation time period. Each mutation time period includes a start time and an end time. The mean signal-to-noise ratio and mean interference power within each mutation period are calculated and compared with the mean signal-to-noise ratio and mean interference power of the stable periods before and after the mutation period. If the difference exceeds the preset difference threshold, the mutation period is confirmed as a period of unstable link quality. All confirmed unstable link quality periods are arranged in chronological order to generate a set of unstable link quality period identifiers. Each entry in the set of unstable link quality period identifiers includes the start time point, end time point, signal-to-noise ratio fluctuation amplitude, and interference power fluctuation amplitude within the unstable period. For each period of unstable link quality, the predicted motion trajectory of the environmental grid model and the underwater target is traced back to identify the physical cause of the instability during that period of unstable link quality, and the physical cause is recorded as an additional attribute in the unstable link quality period identifier.

9. An underwater target cooperative communication guidance system based on multimodal perception, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the underwater target cooperative communication guidance method based on multimodal perception as described in any one of claims 1 to 8 by executing the machine-executable instructions.

10. A computer program product, characterized in that, The computer program product includes machine-executable instructions stored in a computer-readable storage medium. The processor of the multimodal sensing-based underwater target cooperative communication guidance system reads the machine-executable instructions from the computer-readable storage medium and executes the machine-executable instructions, causing the multimodal sensing-based underwater target cooperative communication guidance system to perform the multimodal sensing-based underwater target cooperative communication guidance method according to any one of claims 1 to 8.

Images