Multisource Seismic Full Waveform Inversion of Metal Ore Bodies

: The seismic exploration method could explore deep metal ore bodies (depth > 1000 m). However, it is difﬁcult to describe the geometry of the complex metal ore body accurately. Seismic full waveform inversion is a relatively new method to achieve accurate imaging of subsurface structures, but its success requires better initial models and low-frequency data. The seismic data acquired in the metal mine area is usually difﬁcult to meet the requirements of full waveform inversion. The passive seismic data usually contains good low frequency information. In this paper, we use both passive and active seismic datasets to improve the full waveform inversion results in the metal mining area. The results show that the multisource seismic full waveform inversion could obtain a suitable result for high-resolution seismic imaging of metal ore bodies.


Introduction
As the demand for metal ores has increased in recent years, discovering new shallow deposits has become increasingly difficult. Therefore, the target for the exploration and development of metal mineral resources has gradually shifted to deeper parts of the Earth (>1 km). Potential field and electromagnetic methods are the most common geophysical method in mineral exploration. However, these methods have a limited penetration depth with high resolution. The seismic reflection method is the only exploration method to provide high-resolution images of the subsurface and have the penetration depth required for current mineral planning and exploration for deep mineral resources. There have been many successful cases of the seismic exploration of metal mines [1][2][3][4][5].
However, the metal ore bodies are usually irregular and small. The impedance difference between the ore body and the surrounding rocks is relatively small, and the reflected energy is weak. These factors lead to poor conventional seismic imaging. Full waveform inversion (FWI) is a method that can invert seismic velocity with high precision. It is the most accurate velocity inversion method in seismic exploration and has great potential in metal mine exploration. FWI was initially proposed by Lailly [6] and Tarantola [7] and in the time domain. In the 1990s, Pratt and his colleagues extended the method to the frequency domain and made the industrial application of the method possible [8,9]. Currently, FWI has a wide range of applications in different exploration fields [10][11][12][13][14][15][16].
There are only a few applications of FWI in metal mine exploration. Sun et al. used visibility analysis and energy compensation in FWI to invert the metal ore model and achieved good results [17]. Egorov et al. tested the potential use of FWI in VSP data from hard rock environments [18]. Mao et al. used the direct wave and the adjoint-state source function inversion method to invert the source wavelet accurately and the similarity phenomenon to reconstruct the low-frequency data for the FWI application in deep metal ore bodies [19]. Hlousek et al. used FWI for velocity model building in a tailored workflow for advanced high-resolution seismic imaging of mineral exploration targets [20]. Singh et al. applied FWI for 3D velocity model building in an iron-oxide mining site in Sweden [21].

Synthetic Datasets
Our goal in generating the synthetic data was to verify the effectiveness of our inversion method. The complexity of the synthetic dataset is not comparable with the field data, so we tried our best to make the synthetic datasets close to the filed datasets. We used a velocity model of a complex metal ore body. The model was modified from a geological ore body model in the Luzong Basin, China and is shown in Figure 1. The deposit consists of several ore bodies. The overall ore body is complexly layered and gently lenticular. Its plane projection is elliptical. The ore body is dome-like in spatial expression, with leached poor iron ore dominating in the center and rich and thick ore surrounding it. Only two of them are relatively large in scale, and the rest are small [34].
We used our own finite-difference program based on the MATLAB software to model both active and passive datasets, which were used for this work. The finite-difference program uses a 8th order accuracy and PML (perfect match layer) boundary condition. The modeling parameters are shown in Table 1. The model size was 376 * 126 with a grid size of 10 m for both dx and dz. There were 376 receivers at the model surface as shown in Figure 1, marked with blue inverted triangles. The distance between each receiver is 10 m. The receiver locations are the same for both active and passive data sets. The active source locations marked with red dots are shown in Figure 1. We use a 20 Hz Ricker wavelet as the source wavelet for the active data set, and 75 shots are generated. The distance between each shot is 50 m. The passive source locations marked with red stars are shown in Figure 1. Fifty passives shots were generated. We assumed that the passive sources had a relative uniform distribution and put the passive source locations around the ore bodies. We generated two passive source datasets. First, we used a 10 Hz source wavelet for all the source locations for one passive source dataset. This dataset was used to test whether the method will work in an ideal case. To make the synthetic test more realistic, another passive dataset used different source wavelets for the other source locations. We also added random noise to different datasets. Six different datasets were generated for future tests, as shown in Table 2.
Ricker wavelet as the source wavelet for the active data set, and 75 shots are generated. The distance between each shot is 50 m. The passive source locations marked with red stars are shown in Figure 1. Fifty passives shots were generated. We assumed that the passive sources had a relative uniform distribution and put the passive source locations around the ore bodies. We generated two passive source datasets. First, we used a 10 Hz source wavelet for all the source locations for one passive source dataset. This dataset was used to test whether the method will work in an ideal case. To make the synthetic test more realistic, another passive dataset used different source wavelets for the other source locations. We also added random noise to different datasets. Six different datasets were generated for future tests, as shown in Table 2.
In order to simulate the lack of low-frequency information in the active source data, we performed high-pass filtering with low cut frequency at 5 Hz on the active source datasets. Figure 2 shows the example shot gathers of the active and passive synthetic seismic datasets. Figure 2a is the active seismic shot gather that was shot at the center of the survey, and Figure 2c is the same shot gather with random noise. Figure 2b is an example passive seismic shot gather with a 10 Hz source wavelet, and Figure 2d is the same shot gather with random noise.   In order to simulate the lack of low-frequency information in the active source data, we performed high-pass filtering with low cut frequency at 5 Hz on the active source datasets. Figure 2 shows the example shot gathers of the active and passive synthetic seismic datasets. Figure 2a is the active seismic shot gather that was shot at the center of the survey, and Figure 2c is the same shot gather with random noise. Figure 2b is an example passive seismic shot gather with a 10 Hz source wavelet, and Figure 2d is the same shot gather with random noise.   Table 2. The active and passive seismic datasets were generated for the tests.
Dataset 3: Passive dataset with 10 Hz source wavelet for all the source locations; 4.
Dataset 4: Passive dataset with different source wavelets for different source locations; 5.
Dataset 7: Ten randomly picked pass shot gathers from dataset 4 with real noise, and every 5th receiver channel is used; 8.
Dataset 8: Dataset 1 with real noise and every 5 th receiver channel is used.

Conventional Full Waveform Inversion
FWI is used to find a model that generates a synthetic seismic wavefield that fits the real observed seismic wavefield. The conventional FWI workflow is shown in Figure 3. The objective function E can be defined as:

Conventional Full Waveform Inversion
FWI is used to find a model that generates a synthetic seismic wavefield that fits the real observed seismic wavefield. The conventional FWI workflow is shown in Figure 3. The objective function E can be defined as: where u ij is the synthetic seismic data; d ij is the observed seismic data; i is the source location; and j is the receiver location. This is a nonlinear inversion problem. FWI usually uses a local optimization algorithm to solve it.
where uij is the synthetic seismic data; dij is the observed seismic data; i is the source location; and j is the receiver location. This is a nonlinear inversion problem. FWI usually uses a local optimization algorithm to solve it.

Source Independent Full Waveform Inversion
A successful full waveform inversion needs a good source wavelet estimation of the observed seismic dataset. In real conditions, it is very hard to obtain a good source estimation, especially for the passive seismic dataset. There is one approach called source

Source Independent Full Waveform Inversion
A successful full waveform inversion needs a good source wavelet estimation of the observed seismic dataset. In real conditions, it is very hard to obtain a good source estimation, especially for the passive seismic dataset. There is one approach called source independent FWI that can solve this problem [32].
We could use Green's function and source wavelet to replace the u ij (m) and d ij in Equation (1). The new misfit function is: where G and S are the Green's function and source wavelet; the * represents the convolution process; and superscripts u and d are the modeled and observed seismic data, respectively. The source independent FWI uses a special misfit function that consists of the convolution of the observed wavefields with a reference trace from the modeled wavefield and the convolution of the modeled wavefields with a reference trace from the observed wavefield.
The new misfit function can be written as follows: where u ik and d ik are the reference traces from the modeled and the observed seismic data at the kth receiver position, respectively. The misfit function E can be rewritten with G and S: The effects of the source wavelets are eliminated as the source wavelet of the observed and the modeled wavefields are equally convolved with both terms in the new misfit function.

Multisource FWI Workflow
FWI is based on the Born approximation and can be considered as a local optimization problem. The success of active source FWI depends on the good initial model and lowfrequency data. However, the real seismic data may not contain good low-frequency information or have a very low signal-to-noise ratio (SNR) in the low-frequency band. The passive seismic data usually carries low-frequency information, and when used to compensate for the active seismic data, it has a more physical basis than pure mathematical methods. Although passive seismic data have the advantages above-mentioned, its weak energy, low signal-to-noise ratio, and unknown source locations and source wavelets also limit its application in FWI.
There are many forms of joint FWI of active source and passive source seismic data. The main idea is to compensate the active source FWI by taking advantage of the richer low-frequency information contained in the passive source data. In general, there are three different approaches.

1.
Directly apply FWI to the passive source seismic data to construct an initial model for the active source seismic data FWI; 2.
Use the passive seismic data to compensate for the insufficient illumination of the active seismic data [28,30]; 3.
Merge the active source and passive source seismic data using a specific method to compensate for the low-frequency information; and 4. Use the seismic interferometry method to process the passive source data and generate virtual shot gathers, then directly inverts it to provide an initial model for the active seismic FWI [29].
In this study, we chose the first multisource FWI approach. To further carry out our research, we made the following assumptions. We assumed we could relative correctly locate the passive source locations. This is a strong assumption, but in some cases, we could achieve it. We did not consider the elastic problem. Both forward modeling and inversion were conducted using the acoustic wave equation. The workflow is shown in Figure 4.
The main idea is to compensate the active source FWI by taking advantage of the richer low-frequency information contained in the passive source data. In general, there are three different approaches.
1. Directly apply FWI to the passive source seismic data to construct an initial model for the active source seismic data FWI; 2. Use the passive seismic data to compensate for the insufficient illumination of the active seismic data [28,30]; 3. Merge the active source and passive source seismic data using a specific method to compensate for the low-frequency information; and 4. Use the seismic interferometry method to process the passive source data and generate virtual shot gathers, then directly inverts it to provide an initial model for the active seismic FWI [29].
In this study, we chose the first multisource FWI approach. To further carry out our research, we made the following assumptions. We assumed we could relative correctly locate the passive source locations. This is a strong assumption, but in some cases, we could achieve it. We did not consider the elastic problem. Both forward modeling and inversion were conducted using the acoustic wave equation. The workflow is shown in Figure 4.

Numerical Examples
To verify our multisource FWI method and workflow, we designed three different numerical tests. We used the same initial velocity model shown in Figure 5 for all tests. First, we tested the method in a near-perfect condition. Second, we tested the condition where we did not know the source wavelet for the active and passive seismic datasets. Finally, we tested how the method works with noise in both the active and passive datasets. We conducted five different tests, as shown in Table 3.

Numerical Examples
To verify our multisource FWI method and workflow, we designed three different numerical tests. We used the same initial velocity model shown in Figure 5 for all tests. First, we tested the method in a near-perfect condition. Second, we tested the condition where we did not know the source wavelet for the active and passive seismic datasets. Finally, we tested how the method works with noise in both the active and passive datasets. We conducted five different tests, as shown in Table 3.  We conducted a simple numerical test to compare the multisource FWI and conven-  We conducted a simple numerical test to compare the multisource FWI and conventional active source FWI. We used the synthetic datasets 1 and 3, as shown in Table 2. We assumed that we knew the passive source location and source wavelets for both passive and active sources. First, we applied the multiscale FWI to the active source seismic dataset. The multiscale approach was conducted by filtering the data with different frequency bands (low pass filter with a stop frequency at 24 Hz, 32 Hz, 40 Hz, and 50 Hz). The result of the lower frequency band data was used as the initial model for the higher frequency band data. The upper figure in Figure 6 shows the active source FWI result. Compared with the true velocity model, it provides a relatively correct shallow part of the model. We can roughly see some outlines of the ore body. However, it is difficult to obtain a clear idea of the ore body's shape. The main reason is that the active source survey is at the surface. The gradients were weak in the deep parts of the model. The insufficient illumination will influence the final results of FWI, and it is hard to invert the deep parts of the model. The lack of low-frequency information of the active seismic dataset also limits the accuracy of the inversion results.   Table 2). The middle figure shows that the FWI result only uses sive dataset (dataset 3 in Table 2). The lower figure shows the FWI result when using bot and passive datasets.

Test 2: Numerical Test with Unknown Source Wavelets
The previous simple test showed that our multisource FWI approach work ideal condition. In this test, we weakened the preconditions for our method. In a re the inaccuracy of the source location and wavelets will also affect the inversion resu source wavelets are more difficult to obtain compared with source locations. We lated the condition that we did not know the source wavelets for the passive and  Table 2). The middle figure shows that the FWI result only uses the passive dataset (dataset 3 in Table 2). The lower figure shows the FWI result when using both active and passive datasets.
For the multisource FWI, we used the same approach as the active FWI. However, we used the FWI result of the passive source dataset as the initial model for the active source dataset. The middle figure in Figure 6 shows the FWI result of the passive source dataset. Although the resolution of the shallow layer was not as good as the result of the active source dataset, it provided a much better result of the ore body. This is because the passive sources are located around the ore body, which has a better illumination condition. The rich information of the low-frequency information of the passive dataset also improves the stability and accuracy of the inversion. The lower figure of Figure 6 shows the FWI results of the active source dataset using the passive source dataset that was the result of the initial model. The result velocity model was significantly improved by using both passive and active source data. We could see the details of both the shallow and deep parts of the model. The shapes of the two relatively large-scale and even some small ore bodies are clearly imaged.

Test 2: Numerical Test with Unknown Source Wavelets
The previous simple test showed that our multisource FWI approach works in an ideal condition. In this test, we weakened the preconditions for our method. In a real case, the inaccuracy of the source location and wavelets will also affect the inversion result. The source wavelets are more difficult to obtain compared with source locations. We simulated the condition that we did not know the source wavelets for the passive and active datasets. We used the synthetic datasets 1 and 4 as shown in Table 2. Figure 7 shows the source wavelets used to generate the synthetic datasets. The passive source wavelets at different source locations were generated by using a 10 Hz Ricker wavelet convoluted with a random sequence. We used the source independent FWI algorithm that did not rely on the source wavelet, which can overcome the influence of the inaccuracy of the passive source wavelets on the inversion result. The low-frequency information of the passive source data can be used to obtain the large-scale structure information of the model. The resulting model of the passive source data can be used as the initial model for the FWI of the active source data that is missing low-frequency information. For the active source data FWI, we used the same multiscale approach as in the previous test.     Figure 6. This is because both the passive datasets and FWI approach were different. The lower figure of Figure 8 shows the source independent active FWI result using the passive source FWI result as the initial model. The result shows that we could obtain a similarly good result even when we did not know the correct source wavelets.

Test 3: Numerical Test with Noise and Unknown Source Wavelets
We also added noise to the passive and active source datasets to test the multisource FWI in a more realistic condition. We used synthetic datasets 2 and 5, as shown in Table  2, in this test. A relatively strong noise (30 db Gaussian white noise) was added to both the active and passive source seismic datasets. We chose this noise level because most of the real datasets could reach this level after processing. Figure 9 shows the conventional FWI result of the noisy passive source seismic dataset. The inversion failed to converge due to the noise and unknown source wavelets when using the conventional FWI approach. Then, we tried the source independent FWI approach that was the same as test 2. Figure 10 shows the source independent FWI results of the noise datasets. The result is very similar to the FWI result in test 1, which is shown in Figure 6. This means that the source independent FWI method can handle the noise dataset quite well. However, due to the influence of noise, the final inversion results were slightly worse than test 1.

Test 3: Numerical Test with Noise and Unknown Source Wavelets
We also added noise to the passive and active source datasets to test the multisource FWI in a more realistic condition. We used synthetic datasets 2 and 5, as shown in Table 2, in this test. A relatively strong noise (30 db Gaussian white noise) was added to both the active and passive source seismic datasets. We chose this noise level because most of the real datasets could reach this level after processing. Figure 9 shows the conventional FWI result of the noisy passive source seismic dataset. The inversion failed to converge due to the noise and unknown source wavelets when using the conventional FWI approach. Then, we tried the source independent FWI approach that was the same as test 2. Figure 10 shows the source independent FWI results of the noise datasets. The result is very similar to the FWI result in test 1, which is shown in Figure 6. This means that the source independent FWI method can handle the noise dataset quite well. However, due to the influence of noise, the final inversion results were slightly worse than test 1.
Minerals 2022, 12, x FOR PEER REVIEW 10 of 16 Figure 9. The conventional FWI result of the passive source data with noise. Figure 9. The conventional FWI result of the passive source data with noise.

Test 4: Numerical Test with Noise, Unknown Source Wavelets, and Less Passive Shots
In a real case, we may not get 50 usable passive sources records as we did here. We also tested the condition where we were only able to obtain 10 passive source records. We randomly picked 10 passive shots from dataset 5. We applied the same inversion procedure to datasets 2 and 6. First, we applied source independent FWI to dataset 6. Then, we used the inversion result as the initial model for dataset 2. Figure 11 shows the inversion results. From the results, we could see that even if we only used 10 passive shots, the passive dataset inversion still builds up large-scale information of the velocity structure. This provides a useful initial model for the active dataset inversion. The active dataset inversion result was comparable with the previous one.

Test 4: Numerical Test with Noise, Unknown Source Wavelets, and Less Passive Shots
In a real case, we may not get 50 usable passive sources records as we did here. We also tested the condition where we were only able to obtain 10 passive source records. We randomly picked 10 passive shots from dataset 5. We applied the same inversion procedure to datasets 2 and 6. First, we applied source independent FWI to dataset 6. Then, we used the inversion result as the initial model for dataset 2. Figure 11 shows the inversion results. From the results, we could see that even if we only used 10 passive shots, the passive dataset inversion still builds up large-scale information of the velocity structure. This provides a useful initial model for the active dataset inversion. The active dataset inversion result was comparable with the previous one.  In addition, we also tested our method with real noise data and less dense acquisition geometry. We took the real background noise from a mining area. The noise was directly

Test 5: Numerical Test with Real Noise, Unknown Source Wavelets and Fewer Passive Shots and Fewer Receivers for Both Passive and Active Dataset
In addition, we also tested our method with real noise data and less dense acquisition geometry. We took the real background noise from a mining area. The noise was directly added to the synthetic datasets. We also used only 10 randomly picked passive shots (dataset 7). We took every 5th receiver channel from both passive and active datasets to simulate the less dense acquisition. Figure 12 shows the example shot gathers with real noise and less dense acquisition geometry. The same inversion perdure was applied to datasets 7 and 8. First, we applied source independent FWI to dataset 7. Then, we used the inversion result as the initial model for dataset 8. Figure 13 shows the inversion results. The results show that the method works fine in real noise conditions and sparse acquisition geometry. Figure 11. The source independent multisource FWI results of the noise datasets with less pa shots. The upper figure shows the FWI result only using the passive dataset. The lower figure s the FWI result when using both active and passive datasets.

Test 5: Numerical Test with Real Noise, Unknown Source Wavelets and Fewer Passive Shots and Fewer Receivers for Both Passive and Active Dataset
In addition, we also tested our method with real noise data and less dense acquis geometry. We took the real background noise from a mining area. The noise was dir added to the synthetic datasets. We also used only 10 randomly picked passive shots taset 7). We took every 5th receiver channel from both passive and active datasets to ulate the less dense acquisition. Figure 12 shows the example shot gathers with real n and less dense acquisition geometry. The same inversion perdure was applied to dat 7 and 8. First, we applied source independent FWI to dataset 7. Then, we used the in sion result as the initial model for dataset 8. Figure 13 shows the inversion results results show that the method works fine in real noise conditions and sparse acquis geometry.   Figure 14 shows the final FWI result and the differences between the final FWI and the true model for each test. The ideal test gave the best result. We can clearly iden the geometry of the main ore bodies. However, the velocity of the ore body obtained  Figure 14 shows the final FWI result and the differences between the final FWI test and the true model for each test. The ideal test gave the best result. We can clearly identify the geometry of the main ore bodies. However, the velocity of the ore body obtained by the FWI was about 10% lower than the actual velocity. Figure 15 shows the comparison between the observed active shot gather (generated from the true model) and the synthetic active shot gather (generated from each multisource FWI test). We found that the data fitting was pretty good for all the tests. When adding different noise and making the data set sparse, we could still obtain a reasonable result using our multisource FWI approach.

Conclusions
The active seismic data acquired in the metal mine area usually lacks low-frequency information. It is challenging to use only active source data to obtain good imaging of the ore bodies. The multisource FWI method uses effective low-frequency signals from the original passive source information, proving the potential use of low-frequency passive source information in constructing large-scale velocity models. The passive seismic FWI can provide a good initial model for the active seismic FWI, thus reducing the dependence of the active seismic FWI on the initial model and low-frequency information. We chose a complex velocity model from a mining area and generated serval datasets to test our method. The numerical results of the perfect condition show that the multisource seismic full waveform inversion could obtain a suitable result for high-resolution seismic imaging of metal ore bodies. The method also worked fine when the source wavelets are unknown, and the datasets had a relatively high noise level by introducing the source independent approach. Figure 15. The comparison between the observed shot gather and synthetic shot gather from the different multisource FWI tests. The synthetic shot gather (black) is plotted on top of the observed shot gather (red); (a-e) are the comparison for tests 1-5. There are also some limitations to our current numerical tests. We chose a complex model built from a real mining area, added real noise to the synthetic datasets, and used sparse acquisition geometry. The synthetic datasets were still idealized compared with real datasets. We may also face many other issues in the real datasets such as complex noise characteristics, discontinuity, and complexity of the shot gather signal and elastic effects. However, these were not within our study scope. It is possible to minimize these effects by preprocessing the dataset in some real cases.
Another main assumption was that the location of the seismic sources is known. Numerical tests 4 and 5 showed that we could only use a few passive shots to build a suitable initial model for the active data inversion. We believe that it is possible to obtain a few usable passive shots in some actual cases.
The best way to test our multisource FWI approach is to apply it to an actual dataset. However, there is no suitable real dataset available for testing at present, which can be considered as research content in the future.

Conclusions
The active seismic data acquired in the metal mine area usually lacks low-frequency information. It is challenging to use only active source data to obtain good imaging of the ore bodies. The multisource FWI method uses effective low-frequency signals from the original passive source information, proving the potential use of low-frequency passive source information in constructing large-scale velocity models. The passive seismic FWI can provide a good initial model for the active seismic FWI, thus reducing the dependence of the active seismic FWI on the initial model and low-frequency information. We chose a complex velocity model from a mining area and generated serval datasets to test our method. The numerical results of the perfect condition show that the multisource seismic full waveform inversion could obtain a suitable result for high-resolution seismic imaging of metal ore bodies. The method also worked fine when the source wavelets are unknown, and the datasets had a relatively high noise level by introducing the source independent approach.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.