Subharmonic ocean turbulence phase screen simulation method based on N×N sampling method

Abstract

This invention discloses a subharmonic ocean turbulence phase screen simulation method using an N×N sampling scheme. Step 1: Generate a Gaussian random complex matrix with a mean of 0 and a variance of 1. Step 2: Filter the random complex matrix using the ocean turbulence phase perturbation power spectrum function, and then perform an inverse Fourier transform to obtain a random phase screen. Step 3: Compensate the low-frequency component of the random phase screen generated in Step 2 using N×N subharmonic technology to obtain a compensated low-frequency subphase screen. Step 4: Superimpose the random phase screen obtained in Step 2 and the compensated low-frequency phase screen generated in Step 3 to obtain the final phase screen. This invention uses an N×N sampling subharmonic method to compensate for the insufficient low-frequency components in the power spectrum inversion method. It can improve the compensation effect not only by increasing the harmonic order but also by changing the compensation level and the next sampling region to improve the accuracy of the phase screen simulation. The optimal compensation effect is achieved when N = (q+1) / 2.

Description

Technical Field

[0001] The invention belongs to the technical field of underwater wireless optical communication, and in particular relates to a subharmonic ocean turbulence phase screen simulation method in an N×N sampling mode. Background Art

[0002] Underwater optical communication has the advantages of large bandwidth, high communication rate, strong confidentiality and low latency. When light is transmitted in seawater, it will be affected by ocean turbulence. In practice, due to the limitations of environmental and experimental conditions, the method of establishing a random phase screen is often used to simulate ocean turbulence. Therefore, it is necessary to study ocean turbulence and underwater optical communication by establishing an accurate ocean turbulence phase screen to simulate the ocean turbulence environment.

[0003] At present, the commonly used methods for numerical simulation of ocean turbulence phase screens include spectral inversion based on Fast Fourier Transform (FFT), Zernike polynomial method based on orthogonal polynomials, fractal method based on random midpoint displacement, generation method based on random data element expansion, and multi-scale method based on wavelet analysis. Since the spectral inversion method has the advantages of fast calculation speed and applicability to different forms of power spectrum models, it has become one of the most mainstream phase screen simulation methods, but it has the disadvantage of insufficient low-frequency components, which can be improved by subharmonic compensation technology. The accuracy of the phase screen generated by subharmonic compensation technology is related to the subharmonic order. As the harmonic order increases, the accuracy of the generated phase screen is higher. However, after the harmonic order exceeds a certain value, the accuracy of the phase screen no longer increases significantly. Continuing to increase the harmonic order can only increase the complexity of the algorithm, and even over-compensation may occur. How to generate a more accurate phase screen is an urgent problem to be solved. Summary of the invention

[0004] The purpose of the present invention is to provide a subharmonic ocean turbulence phase screen simulation method with an N×N sampling method, which solves the problem of insufficient low-frequency components in the traditional power spectrum inversion method and achieves the purpose of improving the accuracy of the simulated phase screen by increasing the harmonic order, changing the compensation level and the next sampling area.

[0005] The technical solution adopted by the present invention is: a subharmonic ocean turbulence phase screen simulation method in an N×N sampling mode, which specifically includes the following steps:

[0006] Step 1: Generate a Gaussian random complex matrix with a mean of 0 and a variance of 1;

[0007] Step 2: Use the ocean turbulence phase disturbance power spectrum function to filter the random complex matrix, and then perform inverse Fourier transform on it to obtain the random phase screen;

[0008] Step 3: Compensate the low-frequency part of the random phase screen generated in step 2 by the subharmonic method of variable N×N sampling mode to obtain the compensated low-frequency subphase screen;

[0009] Step 4: The random phase screen obtained in step 2 and the compensated low-frequency phase screen generated in step 3 are superimposed to obtain the final phase screen.

[0010] The power spectrum function of step 2 is:

[0011]

[0012] Where ε is the kinetic energy dissipation rate per unit seawater, and its value range from the depth to the surface of the seawater is [10 -10 ,10 -1 ]m 2 / s 3 ;

[0013] χ T is the mean square temperature dissipation rate of seawater, and its value range from the deep seawater to the seawater surface is [10 -10 ,10 -4 ]K 2 / s;

[0014] w is defined as the ratio of ocean turbulence caused by temperature to that caused by salinity, and its value range is [-5,0], where -5 means that the ocean turbulence is completely caused by temperature, and 0 means that the ocean turbulence is completely caused by salinity;

[0015] η is the internal scale of ocean turbulence,

[0016] A T =1.863×10 -2 ,A s =1.9×10 -4 ,A Ts =9.41×10 -3 .

[0017] Step 2 is implemented according to the following specific steps:

[0018] Step 2.1, based on the power spectrum function of the ocean turbulence phase disturbance, calculate the seawater phase spectrum on any slice perpendicular to the beam propagation direction at this time:

[0019] F φ (k x ,k y )=2πk 2 Δzφ(k x ,k y )(2)

[0020] Where Δz is the thickness of the turbulent layer or the propagation distance of the light beam, k = 2π / λ is the beam wave, and λ is the wavelength of the light beam;

[0021] Step 2.2, use F φ (k x ,k y ) for the Gaussian random complex matrix h(k x ,k y ) is filtered and then inverse Fourier transform is performed to obtain the spatial distribution of the phase screen:

[0022]

[0023] Where C is the constant factor controlling the phase screen variance (Δk x Δk y ) 1 / 2 , in the discretized space domain, x = mΔx, y = nΔy, Δx, Δy are the sampling intervals of the space domain, let Δx = Δy, m, n are integers; i is an imaginary number; x, y are the coordinates of the sampling points in the discretized space domain;

[0024] In the wave number domain k x =m'Δk x ,k y =n'Δk y , m', n' are integers; Δk x , Δk y is the sampling interval in the wavenumber domain, M is the number of grids of the phase screen; k x , k y are the coordinates of the sample points in the wavenumber domain.

[0025] Formula (3) can be written as:

[0026]

[0027] make The expression of the phase screen is obtained as follows:

[0028]

[0029] FFT refers to Fourier Transform.

[0030] The specific steps of using the subharmonic technology of the N×N sampling method to compensate the low-frequency components of the power spectrum inversion method are as follows:

[0031] When the compensation level is q, the low-frequency sampling area is centered at the origin (NΔk x )×(NΔk y ), and then divide the low-frequency sampling area into q×q blocks, (q 2 -N 2) / q 2 The outermost area is the first sampling area, and the middle N 2 / q 2 The area is the next sampling area (1≤N≤q), and the above steps are repeated until the set harmonic order p is reached. (The best compensation effect can be achieved when N=(q+1) / 2);

[0032] The specific steps of subharmonic compensation technology of N×N sampling method are as follows:

[0033] The low-frequency harmonics of the power spectrum can be expressed as:

[0034] F sub (k x ,k y )=(q / N) -2p F(k lx ,k ly ) (6)

[0035] Where P is the harmonic order, q is the compensation level; N is the root mean square of the number of compensation areas in the next compensation, (q / N) -2p is the low frequency harmonic weight, k lx =(q / N) -2p k x ,k ly =(q / N) -2p k y They are low-frequency harmonic frequency, discretized, and improved compensated low-frequency phase screen υ L The expression for (m,n) is:

[0036]

[0037] Where m, m', n, n' are integers, M x and M y are the number of sampling points in the x and y directions respectively, and h(m',n') is a complex Gaussian random number with zero mean and unit variance.

[0038] Verifying the accuracy of ocean turbulence phase screen using phase structure function

[0039] For the discrete phase screen generated by numerical simulation, its phase structure function is obtained by calculating the ensemble average of the phase difference between two points, which is defined as:

[0040]

[0041] where <·> is the ensemble average, is the phase screen matrix, is the phase magnitude at point r, r and r' are two points on the phase screen. In the plane rectangular coordinate system, r can be decomposed into The phase structure function of the theoretical plane wave phase can be expressed as:

[0042]

[0043] Where η is the intrinsic scale factor of the beam. r0 is the coherence radius of the plane wave. The coherence radius of the plane wave here is:

[0044]

[0045] Among them, ε, χ T , ω are ocean turbulence parameters. Under the Rytov approximation, the Nikishov refractive index power spectrum is used. The phase structure function of the phase screen of ocean turbulence is still within the inertial region and conforms to the 3 / 5 law.

[0046] The beneficial effects of the present invention are:

[0047] 1. The present invention adopts the subharmonic method of N×N sampling mode to compensate for the defect of insufficient low-frequency components of the power spectrum inversion method;

[0048] 2. Its characteristics are that it can improve the compensation effect not only by increasing the harmonic order, but also by changing the compensation level and the next sampling area, and improve the operation efficiency by selecting the appropriate compensation level and the next sampling area according to the accuracy requirements;

[0049] 3. When q=9, N=5, the compensation effect error approaches 0.03%. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 It is a schematic diagram of the relationship between the compensation level q=3 and the root mean square number of the next compensation area N=2 sampling points;

[0051] Figure 2 It is a schematic diagram of the relationship between the compensation level q=5 and the root mean square number of the next compensation area N=3 sampling points;

[0052] Figure 3 It is a schematic diagram of the relationship between the compensation level q=7 and the root mean square number of the next compensation area N=4 sampling points;

[0053] Figure 4 It is a schematic diagram of a method for sampling harmonic compensation with a compensation level of q=3 and a root mean square number of the next compensation area of ​​N=2;

[0054] Figure 5 is the FFT ocean turbulence phase screen;

[0055] Figure 6 It is an ocean turbulence phase screen with compensation order q = 3 and the root mean square number of compensation areas N = 2;

[0056] Figure 7 It is an ocean turbulence phase screen with compensation order q = 5 and the root mean square number of compensation areas N = 3;

[0057] Figure 8 It is an ocean turbulence phase screen with compensation order q = 7 and the root mean square number of compensation areas N = 4 in the next compensation;

[0058] Fig. 9 It is a comparison diagram of compensation effects of different compensation levels and different RMS N of the number of next compensation areas when the number of compensation times p=10;

[0059] Fig.10 is the structure function D(r) of the ocean turbulence phase screen with different harmonic orders of q / N=3 / 2;

[0060] Fig.11 is the structure function D(r) of the ocean turbulence phase screen with different harmonic orders of q / N=5 / 3;

[0061] Fig.12 is the structure function D(r) of the ocean turbulence phase screen with different harmonic orders of q / N=7 / 4;

[0062] Fig.13 It is the phase screen structure function D(r) of ocean turbulence with different harmonic orders of q / N=9 / 5. DETAILED DESCRIPTION

[0063] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0064] Example 1

[0065] The simulation parameters are set as follows: the phase screen is a square with a grid number of M = 500 and a length of L = 1 m, and the turbulence parameters are set to η = 0.0001, ε = 0.001 m 2 / s 3 ,ω=-3,χ T =10 -9 K 2 / S, wavelength λ=628nm, transmission distance Δz=30m, compensation order q=3, root mean square number of next compensation areas N=2, harmonic order p=5, 10, 15, 25, 30, and the phase screen structure function is the ensemble average of 1000 times.

[0066] The refractive index fluctuation spectrum of seawater is:

[0067]

[0068] Where qε is the kinetic energy dissipation rate per unit seawater, and its value range from the depth to the surface of the seawater is [10 -10 ,10-1 ]m 2 / s 3 . χ T is the mean square temperature dissipation rate of seawater, and its value range from the deep seawater to the seawater surface is [10 -10 ,10 -4 ]K 2 / s. w is defined as the ratio of ocean turbulence caused by temperature to ocean turbulence caused by salinity, and the value range is [-5,0], where -5 means that ocean turbulence is completely caused by temperature, and 0 means that ocean turbulence is completely caused by salinity. η is the Kolmogorov microscale (the inner scale of ocean turbulence). Its parameters are measured through a large number of experiments and are respectively A T =1.863×10 -2 ,A s =1.9×10 -4 ,A Ts =9.41×10 -3 .

[0069] Step 2 is implemented according to the following specific steps:

[0070] Step 2.1, based on the seawater refractive index fluctuation spectrum, calculate that the seawater phase spectrum on any slice perpendicular to the beam propagation direction (z-axis direction) should be:

[0071] F φ (k x ,k y )=2πk 2 Δzφ(k x ,k y )(2)

[0072] Where Δz is the thickness of the turbulent layer or the propagation distance of the light beam, and k = 2π / λ is the beam wave.

[0073] Step 2.2, use F φ (k x ,k y ) for the Gaussian random complex matrix h(k x ,k y ) is filtered and then inverse Fourier transform is performed to obtain the spatial distribution of the phase screen:

[0074]

[0075] Where C is the constant factor controlling the phase screen variance (Δk x Δk y ) 1 / 2 In the discretized space domain, x = mΔx, y = nΔy, Δx, Δy are the sampling intervals in the space domain. For convenience, Δx = Δy, m and n are integers; in the wave number domain, kx =m'Δk x ,k y =m'Δk y ; Δk x , Δk y The sampling interval M is the number of grids in the phase screen.

[0076] Step 3 is implemented according to the following specific steps:

[0077] Formula (3) can be written as:

[0078]

[0079] make The expression of the phase screen is:

[0080]

[0081] Step 4 is implemented according to the following specific steps:

[0082] The specific steps of using the subharmonic method of N×N sampling mode to compensate for the insufficient low-frequency components of the power spectrum inversion method are as follows:

[0083] In order to fully sample the low-frequency part of the power spectrum, when the compensation level is 2, the low-frequency compensation area is (2Δk x )×(2Δk y ), at this time, the low-frequency compensation area overlaps with the high-frequency phase screen sampling point, and then (2Δk x )×(2Δk y ) area is divided into 3×3 sub-areas, the outermost 5 / 9 sub-areas are the first sampling sub-areas, and the remaining 4 / 9 sub-areas are used for the next level of harmonic compensation; finally, the above steps are repeated until the set harmonic order p is reached.

[0084] The specific steps of using the N×N sampling mode subharmonic method to compensate for the insufficient low-frequency components of the power spectrum inversion method are as follows:

[0085] The construction of the phase screen by the power spectrum inversion method causes the loss of low-frequency components and the superposition of low-frequency harmonics. The low-frequency harmonics of the power spectrum can be expressed as:

[0086] F sub (k x ,k y )=(3 / 2) -2p F(k lx ,k ly )(6)

[0087] Where p is the low-frequency harmonic series, (3 / 2) -2p is the low frequency harmonic weight, klx =(3 / 2) -2p k x ,k ly =(3 / 2) -2p k y are the low-frequency harmonic frequencies, and after discretization, the improved compensation low-frequency harmonic phase screen φ L The expression for (m,n) is:

[0088]

[0089] Where m, m', n, n' are integers, M x and M y are the number of sampling points in the x and y directions respectively, and h(m',n') is a complex Gaussian random number with zero mean and unit variance.

[0090] Step 5: Superimpose the ocean turbulence random phase screen generated by the power spectrum inversion method and the low-frequency component random phase screen compensated by the subharmonic method of the N×N sampling method to obtain the final ocean turbulence random phase screen. The accuracy of the ocean turbulence phase screen is verified by the phase structure function.

[0091] For the discrete phase screen generated by numerical simulation, its phase structure function is obtained by calculating the ensemble average of the phase difference between two points, which is defined as:

[0092]

[0093] where <·> is the ensemble average, is the phase screen matrix, is the phase magnitude at point r, r and r' are two points on the phase screen. In the plane rectangular coordinate system, r can be decomposed into The phase structure function of the theoretical plane wave phase can be expressed as:

[0094]

[0095] Where η is the intrinsic scale factor of the beam. r0 is the coherence radius of the plane wave. The coherence radius of the plane wave here is:

[0096]

[0097] Among them, ε, χ T , ω are ocean turbulence parameters. Under the Rytov approximation, the Nikishov refractive index power spectrum is used. The phase structure function of the phase screen of ocean turbulence is still within the inertial region and conforms to the 3 / 5 law.

[0098] Example 2

[0099] The simulation parameters are set as follows: the phase screen is a square with a grid number of M = 500 and a length of L = 1 m, and the turbulence parameters are set to η = 0.0001, ε = 0.001 m 2 / s 3 ,ω=-3,χ T =10 -9 K 2 / S, wavelength λ = 628nm, transmission distance Δz = 30m, compensation order q = 5, the root mean square of the number of subharmonic compensation areas of the next order N = 3, the number of compensation times p = 5, 10, 15, 20, 25, the phase screen structure function is the ensemble average of 1000 times. The remaining steps are the same as those in Example 1.

[0100] Example 3

[0101] The simulation parameters are set as follows: the phase screen is a square with a grid number of M = 500 and a length of L = 1 m, and the turbulence parameters are set to η = 0.0001, ε = 0.001 m 2 / s 3 ,ω=-3,χ T =10 -9 K 2 / S, wavelength λ=628nm, transmission distance Δz=30m, compensation order q=7, the root mean square number of the next compensation area N=4, the harmonic order p=5, 10, 13, 15, 16, and the phase screen structure function is the ensemble average of 1000 times.

[0102] The remaining steps are the same as those in Example 1.

[0103] Example 4

[0104] The simulation parameters are set as follows: the phase screen is a square with a grid number of M = 500 and a length of L = 1 m, and the turbulence parameters are set to η = 0.0001, ε = 0.001 m 2 / s 3 ,ω=-3,χ T =10 -9 K 2 / S, wavelength λ = 628nm, transmission distance Δz = 30m, compensation order q = 9, the root mean square number of the next compensation area N = 5, harmonic order p = 5, 10, 12, 13, 14, the phase screen structure function is an ensemble average of 1000 times. The remaining steps are the same as those in Example 1.

[0105] like Figure 1-4 Schematic diagram of the sampling method structure under different compensation levels.

[0106] Figure 5This is an uncompensated phase screen. In the figure, the phase fluctuates violently at close distances, indicating that the high-frequency components are abundant, but the phases of two points at longer distances are almost unrelated, indicating that the low-frequency components are insufficient. Figure 6 The number of compensation levels is q = 3, and the number of compensation areas for the next time is N × N = 2 × 2; Figure 7 The number of compensation levels is q = 5, and the number of compensation areas for the next time is N × N = 3 × 3; Figure 8 The number of compensation levels is q = 7, and the number of compensation areas for the next time is N × N = 4 × 4. Figure 6 , 7 8, it can be seen that with the reasonable increase of the number of compensation levels and the number of the next compensation areas, the phase correlation between the two points at a longer distance increases, indicating that with the reasonable increase of the number of compensation levels and the number of the next compensation areas, the low-frequency component also increases.

[0107] Fig. 9 The comparison diagram of the phase structure function and theoretical value of the ocean turbulence phase screen in different next compensation areas under different compensation levels for the harmonic order p = 10. It can be seen that when the harmonic order and the compensation level are fixed, the compensation effect first increases and then decreases. When the root mean square number of the next compensation area N = (q+1) / 2, the low-frequency compensation effect can be achieved at the best.

[0108] Fig.10 , 11 , 12, 13 When q / N=3 / 2 and the harmonic order p reaches 30, the compensation effect basically stops increasing. When q / N=5 / 3, 7 / 4, and 9 / 5, the harmonic order p reaches the critical point of compensation at 25, 16, and 14 respectively. If the harmonic order is increased further, it will be over-compensated. When q / N=5 / 3, p=25, the accuracy can reach 1%; when q / N=7 / 4, p=16, the accuracy can reach 0.08%; when q / a=9 / 5, p=14, the accuracy can reach 0.03%.

Claims

1. A subharmonic ocean turbulence phase screen simulation method using an N×N sampling method, characterized in that: The specific steps include: Step 1: Generate a Gaussian random complex matrix with a mean of 0 and a variance of 1; Step 2: Use the ocean turbulence phase disturbance power spectrum function to filter the random complex matrix, and then perform inverse Fourier transform on it to obtain the random phase screen; specifically, follow the following specific steps: Step 2.1, based on the power spectrum function of the ocean turbulence phase disturbance, calculate the seawater phase spectrum on any slice perpendicular to the beam propagation direction at this time: F φ (k x ,k y )=2πk 2 Δzφ(k x ,k y ) (2) Where Δz is the thickness of the turbulent layer or the propagation distance of the light beam, k = 2π / λ is the beam wave, λ is the wavelength of the light beam; φ(k x ,k y ) is the power spectrum function; Step 2.2, use F φ (k x ,k y ) for the Gaussian random complex matrix h(k x ,k y ) is filtered, and then the inverse Fourier transform is performed to obtain the spatial distribution of the phase screen: Where C is the constant factor controlling the phase screen variance (Δk x Δk y ) 1 / 2 , in the discretized space domain, x = mΔx, y = nΔy, Δx, Δy are the sampling intervals of the space domain, let Δx = Δy, m, n are integers; In the wave number domain k x =m'Δk x ,k y =n'Δk y , m', n' are integers; Δk x , Δk y is the sampling interval in the wavenumber domain, M is the number of grids of the phase screen; i is an imaginary number; Step 3: Compensate the low-frequency part of the random phase screen generated in step 2 by the subharmonic method of N×N sampling mode to obtain the compensated low-frequency subphase screen; the specific steps of compensating the low-frequency component of the power spectrum inversion method by the subharmonic method of N×N sampling mode are as follows: When the compensation level is q, the low-frequency sampling area is centered at the origin (NΔk x )×(NΔk y ), and then divide the low-frequency sampling area into q×q blocks, the outermost (q 2 -N 2 ) / q 2 The area is the first sampling area, the middle N 2 / q 2 The area is the next sampling area (1≤N≤q), and the above steps are repeated until the set harmonic order p is reached. The best compensation effect is achieved when N=(q+1) / 2; The specific steps of the subharmonic method of N×N sampling method are: The low-frequency harmonics of the power spectrum are expressed as: F sub (k x ,k y )=(q / N) -2p F(k lx ,k ly ) (6) Where p is the harmonic order, q is the compensation level; N is the RMS number of the next harmonic compensation area, (q / N) -2p is the low frequency harmonic weight, k lx =(q / N) -2p k x ,k ly =(q / N) -2p k y are the low-frequency harmonic frequency, discretized, and the improved compensated low-frequency phase screen φ L The expression for (m,n) is: Where m, m', n, n' are integers, M x and M y are the number of sampling points in the x and y directions respectively, and h(m',n') is a complex Gaussian random number with zero mean and unit variance; Step 4: The random phase screen obtained in step 2 and the compensated low-frequency phase screen generated in step 3 are superimposed to obtain the final phase screen.

2. The subharmonic ocean turbulence phase screen simulation method of N×N sampling mode according to claim 1 is characterized in that: The power spectrum function in step 2 is: Where ε is the kinetic energy dissipation rate per unit seawater, and its value range from the depth to the surface of the seawater is [10 -10 ,10 -1 ]m 2 / s 3 ; χ T is the mean square temperature dissipation rate of seawater, and its value range from the deep seawater to the seawater surface is [10 -10 ,10 -4 ]K 2 / s; w is defined as the ratio of ocean turbulence caused by temperature to that caused by salinity, and its value range is [-5,0], where -5 means that the ocean turbulence is completely caused by temperature, and 0 means that the ocean turbulence is completely caused by salinity; η is the internal scale of ocean turbulence, A T =1.863×10 -2 ,A s =1.9×10 -4 , A Ts ＝9.41×10 -3 。 3. The subharmonic ocean turbulence phase screen simulation method of N×N sampling mode according to claim 1 is characterized in that: Formula (3) can be written as: make The expression of the phase screen is obtained as follows:

Images