Time-varying tracking and virtual signal combined underwater acoustic channel estimation method

Abstract

The invention belongs to the technical field of underwater acoustic communication, and discloses an underwater acoustic channel estimation method based on combination of time-varying tracking and virtual signals, which is applied to an OFDM (Orthogonal Frequency Division Multiplexing) underwater acoustic communication system adopting a cyclic prefix. The method comprises the following steps: firstly, performing singular value decomposition denoising preprocessing on a received pilot signal, and setting a threshold value based on an energy function to screen an effective path; secondly, path clusters are divided through a clustering algorithm by using historical channel time delay parameters; estimating a time-varying offset parameter of each cluster by adopting an orthogonal matching pursuit algorithm; finally, after a virtual channel signal is constructed, the virtual channel signal and an actual observation signal are jointly input into a time sequence multiple sparse Bayesian learning algorithm for joint channel estimation. According to the method, the noise robustness is enhanced through a dual-signal mechanism, the fast time-varying channel characteristic is accurately tracked by using the cluster-level offset parameter model, and the channel estimation accuracy and stability in the low signal-to-noise ratio and fast time-varying environment are remarkably improved.

Description

Technical Field

[0001] This invention belongs to the field of underwater acoustic communication and is applicable to underwater high-speed communication systems in major marine missions such as deep-sea resource exploration, long-term ecological environment monitoring, and seabed observation networking. It relates to an underwater acoustic channel estimation method that combines time-varying tracking and virtual signals. Background Technology

[0002] Major missions such as marine resource exploration, ecological environment monitoring, seabed observation networking, and underwater security have created an urgent need for high-speed and reliable underwater wireless communication technologies. Underwater acoustic communication technology, as the main technological approach for achieving medium- and long-distance underwater information transmission, is a key support for building a "transparent ocean" and a smart ocean system.

[0003] However, as an extremely complex sound propagation medium, ocean water possesses unique physical characteristics—including spatiotemporally inhomogeneous temperature-salinity stratification, internal fluctuations, turbulence, and the time-varying boundary formed by the sea surface and seabed—which together cause underwater acoustic channels to exhibit strong time-varying, space-varying, and frequency-varying characteristics. Specifically, this manifests as: (1) severe multipath spread, where signals are reflected and scattered through complex interfaces, forming a densely arriving, time-varying diffuse multipath structure; (2) significant Doppler frequency shift and spread, caused by ocean currents, carrier motion, and internal waves, resulting in carrier offset and waveform distortion; and (3) strong non-stationary environmental noise, encompassing noise from shipping, marine life, and wave breaking, which dynamically changes with sea conditions and sea area. These composite interferences introduced by unique ocean physical processes pose severe reliability challenges to underwater communication links, especially in typical marine application scenarios such as single-node communication equipment, AUV collaboration, and seabed sensor networks.

[0004] As the foundational architecture of underwater acoustic communication systems, the transmission performance of single-node communication devices in such complex marine channel environments is highly dependent on the receiver's accuracy in sensing the channel state. Traditional underwater acoustic channel estimation methods generally suffer from the following limitations when dealing with the aforementioned marine characteristics: First, in low signal-to-noise ratio environments, traditional estimators based on least squares and other criteria are extremely sensitive to non-Gaussian noise in the ocean, easily misclassifying strong noise or interference as effective multipath, leading to significant deviations in the estimation results; Second, when identifying effective paths, traditional parameterization methods often rely on isolated experience for threshold settings, lacking consideration of the real-time energy statistical characteristics of the channel, thus making it difficult to track the drastic fluctuations and dynamic evolution of the channel energy and sparse structure caused by internal waves, thermocline changes, seawater movement, etc., often resulting in missed detection of weak effective paths or false inclusion of noisy paths, leading to insufficient estimation robustness; Third, most existing methods process individual snapshots or data blocks in isolation, ignoring the inherent cluster structure of underwater acoustic multipaths (originating from a common reflection interface) and their spatiotemporal correlation between adjacent moments, failing to fully utilize the continuous information of channel state evolution to improve estimation accuracy and tracking capability.

[0005] To address the core challenges posed by complex underwater acoustic channels in the ocean, this invention proposes a high-precision underwater acoustic channel estimation method for marine science and technology applications. This method first employs a dual preprocessing mechanism combining Singular Value Decomposition (SVD) denoising and adaptive energy thresholding to effectively suppress marine environmental noise and extract real underwater acoustic multipath propagation. Next, cluster analysis is used to identify path clusters with similar delay characteristics in the underwater acoustic channel, accurately capturing their multipath propagation features. Based on this, the Orthogonal Matching Pursuit (OMP) algorithm is used to track the gain and delay offset parameters of each path cluster at different times, accurately characterizing the time-varying evolution of the marine channel. Finally, a virtual channel signal is innovatively constructed and input together with the actual observed signal into a Time-Series Multiple Sparse Bayesian Learning (TMSBL) algorithm for joint optimization estimation, significantly improving estimation accuracy and stability. This invention aims to significantly improve the reliability of channel estimation for underwater acoustic communication systems in real marine environments, providing stronger underwater communication technology support for marine resource development, environmental monitoring, and underwater defense activities. Summary of the Invention

[0006] 1. A time-varying tracking and virtual signal combined underwater acoustic channel estimation method, applied to an orthogonal frequency division multiplexing (OFDM) underwater acoustic communication system using a cyclic prefix, wherein the system has a number of subcarriers. pilot subcarrier number Symbol period The channel model has a path count of The sparse time-varying channel, characterized by the underwater acoustic channel estimation method combining time-varying tracking and virtual signals, includes the following steps:

[0007] S1. The receiver performs time synchronization on the OFDM symbol blocks transmitted in the underwater acoustic channel, removes the cyclic prefix, and obtains... Point sampling data, This indicates the number of points in the Fast Fourier Transform operation; then, the following steps are performed. Point FFT is used to obtain the complete frequency domain vector. Finally, based on the index of the system pilot subcarriers, the pilot data at the corresponding position is extracted from the frequency domain vector to form the pilot observation vector. Continuous processing After one OFDM symbol block, the vectors are arranged into a pilot observation matrix. ;for Perform singular value decomposition:

[0008]

[0009] in, , It is a singular value matrix. Take the matrix of A singular value ,in Calculate the average value Keep all greater than The singular values ​​form a new diagonal matrix Reconstruction and Same dimension matrix for:

[0010]

[0011] The denoised receiver pilot matrix is ​​as follows: ;

[0012] S2. Utilize Calculate the frequency domain channel response , It is a diagonal matrix composed of known symbols transmitted on pilot subcarriers; for Perform an IFFT on each column to obtain the contents Time-domain channel impulse response of a delay tap :

[0013]

[0014] in, This represents the inverse fast Fourier transform, and Numerically with Equal; Calculate Energy superposition function Set threshold Based on threshold Perform channel filtering and generate a binary index set. ,Will The index with a median value of 0 is denoted as Then define the dictionary matrix. The specific formula is as follows:

[0015]

[0016] in, The DFT matrix is ​​used; by using the dictionary matrix In and index set Setting the corresponding column elements to zero yields a simplified dictionary matrix. The time-domain channel impulse response The simplified time-domain channel impulse response is obtained by setting the corresponding row to zero. Channel matrix Simplified to ;

[0017] S3. Through Obtain the set of valid paths:

[0018]

[0019] The number of elements in the set is the number of valid paths. ,and Then the simplified channel matrix In the Each path in an OFDM block uses a gain parameter. and delay parameters Parametric representation:

[0020]

[0021] in, Indicates the first One OFDM block in the path The complex amplitude of is expressed as:

[0022]

[0023] Includes the first One OFDM block in this path Amplitude attenuation and by carrier frequency Timely delay The phase shift is jointly determined; For the time-delay-dependent steering vector, its first... The element (corresponding to the ) (a number of subcarriers) is defined as:

[0024]

[0025] in, , The system symbol period; The effective path is represented as follows: For the first OFDM symbol block, the effective path is obtained after executing S1 and S2. Then, the obvious peaks in the obtained time-domain response are extracted as amplitude parameters, and the indices corresponding to the peaks are extracted as time delay parameters. The channel matrix is ​​then processed using the K-means clustering algorithm. The path partitioning in the optimal cluster number And obtain the cluster allocation results. , of which A cluster is a set of path indexes with similar latency characteristics, represented as: , The number of paths in the corresponding data is: and cluster centroid ;

[0026] S4. Based on the division in S3 Clusters and the use of historical blocks The channel parameters are used to calculate the cluster average gain and cluster average delay, which are then used as the first... The representative characteristics of this cluster when it is in block form:

[0027]

[0028]

[0029] Furthermore, the time-varying offset parameter set is iteratively filtered using the orthogonal matching pursuit algorithm. And the channel parameters are predicted as follows:

[0030]

[0031] in This is the path gain offset. This is the time delay offset;

[0032] S5. Time-varying migration parameters estimated based on S4 Construct the first Block 1 The cluster's virtual channel response is:

[0033]

[0034] By combining the predicted responses of all clusters, we obtain the complete predicted simplified channel response for the current block:

[0035]

[0036] All The predicted simplified channel responses of each block are arranged in a matrix. To further improve the accuracy of the estimation results, based on Constructing virtual signals The actual received signal is used as input to the sparse Bayesian learning algorithm for iterative estimation, and the final output is the estimated channel impulse response. .

[0037] Furthermore, the binary index set mentioned in step S2 The calculation process is as follows:

[0038] S2.1: Determine the energy coefficient To preset a fixed value, The value of depends on prior knowledge of the underwater acoustic channel, such as noise level and signal sparsity, and is not dynamically adjusted with the real-time state of the channel during the channel estimation process.

[0039] S2.2: Threshold The following steps are used to set the parameters: First, calculate the energy superposition function of the channel impulse response estimate. :

[0040]

[0041] in, No. Average energy per time-delay tap Represented as In the On the OFDM block The response value at each time delay tap; then the energy coefficient Maximum value of the superposition function with energy Multiply to obtain the threshold for:

[0042]

[0043] in The value of the threshold makes It can distinguish between valid channel taps and spurious taps caused by noise;

[0044] Binary index set for:

[0045] .

[0046] Furthermore, the optimal clustering number mentioned in step S3 The determination process is as follows:

[0047] Set a number of candidate clusters The search range is 1 to For each candidate The values ​​are used to perform K-means clustering and then the sum of squared intra-cluster distances is calculated:

[0048]

[0049] in Indicates the first Cluster , Indicates the first The centroid of a cluster, For delay parameters;

[0050] The results of WCSS are recorded as a decreasing sequence of values. Then calculate the relative rate of change: When the relative rate of change increases sharply, the corresponding The value was determined as the optimal cluster number. .

[0051] Furthermore, the OMP algorithm process described in step S4 is as follows:

[0052] First, the delay offset for each cluster. Define a discrete search grid:

[0053]

[0054] in, The preset maximum time delay offset threshold, The search step size is [value], and the grid size is: Then, for each time delay offset point in the search grid... According to the prediction formula for channel parameters:

[0055]

[0056] Constructing candidate atoms:

[0057]

[0058] Combine all candidate atoms of all clusters to construct a global offset parameter dictionary matrix. :

[0059] ;

[0060] Set initial residuals The observed pilot signal vector for the current block A copy, initializing the active set An empty set:

[0061]

[0062] Calculate the current residual With dictionary Find the most relevant atom by taking the absolute value of the inner product of all atoms in the atom:

[0063]

[0064] in This represents the number of iterations.

[0065] Secondly, select the atom index corresponding clusters Join the active group And record its corresponding time delay offset. :

[0066]

[0067] Secondly, based on the current active set All atoms in the cluster are used to jointly estimate the gain offset of these clusters using the least squares method. To minimize the residuals:

[0068]

[0069] Update the residuals using the estimated gain offset:

[0070]

[0071] Finally, when the active set The size equals At that time, or residual norm Less than the preset threshold When the iteration terminates, the output is... .

[0072] Furthermore, the TMSBL algorithm process described in step S5 is as follows:

[0073] based on Constructing Channel Response , and thus The channel responses of each block are arranged into a channel matrix. ,according to Constructing virtual signals and in conjunction with the actual received signal As input to the algorithm:

[0074]

[0075]

[0076] Set the initial values ​​for the following hyperparameters and variables:

[0077] Channel covariance matrix Initialize to the identity matrix Its diagonal elements Representing the Channel delay taps The initial power;

[0078] Time correlation matrix Initialize as an identity matrix ;

[0079] noise variance Initial value of noise variance , This is the matrix of noise data received during the silent period;

[0080] Set convergence threshold With maximum number of iterations ;

[0081] Then vectorize it:

[0082] Joint observation signal vectorization: ;

[0083] Constructing a block diagonal sensing matrix: ,in Represents the Kronecker product;

[0084] Construct the total prior covariance matrix: , For the first The next iteration;

[0085] Then calculate the posterior covariance matrix and mean vector:

[0086]

[0087]

[0088] posterior mean vector Inverse vectorization to matrix form ;

[0089] Then update the channel covariance matrix. Time correlation matrix :

[0090]

[0091]

[0092]

[0093] in, and Represent the first and second parts of the matrix respectively. Element and the diagonal elements; for No. The transpose of the row represents the first row. Each time-domain tap in all Gain sequence on symbol blocks From the posterior covariance matrix Extracted from, corresponding to the first A matrix with time-domain taps;

[0094] Update noise variance:

[0095]

[0096] in, It is the Frobenius norm. The trace of the matrix;

[0097] Finally, calculate the hyperparameter vectors generated between two adjacent iterations. and The maximum relative change, if it satisfies:

[0098]

[0099] Or the number of iterations Reach the maximum preset value The iteration terminates and the final iteration number is recorded. Otherwise, let Continue iterating; after the algorithm converges, output the posterior mean matrix obtained from the last calculation:

[0100]

[0101] This matrix is ​​the final estimated time-domain channel impulse response matrix, i.e. In this embodiment, ,

[0102] Furthermore, the specific process of the K-means algorithm is as follows:

[0103] For the first centroid, from the time delay parameter vector A time delay value is randomly and uniformly selected as the first initial centroid. ; Calculate each delay value The square of the shortest distance to all selected centroids Then with probability Select the next centroid from the remaining delay values, and repeat until the next centroid is selected. The feature points are used as the initial cluster centroids. ;

[0104] Then each feature point is assigned to the cluster to which the nearest centroid belongs:

[0105]

[0106] Represented as the number of iterations;

[0107] Recalculate the centroid of each cluster:

[0108]

[0109] when The iteration stops when the number of iterations is less than the threshold or the maximum number of iterations is reached. After the iteration is complete, the cluster assignment results are output. Each cluster contains a set of path indexes with similar latency characteristics. , The number of paths in the corresponding data is: and cluster centroid .

[0110] Compared with existing methods, the advantages of this method are as follows:

[0111] (1) A dual filtering mechanism combining SVD denoising preprocessing and parameterized threshold screening based on energy function is adopted to suppress generalized noise and false paths at the signal matrix and channel tap levels respectively, thereby improving the signal purity under low signal-to-noise ratio.

[0112] (2) Introduce a cluster-level time-varying offset parameter model, identify the cluster structure of the path through cluster analysis, and use the OMP algorithm to independently estimate the gain and delay offset of each cluster to accurately characterize the fast time-varying characteristics of the channel.

[0113] (3) The number of effective paths is significantly reduced by threshold screening, and path-level management is transformed into cluster-level parameter estimation by combining clustering. While maintaining temporal correlation, the model complexity is reduced, and the computational load and estimation accuracy are effectively balanced.

[0114] (4) Construct a virtual channel signal and input it together with the actual observed signal into the time-series multiple sparse Bayesian learning algorithm for joint estimation. Make full use of time-series correlation and ensure the stability of the result through the algorithm’s inherent noise suppression mechanism when there is a prediction bias. Attached Figure Description

[0115] Figure 1 This is a flowchart of the present invention;

[0116] Figure 2This invention compares its performance with the LS algorithm, OMP algorithm, and TMSBL algorithm under both fast and slow channel conditions.

[0117] Figure 3 The comparison is based on the mean squared error with different α values ​​and whether or not noise reduction is used as variables;

[0118] Figure 4 The comparison of bit error rates is based on different α values ​​and whether or not noise reduction is used as variables. Detailed Implementation

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

[0120] This embodiment uses the number of system subcarriers. pilot subcarrier number The symbol period is The orthogonal frequency division multiplexing (OFDM) underwater acoustic communication system has a channel model based on path number. sparse time-varying channels.

[0121] like Figure 1 As shown, a method for underwater acoustic channel estimation combining time-varying tracking and virtual signals includes the following steps:

[0122] S1. The receiver performs time synchronization on the OFDM symbol blocks transmitted in the underwater acoustic channel, removes the cyclic prefix, and obtains... Point sampling data, This indicates the number of points in the Fast Fourier Transform operation; then, the following steps are performed. Point FFT is used to obtain the complete frequency domain vector. Finally, based on the index of the system pilot subcarriers, the pilot data at the corresponding position is extracted from the frequency domain vector to form the pilot observation vector. Continuous processing After one OFDM symbol block, the vectors are arranged into a pilot observation matrix. ;for Perform singular value decomposition:

[0123]

[0124] in, , It is a singular value matrix. Take the matrix of A singular value ,in Calculate the average value Keep all greater than The singular values ​​form a new diagonal matrix Reconstruction and Same dimension matrix for:

[0125]

[0126] The denoised receiver pilot matrix is ​​as follows: In this embodiment, ;

[0127] S2. Utilize Calculate the frequency domain channel response , It is a diagonal matrix composed of known symbols transmitted on pilot subcarriers; for Perform an IFFT on each column to obtain the contents Time-domain channel impulse response of a delay tap :

[0128]

[0129] in, This represents the inverse fast Fourier transform, and Numerically with Equal; Calculate Energy superposition function Set threshold Based on threshold Perform channel filtering and generate a binary index set. ,Will The index with a median value of 0 is denoted as Then define the dictionary matrix. The specific formula is as follows:

[0130]

[0131] in, The DFT matrix is ​​used; by using the dictionary matrix In and index set Setting the corresponding column elements to zero yields a simplified dictionary matrix. The time-domain channel impulse response The simplified time-domain channel impulse response is obtained by setting the corresponding row to zero. Channel matrix Simplified to ;

[0132] S3. Through Obtain the set of valid paths:

[0133]

[0134] The number of elements in the set is the number of valid paths. ,and Then the simplified channel matrix In the Each path in an OFDM block uses a gain parameter. and delay parameters Parametric representation:

[0135]

[0136] in, Indicates the first One OFDM block in the path The complex amplitude of is expressed as:

[0137]

[0138] Includes the first One OFDM block in this path Amplitude attenuation and by carrier frequency Timely delay The phase offset is jointly determined; in this embodiment, the carrier frequency... ; For the time-delay-dependent steering vector, its first... The element (corresponding to the ) (a number of subcarriers) is defined as:

[0139]

[0140] in, , The system symbol period; The effective path is represented as follows: For the first OFDM symbol block, the effective path is obtained after executing S1 and S2. Then, the obvious peaks in the obtained time-domain response are extracted as amplitude parameters, and the indices corresponding to the peaks are extracted as time delay parameters. The channel matrix is ​​then processed using the K-means clustering algorithm. The path partitioning in the optimal cluster number And obtain the cluster allocation results. , of which A cluster is a set of path indexes with similar latency characteristics, represented as: , The number of paths in the corresponding data is: and cluster centroid ;

[0141] S4. Based on the division in S3 Clusters and the use of historical blocks The channel parameters are used to calculate the cluster average gain and cluster average delay, which are then used as the first... The representative characteristics of this cluster when it is in block form:

[0142]

[0143]

[0144] Furthermore, the time-varying offset parameter set is iteratively filtered using the orthogonal matching pursuit algorithm. And the channel parameters are predicted as follows:

[0145]

[0146] in This is the path gain offset. This is the time delay offset;

[0147] S5. Time-varying migration parameters estimated based on S4 Construct the first Block 1 The cluster's virtual channel response is:

[0148]

[0149] By combining the predicted responses of all clusters, we obtain the complete predicted simplified channel response for the current block:

[0150]

[0151] All The predicted simplified channel responses of each block are arranged in a matrix. To further improve the accuracy of the estimation results, based on Constructing virtual signals The actual received signal is used as input to the sparse Bayesian learning algorithm for iterative estimation, and the final output is the estimated channel impulse response. .

[0152] Furthermore, the binary index set mentioned in step S2 The calculation process is as follows:

[0153] S2.1: Determine the energy coefficient To preset a fixed value, The value of depends on prior knowledge of the underwater acoustic channel, such as noise level and signal sparsity, and is not dynamically adjusted with the real-time state of the channel during the channel estimation process.

[0154] S2.2: Threshold The following steps are used to set the parameters: First, calculate the energy superposition function of the channel impulse response estimate. :

[0155]

[0156] in, No. Average energy per time-delay tap Represented as In the On the OFDM block The response value at each time delay tap; then the energy coefficient Maximum value of the superposition function with energy Multiply to obtain the threshold for:

[0157]

[0158] in The value of the threshold makes It can distinguish between valid channel taps and spurious taps caused by noise;

[0159] Binary index set for:

[0160]

[0161] Furthermore, the optimal clustering number mentioned in step S3 The determination process is as follows:

[0162] Set a number of candidate clusters The search range is 1 to For each candidate The values ​​are used to perform K-means clustering and then the sum of squared intra-cluster distances is calculated:

[0163]

[0164] in Indicates the first Cluster , Indicates the first The centroid of a cluster, For delay parameters;

[0165] The results of WCSS are recorded as a decreasing sequence of values. Then calculate the relative rate of change: When the relative rate of change increases sharply, the corresponding The value was determined as the optimal cluster number. .

[0166] Furthermore, the OMP algorithm process described in step S4 is as follows:

[0167] First, the delay offset for each cluster. Define a discrete search grid:

[0168]

[0169] in, The preset maximum time delay offset threshold, The search step size is [value], and the grid size is: In this embodiment, Then, for each time delay offset point in the search grid... According to the prediction formula for channel parameters:

[0170]

[0171] Constructing candidate atoms:

[0172]

[0173] Combine all candidate atoms of all clusters to construct a global offset parameter dictionary matrix. :

[0174] ;

[0175] Set initial residuals The observed pilot signal vector for the current block A copy, initializing the active set An empty set:

[0176]

[0177] Calculate the current residual With dictionary Find the most relevant atom by taking the absolute value of the inner product of all atoms in the atom:

[0178]

[0179] in This represents the number of iterations.

[0180] Secondly, select the atom index corresponding clusters Join the active group And record its corresponding time delay offset. :

[0181]

[0182] Secondly, based on the current active set All atoms in the cluster are used to jointly estimate the gain offset of these clusters using the least squares method. To minimize the residuals:

[0183]

[0184] Update the residuals using the estimated gain offset:

[0185]

[0186] Finally, when the active set The size equals At that time, or residual norm Less than the preset threshold When the iteration terminates, the output is... In this embodiment, a preset threshold is used. .

[0187] Furthermore, the TMSBL algorithm process described in step S5 is as follows:

[0188] based on Constructing Channel Response , and thus The channel responses of each block are arranged into a channel matrix. ,according to Constructing virtual signals and in conjunction with the actual received signal As input to the algorithm:

[0189]

[0190]

[0191] Set the initial values ​​for the following hyperparameters and variables:

[0192] Channel covariance matrix Initialize to the identity matrix Its diagonal elements Representing the Channel delay taps The initial power;

[0193] Time correlation matrix Initialize as an identity matrix ;

[0194] noise variance Initial value of noise variance , This is the matrix of noise data received during the silent period;

[0195] Set convergence threshold With maximum number of iterations ;

[0196] Then vectorize it:

[0197] Joint observation signal vectorization: ;

[0198] Constructing a block diagonal sensing matrix: ,in Represents the Kronecker product;

[0199] Construct the total prior covariance matrix: , For the first The next iteration;

[0200] Then calculate the posterior covariance matrix and mean vector:

[0201]

[0202]

[0203] posterior mean vector Inverse vectorization to matrix form ;

[0204] Then update the channel covariance matrix. Time correlation matrix :

[0205]

[0206]

[0207]

[0208] in, and Represent the first and second parts of the matrix respectively. Element and the diagonal elements; for No. The transpose of the row represents the first row. Each time-domain tap in all Gain sequence on symbol blocks From the posterior covariance matrix Extracted from, corresponding to the first A matrix with time-domain taps;

[0209] Update noise variance:

[0210]

[0211] in, It is the Frobenius norm. The trace of the matrix;

[0212] Finally, calculate the hyperparameter vectors generated between two adjacent iterations. and The maximum relative change, if it satisfies:

[0213]

[0214] Or the number of iterations Reach the maximum preset value The iteration terminates and the final iteration number is recorded. Otherwise, let Continue iterating; after the algorithm converges, output the posterior mean matrix obtained from the last calculation:

[0215]

[0216] This matrix is ​​the final estimated time-domain channel impulse response matrix, i.e. In this embodiment, , .

[0217] Furthermore, the specific process of the K-means algorithm is as follows:

[0218] For the first centroid, from the time delay parameter vector A time delay value is randomly and uniformly selected as the first initial centroid. ; Calculate each delay value The square of the shortest distance to all selected centroids Then with probability Select the next centroid from the remaining delay values, and repeat until the next centroid is selected. The feature points are used as the initial cluster centroids. ;

[0219] Then each feature point is assigned to the cluster to which the nearest centroid belongs:

[0220]

[0221] Represented as the number of iterations;

[0222] Recalculate the centroid of each cluster:

[0223]

[0224] when Less than the threshold Or reach the maximum number of iterations Stop when the iteration is complete, and output the cluster allocation results after the iteration is finished. Each cluster contains a set of path indexes with similar latency characteristics. , The number of paths in the corresponding data is: and cluster centroid In this embodiment, .

[0225] This invention proposes a combined time-varying tracking and virtual signal underwater acoustic channel estimation method. By classifying the channel into clusters and iteratively estimating the offset parameters using the OMP algorithm, the time-varying nature of the channel is characterized, achieving good channel estimation performance in rapidly time-varying underwater channel environments. Performance results are as follows: Figure 2 , Figure 3 , Figure 4 As shown, from Figure 2 As can be seen, the method disclosed in this patent yields the best results under both slow and fast time-varying channel conditions; from Figure 3 , Figure 4 It can be seen in different Under the given conditions, the TMSBL algorithm with noise reduction prior performs best.

[0226] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for underwater acoustic channel estimation combining time-varying tracking and virtual signals, applied to an orthogonal frequency division multiplexing (OFDM) underwater acoustic communication system using a cyclic prefix, wherein the system has a number of subcarriers. pilot subcarrier number Symbol period The channel model has a path count of The sparse time-varying channel is characterized by, The underwater acoustic channel estimation method combining time-varying tracking and virtual signals includes the following steps: S1. The receiver performs time synchronization on the OFDM symbol blocks transmitted in the underwater acoustic channel, removes the cyclic prefix, and obtains... Point sampling data, This indicates the number of points in the Fast Fourier Transform operation; then, the following steps are performed. Point FFT is used to obtain the complete frequency domain vector. Finally, based on the index of the system pilot subcarriers, the pilot data at the corresponding position is extracted from the frequency domain vector to form the pilot observation vector. Continuous processing After one OFDM symbol block, the vectors are arranged into a pilot observation matrix. ;for Perform singular value decomposition: in, , It is a singular value matrix. Take the matrix of A singular value ,in Calculate the average value Keep all greater than The singular values ​​form a new diagonal matrix Reconstruction and Same dimension matrix for: The denoised receiver pilot matrix is ​​as follows: ; S2. Utilize Calculate the frequency domain channel response , It is a diagonal matrix composed of known symbols transmitted on pilot subcarriers; for Perform an IFFT on each column to obtain the contents Time-domain channel impulse response of a delay tap : in, This represents the inverse fast Fourier transform, and Numerically with Equal; Calculate Energy superposition function Set threshold Based on threshold Perform channel filtering and generate a binary index set. ,Will The index with a median value of 0 is denoted as Then define the dictionary matrix. The specific formula is as follows: in, The DFT matrix is ​​used; by using the dictionary matrix In and index set Setting the corresponding column elements to zero yields a simplified dictionary matrix. The time-domain channel impulse response The simplified time-domain channel impulse response is obtained by setting the corresponding row to zero. Channel matrix Simplified to ; S3. Through Obtain the set of valid paths: The number of elements in the set is the number of valid paths. ,and Then the simplified channel matrix In the Each path in an OFDM block uses a gain parameter. and delay parameters Parametric representation: in, Indicates the first One OFDM block in the path The complex amplitude of is expressed as: Includes the first One OFDM block in this path Amplitude attenuation and by carrier frequency Timely delay The phase shift is jointly determined; For the time-delay-dependent steering vector, its first... The element that corresponds to the th element The elements of each subcarrier are defined as: in, , The system symbol period; The effective path is represented as follows: For the first OFDM symbol block, the effective path is obtained after executing S1 and S2. Then, the obvious peaks in the obtained time-domain response are extracted as amplitude parameters, and the indices corresponding to the peaks are extracted as time delay parameters. The K-means clustering algorithm is used to cluster the channel matrix. The path partitioning in the optimal cluster number And obtain the cluster allocation results. , of which A cluster is a set of path indexes with similar latency characteristics, represented as: , The number of paths in the corresponding data is: and cluster centroid ; S4. Based on the division in S3 Clusters and the use of historical blocks The channel parameters are used to calculate the cluster average gain and cluster average delay, which are then used as the first... The representative characteristics of this cluster when it is in block form: Furthermore, the time-varying offset parameter set is iteratively filtered using the orthogonal matching pursuit algorithm. And the channel parameters are predicted as follows: in This is the path gain offset. This is the time delay offset; S5. Time-varying migration parameters estimated based on S4 Construct the first Block 1 The cluster's virtual channel response is: By combining the predicted responses of all clusters, we obtain the complete predicted simplified channel response for the current block: All The predicted simplified channel responses of each block are arranged in a matrix. To further improve the accuracy of the estimation results, based on Constructing virtual signals The actual received signal is used as input to the sparse Bayesian learning algorithm for iterative estimation, and the final output is the estimated channel impulse response. .

2. The underwater acoustic channel estimation method combining time-varying tracking and virtual signals according to claim 1, characterized in that, The binary index set mentioned in step S2 The calculation process is as follows: S2.1: Determine the energy coefficient To preset a fixed value, The value of depends on prior knowledge of the underwater acoustic channel, such as noise level and signal sparsity, and is not dynamically adjusted with the real-time state of the channel during the channel estimation process. S2.2: Threshold The following steps are used to set the parameters: First, calculate the energy superposition function of the channel impulse response estimate. : in, No. Average energy per time-delay tap Represented as In the On the OFDM block, the first The response value at each time delay tap; then the energy coefficient Maximum value of the superposition function with energy Multiply to obtain the threshold for: in The value of the threshold makes It can distinguish between valid channel taps and spurious taps caused by noise; Binary index set for: 。 3. The underwater acoustic channel estimation method combining time-varying tracking and virtual signals according to claim 1, characterized in that, The optimal number of clusters mentioned in step S3 The determination process is as follows: Set a number of candidate clusters The search range is 1 to For each candidate The values ​​are used to perform K-means clustering and then the sum of squared intra-cluster distances is calculated: in Indicates the first Cluster , Indicates the first The centroid of a cluster, For delay parameters; The results of WCSS are recorded as a decreasing sequence of values. Then calculate the relative rate of change: When the relative rate of change increases sharply, the corresponding The value was determined as the optimal cluster number. .

4. The underwater acoustic channel estimation method combining time-varying tracking and virtual signals according to claim 1, characterized in that, The OMP algorithm process described in step S4 is as follows: First, the delay offset for each cluster. Define a discrete search grid: in, The preset maximum time delay offset threshold, The search step size is [value], and the grid size is: Then, for each time delay offset point in the search grid... According to the prediction formula for channel parameters: Constructing candidate atoms: Combine all candidate atoms of all clusters to construct a global offset parameter dictionary matrix. : ； Set initial residuals The observed pilot signal vector for the current block A copy, initializing the active set An empty set: Calculate the current residual With dictionary Find the most relevant atom by taking the absolute value of the inner product of all atoms in the atom: in This represents the number of iterations. Then, select the atom index corresponding clusters Join the active group And record its corresponding time delay offset. : Then, based on the current active set All atoms in the cluster are used to jointly estimate the gain offset of these clusters using the least squares method. To minimize the residuals: Update the residuals using the estimated gain offset: Finally, when the active set The size equals At that time, or residual norm Less than the preset threshold When the iteration terminates, the output is... .

5. The underwater acoustic channel estimation method combining time-varying tracking and virtual signals according to claim 1, characterized in that, The TMSBL algorithm process described in step S5 is as follows: based on Constructing Channel Response , and thus The channel responses of each block are arranged into a channel matrix. ,according to Constructing virtual signals and in conjunction with the actual received signal As input to the algorithm: Set the initial values ​​for the following hyperparameters and variables: Channel covariance matrix Initialize to the identity matrix Its diagonal elements Representing the Channel delay taps The initial power; Time correlation matrix Initialize as an identity matrix ; noise variance Initial value of noise variance , This is the matrix of noise data received during the silent period; Set convergence threshold With maximum number of iterations ; Then vectorize it: Joint observation signal vectorization: ; Constructing a block diagonal sensing matrix: ,in Represents the Kronecker product; Construct the total prior covariance matrix: , For the first The next iteration; Then calculate the posterior covariance matrix and mean vector: posterior mean vector Inverse vectorization to matrix form ; Then update the channel covariance matrix. Time correlation matrix : in, and Represent the first and second parts of the matrix respectively. Element and the diagonal elements; for No. The transpose of the row represents the first row. Each time-domain tap in all Gain sequence on symbol blocks From the posterior covariance matrix Extracted from, corresponding to the first A matrix with time-domain taps; Update noise variance: in, It is the Frobenius norm. The trace of the matrix; Finally, calculate the hyperparameter vectors generated between two adjacent iterations. and The maximum relative change, if it satisfies: Or the number of iterations Reach the maximum preset value The iteration terminates and the final iteration number is recorded. Otherwise, let Continue iterating; after the algorithm converges, output the posterior mean matrix obtained from the last calculation: This matrix is ​​the final estimated time-domain channel impulse response matrix, i.e. .

6. The underwater acoustic channel estimation method combining time-varying tracking and virtual signals according to claim 3, characterized in that, The specific process of the K-means algorithm is as follows: For the first centroid, from the time delay parameter vector A time delay value is randomly and uniformly selected as the first initial centroid. ; Calculate each delay value The square of the shortest distance to all selected centroids Then with probability Select the next centroid from the remaining delay values, and repeat until the next centroid is selected. The feature points are used as the initial cluster centroids. ; Then each feature point is assigned to the cluster to which the nearest centroid belongs: Represented as the number of iterations; Recalculate the centroid of each cluster: when The iteration stops when the number of iterations is less than the threshold or the maximum number of iterations is reached. After the iteration is complete, the cluster assignment results are output. Each cluster contains a set of path indexes with similar latency characteristics. , The number of paths in the corresponding data is: and cluster centroid .

Images