Numerical method for solving one bathymetry problem Elizaveta R. Liua, Andrei A. Sushchenkoa,b, Vladimir A. Kana and Alexander Yu. Chebotareva,b a Far Eastern Federal University, 8 Sukhanova St., Vladivostok, 690950, Russia b Institute of Applied Mathematics, Far Eastern Branch, Russian Academy of Sciences, 7 Radio St., Vladivostok, 690041, Russia Abstract Bathymetry problem based on the mathematical model of the acoustic signal propagation is considered. In the case of the single scattering approximation we obtained solution of the direct problem. As a solution to the inverse problem, which consists in determination of the function describing deviation from a reference value, we obtained a non-linear differential equation with some assumptions. A numerical analysis of the bathymetry problem is conducted for real data and the influence of the bottom scattering coefficient on the bathymetric function reconstruction is investigated. Keywords 1 radiation transfer equation; diffuse scattering; side-scan sonar; receive antenna pattern, bottom scattering 1. Introduction Investigation of the ocean floors is still priority for the world community. Many research complexes are developed to solve bathymetry problems. Currently, ocean mapping approach using an autonomous unmanned underwater vehicle, which are equipped with the side-scan sonars, is very relevant and promising. A sonar emits pulses of sound and detects echoes. Gauging water depth using acoustic (sonar) technology involves measuring the time taken for sound waves to travel between the vessel and the seafloor. Acoustic image is formed on the starboard and the port side of the underwater vehicle while the sonar antenna moves. It is worth to note that there are methods of the satellite bathymetry [1], [2]. An approach proposed in the article could be extended for this technology. To describe sound propagation in a fluctuating ocean, mathematical model based on the radiative transfer equation is used (see, e.g., [3], [4], [5]). Also, the kinetic model could be applicable in various scientific fields, e.g. computed tomography [6], [7], visualization of images, optical light transmission in the Earth atmosphere [8], modeling thermal and radiative processes in bio-tissues during sun irradiation [9], laser therapy [10], etc. In the paper, the transport equation is used to simulate the propagation of acoustic waves energy in the case of multiple scattering, and the mathematical model takes into account the sea bottom scattering [11], [12]. 2. Mathematical model Process of acoustic waves propagation on frequencies of the order of tens kilohertz is described by the following integral-differential radiative transfer equation [13]: 1 πœ•πΌ 𝜎 + π’Œ β‹… βˆ‡r 𝐼(𝒓, π’Œ, 𝑑) + πœ‡ 𝐼(𝒓, π’Œ, 𝑑) = ∫ 𝐼(𝒓, π’Œβ€² , 𝑑)π‘‘π’Œβ€² + 𝐽(𝒓, π’Œ, 𝑑). (1) 𝑐 πœ•π‘‘ 2πœ‹ Ξ© Far Eastern Workshop on Computational Technologies and Intelligent Systems, March 2–3, 2021, Khabarovsk, Russia EMAIL: liu.er@dvfu.ru (E.R. Liu); sushchenko.aa@dvfu.ru (A.A. Sushchenko); kan_va@dvfu.ru (V.A. Kan), chebotarev.ayu@dvfu.ru (A. Yu. Chebotarev) ORCID: 0000-0001-9572-0913 (E.R. Liu); 0000-0001-8573-7172 (A.A. Sushchenko); 0000-0001-6122-2000 (A. Yu. Chebotarev) ©️ 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) where 𝒓 ∈ ℝ2 , t ∈ [0, T], the wave vector π’Œ belongs to the unit sphere Ω = {k ∈ R2 ∢ |k| = 1}. The function 𝐼(𝒓, π’Œ, 𝑑) is the wave’s radiation intensity in the moment 𝑑 and the point 𝒓, propagated in the π’Œ direction with the sound speed 𝑐. The coefficients πœ‡ and 𝜎 denote spatially varying coefficients of the attenuation and the scattering, respectively, 𝐽(𝒓, π’Œ, 𝑑) describes the source of the sound field. The process of echo signal propagation is considered in the domain 𝐺: = {𝒓 ∈ ℝ2 : π‘Ÿ2 > βˆ’π‘™ + 𝑒(π‘Ÿ1 )}, which is the upper half-space bounded by a curve πœ•πΊ = 𝛾 = {π’š ∈ ℝ2 : 𝑦2 = βˆ’π‘™ + 𝑒(𝑦1 )}, where the function 𝑒(𝑦1 ) describes the deviation of the sea bottom relief. Assuming there are no sound sources in the medium while 𝑑 < 0, we supplement equation (1) with the initial condition: 𝐼|𝑑<0 = 0, (2) and the boundary condition on the surface Ξ³ which obeys Lambert’s law: 𝐼(𝒓, π’Œ, 𝑑) = 2πœŽπ‘‘ ∫ |𝒏(𝒓) β‹… π’Œβ€²| 𝐼(𝒓, π’Œβ€², 𝑑)π‘‘π’Œβ€² , 𝒓 ∈ 𝛾, π’Œ ∈ Ξ©βˆ’ (𝒓). (3) Ξ©+ (π’š) Here, Ω± (𝐲) = {π’Œ ∈ Ξ©, sgn(π’Œ β‹… 𝒏(π’š)) = Β±1}, πœŽπ‘‘ denotes the constant sea bottom reflection coefficient, 𝒏(π’š) denotes the external normal to πœ•πΊ. Bathymetric function is involved in the boundary condition. Moreover, we introduce the function 𝐽(𝒓, π’Œ, 𝑑) to describe a point of isotropic sound source: 𝐽(𝒓, π’Œ, 𝑑) = 𝐽0 𝛿(𝒓) 𝛿(𝑑). (4) Here, 𝛿 denotes the Dirac delta function and 𝐽0 is the source power. The initial-boundary value problem (1) – (3) is supplemented by the condition on the vehicle: ∫ 𝑆 Β± (π’Œ)𝐼|Π“Β± (𝑢, π’Œ, 𝑑)π‘‘π’Œ = 𝐼 Β± (𝑑), (5) Ξ©+ (π’š) where 𝐼 Β± (𝑑) define the total intensity of the received signal on the starboard and the port side. The functions 𝑆 Β± (π’Œ) characterize the radiation pattern of the receiving antenna on the starboard and the port side of the carrier, respectively. 2.1. Direct problem We will use approximation of non-scattering media (𝜎 = 0), so the solution of the direct problem in the single scattering approximation can be represented in the following form: 2 8 πœŽπ‘‘ 𝐽0 exp(βˆ’πœ‡π‘π‘‘) ( 𝑦1 𝑒′ 𝑦1 + 𝑙 – 𝑒(𝑦1 )) 𝐼 Β± (𝑑) = 𝛸[0,±∞] , (6) 2 𝑐 2 𝑑 3 |𝑒′ 𝑦1 (𝑙 βˆ’ 𝑒(𝑦1 )) βˆ’ 𝑦1 | √1 + (𝑒 𝑦1 ) β€² where π’š, 𝒛 ∈ 𝛾. It is worth to note that the first term in equation (6) corresponds to a single scattered signal (𝐼1 ), and the second is to the double scattering (𝐼2 ). 2.2. Inverse problem Solution of the inverse problem was obtained as a nonlinear differential equation for the function 𝑒 in the single scattering approximation and a narrow radiation pattern of the receiving antenna. Thus, for the numerical solution of the nonlinear differential equation the following scheme was constructed: 1 4 𝐼 Β± (𝑑𝑖 )𝑐 2 𝑑𝑖3 |𝑣0,π‘–βˆ’1 (𝑙 βˆ’ π‘’π‘–βˆ’1 ) βˆ’ 𝑦1,𝑖 | 𝑒𝑖′ = 2 (π‘’π‘–βˆ’1 βˆ’ 𝑙 + √1 + 𝑣0,π‘–βˆ’1 √ ), (7) 𝑦1,𝑖 8πœŽπ‘‘ 𝐽0 exp(βˆ’πœ‡π‘π‘‘π‘– ) 2 1/2 where 𝑑𝑖 = 2𝑐 βˆ’1 (𝑦1,𝑖 + (𝑙 βˆ’ π‘’π‘–βˆ’1 )2 ) , 𝑣0,𝑖 = 𝑒𝑖′ and 𝑒(0) = 0. Value of the 𝑣0 is set as derivative 𝑒′, which is calculated in the previous node. This approach requires a second initial condition, which can be obtained by solving an algebraic equation of forth degree using the initial condition 𝑒 (0) = 0. In the numerical algorithm for relation (7), we used the modified Euler’s method with an accuracy of 1.0𝐸 βˆ’ 10. 3. Numerical experiments In the case of a narrow directivity pattern of the receiving antenna, the problem of remote sensing by the SSS, moving with a constant velocity 𝑽 along the axis π‘Ÿ3 , is reduced to solving the problem (1) – (4) and is solved independently at each probing interval. As said before, the modified Euler’s method is used to conduct computational experiments. The sounding parameters for the computational experiments are presented in Table 1. Table 1 Probing parameters. πœ‡, π‘šβˆ’1 πœŽπ‘‘ 𝑐, m/s 𝐽0 𝑙, m 𝑦1 , m 𝑦3 , m 0.018 1 1500 1 10 [0, 300] [0, 80] The function that describes the topography of the seabed has the following form: 𝑁 𝑁 𝑦1 𝑦3 𝑒(𝑦1 , 𝑦3 ) = βˆ‘ π‘Žπ‘– 𝑠𝑖𝑛 (2πœ‹π‘– 𝐾𝑦1 ) + βˆ‘ 𝑏𝑖 𝑠𝑖𝑛 (2πœ‹π‘– 𝐾𝑦3 )+ 𝑦1 π‘šπ‘Žπ‘₯ 𝑦3 π‘šπ‘Žπ‘₯ 𝑖=1 𝑖=1 𝑁 𝑁 (8) 𝑦1 𝑦3 + βˆ‘ 𝑐𝑖 π‘π‘œπ‘  (2πœ‹π‘– 𝐾𝑦1 ) + βˆ‘ 𝑑𝑖 π‘π‘œπ‘  (2πœ‹π‘– 𝐾𝑦3 ), 𝑦1 π‘šπ‘Žπ‘₯ 𝑦3 π‘šπ‘Žπ‘₯ 𝑖=1 𝑖=1 where π‘Žπ‘– , 𝑏𝑖 , 𝑐𝑖 , 𝑑𝑖 are random coefficients evenly distributed on the interval [0, 1], 𝐾𝑦1 , 𝐾𝑦3 are the number of local extremums on the bottom surface, 𝑁 denotes the number of harmonics. Figure 1 shows the surface which describes the sea bottom relief obtained by the formula (8). The bottom was generated with two local extremums in each direction (𝐾𝑦1 = 𝐾𝑦3 = 2) and ten harmonics (𝑁 = 10). Such approach for the generation of the bottom surface allows to create more realistic relief of the seabed. Figure 2 presents graph of the error βˆ†π‘’ in the case, when the measured signal 𝐼(𝑑) was calculated by formula (6) which takes single scattering signal into account. The maximum error from figure 2 is about 0.5% which is small enough for such an oscillating bottom and it indicates high accuracy of the numerical method. It is worth to note that this experiment was conducted with 100% reflection of the signal from the bottom, i.e. πœŽπ‘‘ = 1. However, in real experiments the bottom scattering coefficient is up to 10% of total reflected signal. Hence, the error can be reduced by a factor of 10. Thus, the aim of the next experiment is to analyze solution of the bathymetry problem with variable coefficient πœŽπ‘‘ ∈ [0.1, 1] in sing scattering approximation. Figure 1: Exact solution 𝑒(𝑦1 , 𝑦2 ). Figure 2: Error π›₯𝑒 in the single scattering approximation. 3.1. Experiments with real data Earlier in paper [14], using the methods of the theory of radiation transfer, an explicit formula was obtained for determining the bottom topography function: 𝑦12 𝑒(𝑦1 , 𝑦2,𝑗 ) = Γ— πœŽπ‘‘ (2𝑙 βˆ’ 𝑙3 /(𝑦12 + 𝑙2 )) Β± 2 2 𝐼 (𝑑)2πœ‹ 𝑦1 (𝑦1 + 𝑙 ) exp (2πœ‡βˆšπ‘¦1 + 𝑙 ) 2 2 (9) πœŽπ‘‘ 𝑙2 Γ—( 2 βˆ’ ), 𝑦1 + 𝑙2 πœŽπ‘‘ 𝑐𝐽𝑗 where 𝑦2,𝑗 = 𝑉𝑑𝑗 , 𝑦1 – the bottom point, 𝑉 – the speed of the source, 𝐽𝑗 – the power of the sound source. Formula (9) was obtained in the single scattering approximation and with the condition that the bottom topography function is weakly varying 𝑒𝑦′ 1 β‰ͺ 1, 𝑒𝑦′ 2 β‰ͺ 1. A numerical analysis of formula (9) was carried out on the basis of synthetic data in paper [15], however, the practical application of the obtained mathematical model is very interesting. Thus, we applied the formula (9) to real data and investigated the effect of bottom scattering. The sounding parameters for the computational experiments are presented in Table 2. Table 2 Probing parameters. πœ‡, π‘šβˆ’1 V, m/s 𝑐, m/s J 𝑙, m 𝑦1 , m 𝑦3 , m 0.01 1 1500 9.156 10 [0, 100] [0, 20] In this case we calculated the bottom scattering coefficient πœŽπ‘‘ in each point, then based on the found values, the function 𝑒 (𝑦1 , 𝑦2 ) corresponding to the seabed profile was calculated. Further, we fixed the coefficient of the bottom reflection πœŽπ‘‘ at each point and, with the constant coefficient, we tried to reconstruct the bottom surface. Figure 3 shows the seabed profile at varying values of the bottom reflection coefficient at each point. Some reverberations of the seabed profile are visible here at 10-40 meters and according to the color scale, we concluded these are hills, and the darkening in the place about 25 and 30 meters after the peaks is explained by the presence of depressions, which’s deep is about 3-5 meters creating shaded areas after. Figure 3: Reconstructed seabed profile by finding the bottom reflection coefficient at each point. Figure 4: Reconstructed seabed profile with constant bottom reflection coefficient Figure 4 shows the reconstructed values of the seabed function 𝑒 (𝑦1 , 𝑦2 ) at a constant coefficient πœŽπ‘‘ . In comparison with Figure 3, where coefficient πœŽπ‘‘ was calculated at each point of the bottom, in that case the function 𝑒 (𝑦1 , 𝑦2 ) is restored better than with constant πœŽπ‘‘ . However, it should be noted that Figure 4 still reflects the main disturbances of the function 𝑒 (𝑦1 , 𝑦2 ), which conveys the main information. Figure 5: Reconstructed seabed profile by finding the bottom reflection coefficient at each point. Figure 6: Reconstructed seabed profile with constant bottom reflection coefficient. In this experiment, the bottom is considered with less changes than in previous one. On Figure 5 over the entire interval we can see the influence of the bottom reflection on the restoration of the function 𝑒 (𝑦1 , 𝑦2 ). Figure 6 shows that neglecting some values of the coefficient πœŽπ‘‘ and considering only one constant value, the image of the reconstructed bottom shows less information than in Figure 5. 4. Conclusion The results support the hypothesis that solution of the bathymetry problem in the single scattering approximation is stable. We simulated the side-scan sonar signal in the single scattering approximation and used it as input parameter for the problem of reconstructing the seabed topography, the solution of which was obtained. The experiments show that bottom scattering significantly affect the restoration of the seabed. Therefore, it is necessary to search for πœŽπ‘‘ for each bottom point. Formula (9) was obtained with conditions of a weakly varying bottom, i.e. formula (9) did not contain the function 𝑒′(𝑦1 , 𝑦2 ), which describes the seabed topography. Such conditions greatly facilitate the adaptation of this mathematical model to real data. It is worth noting that working with real data requires varying many values and selecting the most appropriate value for each of them. So there are difficulties with testing on real data the formula (7) in case of a strongly oscillating bottom. The presence of the derivative 𝑒′(𝑦1 , 𝑦2 ) in formula (7) generates a nonlinear differential equation. 5. Acknowledgments The work was performed as part of the state assignment β„– 075-01095-20-00, funded by Russian Foundation for Basic Research (project number 20-01-00173) and supported by the Ministry of Science and Higher Education of the Russian Federation (contract no. 075-02-2020-1482-1, 21.04.2020). 6. References [1] G. Doxani, M. Papadopoulou, P. Lafazani, C. Pikridas, M. Tsakiri-Strati, Shallow-water bathymetry over variable bottom types using multispectral worldview-2 image, International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences (2012) 159- 164, volume 39. doi:10.5194/isprsarchives-XXXIX-B8-159-2012 [2] Y. Kazama, T. Yamamoto, Shallow water bathymetry correction using sea bottom classification with multispectral satellite imagery, in: Proceedings of SPIE The International Society for Optical Engineering, volume 10422, Warsaw, Poland, 2017. doi:10.1117/12.2293011. [3] C. Poulin, X. Zhang, P. Yang, Y. Huot, Diel variations of the attenuation, backscattering and absorption coefficients of four phytoplankton species and comparison with spherical, coated spherical and hexahedral particle optical models, Journal of Quantitative Spectroscopy and Radiative Transfer (2018) 288 - 304, volume 217. doi:10.1016/j.jqsrt.2018.05.035. [4] I. Prokhorov, A. Sushchenko, Analysis of the impact of volume scattering and radiation pattern on the side-scan sonar images, in: Proceedings of the 5th. Pacific Rim Underwater Acoustics Conference, volume 24, Vladivostok, Russia, 2015. doi:10.1121/2.0000159. [5] J. A. Turner, R. L. Weaver, Radiative transfer of ultrasound, Journal of the Acoustical Society of America (1994) 3654-3674. doi:10.1121/1.410586. [6] D. S. Anikonov, I. V. Prokhorov, The statement and numerical solution of an optimization problem in X-ray tomography, Computational Mathematics and Mathematical Physics (2002) 16-22. doi:10.1134/S0965542506010040. [7] I. V. Prokhorov, I. P. Yarovenko and V. G. Nazarov, Optical tomography problems at layered, Inverse Problems (2008), volume 24. doi:10.1088/0266-5611/24/2/025019. [8] S. M. Prigarin, Monte Carlo simulation of the effects caused by multiple scattering of ground- based and spaceborne lidar pulses in clouds, Atmospheric and Oceanic Optics (2017) 79-83. doi:10.1134/S1024856017010110. [9] I. V. Krasnikov, A. Y. Seteikin, A. P. Popov, Simulation of the effect of photoprotective titanium dioxide (TiO2) and zinc oxide (ZnO) nanoparticles on the thermal response and optical characteristics of skin, Optics and Spectroscopy (2015) 668-673. doi:10.1134/S0030400X15040116. [10] A. E. Kovtanyuk, A. Y. Chebotarev, A. A. Astrakhantseva, A. A. Sushchenko, Optimal Control of Endovenous Laser Ablation, Optics and Spectroscopy (2020) 1508-1516. doi:10.1134/S0030400X20090131. [11] A. A. Sushchenko, V. A. Kan, E. R. Lyu, Double-scattering approximation for the bathymetry problem, in: Proceedings of the 25th. International Symposium on Atmospheric and Ocean Optics: Atmospheric Physics, volume 11208, 2019, pp. 112085Q. doi:10.1117/12.2540982. [12] V. A. Kan, I. V. Prokhorov, Reconstruction of the Lambert Curve in a Scattering Medium by Using Pulsed Sounding, Journal of Applied and Industrial Mathematics (2020) 321-329, volume 14. doi:10.1134/s1990478920020106. [13] S. Y. Zinkov, A. A. Sushchenko, K. V. Sushchenko, Analysis of surface and volume scattering in the problem of seabottom sounding, Siberian Electronic Mathematical Reports (2018) 1361-1377, volume 15. doi:10.17377/semi.2018.15.112. [14] I. V. Prokhorov, A. A. Sushchenko, V. A. Kan, On the Problem of Reconstructing the Floor Topography of a Fluctuating Ocean, Journal of Applied and Industrial Mathematics (2015) 412- 422, volume 9. doi:10.1134/S1990478915030126. [15] I. Prokhorov, A. Sushchenko, V. Kan, E. Kovalenko, Simulation of Sonar Signal Propagation in a Fluctuating Ocean, in: Physics Procedia, volume 70, 2015, pp. 690–694. doi:10.1016/j.phpro.2015.08.087.