A recurrent plot based stochastic nonlinear ray propagation model for underwater signal propagation

A stochastic nonlinear ray propagation model is proposed to carry out an exploration of the nonlinear ray theory in underwater signal propagation. The recurrence plot method is proposed to quantify the ray chaos and stochastics to optimize the model. Based on this method, the distribution function of the control parameter δ is derived. Experiments and simulations indicate that this stochastic nonlinear ray propagation model provides a good explanation and description on the stochastic frequency shift in underwater signal propagation.


Introduction
Our oceans play an important role in daily life and future development of technology and economy [1][2][3]. However, indiscriminate and high-risk exploration is placing great pressure on the ocean ecosystem, leading to pollution, biodiversity declining, and even climate change [4][5][6][7][8]. Therefore, we cannot precisely and efficiently detect and track a specific signal after a long range underwater transformation, especially in the deep ocean. One of the primary difficulties to solve this problem is that underwater propagation principles in complex marine environments are unknown. During resource detection and exploration, the initial signal might be seriously 'damaged' in propagation (as depicted in the upper side of figure 1) leading to wrong decision-making and resource location deviation. Therefore, it is significant to have a precise and effective model of underwater signal propagation.
The complexity of an underwater signal propagation channel is mainly relevant to multi-paths and multi-effects, as depicted in figure 1. The description can be as follows; (1) an acoustic signal has a weak directivity feature, the direction angle raises with the increasing of the distance and the signal transforms like multi-path, see the 'propagation paths' in figure 1 [9]. (2) The ocean ambient noise is affected by multiple elements, such as thermal current, marine traffic, marine organism, etc, which contain many nonlinear components. (3) The propagation trajectories are affected by internal waves and turbulence (IWaT), which can lead to chaos and stochasticity [10].
There have been many ways to model underwater signal propagation in both statistic and nonlinear theories. The underwater signal is primarily modelled as, Figure 1. The ocean ambient noise is composed of thermal current, hydrodynamics, marine traffic, marine organism, sea quake, internal waves, etc. A sent signal is affected by this noise as well as by the multi-path effect in its propagation. Hence, the signal might be damaged by nonlinearities and noises after a long range propagation.
propagating signal. However, they can barely depict and interpret the time-varying, frequency shifting, and nonlinear features during the propagation. Some studies use multiple hydrophone arrays to strengthen the s t (t − x) or to reduce the complex ambient noise b(t, x) [13,14], which is costly and hard launching in the ocean. Time frequency operators [b(t, x) varies in both t and x] and nonlinear filters [b(t, x) varies nonlinearly] were also used to enhance the signal-to-noise ratio [15,16]. However, these studies neglect the multi-paths feature, which may lead to false-alarm and wrong calculation of the detection probability. The multi-path effect has also been widely studied based on the linear ray theory for underwater communication, signal detection and channel estimation [the x in b(t, x) is two-dimensional] [17][18][19][20]. The propagation channel can be estimated precisely under relatively ideal circumstances using Green's function with the assumption of multi-paths effect, and the propagating signal can be enhanced combining the multi-path effect with the time-reversal method. However, these multi-path-based methods are all linear, which neglect the IWaT, which is ubiquitous in the ocean, and can strongly affect the propagation trajectories. Furthermore, nonlinear features can severely affect the intensity and frequency of the propagating signal. To sum up, the propagation channel b(t, x) in (1) is relevant to both t − x nonlinear (temporal nonlinear) and two-dimensional x nonlinear (spatial nonlinear).
The temporal nonlinearity has been widely studied in underwater acoustic propagation in multiple subjects, including relaxing fluids, bubbly liquids, fluids in saturate porous solids, etc [21,22]. They mostly built extended nonlinear acoustic equations based on a fundamental wave equation by adding nonlinear terms to model the temporal nonlinearity. However, the spatial nonlinearity has been comparatively rarely included in the underwater signal processing. The main reason is (1) IWaT are assumed not strong and general enough to affect the propagation channel; (2)the influence of IWaT on the propagation channels is hard to quantify; (3) the generation and variation of IWaT have no regularities. It is important to emphasise that, the hydrological environment is changing along with the increase of human maritime activities, such as marine traffic. Different kinds of IWaT are found in many areas using new techniques and equipment [23][24][25]. IWaT model opposes to the assumption that IWaT are not strong and ubiquitous enough to affect the propagation channel. The underwater propagation theory demonstrates that chaos and stochastic occur in underwater propagation channels [26,27]. To apply the spatial nonlinearity into underwater propagation, we should, on the one hand, find a way to quantify the nonlinear influence of IWaT; and, on the other hand, provide a description of the irregular appearance by a new ray theory propagation model.
In this paper, we focus on two principles to explore the application of the nonlinear ray theory to underwater signal processing. First, to quantify the chaotic and stochastic features of ray trajectories, the nonlinear features of the nonlinear ray theory under different parameters are analyzed by Poincare sections and maximum Lyapunov exponents. Then, based on the recurrence plot analysis, a quantitative method is proposed. Second, for the generation and variation of the irregularity of IWaT, we introduce a stochastic control parameter δ. We get the distribution of δ which is related with two other stochastic parameters' distributions by the recurrence plot quantitative method, and the stochastic control parameter is optimized based on real data by the least square method. Then, the stochastic nonlinear ray propagation (SNRP) model is proposed. Finally, we use the SNRP model to interpret the frequency transform phenomenon in a lake experiment. This paper is organised as follows: section 2 introduces the fundamental theory of the wave equation and derivative process of the ray theory. The stochastic nonlinear ray propagation model is built in section 3. A nonlinear behavior analysis is proposed in section 3.1 by using Poincare section and maximum Lyapunov exponent. The recurrence plot quantitative method is proposed in section 3.2. Using this quantitative method, the control parameter δ distribution function in the model is proposed in section 3.3. In section 4, we provide an application of the stochastic nonlinear ray propagation model in underwater single-frequency signal propagation.

Basic theory: wave equation and ray theory
The propagation of an underwater acoustic signal follows the wave equation derived by the features of the underwater circumstances. If the parameters are constant, the wave equation is stable to express the ideal propagation circumstances. If part of the parameters are not constant, the wave equation can be chaotic or stochastic with a variation of the parameters to express a more realistic underwater environment. The wave equation can be transformed into the ray equation to express the propagation channel geometrically. This section introduces the fundamental wave equations and the ray theory.

Frequency domain plane wave equation
The wave equation in an 'ideal fluid' is based on the following basics, the equation for conservation of mass, Euler's equation, and the adiabatic equation of state, respectively [28].
∂v ∂t and define, where c is the speed of underwater acoustics, ρ the density, v the particle velocity, p the pressure, and the subscript S denotes that the thermodynamic partial derivatives are taken at constant entropy. The subscript 0 of p 0 in (4) denotes a quiescent (time independent) medium, and, then, the pressure and density can be expressed as, p = p 0 + p , ρ = ρ 0 + ρ , with small perturbations p and ρ . To get the fundamental wave equation, (2)-(4) are linearized as, As the changes of the oceanographic time scale is much longer than the acoustic signal propagation, both ρ 0 and c 2 can be assumed of the independent of time, and the density ρ is assumed constant in space [28]. Then the standard wave equation of pressure is derived, Equation (9) can be reduced to three dimensions by the Fourier transformation, called frequency domain wave equation, or Helmholtz equation.
where k(r) is the medium wavenumber, k (r) = ω c (r) (11) where ω is frequency and c(r) is the geometry-dependent sound speed. Under the Cartesian coordinate (x, y, z), the y-axis can be neglected, because the underwater acoustic is a spherical or a cylindrical wave. Hence, the x-axis is equal to the y-axis, as depicted in figure 1. The two dimensional Helmholtz equation is, There are several ways to model the underwater propagation based on equation (12), such as normal mode, parabolic approximation, ray theory, etc. The ray theory has no frequency limitation compared to other approximation methods and can acquire the propagation time as the ray theory is on the geometry perspective [28].

Nonlinear ray theory equation
The ray propagation model expresses the underwater acoustic signal propagation in geometric mode by range and angle. The reliability and performance of such systems depend on the environmental conditions, noises, depth, internal waves, etc [28,29].
The nonlinear ray theory equation is based on the Helmholtz equation, a standard parabolic equation, and a parabolic approximation [26]. The solution of equation (12) can be expressed as, (13) where k = ω/c 0 is the wave number with a constant sound speed c 0 and p (x, z, ω) is the complex-valued field amplitude. Substituting (13) into (12) we get the standard parabolic equation, where If the wavelength λ = 2π/k is small compared with the propagation range, the solution of (14) can be assumed as, where A(x, z) is the amplitude and S(x, z) is the ray eikonal.
A(x, z) and S(x, z) can be expressed via the Hamilton's formalism to describe the trajectories. Each trajectory at distance x can be characterized by the depth z and the generalized momentum g, where g = tan θ, and θ is the ray grazing angle. Then, the fundamental ray theory equations can be formalized by, where H(x, g, z) is the Hamiltonian of the ray system with the form, where k 0 = ω/c 0 , ω is the acoustic frequency. As the purpose of this paper is analyzing the characteristics of the ray chaos, it is reasonable to set the frequency equals to c 0 , i.e. k 0 = 1 [30]. System (16) and (17) is the geometric expression of system (12) and can be solved to model the propagation of the underwater wave. As the circumstance features change in the process of propagation, the propagation channels might be disturbed by IWaT or noises. The disturbance is modelled by the sound speed c(x, z) nonlinear variation.
The underwater sound speed c is affected by temperature, salinity, and pressure. The water is layered by different temperature, salinity, and pressure in different depth, so the gradient of c varies with depth [30]. On the one hand, the underwater sound profile, especially in the ocean, exhibits different kinds of gradients in different areas or seasons. On the other hand, the variation of the sound speed profile is relatively small in a specific area, and the profiles can be calculated according to experimental data and gradient approximation [29,31]. So, the sound speed equation in this paper is presented by the canonical sound speed profile (Munk profile).
According to (14), the canonical sound speed profile can be expressed as, where with the depth z defined positive downwards and η = 2(z − z a )/B, z a is the sound channel axis depth, B the depth scale, c 0 the reference sound speed, δ dimensionless magnitude, and λ the wavelength [32,33].
V(x, z) describes the vibration contributed by the variation of temperature, salinity, and pressure, and can be regarded as an internal-wave-induced sound speed perturbation (the δ can be seen as the strength of the internal wave in this sense) which is one of the main contributions of ray chaos.

Stochastic nonlinear ray propagation model
Nonlinear ray theory equations (16)-(18) represent the ray eikonal in (15) to model the geometry aspect of the acoustic underwater propagation, and V(x, z) in (18) provides a description of chaos or stochastics. According to (16)- (18), the fundamental ray theory equations can be expressed by, where δ is the dimensionless magnitude. The nonlinear part V(x, z) in (18) drives the ray theory equation (19) to fit for nonlinear properties caused by IWaT. The main parameters in V(x, z) to control the nonlinearity are δ and λ, which can be interpreted as internal wave intensity and wave length, respectively. In the meantime, the different initial inputs g 0 and z 0 can also lead to a distinct nonlinear output because of the nonlinear part. In a specific detection or monitoring situation, each target having comparatively the same g as the underwater acoustic wave is spherical; however, the initial depth is different and stochastic among each targets. So, the three parameters δ, λ, and z 0 primarily affect the nonlinear property of the nonlinear ray theory equation (19) output.
As we illustrated above, the appearance of IWaT has no regulation or pattern, i.e. the appearance of nonlinearity or chaos is random in the time and spatial domain. The values of the nonlinear control parameters δ and λ are varying in different circumstances. So, there are three stochastic parameters in the SNRP model, where Δ is the sample space of the stochastic variable δ, Λ is the sample space of the stochastic variable λ, and Z 0 is the sample space of the stochastic variable z 0 which we cannot see in the equation.
The three stochastic parameters can model the random appearance of nonlinearity or chaos. First, to demonstrate the stochastic variables, we should provide definitions of the probability triple: sample space, set of events, and assignment probabilities to the events. Based on the weak nonlinear property of underwater propagation [34], δ belongs to a bounded interval Δ = [0, δ m ], which represent the sample space. A large λ leads to relatively fast nonlinear variations and vice versa, and the IWaT varies in a bounded speed. To model the propagation channel variation caused by IWaT [35], the wavelength Λ should also be bounded, Λ = [λ min , λ max ], where λ min and λ max are distinct in different hydrological environments. The initial depth is also bounded as the depth of ocean is limited, Z = [z min , z max ]. Secondly, the set of events are δ ∈ Δ, λ ∈ Λ, and z ∈ Z, respectively. Thirdly, the probability measures P 1 : Δ → [0, 1], P 2 : Λ → [0, 1], and P 3 : Z → [0, 1], which are countably additive and P 1 (Δ) = P 2 (Λ) = P 3 (Z) = 1. So, we should get three effective distribution functions of three parameters not only conform the features of the probability measure, but also fit for real propagation characters.
There are several methods to get optimal parameters of distribution functions, such as least square methods and machine learning methods. However, on the one hand, to optimize the three distribution functions in the same time need tedious computation. On the other hand, the three parameters are interacted among themselves leading to different kinds of nonlinear outputs. We choose some typical values of each parameter to analyze nonlinear output features using Poincare slices and the maximum Lyapunov exponents to identify the relationship among the parameters.

Nonlinear behavior analysis
Nonlinear features of the ray theory equation (20) are mainly represented as chaotic and/or stochastic behavior in phase space. Poincare sections can reveal the non-normal behaviors and depict them visually. On the other hand, non-normal and normal trajectories diverge exponentially or according to a power law, respectively. Lyapunov exponent quantifies the divergence to verify chaotic or stochastic properties [27]. In this section, we analyze the nonlinear behavior of the ray theory equation (20) by the Poincare sections and the maximum Lyapunov exponents.
According to the canonical sound profile (18), we set z a = 1.3 km, B = 1.3 km, = 7.4 × 10 −3 in our simulation. The initial angles are set from 5 • to 50 • , and the step length is 1 • . If the initial angle was larger than a threshold, the propagating acoustic signal reflected many times that the corresponding ray energy decreased quickly. The angle threshold increases with depth, and under the canonical sound profile, we set a relatively large angle scope, from 5 • to 50 • .
The representative parameters are selected to show and analyze chaotic and stochastic behaviors in different circumstances. The wavelength λ are set as 1 km, 4 km, 7 km, and 10 km. It is related with the sound profile. Under the canonical sound profile, the center wavelength which can lead to chaos is from 4 km to 7 km. We select a larger range in the simulations from 1 km to 10 km. The initial depths z(0) are set as 0.01 km, 0.1 km, 1 km, and 2 km. The underwater targets are mostly active from 0.01 km to 2 km, where 0.01 km represents the water surface, 0.1 km the shallow water, 1 km the depth of the eigen-ray, and 2 km the deep water. The internal wave strengths are set as 0.001, 0.005, 0.01, and 0.05. It is also related with the sound profile and the reason to choose these values is the same as the wavelength.

Poincare section
The ray trajectories lie in a bounded three-dimensional space (g, z, x). If the system is integrable, there is a two-dimensional slice which all trajectories lie on as a constant. The slice becomes changeable when the system is non-integrable. It is called the Poincare section which can distinguish the chaotic from periodic behavior of ray trajectories. Poincare sections are extracted by stroboscopically viewing trajectories [36].
where n = 0, 1, 2, . . . , and λ is the wavelength. Figure 2 shows a set of Poincare sections of four wavelength values λ = 1 km, 4 km, 7 km, 10 km, with initial depth z 0 = 1 km, and internal wave strength δ = 0.01. For these parameter circumstances, the ray trajectories fill more speckled regions (non-normal islands) as λ increases. Different wavelengths of the internal waves lead to different non-normal behaviors. With the increase of λ the normal regions in the center (curves without islands) are shrinking. Furthermore, with the rise of the angle, the Poincare sections in figures 2(b)-(d) depict no islands, called stochastic region, and with the raise of λ, the stochastic region becoming larger. That is, with the increasing of the interference, the ray trajectories express increasing non-normal features. Figure 3 shows a set of Poincare sections for four initial depth values, z(0) = 0.01 km, 0.1 km, 1 km, 2 km. As the variation of the initial depth value, both the number of normal trajectories and the type of speckled regions are different. With the increase of the initial depth from 0.1 km to 1 km, the center blank area strongly shrinks. The phase space have both a chaotic area and a stochastic area, and the boundaries between normal and non-normal trajectories are not clear. But we can see an increasing trend of the normal trajectories as z(0) becomes larger. Figure 4 shows a set of Poincare sections for four internal wave strength values, δ = 0.001, 0.005, 0.01, 0.05, as the strengths in this region behave mostly non-normal with the internal wavelength λ = 5 km. The normal region becomes dominating as the wave strength decrease, and with the rise of the wave strength, the non-normal region (the chaotic area and stochastic region) is getting bigger.
These Poincare sections illustrate the appearance of normal and non-normal trajectories relatively related with the three parameters in equation (19), and depict the changing trend for different parameters. However in some circumstances, the variation is not very clear. Here, we should quantify these islands and stochastic regions to verify its chaotic and stochastic nature.

Maximum Lyapunov exponent
Chaotic trajectories stretch, fold and squeeze toward one another with the growth of range (or time). The uncertainty of these variations causes two manifolds to squeeze closer to one another with uncertainty fluctuations. The maximum Lyapunov exponent expresses the exponential variation of different trajectories quantitatively [37].
Here d(x) is a measure of a pair neighbor trajectories in phase space. The measure expresses the separation of two chosen trajectories,  Figure 2(c) shows the transformation from a normal region to a stochastic one; figure 3(b) shows the transformation from a normal region to a chaotic one and then to a stochastic one; figure 4(a) shows a normal region. The Lyapunov exponent (22) should be calculated in a sufficiently long-range to reflect structures of the system, the data length expands to 20 000 points. Six angles are chosen to represent the maximum Lyapunov exponents' variation with the increase of the angle under each circumstance, and two of them are picked to analyze.
According to figures 5-7, all the maximum Lyapunov exponents of different trajectories are positive. From the Lyapunov exponent point of view, the normal trajectories are not 'real periodic', as the value is not zero. The maximum Lyapunov exponents' value decrease with periodic vibration illustrates that each distance in phase space is not constant, called 'gradual change periodic'. The normal and non-normal trajectories can not be distinguished just by the sign of the Lyapunov exponents' values. Furthermore, the value changes less regular with the increase of the initial angle. The non-normal trajectories defined by Poincare slices are chaotic in the sense of maximum Lyapunov exponents. The maximum Lyapunov exponents can express the differences between normal trajectories and non-normal trajectories by the distinct variation as depicted in figure 5.
From the analysis above, the appearance of chaotic and stochastic trajectories is related to δ, λ, and z(0), respectively. The differences between normal and non-normal trajectories can be recognized via maximum Lyapunov exponents. But the geometric characteristics are neglected, and one threshold is needed to express the differences.

Recurrence plot
The recurrence plot (RP) visualizes the trajectories recurrences in phase space. It transforms the trajectory trend to a 0-1 matrix, which can not only be transferred into a two-dimension graph to visualize the chaos characteristic, but also to provide quantitative parameters. Both the graph and the parameters behave differently among chaotic, stochastic and periodic regimes. They can be transformed into quantification in both one and two dimensions to satisfy the model demand of signal processing [38].
The RP matrix can be expressed as, where N is the number of measured points, is a threshold, · is metric (Euclidean distance, for example), and Θ(·) the Heaviside function, Under different choices of , the definition equals to a notion, The threshold is a crucial parameter for RP. If is too small, the RP has almost no recurrence points which can not express any recurrence structure. On the other hand, if the is too large, most points are neighbors, hence the 'real structure' will be hidden in artifacts [38]. This paper uses the position of phase space points to determine . On the one hand, from the Poincare slices in figures 2-4, the ray chaos appears  in the form of 'small islands'. So the differences between normal and non-normal trajectories in the phase space are on a small scale. For another, with the increase of the initial angle, the range of the phase space is getting bigger. According to our simulations, is set as 1% of the maximum phase space diameter [39,40].
Reconstitution of the time series is another important step of RP. Reconstruction methods are based on Takens' theorem [41], using time-delay embedding with delay τ and embedding dimension m. According to our simulations and the chaotic property in the Poincare slices, the delay τ is set as the wavelength λ in the propagation system (16) and (17), and m = 2.
where m is the embedding dimension and λ is the delay. RP is depicted by the matrix (24) to plot a black dot at the coordinates (i, j), if R i,j = 1, and a white dot if R i,j = 0. Both axis of the RP in the ray system are the range, and the RP lines show rightwards and upwards. The RP is symmetric by definition with respect to the main diagonal.  Recurrence points occur when the points satisfy the definition notion (25), and the Euclidean distance is chosen, In our simulations, the parameters of the propagation system (16) and (17) are set the same as in section 3.1.2, situation 1 to 3. Two representative RPs are chosen from each situation to show differences between normal trajectories and non-normal ones. Figures 8-10 show that the normal trajectories have RPs which are diagonal oriented. But the density of each diagonal is different and irregular. The non-stable density of diagonals has once again proven that the normal trajectories are not really periodic. On the other hand, the non-normal trajectories have RPs with isolated recurrence points and some proportional short diagonals. According to the distribution of the recurrence points, a quantification can be based on the density of the points to distinguish normal trajectories for non-normal ones.
The parameter recurrence rate (RR) can be expressed as, where n k is the number of trajectories. According to (28), the RR of the ray propagation system is depicted in figure 10.
From the RR estimation, we find a joint threshold, = 0.01, to distinguish normal and non-normal trajectories. The recurrence rate points which are larger than 0.01 correspond to the normal trajectories depicted in the Poincare slices, and the lesser ones corresponding to non-normal trajectories. Furthermore, figure 11(c) shows that the RR have fluctuations around 0.01 when θ > 25 • . The reason is that, when δ = 0.001, the outputs have relatively weak chaotic or stochastic features, as depicted in figure 4. The vibrations also meet with trajectories property in the Poincare section. So, based on the quantitative RP method, we can get the relationship among δ, λ, and z.

Distribution feature estimation of the SNRP equation
As shown, the wave strength δ, wave length λ, and depth z can lead to non-normal or normal trajectory outputs. We define n t as the nth trajectory that transforms a normal into a non-normal one, called transform trajectory. The relationship between δ, λ, z and n t can be depicted based on the quantitative method above (figure 12).
We find from figures 12(a) and (b) that, with the increase of δ and λ, n t monotonically decreases, respectively. Using the first degree polynomial equation to fit the data points, we have, where a δ , b δ , a λ , and b λ are polynomial fit coefficients.  From figure 12(c), the relation between n t and z exhibits a peak value around the axial sound channel, z a , which is the underwater acoustic convergence depth. We use a piecewise linear function to fit, where a z , b z , c z , and d z are polynomial fit coefficients.
According to (29)- (31), the relationship between δ and λ, z relatively are, A more precise relation can be derived with a higher degree polynomial equation, but that is not the focus of this article. The distribution function can also be derived based on real hydrology circumstances and some assumptions. On the one hand, λ is assumed as the internal wave length and bounded as analyzed above, but we do not really know the appearance location and time of the internal wave. So, we simply present the distribution of λ as a continuous uniform distribution, On the other hand, in underwater passive signal detection, the depth of the target is unknown and bounded as analyzed above. The depth even sometimes changes with time. According to the relationship (31), we define the distribution as a piecewise uniform distribution, Furthermore, a more precise distribution density function of z can be provided by analyzing the features of specific objects, and so we can get a more precise distribution of λ through multiple internal searches [42,43]. But these are not the focus of this article. The distribution of (34) and (35) is independent of the physical property as wave length and depth. We can finally get the distribution of the wave strength δ, which is the joint probability distribution of λ and z, Optimizing the distribution parameter of δ is also relevant to the affect on λ and z. The initial depth can be set as the axial sound channel depth z a , and the wave length set as the average value λ a of the canonical sound speed profile (18).
Finally, we get the SNRP model: where Δ ∼ f Δ (δ). The distribution parameters can be calculated by the least square method.

Application: an explanation of frequency shift
We did a signal propagation experiment to test the application of the SNRP model (37). The parameters of the model are optimized based on the experiment data. The model can provide a better explanation on frequency shift than classic linear ray theory. The experimental background is as follows, • Location: Danjiangkou Reservoir, China; • Total area: 1022.75 square kilometres; • Deepest water level: 170 m; • Devices: smart-hydrophone, signal emitter, hand-held GPS, two trial vessels (depicted in figure 13); • Signal feature: 10 kHz single-frequency signal; • Sample frequency: 100 kHz.
• Experiment arrangement and process: dip the smart-hydrophone into 10 m depth from the ship one broadsides, and dip the signal emitter into 10 m depth from the ship two broadside. Anchoring ship one in the middle of the reservoir, and drive ship two to different distances (use the hand-held GPS to confirm distance value), anchor ship two in each distance and send a signal.
There are a few points that should be mentioned about the experiment. First, we calculated the axial sound channel of the reservoir. However, the bottom varies in different locations, and the smart-hydrophone posture is not stable. So we set the depth as 10 m. Second, the scale of the reservoir can roughly imitate ocean circumstance, but nonlinear affects are mostly decreased. Third, the sound speed profile in (37) should be changed into the reservoir one. Forth, the reservoir background noise is analyzed, as depicted in figure 14, and use this kind of noise in the linear ray model as a comparison [28]. Fifth, during the signal collection both ships are anchored, so there is no Doppler effect in this experiment data ( figure 14).
We use four ray lines in the SNRP model to estimate the received signal frequency. The results are depicted in figure 15. On the one hand, the SNRP model can fit for the close distance(nearly linear) pattern, as depicted in (a) and (b) in line one in figure 15, as well as can provide a good description of the long  distance(weak nonlinear) pattern, as shown in (a) and (b) in line two and three in figure 15. In close distance, the frequency transform around 0.03 s, 0.06 s, and 0.12 s are estimated by the SNRP model. In long distance, the real data have a frequency shift from 1 to 2 kHz, but the appearance of the shift is irregular. The SNRP model successfully estimates and simulates the frequency shift and the appearance of a stochastic mode. More ray lines could be used to get a more accurate and effective result. On the other hand, the linear ray theory, as depicted in figure 15(c), can only represent the background noise. All three distances are represented as (c) in line one, as the signal intensity is large. To make a further explanation of the noise effect, we simulate two much longer distances, as depicted in (c) line two and line three. These simulations indicate that, in the linear ray theory, the frequency shift can not be simulated even in a very low signal-noise ratio, and the noise can just lead to signal lose.

Conclusion
In this paper, we have presented a stochastic nonlinear ray propagation model to simulate underwater signal propagation. The Poincare section analysis indicates that the nonlinear ray theory can simulate the non-normal trajectories by changing three characteristic parameters δ, λ, and z. The maximum Lyapunov exponent analysis exhibits that the non-normal trajectories caused by the three parameters have chaotic and stochastic features. We have used the recurrent plot method to quantify the non-normal feature. Then the relationship between δ and λ, z are inferred, respectively. The distribution function of the control parameter δ is provided according to its relationship with λ, z and the real data. Then the stochastic nonlinear ray propagation model is proposed to provide a better description of the underwater signal propagation.
Our experiment demonstrates that the SNRP model can offer a good explanation on the signal frequency stochastic shift in reality, which the linear ray model can not. Experiments in ocean will be done in the near future to get a better optimized SNRP model.
Our main contributions are as follows: a new model is proposed for the underwater acoustic signal propagation combining the recurrence plot method and stochastic parameter design with the nonlinear ray model. The proposed model applies the spatial nonlinearity into underwater propagation and quantifies the nonlinear influence of the internal waves and turbulence. The model is available for the multi-path and multi-effects complex underwater environments and provides a description of the irregular appearance of noise in acoustic signal propagation. The stochastic parameters in the model can effectively adapt different underwater circumstances. We demonstrate that the model has theoretical and practical value to the development of underwater signal processing.