Infrasonic source altitude localization based on an infrasound ray tracing propagation model

Since the 1990s, according to the requirements of the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) of the United Nations, there has been an increasing demand for monitoring nuclear tests worldwide. Due to its low frequencies and low atmospheric attenuation, infrasound has gradually become one of the significant approaches used to monitor nuclear tests. To this end, a series of studies has been developed, including infrasound signal processing methods and propagation theory algorithms. Based on existing seismic detection techniques, Cansi¹1. Y. Cansi, “ An automated seismic event processing for detection and location: The P.M.C.C. method,” J. Geophys. Res. Lett. 22, 1021–1024, https://doi.org/10.1029/95GL00468 (1995). proposed the progressive multichannel correlation (PMCC) method, which optimized a monitoring array and obtained the azimuth and apparent velocity of an infrasound wave reaching the array in different frequency bands. At present, researchers in CTBTO countries widely use the azimuth measured by two or more arrays and their geographical coordinates to locate the latitude and longitude of an infrasound source by mapping to the Earth coordinate system. However, this method does not make full use of the information of the apparent velocity and ignores the propagation process of the sound wave in the atmosphere, resulting in the inability to obtain the altitude of the infrasound source.

Ray tracing is an important method for studying the process of infrasound propagation in the atmosphere. As early as 1935, Whipple²2. F. J. W. Whipple, Sc. D., F. Inst. P., “ The propagation of sound to great distances,” Q. J. R. Meteorol. Soc. 61, 285–308 (1935). https://doi.org/10.1002/qj.49706126102 conducted a preliminary study on sound rays and found not only that sound wave propagation paths change with different seasons but also that temperature and wind direction are main influencing factors. Groves³3. G. V. Groves, “ Geometrical theory of sound propagation in the atmosphere,” J. Atmos. Terr. Phys. 7, 113–127 (1955). https://doi.org/10.1016/0021-9169(55)90119-7 derived a theory of acoustic rays by an analogy with geometric optics. Assuming that both the sound speed and wind velocity are functions of altitude, he verified that sound waves can reenter the atmosphere. Barry⁴4. G. Barry, “ Ray tracings of acoustic waves in the upper atmosphere,” J. Atmos. Terr. Phys. 25, 621–629 (1963). https://doi.org/10.1016/0021-9169(63)90156-9 combined an atmospheric layering model to study the frequency and amplitude constraints using ray propagation theory. He simulated the propagation trajectories of rays in a layered medium where the air density and temperature were functions of altitude without considering propagation absorption. Based on the Sommerfeld-Runge law and the law of conservation of peaks, Haselgrove and Haselgrove^5,65. C. B. Haselgrove and J. Haselgrove, “ Twisted ray paths in the ionosphere,” J. Proc. Phys. Soc. 75, 357–363 (1960). https://doi.org/10.1088/0370-1328/75/3/3046. J. Haselgrove, “ The Hamiltonian ray path equations,” J. Atmos. Terr. Phys. 25, 397–399 (1963). https://doi.org/10.1016/0021-9169(63)90173-9 derived Hamilton's equations in general curve coordinates and created another solution to the ray equation by combining the equations with the acoustic dispersion equation. Jones and Stephenson⁷7. R. M. Jones and J. J. Stephenson, “ A versatile three-dimensional ray tracing computer program for radio waves in the ionosphere,” Rural Utilities Service Government Printing Office, Washington, D.C., 1975. proposed a numerical integration method for solving the Hamilton's equations in three dimensions, which could be applied not only to the propagation of electromagnetic waves in the ionosphere but also to long-distance propagation of atmospheric infrasound waves.⁸8. R. M. Jones, E. S. Gu, and A. J. Bedard, Jr., “ Infrasonic atmospheric propagation studies using a 3-D ray trace model,” in Proceedings 22nd Conference on Severe Local Storm, 4–8 October 2004. Jones et al.⁹9. R. M. Jones, J. P. Riley, and T. M. Georges, “ A versatile three-dimensional Hamiltonian ray-tracing program for acoustic waves in the atmosphere above irregular terrain,” Wave Propagation Laboratory, Boulder, CO, 1986. proposed an algorithm and a program for the numerical integral Hamilton method applied to the atmosphere named HARPA. Virieux et al.¹⁰10. J. Virieux, N. Garnier, E. Blanc, and J.-X. Dessa, “ Paraxial ray tracing for atmospheric wave propagation,” J. Geophys. Res. Lett. 31, 349–371, https://doi.org/10.1029/2004GL020514 (2004). modified the Hamiltonian and established a ray tracing model in a two-dimensional Cartesian coordinate system, in which the sound pressure amplitude and propagation time estimation were more accurate due to the addition of layered media and a wind velocity correction. Dessa et al.¹¹11. J.-X. Dessa, J. Virieux, and S. Lambotte, “ Infrasound modeling in a spherical heterogeneous atmosphere,” J. Geophys. Res. Lett. 32, 250–258, https://doi.org/10.1029/2005GL022867 (2005). further derived Hamilton's equations in three-dimensional spherical coordinates to establish a ray propagation model under the curvature of the Earth. These studies were limited to the study of ray propagation models for different coordinate systems and layered media. The ray propagation model in three-dimensional spherical coordinates is more accurate than the other models, while a Cartesian coordinate system will require an increase in computational complexity to increase the accuracy due to the large infrasound propagation distance, which even reaches the same order of magnitude as the Earth's radius. In this paper, a ray tracing method in spherical coordinates is used to solve the problem of the PMCC method not being able to locate the altitude of an infrasound source. A back-projection model of the ray propagation path is combined with the azimuth and apparent velocity calculated by the PMCC method to reconstruct the atmospheric propagation path of the infrasound wave and locate the infrasound source, which is of great significance for infrasound monitoring.

Based on a three-dimensional Hamiltonian ray tracing program⁹9. R. M. Jones, J. P. Riley, and T. M. Georges, “ A versatile three-dimensional Hamiltonian ray-tracing program for acoustic waves in the atmosphere above irregular terrain,” Wave Propagation Laboratory, Boulder, CO, 1986. in spherical coordinates, this paper establishes a back-projection reconstruction theoretical model of an infrasound propagation trajectory in the Earth's atmosphere by determining the main factors affecting long-distance infrasound propagation. Adopting the atmospheric empirical models, the mass spectrometer incoherent scatter radar model (MSISE00) and horizontal wind model (HWM93), the effective sound speed and infrasound propagation paths are calculated. Some paths are used to verify the established back-projection model. Finally, according to existing experimental data, ray trajectories are reconstructed, and the altitude of the source is located accurately, which demonstrates the feasibility of the back-projection model.

Fermat's principle states that the path taken between two points by a ray is the path for which a wave crest traverses the path in the least time (that is, the phase of the wave is a minimum when varying the path). The ray tracing program numerically integrates Hamilton's equations to trace three-dimensional paths of acoustic rays through model atmospheres, which are a differential expression of Fermat's principle. In the high-frequency limit, waves behave like particles and travel along rays according to equations that exactly parallel the equations governing changes in the position and momentum of mechanical systems.⁸8. R. M. Jones, E. S. Gu, and A. J. Bedard, Jr., “ Infrasonic atmospheric propagation studies using a 3-D ray trace model,” in Proceedings 22nd Conference on Severe Local Storm, 4–8 October 2004. A ray bends in the direction where the speed of sound determined by sound speed and wind velocity is low.

In this paper, the azimuth is 0° to the north and positive in the clockwise direction. The value ranges from 0° to 360°. The elevation is defined as the angle with the horizontal ground and ranges from −90° to 90°. Rays with a elevation close to 90° will reach a higher altitude and terminate at the thermosphere eventually due to the sound absorption increasing sharply and exceeding 200 dB, so the elevation range in Fig. 3(b) and Fig. 4(b) is 0°–44°. Affected by the summer westward monsoon, the effective sound speed in Fig. 3(a) has two distinct minimum values, and the propagation trajectories corresponding to the southwest direction in Fig. 3(b) have two long-distance propagation modes called stratospheric and thermospheric returns. However, the effective sound speed in Fig. 3(c) has only one evident minimum value, and the propagation trajectories to the northeast in Fig. 3(b) exhibit only a thermospheric return. Comparing the two propagation directions in Fig. 3(b), the stratospheric monsoon is the main reason for the formation of the stratospheric return. It will exist only when the propagation direction follows the wind and will not appear against the wind. The thermospheric return always exists as the altitude increases. The wind vector in the back-projection model takes the opposite value, so the wind direction turns eastward. Therefore, there are stratospheric and thermospheric returns northeastward and only a thermospheric return southwestward in Fig. 4(b), which corresponds to the conventional model. Figures 3 and 4 illustrate that the downwind direction facilitates the formation of the return and that the windward direction suppresses the formation of the return. After the altitude reaches a certain height, the effective sound velocity will gradually increase, and the thermospheric return will certainly exist.

Figure 5(a) and Fig. 5(b) show that the stratospheric return and the back-projection paths substantially coincide and that the error of the whole process is within 0.1 km. Under the condition of a horizontal propagation distance of 700 km, the latitude and longitude errors of the back-projection paths gradually increases but are still within 0.006° in Figs. 5(c) and 5(d), which are approximately equal to the horizontal distance of 600 m.

Figure 6(a) shows that the thermospheric return and the back-projection paths also substantially coincide. The altitude error fluctuates greatly but is still within 1.5 km in Fig. 6(b). The latitude and longitude errors are both within 0.005° in Figs. 6(c) and 6(d).

All above results show that the back-projection trajectory of a single ray is consistent with the propagation trajectory regardless of whether it is in the stratosphere or thermosphere. For a horizontal propagation distance between 700 to 800 km, the back-projection model has high precision, with an altitude error within 1.5 km and latitude and longitude errors both within 0.008°.

According to NASA, a meteorite exploded over Xilinguole League in Inner Mongolia at 04:13:30 on November 5, 2014, Beijing time. The explosion occurred at 43.1°N, 115.8°E, and the altitude was 22.2 km. The infrasound signal was collected through a three-element array, and the spacings between the elements were 1.39, 1.19, and 0.26 km, as shown in Fig. 7. The average of the element coordinates is taken to obtain the array center coordinates, which are in the southeast direction of the infrasound source.

There are two methods to locate the infrasound source altitude. One is to use signals received by two or more arrays to reconstruct the infrasound ray trajectories, where the intersection or dense point of the ray corresponds to the position of the infrasound source, including latitude, longitude and altitude. The other method relies on knowing any extra information like the latitude and longitude of the infrasound source or the time it takes for sound waves from generation to reception to reconstruct the infrasound ray trajectories. When the horizontal distance of the ray equals the horizontal distance from the infrasound source to the center of the array or when the propagation time of the ray equals the time it takes from generation to reception, the corresponding altitude on the ray is the estimate of infrasound source altitude. Since we only have one array receiving the explosion data, the second method was used.

In November, the northern hemisphere is in winter, and the eastward wind prevails. Southeastward infrasound propagation should have stratospheric returns. However, the infrasound source is 22.2 km above sea level, which is in the position of the minimum effective sound speed of the stratosphere, as shown in Fig. 8(a). As a result, in Fig. 8(b), the infrasound rays will periodically be reflected between 10 and 40 km in the stratosphere. This part of the signal cannot reach the ground without considering special airflow.

The sound speed $c_{\text{true}}$ of the center of the array on that day was approximately 0.337 km/s. Combined with the apparent velocity $c_{\text{apparent}}$, the elevation can be obtained by $\text{arccos}(c_{\text{true}}/c_{\text{apparent}})$. The data acquisition result may contain interference noise. The PMCC result in Fig. 9(b) needs to be filtered according to the actual conditions: the apparent speed of the sound wave needs to be greater than the local sound speed, and because the latitude and longitude of the infrasound source are known, the azimuth can be limited to a certain range. In Fig. 9(a), each array element receives two pulse signals near 04:37:30 and 04:39:20. The pulse signal of 04:37:30 has a larger apparent velocity and higher elevation as seen in Fig. 9(b). Since there is one single array in this explosion, it is necessary to get additional information for source localization. Assuming that the time of the meteorite explosion has been known, the back-projection model is used to reconstruct the ray trajectories for these two pulse signals.

As shown in Fig. 10, the horizontal distance of each localization point on the trajectories of the back projection reconstruction does not exceed 450 km, indicating that this signal does not originate from the final explosion provided by NASA. Considering that the meteorite explosion is a complex process and the signal amplitude at 04:37:30 is less than 04:39:20, it is speculated that the meteorite had a small explosion before the final explosion. This hypothetical explosion was very close to the final explosion. The falling speed of the meteorite could reach 17 km/s, and the time between the two explosions was within 3 s. Although the meteorite made two explosions in the last 3 s, but the meteorite falling speed was much greater than the sound speed, two received signals have a delay of 110 s. Not all of the rays converge at one point because the sound diffusion also shifted the ray during propagation. Most of the localization points are concentrated in a circle with horizontal distance of 433 km and altitude of 35.8 km, with a radius of 4.375 km. The average latitude and longitude of this part is 43.1123°N, 116.2771°E, the horizontal average distance to the array is 433.969 km, and the average altitude is 35.8204 km. So this hypothetical explosion is closer to the array, and the resulting signal is received earlier.

Most of the localization points are distributed near the explosion location provided by NASA in Fig. 11, which means that this signal originates from the final explosion. The back-projection trajectories have only thermospheric returns, which proves that the array collected only the thermospheric infrasound signal on that day. The back-projection trajectories consist of two paths: one path that is reflected by the ground and the other not. According to the simulation in Fig. 8, some rays with a positive elevation directly reach the thermal layer from the infrasound source and finally reach the array, and other rays with a negative elevation first reflect on the ground, then reach the thermal layer and finally reach the array. Although the ray trajectories reflecting on the ground are longer, there are more paths at a lower altitude, which have a greater effective sound speed, and the time to reach the array is approximately the same.

Figure 12(a) shows that most of the location points for the infrasound source have absolute errors of less than 12 km. Except for two rays, the absolute errors of the latitude are below 0.26° in Fig. 12(b) and the absolute errors of the longitude are below 0.31° in Fig. 12(c). The average latitude and longitude of these location points is 43.2362°N, 115.9699°E. The horizontal average distance to the array is 465.645 km, and a relative error of approximately 0.34%. The average altitude is 28.2281 km, and a relative error of approximately 27.15%.

All of the above are carried out under the assumption that the time of the meteorite explosion has been known. If the latitude and longitude of the meteorite explosion are known, the infrasound source altitude can also be calculated.

Through the latitude and longitude of the meteorite explosion and the array, the horizontal distance from the center of the array to the meteorite explosion position is calculated to be approximately 467.22 km. As shown in Fig. 13, when the horizontal distance reaches 467.22 km, the corresponding ray altitudes range from 6.24 to 50.29 km. Among these altitudes, the radiation at 17.42–36.23 km is dense, and the altitude of the infrasound source should be within this range.

Figure 14(a) shows that most of the location points on the back-projection trajectories for the infrasound source have absolute errors of less than 13 km, an average altitude of 27.7415 km, and a relative error of approximately 25.53%. Except for one ray, the absolute errors of the latitude and longitude are below 0.25° in Fig. 14(b) and Fig. 14(c), and the corresponding average is 43.2472°N, 115.9488°E, which is approximately equivalent to the horizontal distance of 10 km.

For the location result, some partial trajectories deviate from the infrasound source, which may be caused by scattering, atmospheric turbulence or noise. In addition to these paths, other rays can accurately locate the latitude, longitude and altitude of the infrasound source. The average result is consistent with the data given by NASA.

Based on a three-dimensional Hamiltonian ray tracing program in spherical coordinates, this paper proposes a back-projection method to locate the altitude of an infrasound source. Combined with the MSISE00 and HWM93 models, infrasound propagation and back-projection paths are simulated and analyzed. According to Hamilton's canonical equations, the effective sound speed determines the trajectory of the ray when the initial position of the ray and initial wavevector are known. The theoretical model of the back-projection trajectory is derived by transforming the propagation direction along the trajectory. The effective sound speed whose minimum values are the main reason for atmospheric returns is a superposition of the sound speed and the inner product of the wind vector and the wavefront normal unit vector. The sound speed is determined by the temperature and will increase sharply at approximately 120 km, forming the thermospheric return that always exists regardless of the season, latitude, longitude, and propagation direction. At an altitude of 100 km or less, there is almost no meridional wind, and the zonal wind plays a major role. The wind has seasonality and directionality, which are the main reasons for the stratospheric return. Tailwind promotes the return, and dead wind suppresses the return. Starting from the numerical simulation results, the back-projection model exhibits good accuracy for the altitude, latitude, and longitude of the back-projection stratospheric and thermospheric paths in the conventional model. The back-projection propagation trajectories are basically consistent with the true values for the meteorite explosion acquired by means of the PMCC method, further explaining the accuracy of the back-projection algorithm.

For long-distance propagation, acoustic scattering can cause path shifts. Atmospheric turbulence affects the angle of incidence of the array. MSISE and HWM are not real-time models, so the atmospheric state on a given day may be inconsistent with that from an empirical model. These factors lead to some errors in the back-projection model. If a more accurate real-time atmospheric model is obtained, the positioning results will be more accurate.