A Fast Analysis Method for Blue-Green Laser Transmission through the Sea Surface

The fast estimation of blue-green laser transmission characteristics through the fluctuating sea surface, such as refraction angles and transmittance, is very important to correct operating parameters, detection depth and anti-detection warning in airborne Light Detection and Ranging (LiDAR) applications. However, the geometry of the sea surface is changed by complex environment factors, such as wind and wave, which significantly affect the rapid acquisition of the blue-green laser transmission characteristics. To address this problem, a fast analysis method is provided to rapidly compute the blue-green laser transmittance and refraction angles through the fluctuating sea surface driven by different wind directions and speeds. In the method, a three-dimensional wave model driven by the wind was built to describe the wave spatial distribution varying with time. Using the wave model, the propagation path of the scanning laser footprint was analyzed using the proposed meshing method, thus the transmittance and refraction angles of the optical path can be fast obtained by using parallel computing. The simulation results imply that the proposed method can reduce the time consumption by 70% compared with the traditional analytical method with sequential computing. This paper provides some statistical laws of refraction angles and transmittance through the fluctuating sea surface under different wind conditions, which may serve as a basic for fast computation of airborne LiDAR transmission characteristics in complex environments.


Introduction
Airborne Light Detection And Ranging (LiDAR) is an effective active detection technique for underwater targets [1]. The airborne LiDAR emits high-power blue-green laser pulses to the underwater target then collects the laser pulses scattered by the target. The position, geometric structures, and physical properties of the target can be calculated by processing the collected pulses [2]. The transmission characteristics of the blue-green laser through the sea surface are the most important factors of detection range and accuracy, which attenuate the energy and change the direction of the blue-green laser pulse. The geometry of the sea surface is changed by the effects of complex environmental factors, such as the fluctuations caused by the dynamic wind. As a result, part of the laser energy will be lost and the propagation path of the laser will deviate during the transmission through the sea surface, significantly affecting the detection range and accuracy. Therefore, it is Sensors 2020, 20, 1758 3 of 21

Modeling Three-Dimensional Dynamic Wave
Due to the dynamic influence of sea wind, the wave inclination of the undulating and rough sea surface changes continuously, which directly affects the transmission process of laser beam on the sea surface, and causes the complicated variation of transmission energy on the interface. In order to simulate the transmission process of the laser beam on the rough interface, the height and normal of the wave surface must be obtained. In this section, the wave spectrum theory and linear superposition method were used to build a three-dimensional dynamic wave model. Next, we divided the simulated wave surface into meshes, then obtained the heights and normal of every mesh nodes by calculating the wave model.

Wave Model
In oceanography, wave spectrum is the main way to describe the wave motion, a complex random process. The wave spectrum is obtained by using the spectrum analysis method to analyze long-term observed wave elements such as wave height and period. It describes the frequency distribution of wave energy as a random process relative to its constituent waves.
Multiple kinds of wave spectrum have been proposed. The one we use to simulate wind-driven dynamic wave is the Pierson-Moscowitz (PM) spectrum, which is widely used because of its better data base and higher accuracy. The spectrum function is denoted as S PM (ω) and described by [11]: where ω is the angular frequency of the wave, g is the gravitational acceleration, u is the wind speed at 19.8 meters above the sea, and the spectral peak frequency is denoted as ω m and described as ω m = 8.565/u (2) The actual wave surface is three-dimensional, which means that its energy is not only distributed in a certain frequency range, but also in a certain direction range. However, the PM spectrum is a kind of one-dimensional spectrum that can only describe the distribution of energy with frequency. Therefore, a directional spreading function D(θ) was defined to more accurately describe the energy distribution of three-dimensional waves, which is denoted as S(ω, θ) and shown as follows.
S(ω, θ) = S PM (ω)·D(θ) (3) In this paper, the directional spreading function proposed by the International Towing Tank Conference (ITTC) [12] is used and shown as: where θ denotes the angle between the wave and the wind direction above the sea, which is defined as the direction angle of the wave. Substituting Equations (1) and (4) into Equation (3), the complete energy distribution spectrum can be obtained as: S(ω, θ) = 8.1 × 10 −3 g 2 ω 5 exp −0.74 g uω 4 · 2 π cos 2 (θ), |θ| ≤ π According to Equation (5), energy distribution of waves at different wind speeds were drawn and are shown in Figure 1. It can be seen that most of the energy of waves is concentrated in a Sensors 2020, 20, 1758 4 of 21 certain frequency range, and the range moves to the low frequency part with the increase of the wind speed. Meanwhile, in terms of direction angle, most of the energy of waves is approximately concentrated in the angle range [−90 • , 90 • ]. According to the linear wave theory, the wave can be regarded as a stationary random process composed of multiple cosine waves with different periods, frequencies, and directions. Therefore, we only needed to apply the linear superposition method in the range of wave energy concentration, which can effectively represent all waves to quickly simulate the three-dimensional wave.
Sensors 2020, 20, x FOR PEER REVIEW 4 of 21 and directions. Therefore, we only needed to apply the linear superposition method in the range of wave energy concentration, which can effectively represent all waves to quickly simulate the threedimensional wave. The geometry of the sea surface is described in a right-hand three-dimensional coordinate frame. The − plane denotes the horizontal datum sea level, and the axis is defined as the zenith direction in the frame. The wave superposition model which defines the geometry of the threedimensional sea surface can be described as follows [13]: where = ( , , ) is the wave height of point ( , ) at time ; is the division number of the angular frequency range; is the division number of the direction angle range. For a component The geometry of the sea surface is described in a right-hand three-dimensional coordinate frame. The x − y plane denotes the horizontal datum sea level, and the z axis is defined as the zenith direction in the frame. The wave superposition model which defines the geometry of the three-dimensional sea surface can be described as follows [13]: where z = H(x, y, t) is the wave height of point (x, y) at time t; m is the division number of the angular frequency range; n is the division number of the direction angle range. For a component wave, a ij is its amplitude, ω i is its angular frequency, k i is the corresponding wave number, θ j is its direction angle, and ε ij is the random initial phase (0 ≤ ε ij ≤ 2π). When we get the parameters of the component waves, we can obtain the shape of the wave under different wind conditions. For a specific wind condition, the spectral peak frequency ω m is firstly obtained according to Equation (2), then we can select an angular frequency range [ω l , ω h ], ω l < ω m < ω h correspondingly, and equally divide the frequency range to get the angular frequency of each component wave: where ω i denotes the angular frequency of a component wave and ∆ω is the division increment of angular frequencies which satisfies m·∆ω = ω h − ω l . Meanwhile, wave number k i and angular frequency ω i satisfy the following relation [14]: Although the theoretical variation range of the direction angle is [−π, π], the actual wave energy is mostly distributed within ± π 2 on both sides of the wind direction. Representing the wind direction by U, we can equally divide the range U − π 2 , U + π 2 to get the direction angle of each component wave: where θ j denotes the direction angle of a component wave and ∆θ is the division increment of direction angles which satisfies n·∆θ = π. According to Equation (5), a ij can be obtained as [13]: Finally, the heights of every point on the wave surface can be calculated by combining Equations (6-10), so we can obtain the shape of the wave at any time under different wind conditions.
After investigation and comparison with typical actual wind conditions given by Beaufort wind force scale, we determined the ranges of angular frequency, the division increments and division numbers correspondingly, as shown in Table 1.
where → N(x 0 , y 0 , z 0 ) denotes the normal vector; For the sea surface S : z = H(x, y, t) defined by Equation (6), the following equation can be derived based on Equation (11): Equations (12) and (13) shows that it is unavoidable to superimpose waves two times to calculate the two partial derivatives in order to obtain → N(x 0 , y 0 , z 0 ). Moreover, the calculation of the two partial derivatives requires 2 × m × n additional multiplication of k i and sin θ j (or cos θ j ). Generally, the extra computation consumption is very considerable especially when m and n are relatively large. To avoid this problem, we obtained the normal vector of a point on the wave surface by calculating the sum of the normal of the adjacent meshes rather than by using the analytical computation given by Equations (11)- (13).
The x − y plane region of a certain size is divided into M × N uniform regular mesh by the step size q. For each mesh node x i , y j , i = 1, 2, . . . , M; j = 1, 2, . . . , N, the wave height z i,j = H x i , y j , t can be calculated according to Equation (6). In this way, an approximate sea surface with M × N mesh nodes is determined. The larger the M and N are, the denser the mesh is, and the smoother the surface is. As Figure 2 shows, 1000 × 1000 mesh of the waves at typical wind speeds is drawn in the x − y − z three-dimensional coordinate system based on the parameters shown in Table 1. The − plane region of a certain size is divided into × uniform regular mesh by the step size . For each mesh node , , = 1,2, … , ; = 1,2, … , , the wave height , = ( , , ) can be calculated according to Equation (6). In this way, an approximate sea surface with × mesh nodes is determined. The larger the and are, the denser the mesh is, and the smoother the surface is. As Figure 2 shows, 1000 × 1000 mesh of the waves at typical wind speeds is drawn in the − − three-dimensional coordinate system based on the parameters shown in Table 1. Each mesh node x i , y j , z i,j on the wave surface is adjacent to four mesh nodes: x i−1 , y j , z i−1, j , x i , y j+1 , z i,j+1 , x i+1 , y j , z i+1, j , x i , y j−1 , z i,j−1 , and x i , y j , z i,j is the intersection of the four surrounding mesh elements, each of which has an independent normal vector, as shown in Figure 3. Each mesh node ( , , , ) on the wave surface is adjacent to four mesh nodes: ( −1 , , −1, ), ( , +1 , , +1 ), ( +1 , , +1, ), ( , −1 , , −1 ), and ( , , , ) is the intersection of the four surrounding mesh elements, each of which has an independent normal vector, as shown in Figure 3.
The normal vectors of the four surrounding mesh elements are cross-products of In the above calculation process, only one wave superposition is needed for each mesh node to get the wave heights, and the normal vector can be calculated by using the height values of adjacent mesh nodes, which avoids a lot of additional computation. The denser the mesh is, the closer the normal vectors obtained by this method is to those obtained by the analytical computation.

Incidence Angle Modeling
For a typical LiDAR detection scene, an aircraft is generally used as the air platform to carry the blue-green laser emitter and detector. The transmission of the blue-green laser on the sea surface includes two parts: the down channel where laser pulse signals are transmitted from the air to the underwater, and the up channel where the reflected pulse echoes are transmitted in the opposite direction.
During the continuous flight of the aircraft, the emitter sent the laser pulses along the down channel in the form of circular scanning. Due to the existence of inherent divergence angle of the laser, a pulse was projected on the sea surface to form a spot in the form of laser beam. Therefore, it was necessary for us to focus on the spot to analyze and calculate the transmittance. After the laser beam was transmitted through the sea surface, it hits the underwater target surface by which the reflected pulse echo was formed. Next, the reflected pulse echo was received by the detector along the up channel through the surface. Figure 4 shows the schematic diagram of a typical scene of LiDAR scanning detection.
As shown in Figure 4, the motion and scanning of airborne LiDAR over the sea surface were quantitatively described in a right-hand three-dimensional coordinate frame which was consistent with the previous one used to describe the three-dimensional sea surface. The x − y plane denotes the horizontal datum sea level and the vertical z axis faced upward; O indicates the initial horizontal position of the aircraft and h represents the fixed flight altitude of the aircraft. In Figure 4, the scanning trajectory was circular, while the actual scanning trajectory on the sea surface was helix due to the move of the aircraft; δ denotes the laser scanning angle, which is the angle between the laser emission direction and the z axis and φ is the laser beam divergence angle. The laser pulses were emitted every ∆t at the frequency of f emit (∆t = 1/ f emit ), and each laser pulse was projected on the sea surface to form a spot. Because φ is generally very small, the pulse spot could be approximately regarded as a solid circle, and its diameter d could be approximately obtained by the following formula: Sensors 2020, 20, 1758

of 21
Sensors 2020, 20, x FOR PEER REVIEW 8 of 21 channel through the surface. Figure 4 shows the schematic diagram of a typical scene of LiDAR scanning detection. As shown in Figure 4, the motion and scanning of airborne LiDAR over the sea surface were quantitatively described in a right-hand three-dimensional coordinate frame which was consistent with the previous one used to describe the three-dimensional sea surface. The − plane denotes the horizontal datum sea level and the vertical axis faced upward; indicates the initial horizontal position of the aircraft and ℎ represents the fixed flight altitude of the aircraft. In Figure  4, the scanning trajectory was circular, while the actual scanning trajectory on the sea surface was helix due to the move of the aircraft; denotes the laser scanning angle, which is the angle between the laser emission direction and the axis and is the laser beam divergence angle. The laser pulses were emitted every Δ at the frequency of (Δ = 1/ ), and each laser pulse was projected on the sea surface to form a spot. Because is generally very small, the pulse spot could be approximately regarded as a solid circle, and its diameter could be approximately obtained by the following formula: The radius of the scanning trajectory can be obtained as: Figure 5 shows the laser scanning process from the top view. For the convenience of analysis, the plane was set to fly in the positive direction of the axis with a flight speed of ; The LiDAR was scanning in the counterclockwise direction with a scanning frequency of in revolutions per second, which indicates the number of scanning trajectory circles per second. The radius of the scanning trajectory R can be obtained as:  In Figure 5, the coordinate of a pulse spot center is used to represent the pulse spot position; Δ denotes the scanning azimuth step between the next pulse spot and the previous pulse spot can be obtained as: We set the aircraft coordinate at the initial time 0 as ( 0 , 0 , ℎ) , and let the pulse spot coordinate be ( 0 , 0 , 0 ) correspondingly. There was an obvious relationship as follows. In Figure 5, the coordinate of a pulse spot center is used to represent the pulse spot position; ∆γ denotes the scanning azimuth step between the next pulse spot and the previous pulse spot can be obtained as: Sensors 2020, 20, 1758 10 of 21 We set the aircraft coordinate at the initial time t 0 as X c t0 , Y c t0 , h , and let the pulse spot coordinate be (x t0 , y t0 , z t0 ) correspondingly. There was an obvious relationship as follows.
We let the successive aircraft coordinate at the time t i (t i = t 0 + i·∆t) be X c ti , Y c ti , h , and let the pulse spot coordinate be (x ti , y ti , z ti ) correspondingly. Because the aircraft was set to fly in the positive direction of the x axis, the following relationship can be easily obtained as: According to the angle relationship, considering the displacement caused by the aircraft flight, (x ti , y ti , z ti ) can be derived as follows: Substituting Equations (18) and (19) into Equation (21), the following expression can be obtained as: After obtaining the coordinate of the pulse spot, we focused on the local wave mesh at the pulse spot position. As Figure 6 shows, multiple mesh nodes are covered by the pulse spot at time t i . Next, we began to analyze the transmission characteristics of the mesh nodes inside the pulse spot. The mesh nodes inside the pulse spot at time form a point set and shown as: where , , , indicates the coordinates of these internal nodes; is the spot diameter given by Equation (16); and is the step size of the wave mesh mentioned in Section 2.1.2. Figure 7 shows the incident angles of the laser light in a laser beam. Naturally, for each laser light line in the laser beam, the falling point of them were considered different. However, only those  The mesh nodes inside the pulse spot at time t i form a point set I ti and shown as: where x p , y q , z p,q indicates the coordinates of these internal nodes; d is the spot diameter given by Equation (16); and q is the step size of the wave mesh mentioned in Section 2.1.2. Figure 7 shows the incident angles of the laser light in a laser beam. Naturally, for each laser light line in the laser beam, the falling point of them were considered different. However, only those light lines falling on the mesh nodes were considered.
where ( , , , ) indicates the coordinates of these internal nodes; is the spot diameter given by Equation (16); and is the step size of the wave mesh mentioned in Section 2.1.2. Figure 7 shows the incident angles of the laser light in a laser beam. Naturally, for each laser light line in the laser beam, the falling point of them were considered different. However, only those light lines falling on the mesh nodes were considered. ⃗ ⃗ denotes the vector of the laser light line falling on a node; ⃗⃗ is the normal vector of the node; and is the incident angle of the laser light line.
Combining with the normal vector of the mesh node given by Equation (15) For every mesh node in point set I ti , the vector of the light line falling on the node can be obtained as: Combining with the normal vector of the mesh node given by Equation (15), the incident angle of the light line falling on x p , y q , z p,q can be obtained as:

Transmittance Modeling
Because the wavelength of the blue-green laser (470-580 nm) was far less than the millimeter wave band where the wave characteristic wavelength was located [15], the geometrical optics theory based on the Fresnel formula was applicable to the analysis of the transmission characteristics on the interface.
The refraction occurred when the laser was transmitted through the sea surface twice along the down channel and the up channel respectively. The wave motion may change the light path in this short time between the two transmissions. For the down channel, the transmission medium of the refracted light was sea water. We let the length of the underwater transmission path between the incident point and the target be L. The physical definition of the optical path: where L r denotes the optical path of the refracted light in sea water and n SA is the refractive index of the sea water relative to the air (n SA = 1.33). At present, the depth of the underwater target is within 100m. We set L = 100 m and the light speed in vacuum c = 3 × 10 8 m/s. The time difference ∆t between the two transmissions through the sea surface can be calculated as: The phase difference ∆ϕ of the wave between the two transmissions can be obtained as: where T s is the average period of the wave motion, generally at 1-9 s according to the observation and statistics [16][17][18][19]. Substituting the range of T s into Equation (28), the corresponding range of ∆ϕ can be obtained as: It can be seen that ∆ϕ was extremely small, so the motion of the wave in this period can be ignored. It is considered that the wave was practically stationary in this time between the two transmissions. As a result, the light paths of the down channel and the up channel were reversed, and the laser light lines along the two channels passed through the same incident point on the rough sea surface, which is shown in Figure 8.
down channel and the up channel respectively. The wave motion may change the light path in this short time between the two transmissions. For the down channel, the transmission medium of the refracted light was sea water. We let the length of the underwater transmission path between the incident point and the target be . The physical definition of the optical path: where denotes the optical path of the refracted light in sea water and is the refractive index of the sea water relative to the air ( = 1.33). At present, the depth of the underwater target is within 100m. We set = 100m and the light speed in vacuum = 3 × 10 8 m/s. The time difference between the two transmissions through the sea surface can be calculated as: The phase difference Δ of the wave between the two transmissions can be obtained as: where is the average period of the wave motion, generally at 1-9 s according to the observation and statistics [16][17][18][19]. Substituting the range of into Equation (28), the corresponding range of Δ can be obtained as: 1.98 × 10 −7 < < 1.78 × 10 −6 (29) It can be seen that Δ was extremely small, so the motion of the wave in this period can be ignored. It is considered that the wave was practically stationary in this time between the two transmissions. As a result, the light paths of the down channel and the up channel were reversed, and the laser light lines along the two channels passed through the same incident point on the rough sea surface, which is shown in Figure 8. In Figure 8, ( , , , ) is a mesh node inside the pulse spot on the interface, as the same incident point of the laser light lines along the down channel and up channel respectively; For the down channel, represents the angle of incidence, and ′ represents the angle of refraction; On the In Figure 8, x p , y q , z p,q is a mesh node inside the pulse spot on the interface, as the same incident point of the laser light lines along the down channel and up channel respectively; For the down channel, α represents the angle of incidence, and α represents the angle of refraction; On the contrary, α represents the angle of refraction and α represents the angle of incidence for the up channel. We can analysis the transmittance based on α which can be obtained by Equation (25).
The Snell Law gives a relationship between α and α as follows: When the laser is transmitted from the first medium to the second medium, it could be decomposed into two mutually perpendicular components, the S wave and P wave. The S wave was the light component perpendicular to the incident plane and the P wave was the one parallel to the incident Sensors 2020, 20, 1758 13 of 21 plane. The amplitude transmission coefficients of the S wave and P wave were given by the Fresnel formula as follows: where t s and t p are the amplitude transmission coefficients of the S wave and P wave respectively; θ x is the incident angle in the first medium, and θ y is the refraction angle in the second medium.
Considering the polarization of the laser light, and according to the Malus Law, the transmittance can be calculated as follows: where ρ is the transmittance of the refracted light energy to the incident light energy; n 1 and n 2 are the refractive index of the first medium and the second medium respectively; and β represents the polarization azimuth of the incident light relative to the incident plane. β represents the polarization azimuth of the refracted light relative to the incident plane. According to Figure 8, for the down channel, the following equation can be obviously obtained: In the case that the reflection characteristics of the target surface was not available, it was assumed that the specular reflection occurred on the surface of the underwater target. As a result, the polarization characteristics of reflected light were not changed. Therefore, for the up channel, the corresponding equation is: Considering Equations (30)-(35) comprehensively, the transmittance of the laser light line passing through the mesh node x p , y q , z p,q along the down channel and the up channel could be calculated respectively. Furthermore, the total transmittance at x p , y q , z p,q can be obtained by multiplying them as follows. ρ x p , y q , z p,q = ρ d x p , y q , z p,q ·ρ u x p , y q , z p,q where ρ x p , y q , z p,q denotes the total transmittance; ρ d x p , y q , z p,q denotes the transmittance of the down light line, and ρ u x p , y q , z p,q denotes the transmittance of the up laser light line.
For the whole pulse spot on the rough sea surface, the average value of the transmittance at all mesh nodes inside the spot was taken as the transmittance of the spot, as follows: where ρ spot denotes the transmittance of the pulse spot and N denotes the number of the mesh nodes inside the pulse spot.

Simulations and Results
Based on the models and simulation methods presented in Section 2, several simulations for the pulse spot transmittance and refraction angles under different wind conditions were performed. Section 3.1 described simulation conditions and parameters. Section 3.2 compared the simulation speeds of the proposed method with parallel computing with those of the traditional analytical method with sequential computing. Section 3.2 provided the statistics of the pulse spot refraction angles and transmittance.

Simulation Conditions
We set the coordinate of the aircraft at initial time as (0, 0, h); The step size of the mesh is set as q = 0.01 m; In the first 100 s, we took a pulse spot every 0.1 s, and computed the transmittance of the 1000 pulse spots in different typical wind speeds and wind directions. The simulation environment is shown in Table 2. We took the working condition parameters of mainstream airborne LiDAR systems as the simulation parameters, as shown in Table 3.

Simulation Time
In the simulation, we first used the traditional differential analytic method to calculate the normal vector of the sea surface, and simulated the transmission of 1000 pulse spots with sequential computing. To illustrate the rapidity of the simulation method, we then used the mesh division method used in this paper to calculate the normal vector, and simulated the transmission of the pulse spots with parallel computing. The simulation time comparison between the two methods is shown in the Table 4. In simulation, typical wind speeds were taken from Beaufort wind force scale to simulate various wind conditions (include strong wind), which may be encountered in practical application.
From Table 4, it can be seen that the proposed meshing method with parallel computing used in this paper can reduce the simulation time consumption by about 70% compared with the traditional analytical method with sequential computing.

Statistics of Refraction Angles
The refraction angles on the mesh points in a pulse spot were averaged as the refraction angle of the pulse spot. The maximum, minimum, and average of the refraction angles of 1000 pulse spots were obtained, as shown in Table 5.
The histogram of refraction angle frequency distribution of the 1000 pulse spots is shown in Figure 9.
From the statistical results of the refraction angles of pulse spots shown in Table 5, we can see that the refraction angles under each wind condition were concentrated between 0.14 rad and 0.38 rad; Furthermore, Figure 9 shows that the distributions of refraction angles under each wind condition were consistent with the Gaussian distribution. For each wind condition, more than 90% of the pulse spot refraction angles were distributed in three bins of [0.1848 rad, 0.3234 rad], and the maximum frequency was in the bin of [0.2310 rad, 0.2772 rad].
From the results, we can see the dynamic range of refraction angle was relatively small. However, the planimetric effects and depth coordinate errors caused by the refraction angle were considerable to detect the underwater targets at tens or even hundreds of meters [20]. Therefore, it is necessary to quickly calculate the refraction angle under different wind conditions to correct the detection results. The histogram of refraction angle frequency distribution of the 1000 pulse spots is shown in Figure 9. From the statistical results of the refraction angles of pulse spots shown in Table 5, we can see that the refraction angles under each wind condition were concentrated between 0.14 rad and 0.38 rad; Furthermore, Figure 9 shows that the distributions of refraction angles under each wind condition were consistent with the Gaussian distribution. For each wind condition, more than 90% of the pulse spot refraction angles were distributed in three bins of [0.1848 rad, 0.3234 rad], and the maximum frequency was in the bin of [0.2310 rad, 0.2772 rad].
From the results, we can see the dynamic range of refraction angle was relatively small. However, the planimetric effects and depth coordinate errors caused by the refraction angle were considerable to detect the underwater targets at tens or even hundreds of meters [20]. Therefore, it is

Statistics of the Transmittance
In the practical application of airborne LiDAR, the minimum value of transmittance should be considered in the detection, which decides the maximum detection depth. The maximum value of transmittance should be considered in the anti-detection, which decides the minimum safe depth. The maximum, minimum and average of the transmittance of 1000 pulse spots were obtained, as shown in Table 6. The scatter diagram of the transmittance of the 1000 pulse spots are shown in Figure 10.
From the statistical results of the pulse spot transmittance shown in Table 6, we can see that the pulse spot transmittance at different positions under each wind condition had an average value of about 95.9%; Furthermore, as Figure 10 shows, it is apparent that the distribution of the pulse spot transmittance under different wind conditions was similar, and the transmittance values were concentrated between 95.7% and 96%.  Table 6, we can see that the pulse spot transmittance at different positions under each wind condition had an average value of about 95.9%; Furthermore, as Figure 10 shows, it is apparent that the distribution of the pulse spot transmittance under different wind conditions was similar, and the transmittance values were concentrated between 95.7% and 96%.

Discussion
In this paper, we analyzed and computed laser transmittance through the rough wind-driven sea surface and we obtained some statistical properties of the transmittance of pulse spots under different wind conditions. As the statistical results shows, the refraction angles of pulse spots were concentrated between a small range, and the transmittance of pulse spots was generally high and stable under different wind conditions. The results may be explained as follows.
Firstly, we can obtain the relationship between the transmittance and the incident angle of the down channel based on the computational process given by Equations (30)-(36), as shown in Figure 11.

Discussion
In this paper, we analyzed and computed laser transmittance through the rough wind-driven sea surface and we obtained some statistical properties of the transmittance of pulse spots under different wind conditions. As the statistical results shows, the refraction angles of pulse spots were concentrated between a small range, and the transmittance of pulse spots was generally high and stable under different wind conditions. The results may be explained as follows.
Firstly, we can obtain the relationship between the transmittance and the incident angle of the down channel based on the computational process given by Equations (30)-(36), as shown in Figure 11.
It can be seen that the transmittance was above 95% and was almost unchanged when the incident angle changed in the range of about 0 •~4 0 • . Meanwhile, when the incident angle increased to 57 • , the transmittance slowly decreases to 90%; When the incident angle increased in the range of 60 •~9 0 • , the transmittance decreased rapidly to 0. According to our calculation results, the incident angles of the pulse spots were concentrated in the range of 7.9412 • to 29.1177 • , which meant that the transmittance and the refraction angles are correspondingly distributed in a small range. Sensors 2020, 20, x FOR PEER REVIEW 18 of 21 Figure 11. The relationship between transmittance and the incident angle of the down channel.
It can be seen that the transmittance was above 95% and was almost unchanged when the incident angle changed in the range of about 0°~40°. Meanwhile, when the incident angle increased to 57°, the transmittance slowly decreases to 90%; When the incident angle increased in the range of 60°~90°, the transmittance decreased rapidly to 0. According to our calculation results, the incident angles of the pulse spots were concentrated in the range of 7.9412° to 29.1177°, which meant that the transmittance and the refraction angles are correspondingly distributed in a small range.
Furthermore, under the set simulation conditions, according to the actual working conditions of the LiDAR detection, the size of the pulse spot on the sea surface was about 26 cm in diameter, which was much smaller than the scale of the wave. There was little difference in the incident angles of the light lines falling on every mesh node inside the pulse spot. Therefore, the transmittance of a pulse spot was quite close to that of the mesh nodes inside it, which depended on the incident angle. In addition, because of the relatively small rotation scanning angle, the incident angle was almost determined by the inclination of the fluctuating sea surface.
First of all, under low wind speeds, the sea surface was gentle where there was basically no fluctuation, so the inclination of the sea surface was relatively small. Furthermore, with the increase of the wind speed above the sea surface, the larger wave energy described based on the PM spectrum was concentrated in a smaller frequency band. Consequently, the fluctuation of the rough sea surface became much larger, which meant that the wave height in the vertical direction and the wave scale in the horizontal direction became larger at the same time, so the inclination of the sea surface remained relatively small. Moreover, because the wave fluctuates periodically throughout the horizontal plane, the overall range of the sea surface inclination hardly changeed with the wind direction which only affected the distribution of the sea surface inclination at different positions. Therefore, in the description of the PM wave spectrum, the wind speed and directions had little effect on the inclination of the sea surface.
As a result, the incident angles of the down channel in the scanning process were small. Therefore, the simulation results were reasonable and consistent with the conclusion of the research of Dong et al. [21]. The effects of low wind speeds from 2 to 5.25 m/s on the refraction angle were analyzed in Reference [8]. From Figure 10 in Reference [8], it can be seen that no correlation was observed between the means of the laser beam center (refraction angles) and wind speeds. The experimental results are coincident with those in our paper.
However, it should be noted that the outcome of this study cannot be taken as evidence that complicated wind conditions will never affect the transmittance and refraction angle on the sea surface, because we did not consider other relevant factors on the dynamic complex sea surface, such as foams, tides, other random fluctuations, and special wave patterns beyond the description range of the PM spectrum. Specifically speaking, foams covered on the sea surface can significantly enhance surface scattering with the increase of wind speed, which has a great influence on the transmittance [9]. The tides and random fluctuations can increase the geometry uncertainty of the sea surface Furthermore, under the set simulation conditions, according to the actual working conditions of the LiDAR detection, the size of the pulse spot on the sea surface was about 26 cm in diameter, which was much smaller than the scale of the wave. There was little difference in the incident angles of the light lines falling on every mesh node inside the pulse spot. Therefore, the transmittance of a pulse spot was quite close to that of the mesh nodes inside it, which depended on the incident angle. In addition, because of the relatively small rotation scanning angle, the incident angle was almost determined by the inclination of the fluctuating sea surface.
First of all, under low wind speeds, the sea surface was gentle where there was basically no fluctuation, so the inclination of the sea surface was relatively small. Furthermore, with the increase of the wind speed above the sea surface, the larger wave energy described based on the PM spectrum was concentrated in a smaller frequency band. Consequently, the fluctuation of the rough sea surface became much larger, which meant that the wave height in the vertical direction and the wave scale in the horizontal direction became larger at the same time, so the inclination of the sea surface remained relatively small. Moreover, because the wave fluctuates periodically throughout the horizontal plane, the overall range of the sea surface inclination hardly changeed with the wind direction which only affected the distribution of the sea surface inclination at different positions. Therefore, in the description of the PM wave spectrum, the wind speed and directions had little effect on the inclination of the sea surface.
As a result, the incident angles of the down channel in the scanning process were small. Therefore, the simulation results were reasonable and consistent with the conclusion of the research of Dong et al. [21]. The effects of low wind speeds from 2 to 5.25 m/s on the refraction angle were analyzed in Reference [8]. From Figure 10 in Reference [8], it can be seen that no correlation was observed between the means of the laser beam center (refraction angles) and wind speeds. The experimental results are coincident with those in our paper.
However, it should be noted that the outcome of this study cannot be taken as evidence that complicated wind conditions will never affect the transmittance and refraction angle on the sea surface, because we did not consider other relevant factors on the dynamic complex sea surface, such as foams, tides, other random fluctuations, and special wave patterns beyond the description range of the PM spectrum. Specifically speaking, foams covered on the sea surface can significantly enhance surface scattering with the increase of wind speed, which has a great influence on the transmittance [9]. The tides and random fluctuations can increase the geometry uncertainty of the sea surface [10,20], which may greatly affect the transmission path of the laser to influence the transmittance and refraction angles. Therefore, there is abundant room for further progress in establishing a more unified and comprehensive model for more complete and complex environmental conditions by introducing other environmental factors. Previous relevant research may serve as a necessary basis for further study [9,10,22,23].
This finding provides some statistical laws of refraction angles and transmittance through the fluctuating sea surface under complex wind conditions, which may serve as a basis for future fast computation of airborne LiDAR transmission in complex environments. Moreover, this preliminary finding has a certain reference significance for the correlative development, verification, and practical application of the LiDAR. From the perspective of detecting, benefiting from the high concentration of the energy inside the pulse spot, the echo strength, and the detection depth of the LiDAR are almost unaffected by the environment on the sea surface. From the perspective of anti-detecting, the windy weather hardly reduces the risk detected by the LiDAR on account of the invariably guaranteed detection effect. With this fast analysis method, the transmission characteristic in the current complex environment can be quickly found out in order to undertake corresponding correction and anti-detection.

Conclusions
This paper proposes a fast analysis method for blue-green laser transmission through a wind-driven dynamic sea surface. We obtained the transmission characteristic under different wind speeds and directions quickly with the meshing method and parallel calculation method. The simulation results show that our method can reduce the time consumption by 70% compared with the traditional method. The statistical results suggest that the refraction angles and transmittance of pulse spots on the fluctuating rough sea surface is generally stable under different wind conditions, which may provide a basis for the future fast computation of airborne LiDAR transmission in complex environments. Meanwhile, this study provided a fresh and useful mean for the fast analysis of the transmission process of the LiDAR detection, which may have a certain reference significance for the development, verification, and practical application of the LiDAR.
However, one possible limitation involves the actual transmission process, which is extremely complicated and, is not affected only by the sea wind, but rather a variety of coupled environmental factors. One direction of the future work is to take other environmental factors, such as foams, tides, and other random fluctuations, into account in order to analyze more comprehensive transmission mechanisms on the complex sea surface. Furthermore, our study only focuses on the simulation of parameters of mainstream airborne LiDAR systems, and gives some reference conclusions. In future research, it might be possible to conduct more investigations for the relationship of more LiDAR system parameters and their random and systematic effects on the detection results.