CN109039506B - A kind of underwater mobile channel emulation mode - Google Patents
A kind of underwater mobile channel emulation mode Download PDFInfo
- Publication number
- CN109039506B CN109039506B CN201810796972.4A CN201810796972A CN109039506B CN 109039506 B CN109039506 B CN 109039506B CN 201810796972 A CN201810796972 A CN 201810796972A CN 109039506 B CN109039506 B CN 109039506B
- Authority
- CN
- China
- Prior art keywords
- channel
- time
- signal
- platform
- impulse response
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04B—TRANSMISSION
- H04B17/00—Monitoring; Testing
- H04B17/30—Monitoring; Testing of propagation channels
- H04B17/391—Modelling the propagation channel
- H04B17/3912—Simulation models, e.g. distribution of spectral power density or received signal strength indicator [RSSI] for a given geographic region
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04B—TRANSMISSION
- H04B13/00—Transmission systems characterised by the medium used for transmission, not provided for in groups H04B3/00 - H04B11/00
- H04B13/02—Transmission systems in which the medium consists of the earth or a large mass of water thereon, e.g. earth telegraphy
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Signal Processing (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Electromagnetism (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
The invention discloses a kind of underwater mobile channel emulation modes, the described method includes: step 1) sets the environmental parameter and signal parameter of three-dimensional hydrospace communication scenes, the motion model and sea level height model of source emission platform and receiving platform are established；Step 2) obtains the channel snap of different moments a series of, and calculate time-varying channel impulse response accordingly by repeatedly calling ray tracing model to carry out real-time sampling to the state of channel；Step 3) is responded using time-varying channel impulse and carries out time domain dynamic filter to the transmitting signal of source emission platform, obtains receiving signal.Method of the invention guarantees the accuracy and precision of underwater mobile telecommunication channel emulation to a certain extent, meanwhile, it can effectively reduce the complexity of the Realization of Simulation, reduce memory requirements, improve computer operational efficiency.
Description
Technical field
The present invention relates to field of underwater acoustic communication, and in particular to a kind of underwater mobile channel emulation mode.
Background technique
In underwater sound communication, the movement of communications platform and the movement of seawater itself can all bring strong channel time-varying and more
Doppler spread seriously affects communication quality.It is communicated under the movement of subsurface communication platform or propagation medium motion conditions effectively to assess
The performance of system or system module needs to emulate underwater mobile telecommunication channel.
For shallow sea underwater sound mobile communication scene, a kind of fairly simple channel simulator implementation method is: first by penetrating
Line model calculates the relevant information of each intrinsic sound ray of static channel, such as propagation delay, decaying and direction of arrival；Then basis
The speed and direction of arrival of transmit-receive platform calculate the Doppler factor in each intrinsic path；Further according to the time delay of intrinsic sound ray and more
Pu Le successively carries out delay and scaling operation to transmitting signal, the signal copy in different paths is obtained, finally by these signal pairs
This superposition obtains receiving signal.This method only needs to carry out an acoustic ray tracing, and complexity is very low, but application scenarios ten
Branch office's limit, is slowed by, the channel scenario that medium level dependence is weak when being only applicable to.And, platform big for sea, submarine relief
Situations such as vertical movement and platform movement velocity are very fast, the above-mentioned quasi-static emulation mode based on an acoustic ray tracing are just incompetent
For power, because it is unable to get channel multipara structure and intensity changes with time situation.At this time, it is necessary to repeatedly be believed
Road acoustic ray tracing obtains the situation of change in accurate channel Eigen path.
Channel simulator based on multiple acoustic ray tracing usually samples channel status with a certain fixed frequency, wherein
Sampling is known as a channel snap each time.At each channel snap moment, according to the position of current flat pad, battle array is received
Position and sea level height carry out acoustic ray tracing to channel, obtain the relevant information of each intrinsic sound ray, and be calculated work as accordingly
The channel impulse response at preceding moment；The channel impulse response at each moment is combined again to obtain time-varying channel impulse response；
Transmitting signal is responded with time-varying channel impulse finally and carries out time-varying convolution, obtains each reception signal for receiving array element.
The existing underwater mobile channel emulation mode based on an acoustic ray tracing only considers the first-order dynamic characteristic of channel,
Can not in the case that accurate simulation channel time-varying is significant, medium level dependence is strong channel multipara structure variation, the scope of application
Extremely limit to.And the existing emulation mode based on multiple acoustic ray tracing would generally preset the fast beat frequency of fixed channel.
Often there is biggish randomness in this way, lack a specific criterion, and parameter selection mechanism is inflexible, cannot
The sampling frequency of channel status is adjusted according to the speed of real channel variance.Once channel snap set of frequency is unreasonable, whole
The accuracy or efficiency of a channel simulator just will receive strong influence: the fast beat frequency of channel is too low, and simulation result then may nothing
Method keeps up with the quick variation of actual channel structure；And the fast beat frequency of channel is excessive, then will lead to the increase and calculating of operation time
The waste of resource.
Summary of the invention
It is an object of the invention to solve above-mentioned technological deficiency, on the basis of existing multiple acoustic ray tracing emulation mode,
Emphasis improves the selection mechanism of the fast beat frequency of channel, can be automatically adjusted according to the speed of real channel variance
It is real to reduce emulation to be reduced as far as channel snap number under the premise of guaranteeing accuracy of simulation for the fast beat frequency of channel
Existing complexity improves operational efficiency.
To achieve the goals above, the invention proposes a kind of underwater mobile channel emulation modes, which comprises
Step 1) sets the environmental parameter and signal parameter of three-dimensional hydrospace communication scenes, establishes source emission
The motion model and sea level height model of platform and receiving platform；
Step 2) obtains a series of differences by repeatedly calling ray tracing model to carry out real-time sampling to the state of channel
The channel snap at moment, and time-varying channel impulse response is calculated accordingly；
Step 3) is responded using time-varying channel impulse and carries out time domain dynamic filter to the transmitting signal of source emission platform, is obtained
To reception signal.
As a kind of improvement of the above method, the step 1) is specifically included:
Step 1-1) setting underwater 3 D space range: horizontal transverse and longitudinal coordinate range and sea area depth capacity；Import three
Velocity of sound distributed data and submarine topography data in water in dimension space；Set substrate parameter: compression velocity of wave, density and attenuation rate；
Step 1-2) import the required signal of communication emittedSampled data, and Setting signal length T, centre frequency
fc, sample rate fsWith signal bandwidth B；
Step 1-3) establish the motion model of source emission platform and receiving platform:
Three-dimensional cartesian coordinate system xyOz is established in space under water；Set the movement speed of source emission platform and receiving platform
Spending vector is respectively vs(t)=(vsx(t), vsy(t), vszAnd v (t))r(t)=(vrx(t), vry(t), vrz(t))；
Set channel maximum delay extended by taumaxWith signal transmission time Tmax；Source emission platform motion profile is calculated
Ps(t):
Ps(t)=(psx(t), psy(t), psz(t)), psx(t), psy(t), psz(t) P is respectively indicateds(t) x coordinate, y is sat
Mark and z coordinate；Ps0=(xs0, ys0, zs0) be source emission platform initial value；T ∈ [0, Tmax]；
Calculate receiving platform motion profile:
Pr(t)=(prx(t), pry(t), prz(t)), prx(t), pry(t), prz(t) P is respectively indicatedr(t) x coordinate, y is sat
Mark and z coordinate；Pr0=(xr0, yr0, zr0) be receiving platform initial value；
Step 1-4) establish sea level height model:
Set ocean surface wind speed and wind direction；If wind speed is 0, sea does not have stormy waves, then it is assumed that sea level height is always 0 everywhere；
If wind speed is not 0, according to preset wave spectrum model, the sea level height function h changed over time is calculatedati(x, y's, t)
Sampling.
As a kind of improvement of the above method, the signal transmission time should be not less than signal duration T, initial transmission
The sum of delay and channel maximum delay extension three, it may be assumed that
Tmax≥T+|Ps0-Pr0|/c0+τmax. (1)
Wherein, c0Indicate the bulk sound velocity of channel.
As a kind of improvement of the above method, the step 2) is specifically included:
Step 2-1) the setting fast beat frequency of initial channelAnd initial channel number of snapshotsOrder changes
Generation number q=0, time sampling serial number n=0；
Step 2-2) using the source emission platform of step 1-3) and the motion model of receiving platform, when calculating current snap
It carvesThe spatial position coordinate P of lower source emission platforms(tn) and receiving platform spatial position coordinate Pr(tn)；
Current time t is obtained by interpolationnSea level height section in perpendicular locating for lower source emission platform and receiving platform
Data and sounding survey cross-sectional data；
Step 2-3) according to the position of current transmit-receive platform and sea, sea-bottom profile data, utilize acoustic ray tracing model
Obtain current snap moment tnThe path of lower channel reaches wave parameter, amplitude, phase and propagation delay including all propagation paths；
Step 2-4) according to current path wave reach parameter, current snap moment t is calculatednUnder discrete channel impulse ring
It answersM indicates time delay sampling sequence number；
Step 2-5) judge whether time serial number n is less thanIf a determination be made that certainly, then enable n increase by 1,
It is transferred to step 2-2)；Otherwise, by the channel impulse response at each snap moment Temporally first
It sequentially combines afterwards, obtains discrete time-varying channel impulse response h(q)[n, m]；
Step 2-6) judge whether the number of iterations q has reached upper limit Q, if a determination be made that certainly, it is transferred to step 2-
12), otherwise, it is transferred to step 2-7)；
Step 2-7) h is responded to discrete time-varying channel impulse(q)[n, m] does Fast Fourier Transform (FFT) along snap time dimension, obtains
To Discrete Time-Delay-doppler spread function u(q)[k, m], k=0,1 ..., K-1 indicate that Doppler domain sampling sequence number, K are FFT point
Number；And the aliasing degree ρ of doppler spread spectrum is calculated on this basis；
Step 2-8) judge whether doppler spread spectrum aliasing degree ρ is less than or equal to pre-determined threshold ρalias, if it is determined that knot
Fruit is affirmative, is transferred to step 2-12)；Otherwise, it is transferred to step 2-9)；
Step 2-9) the number of iterations q increase by 1, enable snap doubling frequency:Channel number of snapshots are setEnable time serial number n=0；
Step 2-10) judge whether time serial number n is even number, if a determination be made that certainly, then directly according to previous
Secondary iteration result provides the current snap momentUnder discrete channel impulse response It is no
Then, step 2-2 is executed) to 2-4), the discrete channel impulse response inscribed when current snap is calculated
Step 2-11) judge whether time serial number n is less thanIf a determination be made that certainly, then enable n increase
1, it is transferred to step 2-10)；Otherwise, by the channel impulse response at each snap moment On time
Between sequencing combine, obtain discrete time-varying channel impulse response h(q)[n, m], is transferred to step 2-6)；
Step 2-12) iteration is terminated, obtain final channel dispersion time-varying impulse response h [n, m]=h(q)[n, m].
As a kind of improvement of the above method, the step 2-1) in will initial fast beat frequencyIt is set as 2 times most
Maximum Doppler frequency shift, it may be assumed that
Wherein, | vd| indicate the maximum relative motion in acoustical signal transmission process between source emission platform and receiving platform
Velocity amplitude.
As a kind of improvement of the above method, the step 2-4) it specifically includes:
Step 2-4-1) set current time tnObtained channel path is combined into up to wave parameter collectionWherein Npa
For the channel propagation paths sum at current time, ap、φpAnd τpWhen being expressed as the amplitude, phase and propagation of pth paths
Prolong, thus obtain the channel transfer functions H (f) at current time:
Wherein, τ0For with reference to time delay, f ∈ [fc- B/2, fc+B/2]；
Step 2-4-2) give fixed response time sample rate fτ, set channel impulse response lengthH (f) is carried out
Sampling, obtains the channel transfer functions H [l] of discretization:
Step 2-4-3) inverse fast fourier transform is done to H [l], obtain the punching of discrete channel corresponding to the current snap moment
Swash response:
Wherein, m indicates time delay sampling sequence number, m=0,1 ..., M-1.
As a kind of improvement of the above method, the step 2-7) it specifically includes:
Step 2-7-1) calculate channel Discrete Time-Delay-doppler spread function u(q)[k, m]:
Step 2-7-2) calculate channel doppler spread power spectrum D [k]:
Step 2-7-3) calculate the aliasing degree that channel Doppler spread is composed:
As a kind of improvement of the above method, the step 3) is specifically included:
Step 3-1) setting time series hitsWherein fsIt indicates signal passband sample rate, uses ht(τ)
It indicates transient channel impulse response corresponding to moment t, enables time sampling serial number n '=0；
Step 3-2), to discrete time-varying channel impulse response h [n, m] acquired in time point tn'=n '/fsPlace
Cubic spline interpolation is carried out, transient channel impulse response sample sequence h corresponding to the moment is obtainedn' [m]=
Step 3-3) with sample rate fsTo hn' [m] carries out linear interpolation, obtains time delay sample rate with emitting signal sampling rate
Consistent transient channel impulse response
Step 3-4) use hn' [m '] penetrates signal x [n] to base band recurrence and carries out time-varying FIR filtering, and it is a to obtain receiving end the n-th '
Baseband output signal samples y [n ']:
Step 3-5) judge whether time sampling serial number n ' is less than N-1, if a determination be made that certainly, then enable n ' increase
1, it is transferred to step 3-2)；Otherwise, output baseband receiving signals sampling sequence y [n '], n '=0,1 ..., N-1；
Step 3-6) give baseband receiving signals y [n '] to be multiplied by carrier wave, real part is taken, then be superimposed random noise, obtained final
Passband receives signal sample sequence:
Wherein,For passband random noise.
Present invention has an advantage that
1, method of the invention can adjust the fast beat frequency of channel according to the speed of real channel variance, with existing routine
Method is compared, and parameter selection is more reasonable, also more flexible；
2, method of the invention can largely avoid because channel snap set of frequency it is unreasonable caused by simulation accuracy
Decline or computing resource waste；
3, method of the invention guarantees the accuracy and precision of underwater mobile telecommunication channel emulation to a certain extent, together
When, it can effectively reduce the complexity of the Realization of Simulation, reduce memory requirements, improve computer operational efficiency.
Detailed description of the invention
Fig. 1 is underwater mobile telecommunication channel the Realization of Simulation block diagram of the invention；
Fig. 2 is underwater mobile telecommunication channel simulating scenes schematic diagram of the invention；
Fig. 3 (a) is initial time channel sound velocity profile；
Fig. 3 (b) is initial time channel static state acoustic propagation path structure figure；
Fig. 4 (a) is initial time sea level height schematic diagram；
Fig. 4 (b) is that wave of the sea composes schematic diagram；
Fig. 5 is sea-floor relief and transmit-receive platform initial position schematic diagram；
Fig. 6 (a) is flat pad movement locus schematic diagram；
Fig. 6 (b) is receiving platform movement locus schematic diagram；
Fig. 7 (a) is that time-varying channel impulse corresponding to the receiving channel 1 of simulation data responds schematic diagram；
Fig. 7 (b) is that time-varying channel impulse corresponding to the receiving channel 2 of simulation data responds schematic diagram；
Fig. 8 (a) is channel delay-doppler spread function schematic diagram corresponding to the receiving channel 1 of simulation data；
Fig. 8 (b) is channel delay-doppler spread function schematic diagram corresponding to the receiving channel 2 of simulation data；
Fig. 9 (a) is transmitting signal figure compared with the time domain waveform of simulation data signal；
Fig. 9 (b) is transmitting signal figure compared with the Time Domain Spectrum of simulation data signal.
Specific embodiment
The present invention will be described in detail in the following with reference to the drawings and specific embodiments
As shown in Figure 1, firstly, set in three-dimensional hydrospace to communication scenes, mainly include static environment and
The dynamic change modeling of the motion state modeling of the setting of signal of communication, sound source and receiving array and sea level height at any time.
Then, by repeatedly calling ray tracing model to carry out real-time sampling to the state of channel, the letter of different moments a series of is obtained
Road " snap ", and time-varying channel impulse response is calculated accordingly.When time-varying channel impulse response being recycled to carry out transmitting signal
Domain dynamic filter obtains receiving signal.Finally, display output time varying channel response and reception signal waveform.
1, communication scenes are set
Three-dimensional cartesian coordinate system xyOz is established in space under water.Assuming that the position coordinates of flat pad and receiving platform and
Movement velocity changes over time, and therefore, they can be expressed as to the form of the function about time t.Signal transmission time
Interior, the movement velocity vector of transmitting and receiving platform can be expressed as vs(t)=(vsx(t), vsy(t), vszAnd v (t))r(t)
=(vrx(t), vry(t), vrz(t)), emit and the three-dimensional coordinate of receiving platform position is respectively Ps(t)=(psx(t), psy(t),
pszAnd P (t))r(t)=(prx(t), pry(t), prz(t)).Channel simulator scene is as shown in Figure 2.The process of communication scenes setting
It is as follows.
(1) static channel environment is arranged: the range in setting underwater 3 D space, including horizontal coordinate range xmin~xmax、
ymin~ymaxWith sea area depth capacity zmax；Import in three-dimensional space velocity of sound distributed data and submarine topography data (depth measurement number in water
According to)；Set substrate parameter (compression velocity of wave, density and attenuation rate)；Set communications platform initial position Ps0=(xs0, ys0, zs0)、
Pr0=(xr0, yr0, zr0)。
(2) signal of communication is arranged: the signal of communication emitted needed for importingSampled data, and Setting signal length T,
Centre frequency fc, sample rate fsWith signal bandwidth B.
(3) dynamic is set: setting channel maximum delay extended by tau firstmaxWith signal transmission time Tmax.Signal transmission time
It should be not less than the sum of signal duration, initial transmission delay and channel maximum delay extension three, i.e.,
Tmax≥T+|Ps0-Pr0|/c0+τmax. (1)
Then, function v of the flat pad speed about the time is givens(t) and receiving platform velocity function vr(t), and accordingly
Flat pad motion profile is calculated
And receiving platform motion profile
Wherein, t ∈ [0, Tmax].Reset ocean surface wind speed and wind direction.If wind speed is 0, sea does not have stormy waves, then it is assumed that each
Locating sea level height is always 0；If wind speed is not 0, need to model pneumatic wave, according to preset wave spectrum model,
The sea level height function h changed over time is calculatedatiThe sampling of (x, y, t).Wherein level sampling grid can be according to practical need
It is set, Temporal sampling is consistent with signal bandwidth.
2, linear time-variant channel is tracked
Time varying channel is indicated with time-varying channel impulse response function h (t, τ).Since continuous letter can not directly be calculated
The analytical expression of number h (t, τ) needs to carry out discretization to time shaft, to multiple sampling instants in signal transmission time
Channel state information is tracked, and calculates the response of discrete channel impulse corresponding to these moment, finally with these channels
The combination of impulse response approaches h (t, τ).
The tracking of above-mentioned channel state information refers to: in a series of regular time points, by channel interface, propagation medium and
Communications platform all fixes, and obtains the sound field information of these time points reception position by Static Water propagation model.It will be different
Static sound field information corresponding to moment combines sequentially in time, and the dynamic change of approximate reflection channel is carried out with this.I
The process for obtaining channel state information in certain timing node is known as a secondary channel snap.Conjunction is determined by the means of iteration
The fast beat frequency of reason.One initial fast beat frequency is set first, is calculated to channel impulse response corresponding to all snap moment
After completion, judge whether fast beat frequency has kept up with the rate of change of channel impulse response by certain criterion；If cannot,
Fast beat frequency is doubled, is inserted into new snap between the every two adjacent snap moment, and calculate corresponding to these snaps
Channel impulse response, then judge whether new fast beat frequency meets the requirements；If cannot still meet the requirements, continue to repeat above-mentioned
Process, until fast beat frequency is met the requirements or the number of iterations reaches the upper limit.The specific implementation process of linear time-variant channel tracking is such as
Under.
Setting the number of iterations upper limit Q simultaneously enables the number of iterations q=0.The fast beat frequency f of initial channel is setsstAnd channel Doppler
The thresholding ρ of spread spectrum aliasing degreealias.In order to reduce the number of iterations as far as possible, by initial snap set of frequency be 2 times most mostly
General Le frequency displacement, i.e.,
Wherein, | vd| it indicates in acoustical signal transmission process, the maximum speed of related movement between sound source and receiving platform
Size, c0Indicate the bulk sound velocity of channel, fcThe centre frequency of signal is penetrated, B emits the bandwidth of signal.Thresholding ρalias∈ (0,1)
For judging whether current fast beat frequency can keep up with the rate of channel variation.It is required that ρaliasClose to 0.
(2) according to fast beat frequency, channel number of snapshots are setWith snap moment tn=n/fsst, n=0,
1 ..., Nsst-1.Then, to the channel path of each receiving channel corresponding to all newly-installed snap moment up to wave information into
Row tracking.With n-th of snap moment tnFor, according to the motion profile of the transmit-receive platform found out, find out current time hair
Penetrate the spatial position coordinate P of platforms(tn) and receiving platform spatial position coordinate Pr(tn).Using interpolation, current time is obtained
Sea level height cross-sectional data and sounding survey cross-sectional data on perpendicular locating for transmit-receive platform.It is set according to current environment
It sets, runs acoustic ray tracing program, the path for obtaining current time channel reaches wave parameter, amplitude, phase including all propagation paths
Position and propagation delay etc..Finally obtain all NsstThe channel path at a snap moment reaches wave parameter set.
(3) according to the channel path parameter being previously obtained, each corresponding to all newly-installed snap moment connect is calculated
Receive the discrete channel impulse response in channel.By taking n-th of snap moment tn as an example, if the channel path in certain channel that the moment obtains
It is combined into up to wave parameter collectionWherein NPaFor the channel propagation paths sum at current time, αp、φpAnd τpRespectively
It is expressed as the amplitude, phase and propagation delay of pth paths.The then channel transfer functions at the channel current time are as follows:
Wherein, τ0For with reference to time delay, f ∈ [fc- B/2, fc+B/2].Frequency domain discretization is carried out to channel transfer functions, then right
It does inverse fast fourier transform (IFFT), obtains the response of discrete channel impulse corresponding to current snap moment hn[m]=h
(tn, m/fτ), m=0,1 ..., M-1, wherein m indicates time delay sampling sequence number, fτFor time delay sample rate, time delay sampling sumM=0, which corresponds to, refers to delay, τ0.It should be noted that the considerations of for operation efficiency and memory, channel
The sample rate of impulse response should not be too high, need to only guarantee fτ>=2B.Finally obtain all NsstA snap moment corresponding letter
The set of channel shock responseTo these channel impulse responses according to the sequencing permutation and combination at snap moment
Get up, discrete time-varying channel impulse response h [n, m] can be obtained.
(4) to acquire discrete time-varying channel impulse response do Fast Fourier Transform (FFT) (FFT), calculate channel it is discrete when
Prolong-doppler spread function:
Wherein, k=0,1 ..., K-1 indicate that Doppler domain sampling sequence number, K are FFT points.The Doppler for defining channel expands
Open up power spectrum:
The aliasing degree of channel Doppler spread spectrum is defined on this basis:
(5) if Doppler aliasing degree ρ is greater than pre-determined threshold ρalias, wherein ρaliasFor one close to 0 positive real number, then
Think that Doppler's aliasing degree of discrete time varying channel counted at this time is higher, illustrates that the fast beat frequency of channel is too low, cannot keep up with
Channel changes with time.In this case, if the number of iterations is not up to the upper limit, fast beat frequency is doubled, in original
The mid-point at first all adjacent snap moment is inserted into a new snap moment, then repeatedly step (2)~(4), into
Row next iteration；If the number of iterations reaches the upper limit, iteration is terminated, exports discrete time-varying channel impulse response h [n, m].If ρ
No more than pre-determined threshold ρalias, then it is assumed that Doppler's aliasing degree of discrete time varying channel counted at this time is lower, illustrates channel
Fast beat frequency has been able to change with time with upper signal channel, then terminates iteration, export discrete time-varying channel impulse response h [n,
m]。
3, signal is received to calculate
Signal time-varying convolution is responded and emitted by time-varying channel impulse, can be received the signal at end.Time-varying convolution
Operation can be realized by time-varying FIR filter.Time-varying FIR filter requires the update of the sample rate and tap of filter tap
Frequency is identical with the sample rate of input signal.For the application of medium-high frequency shallow-sea underwater acoustic communication, what is be the previously calculated is discrete
Time-varying channel impulse responds the Temporal sampling f of h [n, m]sstWith time delay sample rate fτSignificantly less than the sample rate of transmitting signal
fs, it is thus impossible to which directly h [n, m] is filtered transmitting signal as time-varying FIR filter.Here, consider to h [n, m]
The interpolation in two dimensions is carried out, its sample rate and transmitting signal sampling rate exact matching on time shaft and time delay axis is made.
Directly carrying out two-dimensional interpolation in time dimension and time delay dimension to h [n, m] to emit signal sampling interval will cause huge memory and disappears
Consumption.Here it provides a kind of more feasible implementation method: when calculating the reception signal sampling value at a certain moment, passing through three first
Secondary spline interpolation obtains inscribing time delay sample rate when current time sample being fτDiscrete channel impulse response, then to the impulse
Response carries out linear interpolation, and obtaining sample rate is fsDiscrete impulse response, finally with this discrete impulse response to transmitting believe
Sample sequence in number correlation time block carries out FIR filtering, obtains the sampled value of the time sampling reception signal.
If baseband receiving signals sample sequence y [n] length is N, with the n-th ' a sampling instant tn′=n '/fsFor.First
By carrying out cubic spline interpolation to h [n, m], transient channel impulse response corresponding to the moment is obtained Again with sample rate fsTo hn′[m] carries out linear interpolation, obtains time delay sample rate with hair
Penetrate the consistent transient channel impulse response of signal sampling rate Utilize hn′[m '] penetrates signal to baseband transmission recurrence and carries out time-varying FIR filtering, obtains:
Then baseband output signal is multiplied by carrier wave, takes real part, then be superimposed with the random noise of certain power, obtained final
Passband receive signal sample sequence:
Wherein,For passband random noise.Noise spectrum level and distribution can independently be set as needed.
4, simulated example
If waters depth capacity is 120m.The velocity of sound is uniformly distributed in the horizontal direction, and vertical section is arranged such as Fig. 3 (a) institute
Show.Ocean surface wind speed 4m/s, sea out-of-flatness, the sea level height of each position change at random at any time.Using JONSWAP wave spectrum
Model models random sea level height, does not consider that wind direction bring angle extends.Initial time sea shape and ocean wave spectrum are such as
Shown in Fig. 4 (a) and Fig. 4 (b).Horizontal distance between flat pad and receiving platform is 500m.Emission depth is 80m.It receives
Platform is disposed vertically two hydrophones, and receiving depth is respectively 30m and 35m, corresponds respectively to receiving channel 1 and receiving channel 2.
If there are significant changes with the change of horizontal position in seabed depth.The initial placement and sea-floor relief of sound source and receiving array
As shown in Figure 5.If sound source radiation directive property angle of release is 150 °.Doppler's aliasing degree Doppler aliasing degree thresholding ρ is setalias=
0.2, maximum number of iterations Q=4.The initial time acoustic propagation path structure such as Fig. 3 obtained with BELLHOP acoustic ray tracing model
(b) shown in.The linear motion if flat pad remains a constant speed, speed vs=(0, -2,1) m/s；Receiving platform also remains a constant speed straight line
Movement, speed vr=(0,1,0) m/s.Flat pad and receiving platform are located remotely from each other.In signals transmission, communications platform
Motion profile as shown in Fig. 6 (a) and Fig. 6 (b) shown in.Fig. 7 (a) and Fig. 7 (b) are shown corresponding to two channels of simulation data
Time-varying channel impulse response.It can be seen that the channel in two channels all has very strong time variation.Communications platform level away from
From changing with time, each path of channel is caused to become larger at any time up to the propagation delay of wave.Meanwhile the movement of platform, companion
With the random variation and the violent fluctuating of sea-floor relief of wave of the sea, so that sound ray is in the reflection point position of interface and anti-
Firing angle degree constantly changes, to cause up to wave sum and change with time up to wave time delay and amplitude.Believed according to time-varying
Channel shock response function, acquire the delay-Doppler spread function in two channels as shown in Fig. 8 (a) and Fig. 8 (b) shown in.Fig. 9
It (a) is the comparison for receiving time domain plethysmographic signal and emitting signal waveform of simulation data.Fig. 9 (b) is that two of simulation data are logical
The comparison of the time-frequency spectrum and transmitting signal time-frequency spectrum of channel receiving signal.
It should be noted last that the above examples are only used to illustrate the technical scheme of the present invention and are not limiting.Although ginseng
It is described the invention in detail according to embodiment, those skilled in the art should understand that, to technical side of the invention
Case is modified or replaced equivalently, and without departure from the spirit and scope of technical solution of the present invention, should all be covered in the present invention
Scope of the claims in.
Claims (1)
1. a kind of underwater mobile channel emulation mode, which comprises
Step 1) sets the environmental parameter and signal parameter of three-dimensional hydrospace communication scenes, establishes source emission platform
With the motion model and sea level height model of receiving platform；
Step 2) obtains a series of different moments by repeatedly calling ray tracing model to carry out real-time sampling to the state of channel
Channel snap, and calculate accordingly time-varying channel impulse response；
Step 3) is responded using time-varying channel impulse and carries out time domain dynamic filter to the transmitting signal of source emission platform, is connect
The collection of letters number；
The step 1) specifically includes:
Step 1-1) setting underwater 3 D space range: horizontal transverse and longitudinal coordinate range and sea area depth capacity；Import three-dimensional space
Velocity of sound distributed data and submarine topography data in interior water；Set substrate parameter: compression velocity of wave, density and attenuation rate；
Step 1-2) import the required signal of communication emittedSampled data, and Setting signal length T, centre frequency fc, adopt
Sample rate fsWith signal bandwidth B；
Step 1-3) establish the motion model of source emission platform and receiving platform:
Three-dimensional cartesian coordinate system xyOz is established in space under water；Set the movement velocity arrow of source emission platform and receiving platform
Amount is respectively vs(t)=(vsx(t),vsy(t),vszAnd v (t))r(t)=(vrx(t),vry(t),vrz(t))；
Set channel maximum delay extended by taumaxWith signal transmission time Tmax；Source emission platform motion profile P is calculateds
(t):
Ps(t)=(psx(t),psy(t),psz(t)), psx(t),psy(t),psz(t) P is respectively indicateds(t) x coordinate, y-coordinate and
Z coordinate；Ps0=(xs0,ys0,zs0) be source emission platform initial value；t∈[0,Tmax]；
Calculate receiving platform motion profile:
Pr(t)=(prx(t),pry(t),prz(t)), prx(t),pry(t),prz(t) P is respectively indicatedr(t) x coordinate, y-coordinate and
Z coordinate；Pr0=(xr0,yr0,zr0) be receiving platform initial value；
Step 1-4) establish sea level height model:
Set ocean surface wind speed and wind direction；If wind speed is 0, sea does not have stormy waves, then it is assumed that sea level height is always 0 everywhere；If wind
Speed is not 0, and according to preset wave spectrum model, the sea level height function h changed over time is calculatedati(x, y's, t) adopts
Sample；
The signal transmission time should be not less than signal duration T, initial transmission delay and channel maximum delay extension three
The sum of person, it may be assumed that
Tmax≥T+|Ps0-Pr0|/c0+τmax. (1)
Wherein, c0Indicate the bulk sound velocity of channel；
The step 2) specifically includes:
Step 2-1) the setting fast beat frequency of initial channelAnd initial channel number of snapshotsEnable iteration time
Number q=0, time sampling serial number n=0；
Step 2-2) using the source emission platform of step 1-3) and the motion model of receiving platform, calculate the current snap momentThe spatial position coordinate P of lower source emission platforms(tn) and receiving platform spatial position coordinate Pr(tn)；It is logical
It crosses interpolation and obtains current time tnSea level height section number in perpendicular locating for lower source emission platform and receiving platform
According to sounding survey cross-sectional data；
Step 2-3) according to the position of current transmit-receive platform and sea, sea-bottom profile data, it is obtained using acoustic ray tracing model
Current snap moment tnThe path of lower channel reaches wave parameter, amplitude, phase and propagation delay including all propagation paths；
Step 2-4) according to current path wave reach parameter, current snap moment t is calculatednUnder discrete channel impulse responseM indicates time delay sampling sequence number；
Step 2-5) judge whether time serial number n is less thanIf a determination be made that certainly, then it enables n increase by 1, is transferred to
Step 2-2)；Otherwise, by the channel impulse response at each snap moment It is in chronological sequence suitable
Sequence combines, and obtains discrete time-varying channel impulse response h(q)[n,m]；
Step 2-6) judge whether the number of iterations q has reached upper limit Q, if a determination be made that certainly, it is transferred to step 2-12),
Otherwise, it is transferred to step 2-7)；
Step 2-7) h is responded to discrete time-varying channel impulse(q)[n, m] does Fast Fourier Transform (FFT) along snap time dimension, obtain from
Dissipate delay-Doppler spread function u(q)[k, m], k=0,1 ..., K-1 indicate that Doppler domain sampling sequence number, K are FFT points；And
The aliasing degree ρ of doppler spread spectrum is calculated on this basis；
Step 2-8) judge whether doppler spread spectrum aliasing degree ρ is less than or equal to pre-determined threshold ρalias, if a determination be made that
Certainly, it is transferred to step 2-12)；Otherwise, it is transferred to step 2-9)；
Step 2-9) the number of iterations q increase by 1, enable snap doubling frequency:Channel number of snapshots are setEnable time serial number n=0；
Step 2-10) judge whether time serial number n is even number, if a determination be made that certainly, then it is directly once changed according to preceding
The current snap moment is provided for resultUnder discrete channel impulse response Otherwise,
Execute step 2-2) to 2-4), the discrete channel impulse response inscribed when current snap is calculated
Step 2-11) judge whether time serial number n is less thanIf a determination be made that certainly, then it enables n increase by 1, turns
Enter step 2-10)；Otherwise, by the channel impulse response at each snap moment Temporally first
It sequentially combines afterwards, obtains discrete time-varying channel impulse response h(q)[n, m], is transferred to step 2-6)；
Step 2-12) iteration is terminated, obtain final channel dispersion time-varying impulse response h [n, m]=h(q)[n,m]；
The step 2-1) in will initial fast beat frequencyIt is set as 2 times of maximum doppler frequency, it may be assumed that
Wherein, | vd| indicate the maximum speed of related movement in acoustical signal transmission process between source emission platform and receiving platform
Value；
The step 2-4) it specifically includes:
Step 2-4-1) set current time tnObtained channel path is combined into up to wave parameter collectionWherein NpaTo work as
The channel propagation paths sum at preceding moment, ap、φpAnd τpIt is expressed as the amplitude, phase and propagation delay of pth paths, by
This obtains the channel transfer functions H (f) at current time:
Wherein, τ0For with reference to time delay, f ∈ [fc-B/2,fc+B/2]；
Step 2-4-2) give fixed response time sample rate fτ, set channel impulse response lengthH (f) is sampled,
Obtain the channel transfer functions H [l] of discretization:
Step 2-4-3) inverse fast fourier transform is done to H [l], it obtains discrete channel impulse corresponding to the current snap moment and rings
It answers:
Wherein, m indicates time delay sampling sequence number, m=0,1 ..., M-1；
The step 2-7) it specifically includes:
Step 2-7-1) calculate channel Discrete Time-Delay-doppler spread function u(q)[k, m]:
Step 2-7-2) calculate channel doppler spread power spectrum D [k]:
Step 2-7-3) calculate the aliasing degree that channel Doppler spread is composed:
The step 3) specifically includes:
Step 3-1) setting time series hitsWherein fsIt indicates signal passband sample rate, uses ht(τ) is indicated
Transient channel impulse response corresponding to moment t enables time sampling serial number n '=0；
Step 3-2), to discrete time-varying channel impulse response h [n, m] acquired in time point tn′=n '/fsPlace carries out
Cubic spline interpolation obtains transient channel impulse response sample sequence corresponding to the moment
Step 3-3) with sample rate fsTo hn′[m] carries out linear interpolation, and it is consistent with transmitting signal sampling rate to obtain time delay sample rate
Transient channel impulse response
Step 3-4) use hn′[m '] penetrates signal x [n] to base band recurrence and carries out time-varying FIR filtering, obtains a base band in receiving end the n-th '
Output signal samples y [n ']:
Step 3-5) judge whether time sampling serial number n ' is less than N-1, if a determination be made that certainly, then n ' increase by 1 is enabled, is turned
Enter step 3-2)；Otherwise, output baseband receiving signals sampling sequence y [n '], n '=0,1 ..., N-1；
Step 3-6) give baseband receiving signals y [n '] to be multiplied by carrier wave, real part is taken, then be superimposed random noise, obtains final passband
Receive signal sample sequence:
Wherein,For passband random noise.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810796972.4A CN109039506B (en) | 2018-07-19 | 2018-07-19 | A kind of underwater mobile channel emulation mode |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810796972.4A CN109039506B (en) | 2018-07-19 | 2018-07-19 | A kind of underwater mobile channel emulation mode |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109039506A CN109039506A (en) | 2018-12-18 |
| CN109039506B true CN109039506B (en) | 2019-09-06 |
Family
ID=64644226
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810796972.4A Active CN109039506B (en) | 2018-07-19 | 2018-07-19 | A kind of underwater mobile channel emulation mode |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109039506B (en) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109828264B (en) * | 2019-02-28 | 2022-07-19 | 东南大学 | Underwater acoustic channel phase-frequency response correction method based on RAM model |
| CN110138467B (en) * | 2019-05-15 | 2020-06-23 | 北京大学 | Underwater electric field communication analysis method |
| CN110932796A (en) * | 2019-12-20 | 2020-03-27 | 上海创远仪器技术股份有限公司 | Method for realizing Doppler spectrum generation processing applied to channel simulator platform |
| CN112287527A (en) * | 2020-10-15 | 2021-01-29 | 西北工业大学 | C # based active sonar signal detection simulation method under strong reverberation background |
| CN112684407B (en) * | 2020-12-31 | 2023-10-27 | 深圳市智慧海洋科技有限公司 | Underwater sound signal direction of arrival estimation method and device, underwater sound equipment and storage medium |
| CN115038165B (en) * | 2022-05-17 | 2023-05-12 | 上海船舶运输科学研究所有限公司 | Combined estimation method for target position and environment propagation parameters of underwater wireless sensor network |
| CN115085840B (en) * | 2022-06-15 | 2023-10-31 | 南京捷希科技有限公司 | Complicated mobile time-varying wireless channel simulation method based on ray tracing |
| CN117169862B (en) * | 2023-08-09 | 2024-03-26 | 中国科学院声学研究所 | Deep sea broadband signal waveform rapid simulation method and system based on ray acoustics |
| CN117375731B (en) * | 2023-10-27 | 2024-05-28 | 中国科学院声学研究所南海研究站 | Quantitative analysis method and system for time-varying characteristics of underwater acoustic communication multi-path channel |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101551452A (en) * | 2008-04-01 | 2009-10-07 | 中国科学院声学研究所 | Compensating method with associated movement of synthetic aperture sonar and system |
| CN105144598A (en) * | 2012-12-03 | 2015-12-09 | 数字无线功率有限公司 | Systems and methods for advanced iterative decoding and channel estimation of concatenated coding systems |
| CN106411438A (en) * | 2016-11-02 | 2017-02-15 | 东北农业大学 | Shallow water time-varying multi-path underwater acoustic channel modeling method |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20110075649A1 (en) * | 2009-09-13 | 2011-03-31 | Research Institute Of Tsinghua University In Shenzhen | Method and system of frequency division multiplexing |
| US9172476B2 (en) * | 2012-03-09 | 2015-10-27 | The United States Of America As Represented By The Secretary Of The Army | Method and system for removal of noise in signal |
-
2018
- 2018-07-19 CN CN201810796972.4A patent/CN109039506B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101551452A (en) * | 2008-04-01 | 2009-10-07 | 中国科学院声学研究所 | Compensating method with associated movement of synthetic aperture sonar and system |
| CN105144598A (en) * | 2012-12-03 | 2015-12-09 | 数字无线功率有限公司 | Systems and methods for advanced iterative decoding and channel estimation of concatenated coding systems |
| CN106411438A (en) * | 2016-11-02 | 2017-02-15 | 东北农业大学 | Shallow water time-varying multi-path underwater acoustic channel modeling method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109039506A (en) | 2018-12-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109039506B (en) | A kind of underwater mobile channel emulation mode | |
| Morozs et al. | Channel modeling for underwater acoustic network simulation | |
| Liu et al. | Biologically inspired covert underwater acoustic communication by mimicking dolphin whistles | |
| CN106254010B (en) | A kind of time-varying ocean channel modeling method | |
| Peterson et al. | Ray/beam tracing for modeling the effects of ocean and platform dynamics | |
| CN110535537B (en) | Underwater communication and detection integrated method | |
| CN111580110B (en) | Composite code underwater acoustic ranging method based on shallow sea multipath time delay | |
| Bjerrum-Niese et al. | A simulation tool for high data-rate acoustic communication in a shallow-water, time-varying channel | |
| CN107769862A (en) | A kind of bionical low communication interception method | |
| CN109347568A (en) | A kind of polynary frequency modulation(PFM) underwater acoustic communication method of imitative dolphin whistle continuous phase | |
| CN106411438A (en) | Shallow water time-varying multi-path underwater acoustic channel modeling method | |
| Zhou et al. | Study of propagation channel characteristics for underwater acoustic communication environments | |
| Wolff et al. | Acoustic underwater channel and network simulator | |
| CN110850396B (en) | Electric simulator applied to deep sea black box search and exploration positioning system and track generation method thereof | |
| Lou et al. | Basic principles of underwater acoustic communication | |
| CN116232478A (en) | Underwater non-fixed node communication method based on deep learning and deep migration learning | |
| Roudsari et al. | Channel model for wideband time-varying underwater acoustic systems | |
| Jia et al. | The study on time-variant characteristics of under water acoustic channels | |
| Chen et al. | Signal transmission using underwater acoustic vector transducers | |
| CN115021829B (en) | Digital pulse interval modulation underwater acoustic communication method based on marine environmental noise | |
| Xiang-ping et al. | Analyzing the performance of channel in Underwater Wireless Sensor Networks (UWSN) | |
| CN108494513A (en) | Shallow Water Acoustic Channels model foundation and its computational methods | |
| CN111162845A (en) | Sea area sight distance channel generation method | |
| Lasota et al. | Transmission parameters of underwater communication channels | |
| Murano et al. | Doppler resilience evaluation of different encoding techniques for underwater acoustic ranging systems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20201231 Address after: 311800 Room 101, 1st floor, 22 Zhongjie building, 78 Zhancheng Avenue, Taozhu street, Zhuji City, Shaoxing City, Zhejiang Province Patentee after: Zhejiang wanghaichao Technology Co.,Ltd. Address before: 100190, No. 21 West Fourth Ring Road, Beijing, Haidian District Patentee before: INSTITUTE OF ACOUSTICS, CHINESE ACADEMY OF SCIENCES |
|
| TR01 | Transfer of patent right |