The UXO for the simulation modeling derives from projectiles recovered from the site that have been tested safely by professional authorities. The retrieved unexploded ordnance is measured by hand and modeled with Auto CAD-2018, then imported into COMSOL. The sphere and circular transmitting coil have a better symmetry, and the 2D axial-symmetry of AC/DC module in COMSOL is chosen for modeling, which can reduce the calculation and the simulation time. The analytical and finite element solutions were verified by calculating the spheres created under the identical conditions, owing to the fact that the analytical equations for other 3D conducting bodies could not be given accurately. The variation of the decay curves of the analytic solution of the secondary induction potential and the numerical solution of the secondary induction field in the homogeneous full space is shown in Fig. 2. The simulation model was set up as a homogeneous full-space geoelectric model. The air conductivity is 1E-6 S/m, the relative permeability and relative permittivity are both set to 1, the radius of the circular transmitter loop is 0.26m, the current magnitude is 5A and the number of turns is 30. The target body is located directly below the circular transmitting coil with the parameters σ = 1.12e7 S/m, a = 0.032m and h = 0.4m respectively. The analytical and numerical solutions calculated by Eq. (1) are selected for comparison to verify the reliability of the simulations under the identical conditions. The analysis of the numerical solution compared with the theoretical solution shows that the curve patterns have been basically consistent. The induced voltage arising from the secondary field was chosen as a standard of measurement because it is a key parameter that can be used to discriminate the target body during its excitation by the primary field under the ground. The induced voltages are calculated in COMSOL and will be used for comparison and analysis shown in the subsequent diagram. When the receiving point is positioned at the center of the transmitting coil, the analytical solution for the central point of a circular loop with equal area can be utilized to verify its numerical accuracy. In the subsequent data processing, we normalize the induced voltage and convert it to \({{V(t)} \mathord{\left/ {\vphantom {{V(t)} q}} \right. \kern-0pt} q}\). Where, *q* is the effective area of the receiving coil and *I* is the transmitting current.

Figure 2 shows the variation trend of the induced voltage received by the receiving coil with time, and the curve is plotted on the logarithmic coordinate. The curves are essentially the identical for the late conditions at\(\text{t}>>\tau\)which can be seen from the Fig. 2. The error in the numerical simulation for the late conditions is approximately 4%, indicating that it is feasible to carry out the transient electromagnetic method forward numerical simulation utilizing COMSOL MULTIPHYSICS software. Moreover, the numerical solution is larger than the analytical solution in the early response, indicating that more target body characteristics can be obtained by the numerical solution and the signal-to-noise ratio of the signal can be greatly improved.

Numerical Solution Modeling

In the COMSOL modeling space, each model is centered on a\(25 \times 25 \times 25{m^3}\)cube, which represents the geologic layer. This "ground" block is located below another block of \(25 \times 25 \times 25{m^3}\), which includes the height of the air and probe sensors, where z = 0 represents the earth surface. The set-up of each cubic block has followed the equation relationship shown in (4) and (6), and the relative permeability of each part is 1. The rest of the boundary conditions are selected by default. A very critical factor in finite element calculations is determining the size of the discretized mesh. The mesh discretization in finite elements is the partitioning of the geometric model into simple small units by the idea of discretization. In COMSOL, the geometric model can be partitioned into tetrahedral mesh cells, with every domain consisting of tetrahedral cells. However, the size of the partitioned mesh cells varies between domains. The basis for this is that, considering the depth, attitude, and location of the UXO model, the distribution of the response of buried ferrous material to magnetic fields is different across the ground. The coils on the ground need to receive the secondary fields generated by the UXO, and the magnetic field variations in the upper half of the ground block are of more interest to us compared to the lower half. Therefore, smaller cell meshes are used for dissection near the ground surface and in the UXO region, while larger mesh dissections are applied to regions further from the ground surface. A significant consideration is that in COMSOL numerical simulations, the precision of the solution of a geometric model is related to the quality of the cell mesh. A better cell mesh quality provides a more accurate numerical resolution, but at the cost of longer time and a higher computer configuration. We have optimized the mesh size of our geometric models through computer configuration while ensuring that the modeling produces correct and reliable results when compared to the analytical solution14. In order to take into account both efficiency and calculation accuracy, grid encryption is carried out at the transmitting loop and receiving point, and local grid encryption is also carried out near the surface and target body. We used the above meshing scheme to perform the grid dissection of the modelling region and obtained an average grid cell mass of 0.75 and a minimum grid cell mass of 0.2, where the optimal grid cell mass is 1. The Degrees of Freedom (DOF) within the entire model is 979266. The DOF is related to the number of nodes, the computational points that define the shape of each finite element. The average cell mass of the encrypted local grids of the target body, the transmitting source, the receiving point and the part near the ground surface is 0.95 and the minimum cell mass is 0.8. We perform mesh convergence tests on the forward model discussed in this paper. We have simulated the grid thickness of the air domain and the ground block away from the ground surface and found that the grid size has a negligible effect on the aspect of the secondary induction field, and there is no significant difference in the induced voltage from the fine-grid to the coarse-grid. The numerical simulation results are compared with the analytical solutions, and the results are basically consistent. Therefore, we have chosen to simulate a mesh grid that does not affect the accuracy of the induced voltage calculation and run as stably as the computer configuration permits.

Figure 3(a) shows a schematic of the meshing of the 82mm shell and other geometric regions in COMSOL modelling. The geometric center of gravity of the UXO is used as the parameter for its location and the depth of burial of the UXO is 0.5m in the modelling. Figure 3(b) shows the spatial coordinate system used in COMSOL modeling and the height of the sensor from the soil surface, respectively. The resistivity of the soil is\({\rho _1}=100\Omega \cdot m\). In Fig. 3(b), θ is the angle between the main axis of the axisymmetric target body (red dashed arrow) and the z-axis, and \(\varphi\) is the angle between the projection of the main axis in the xoy-plane and the x-axis.

Numerical Solution Analysis

The In the modeling, air conductivity air σ = 1E -6 S/m, soil resistivity\({\rho _1}=100\Omega \cdot m\), the material of the UXO is iron, and its conductivity is σ = 1.12e7 S/m. The parameters of the sensor and the buried depth of the UXO are consistent with Fig. 3. We investigated the influence of the attitude of UXO on the time-domain response by modeling UXO at different angles. The time-domain three-component response at the receiver point of the homogeneous half-space model was calculated in the observed time range. Figure 4 shows the time domain response of the 82 mm mortar projectile (length 31.6 cm and diameter 82 mm) for three components (Vx,Vy,Vz) and the target body response for inclination angles \(\theta\) of \({0^ \circ },{30^ \circ },{60^ \circ }\) and \({90^ \circ }\), respectively. Although the decay trend of the three-component time domain response is almost identical, the amplitude of the induced voltage response for Vx and Vy is much smaller than that of the induced voltage response for Vz. The UXO has an axisymmetric properties with a similar structure in the x and y directions, and it can be found that Vx and Vy almost coincide. The difference between the decay curves of A and B at the late stage is mainly caused by the asymmetric distribution of fins in the tail of the unexploded ordnance. Therefore, only its two component induced voltages need to be known to obtain the target properties.

It can be seen from the Fig. 4(a) that the long axis Vz of the induced voltage response is much larger than the short axis Vx (or Vy), and that Vz decays more slowly than Vx (or Vy) and takes longer to decay. The time-domain Z-component response at receiver 0 for the UXO model is calculated separately for different deflection angles *θ* in Fig. 4(b). In the Z-component response of different inclination angles, the initial amplitude of the attenuation signal increases slightly with the increase of the angle (*θ*), but the overall trend is basically the same. The long axis response of the target is not only greater than the short axis response, but also the attenuation rate is slow. When *θ* = 0°, the long axis is excited by the primary field and produces a strong secondary field response, which decays at a slower rate. As shown in Fig. 4 (b), the late decay rate of the target body Z-component response accelerates as the target body deflection angle increases. The initial amplitude and decay rate of the electromagnetic response of the target body is related to the inclination of the target body. The effect of the long and short axes on the time domain decay signal is mainly in the decay rate as the inclination of the target body varies. The time domain attenuation signal provides a simple identification of the electromagnetic properties of the target body.

In order to investigate the time-domain electromagnetic response characteristics of the 82 mm UXO calculated by the 3D forward model, field experiments and analyses were carried out. Meanwhile, the UXO was placed vertically in a deep pit (head down tail up). Head down, tail up is the typical attitude of UXO in soil. The target body was placed and backfilled with soil as shown in Fig. 5. We chose a low EMI environment for the field tests, where the distances to high voltage lines and roads were more than 1km and 10m respectively. The sensor is placed horizontally on the ground, the distance between each measuring point is 0.5m and the length of each measuring line is 4m (As shown in Fig. 5). The transmitting and receiving coils are 1m *1m square coincident loop devices with an excitation emission frequency of 75Hz. The target body was buried at a depth of 1m, and a single measuring point was selected for fixed-point detection. The transient electromagnetic sensor is placed directly above the target, and each measuring point is detected four times to eliminate random errors. Finally, we superimposed and averaged the four detection results to obtain the time-domain electromagnetic response curve as shown in Fig. 6.

The results of the measured time-domain electromagnetic response of the 82mm unexploded ordnance in the field are shown in Fig. 6. It is easy to find that the three-dimensional forward numerical simulation results are basically consistent with the trend of the measured results. In this section, the reliability of the finite element 3D forward geoelectric simulation proposed in this paper is verified by comparing the numerical simulation results of UXO in different attitudes with the measured data. Next, we perform the subsequent transient electromagnetic response characterization by adopting the experimentally verified and reliable 3D geoelectric forward model.

The time domain decay signal can be a simple identification of the detection target electromagnetic characteristics, the exploration of further identification of the target also need to contain more information amount of line and plane response 15,16. The measurement line response of the target is the series of transient electromagnetic responses obtained by moving the detection system along a fixed line with a specified step size (shown in Fig. 7). As shown in Fig. 3(c), the center of the UXO as the origin of the coordinate system, the target body projection and the x-axis angle *φ* = 0, the target body long axis in the xoz-plane. The detection target position can be expressed in coordinates and the inclination angle with z-axis is θ. The probe system moves along the line of measurement in the x-axis direction, and more target characteristics can be obtained by calculating the time domain response of both Vx and Vz components. The target body measurement line response can be obtained when the detection system moves to \(({x_0},0,0)\)as shown in Fig. 7.

As shown in Fig. 7, three mortar shells were placed at different attitudes at \(h \geqslant 0.5m\) directly below the measurement line. For the purpose of achieves lightweight and high localization capability while requiring sufficient detection data for subsequent inversions, we designed and proposed one-transmission, three-receiver sensor. The sensor contains three receivers all in the same straight line, located at the center of the transmitting coil at Receiver 0, Receiver 1 at a distance of -D and Receiver 2 at a distance of D (where, D = 0.52m is the diameter of the transmitting coil). The probe system in the − 5–5m measuring line direction with 20cm as the distance, a total of 51 points measured, to get the target body inclination and depth of the line response results for \(( - {30^ \circ }, - 0.5m)\),\(({45^ \circ }, - 0.5m)\) and \(({0^ \circ }, - 0.8m)\) respectively. The detection results for t = 0.1ms and t = 1ms under the z-component response of the receiver coils R1 and R2 are shown in Fig. 7. As shown in Fig. 7, the target bodies in different attitudes are precisely located by receivers 1 and 2. There are many intersection points in the measured line response results of receivers 1 and 2, however, we choose the intersection point of the curve at the site of the maximum decline in receiver 1 results and the site of the maximum rise in receiver 2 as the localization point of the target body. The signal received by the receiver continues to increase as the sensor moves closer and closer to the pre-embedded target body, and continues to decrease as it moves away from the target body. The intersection of receivers 1 and 2 provides precise positioning information of the pre-built target body.

We can mark suspicious locations with the continuous data- sampling mode of the sensing system. The continuous acquisition mode is often used for data acquisition in field exploration. This exploration mode not only improves the work efficiency, but also can obtain enough sampling point data. Fixed-point sampling model will be utilized for further identification of the target body once these positions are marked. The fixed-point sampling model provides additional target body information, including attitude, depth, and category identification. More information on the inclination of the target body is shown in Fig. 8(a) and (b). Figure 8(a) and (b) show the measured line response for UXOs with inclinations θ = 0°, 30°,60° and 90° respectively, where the UXOs are buried at a depth of 0.5m and are located in the center of the line. As shown in Fig. 8, the x-component measurement line responses are all point-symmetric, with the response amplitude at the center point being 0. The Z-component responses are axisymmetric, with the maximum response amplitude at the axis of symmetry. As the inclination angle increases from 0° to 90°, the magnitude of both the x- and z-components change significantly, which allows good discrimination of the inclination information of the target body. The response amplitudes in Fig. 9 all decay rapidly with increasing time. The early signal amplitude is hundreds or even thousands of times that of the late signal, so the use of an effective early response can greatly improve the signal-to-noise ratio, and thus improve the detection depth (As shown in Fig. 9).

The measured line response results for 82 mm mortar rounds with a depth of 0.5-2.5m and a target body inclination of \(\theta ={0^ \circ }\) are calculated in Fig. 10(a) and (b). Figure 10(a) and (b) show the normalized x-component response and z-component response of the target body, respectively. The response amplitude decreases rapidly when the depth of the target body increases, whether it is the x-component measurement line response or the z-component measurement line response. The response amplitude decreases by about 50% for each 0.5m depth from 0.5m to 2.5m depth. The tendency of the response results will not vary significantly with increasing depth for either the x-component response or the z-component response. The magnitude of the electromagnetic response of the unexploded bomb decreases rapidly with increasing depth, and the x and z component response trends are almost invariable. Therefore, the desire to double the depth of detection requires a tens of times increase in the peak magnetic moment of the system and a reduction in the system noise level.

The 2D plane response of the target body is calculated as further identification of the target body attitude parameter variations16,17. After the position of the target body is localized, monostatic measurements are applied to calculate the component response at different attitudes. The x-component response and the horizontal component response are calculated in Fig. 11 and Fig. 13 when phi = 0°,30°,60°,theta = 0°,30°and 60°, respectively. The horizontal component response is calculated by combining the x-component and y-component responses. The target body is located at -0.5m and the plane responds to a square with a side length of 3m. As shown in Fig. 11, when *θ* = 0°, the x-component response has axisymmetric characteristics and the magnitude is symmetric about the y-axis, where the y-axis position is the target body position. The x-component response is not varied with φ when θ = 0°, as shown in Fig. 11 (a), (d) and (g), and always has the symmetry axis at the location of the target body center. As shown in Fig. 11(b), (c) or (e), (f) or (h), (i), the entire x-component response is shifted to the right and the symmetry is broken as θ continues to increase. More importantly, the maximum and minimum values are equal. The corresponding amplitude of the target body reaches a maximum when θ and phi = 0°, and the response amplitude decreases as the inclination angle increases. Moreover, the values of the maximum and minimum values are not the same. The pattern allows a rough identification of the attitude distribution of the underground target body.

As shown in Fig. 12, when theta = 0°, the y-component response moves upward as phi increases. However, it has always been symmetric about x-axis. At the same time, the y- component response has different values for the maximum and minimum values. At phi = 0°, the amplitude of the y- component response decreases as theta increases and the overall distribution shifts to the right. The y-component response is essentially similar to the x-component response, with a good symmetry in the target body x (or y)-component response for either theta = 0 or phi = 0. Furthermore, the values of the magnitude of the maximum and minimum values can also identify the attitude of the unexploded ordnance. The maximum value always corresponds to the direction of the head of the unexploded ordnance.

The horizontal component response of the target body at different attitudes is calculated in Fig. 13, where the target body is located at -0.5m. The three rows of the horizontal component response represent theta = 0°,30° and 60°, and the three columns represent phi = 0°,30° and 60°, respectively. As shown in Fig. 13, the horizontal component response amplitude has circular symmetry when \(\theta\) and \(\varphi ={0^ \circ }\). The position of the target body is its center of symmetry, and the horizontal component response has a minimal value point above the target body. The horizontal component response at the location of the target body is the smallest, and the outward response amplitude increases first and then decreases. The horizontal component response has great and small values on both sides of the long axis position of the target body when θ = 30° and 60°. The target body horizontal component response varies with phi when theta > 0°. The minimum of the horizontal component shifts to the right and the maximum to the left, varying with θ. The response continues to be distributed in circular symmetry after moving away from the target body. The horizontal component of the target body is minimum at θ and phi = 60°, and maintains a short axis symmetric distribution around the target body. As theta increases, the horizontal response moves to the right and the amplitude decreases. The axisymmetric characteristic is gradually broken, and the extreme value on the right side hardly varies corresponding to the tail of the shell, and the extreme value on the left side first decreases and then increases corresponding to the head of the shell when θ increases from 0° to 60°. It can be seen that the center of the line connecting the maximum and minimum of the horizontal component is the location of the long axis of the target body. The rotation angle of the target body can be determined from the symmetry axis of the horizontal response component of the target body. It can be seen that the horizontal component of the planar response contains both the target body inclination information and the position characteristics.

As shown in Fig. 14, the z-component response of the target body is much higher in magnitude than the x-component response, the y-component response and the horizontal component response. When theta = 0°, the z-component response decreases as phi increases and the overall distribution is always symmetric about the y-axis. When theta = 0°, there is always circular symmetry with the position of the target body, and the z-component response is greatest at the location of the target body, decreasing in amplitude outwards. The maximum value of the z-component response of the target body follows the changes in theta and phi. The amplitude of the z-component response of the target body is much larger than the horizontal component response, which is beneficial for improving the signal-to-noise ratio. However, the z-component response identifies the attitude of the target body to a much lesser extent than the horizontal component response. Therefore, a three-component receiver coil needs to be designed to improve the detection performance of the system while obtaining the z-component and horizontal component.

Summarize

We investigated the electromagnetic response model of unexploded ordnance by COMSOL MULTIPHYSICS based on the finite element method. The results of the theoretical response of a solid nonmagnetic iron sphere and the response of a finite element model of the identical geometry are first compared. After establishing confidence in the finite element procedure, we extended the finite element modelling to the complex geometry of real UXO. Next, we adopted a homogeneous half-space for modelling the 3D geological body in order to make the UXO buried in a more realistic geographical environment. The design of the magnetic coupling source and the related parameters are derived and analyzed in detail. The electromagnetic response characteristics of UXO with transient electromagnetic excitation are explored. We conclude that the attitude and position of the target body cannot be identified by the time-domain response of the z-component in the single-transmitter-single-receiver sensor operation mode. The transient EM sensing system of single-transmitter and triple-receiver was established, which can not only realize the rapid localization of the target body, but also obtain richer EM response data to realize the inversion. The multiple attitudes and depths of the target body are investigated in detail through the line and plane responses. The combination of line and plane response provides multi-attitude information about the target body including the principal axis deflection angle and depth of the UXO. The x component is symmetric about x = 0 at θ = 0° in the measured line response. The amplitude of the x-component response on the left side of the symmetry axis is greater than that on the right side indicates θ > 0° and conversely, θ < 0°. Extremely large and small values arise in the two-dimensional plane response, as θ and \(\varphi\) vary. This is related to the shape of the front and rear ends of the UXO. The direction of the main axis of the UXO can be determined by connecting the maximum and minimum values. Reasonable geoelectric modelling combined with field experiments can not only improve the exploration efficiency but also provide assistance for the selection of instrument parameters and data interpretation. The numerical results show that optimized three-component receiver coils placed at suitable locations can greatly improve the detection efficiency and provide a reference for the preparation of sensing systems.