*Step 2: Fast non-dominated sorting procedure*

The algorithm sorts the combined population *R<sup>t</sup>* according to non-domination to obtain different non-dominated fronts *F<sup>i</sup>* .

#### *Step 3: Crowding distance estimation*

To maintain the population diversity in each non-dominated front, the crowding distance is computed; that is, the average distance between two points on either side of a particular solution. Then, repeat the process with other objective functions.

#### *Step 4: Binary tournament selection*

Individuals are selected using a binary tournament selection with a crowded comparison operator to create a mating pool: The best solutions in the combined population are those belonging to the lower-level non-dominated set *F1*. When the two solutions have an equal non-domination level, the solution with a higher crowding distance is preferred to preserve diversity.

#### *Step 5: Crossover and mutation*

The algorithm creates an offspring population from the mating pool by applying the simulated binary crossover operator and the polynomial mutation operator.

#### *Step 6: Recombination and selection*

The last stage of the NSGA-II involves combining the current population with the offspring and then selecting the best solutions from this combined population to create a new population for the next generation. The process is iterated for the desired generations to achieve the best multi-objective solution in the last generation.

#### **3. Results and Discussion**

We compiled several MATLAB functions and scripts to predict geological targets in 3D, such as 3D lithotype estimation and 3D Au grade estimation. We used unevenly distributed borehole datasets, including the lithotype data (Figure 5a) and the Au concentrations in g/t (Figure 5d). Borehole datasets show that the hydrothermal alteration has destroyed magnetite and replaced it with pyrite in the felsic volcanic rocks that are the host rocks of epithermal Au/Ag mineralization [27–30]. Therefore, the high-Au grades are accompanied by low-magnetic anomalies of the felsic rocks. This correlation made it possible to separate the high magnetics of the mafic volcanic rocks and dioritic intrusions from low-magnetic felsic volcanic rock and porphyritic intrusions based on the 3D cooperatively recovered magnetic susceptibility model.

**Figure 5.** Interpolation of geological targets with conventional methods. (**a**) Lithotype data (mafic/felsic). (**b**) Nearest-neighbor interpolation of lithotypes. (**c**) Kriging of lithotypes. (**d**) Au grade data. (**e**) Nearest-neighbor interpolation of Au grades. (**f**) Kriging of Au grades [27]. **Figure 5.** Interpolation of geological targets with conventional methods. (**a**) Lithotype data (mafic/felsic). (**b**) Nearest-neighbor interpolation of lithotypes. (**c**) Kriging of lithotypes. (**d**) Au grade data. (**e**) Nearest-neighbor interpolation of Au grades. (**f**) Kriging of Au grades [27].

Consequently, we assigned three lithotype codes to quantify the lithological variations: Code 1 for felsic volcanic rocks, code 2 for mafic volcanic rocks, and code 3 for dioritic intrusions. The lithotype codes are matched and sorted according to their magnetic properties, from low magnetics of felsic volcanic rocks and intrusions to high magnetics of mafic volcanic rocks and diorites. The datasets are interpolated with the nearest neighbor (direct gridding) and kriging methods for comparison. Figure 5 shows that traditional interpolation methods cannot add any new geological information in places without borehole data. This study aims to solve this 3D interpolation problem by reconstructing lithological and Au grade patterns from 3D geophysical images. Consequently, we assigned three lithotype codes to quantify the lithological variations: Code 1 for felsic volcanic rocks, code 2 for mafic volcanic rocks, and code 3 for dioritic intrusions. The lithotype codes are matched and sorted according to their magnetic properties, from low magnetics of felsic volcanic rocks and intrusions to high magnetics of mafic volcanic rocks and diorites. The datasets are interpolated with the nearest neighbor (direct gridding) and kriging methods for comparison. Figure 5 shows that traditional interpolation methods cannot add any new geological information in places without borehole data. This study aims to solve this 3D interpolation problem by reconstructing lithological and Au grade patterns from 3D geophysical images.

The 3D physical property images are derived from cooperative inversion [27,31] of magnetic and DC/IP datasets (Figure 2b–d). A preliminary Fast-ICA separates the latent features inside the multiple physical property images in the form of three negentropymaximized ICs (spatial feature extraction). Then, three ICs are decomposed to form a set of raw spectral features through a continuous wavelet transform (CWT). The 3D physical property images are derived from cooperative inversion [27,31] of magnetic and DC/IP datasets (Figure 2b–d). A preliminary Fast-ICA separates the latent features inside the multiple physical property images in the form of three negentropymaximized ICs (spatial feature extraction). Then, three ICs are decomposed to form a set of raw spectral features through a continuous wavelet transform (CWT).

We calculated the 2D wavelet coefficients with a gaussian mother wavelet for three scales in 72 directions (every 5◦ ). We used only three scales because we observed that for more than three scales, the changes in low-frequency features (high scales) would disappear, and there will not be any difference in the raw spectral features for scales more than three. The 3D spectral decomposition is iterated for the three initial ICs and has produced 648 (3 × 72 × 3) features containing several frequency-dependent raw features that need further extraction. Figure 6 demonstrates the raw spectral features sliced at an elevation of 1000 m. X and Y UTM coordinates are removed to save space for display. We calculated the 2D wavelet coefficients with a gaussian mother wavelet for three scales in 72 directions (every 5°). We used only three scales because we observed that for more than three scales, the changes in low-frequency features (high scales) would disappear, and there will not be any difference in the raw spectral features for scales more than three. The 3D spectral decomposition is iterated for the three initial ICs and has produced 648 (3 × 72 × 3) features containing several frequency-dependent raw features that need further extraction. Figure 6 demonstrates the raw spectral features sliced at an elevation of 1000 m. X and Y UTM coordinates are removed to save space for display.

**Figure 6.** CWT with a gaussian mother wavelet produced raw 3D spectral features with three scales in 72 directions (0°, 5°, 10°, …, 355°). The extracted features are sliced at an elevation of 1000 m. **Figure 6.** CWT with a gaussian mother wavelet produced raw 3D spectral features with three scales in 72 directions (0◦ , 5◦ , 10◦ , . . . , 355◦ ). The extracted features are sliced at an elevation of 1000 m.

Fast-ICA, through negentropy maximization, separates the 648 raw features to produce the 3D independent spectral inputs necessary for feature selection. The Fast-ICA algorithm can also reduce the dimensionality of the raw features. The best strategy to reduce the dimensionality of features depends on two factors: computer hardware resources and the stability of the SFSS output. Low reductions result in high computational costs, and severe reductions lead to loss of valuable information and poor predictions. There is always an optimum way to moderately reduce the dimensionality of features through several testing and running of the algorithm by the user. In cases with no hardware limitation, this trial-and-error strategy is unnecessary. In this study, we reduced the dimensionality of the raw spectral features by 35 percent to 227 extracted spectral features (Figure 7). Fast-ICA, through negentropy maximization, separates the 648 raw features to produce the 3D independent spectral inputs necessary for feature selection. The Fast-ICA algorithm can also reduce the dimensionality of the raw features. The best strategy to reduce the dimensionality of features depends on two factors: computer hardware resources and the stability of the SFSS output. Low reductions result in high computational costs, and severe reductions lead to loss of valuable information and poor predictions. There is always an optimum way to moderately reduce the dimensionality of features through several testing and running of the algorithm by the user. In cases with no hardware limitation, this trial-and-error strategy is unnecessary. In this study, we reduced the dimensionality of the raw spectral features by 35 percent to 227 extracted spectral features (Figure 7).

**Figure 7.** Total 227 number of separated 3D spectral features after Fast-ICA dimensionality reduction for lithological feature selection. The extracted features are sliced at an elevation of 1000 m. **Figure 7.** Total 227 number of separated 3D spectral features after Fast-ICA dimensionality reduction for lithological feature selection. The extracted features are sliced at an elevation of 1000 m.

We ran the SFSS algorithm two times in parallel to predict lithotypes and Au grades. For feature selection, we used an MLP network with independent spectral features as inputs, 20 neurons in the hidden layer, a maximum of 20 iterations, and one output. Of the borehole data, 70% are used for training and the rest for validation (10%) and testing (20%). The GA is iterated in 20 generations with 20 populations, 70% crossover, 40% of mutation, and a mutation rate of 0.05. We ran the SFSS algorithm two times in parallel to predict lithotypes and Au grades. For feature selection, we used an MLP network with independent spectral features as inputs, 20 neurons in the hidden layer, a maximum of 20 iterations, and one output. Of the borehole data, 70% are used for training and the rest for validation (10%) and testing (20%). The GA is iterated in 20 generations with 20 populations, 70% crossover, 40% of mutation, and a mutation rate of 0.05.

For the lithotypes model reconstruction, evaluation results are compared for both SRL optimized with the Levenberg–Marquardt method and the SFSS algorithm optimized with GA (Figure 8a). As can be seen, the SFSS algorithm gives better validation and test results (98% and 91% R-squares, relatively) compared to the poor results of the conventional SRL method (with 90% and 27% R-squares, relatively). For the lithotypes model reconstruction, evaluation results are compared for both SRL optimized with the Levenberg–Marquardt method and the SFSS algorithm optimized with GA (Figure 8a). As can be seen, the SFSS algorithm gives better validation and test results (98% and 91% R-squares, relatively) compared to the poor results of the conventional SRL method (with 90% and 27% R-squares, relatively).

For Au grade prediction, the same network parameters are used, and the Fast-ICA dimensionality reduction has reduced the number of independent spectral features to 227 (Figure 7). Evaluation results are compared for both SRL optimized with the Levenberg– Marquardt method and the SFSS algorithm optimized with GA (Figure 8b). As can be seen, the SFSS algorithm gives much better validation and test results (with 97% and 88% R-squares, relatively) compared to the poor results of the conventional SRL method (with 95% and 56% R-squares, relatively).

The error landscapes of the lithotype and Au grade prediction by SFSS are also shown in Figure 9. The best results come after the last generation of GA and for the least global error indicated by the dark blue color. For lithotype prediction, the SFSS algorithm selects 200 out of 227 features for spectral learning (Figure 10). As can be seen, the SFSS can recover a much better lithological model compatible with prior geological information. The selected features for lithotype prediction are also shown in Figure 11. In this case, SFSS detects 27 redundant features in the prediction of lithological variations that helps to facilitate the machine learning predictions.

**Figure 8.** Evaluation of lithotype (**a**) and Au (**b**) predictive modeling for SRL and SFSS algorithms. **Figure 8.** Evaluation of lithotype (**a**) and Au (**b**) predictive modeling for SRL and SFSS algorithms. R-sq denotes the R-squared values. *Minerals* **2022**, *12*, x 14 of 19

R-sq denotes the R-squared values.

**Figure 9.** Bi-objective error landscape of SFSS algorithm during lithotype (**a**) and Au (**b**) estimation. **Figure 9.** Bi-objective error landscape of SFSS algorithm during lithotype (**a**) and Au (**b**) estimation.

**Figure 10.** Results of SRL lithotype predictions (**a**) versus SFSS lithotype predictions (**b**,**c**). SFSS resulted in the better reconstruction of lithotypes with a smaller number of input spectral features (Nf

= 200). Results are sliced horizontally at the elevation of 1000 m.

**Figure 10.** Results of SRL lithotype predictions (**a**) versus SFSS lithotype predictions (**b**,**c**). SFSS resulted in the better reconstruction of lithotypes with a smaller number of input spectral features (Nf = 200). Results are sliced horizontally at the elevation of 1000 m. **Figure 10.** Results of SRL lithotype predictions (**a**) versus SFSS lithotype predictions (**b**,**c**). SFSS resulted in the better reconstruction of lithotypes with a smaller number of input spectral features (Nf = 200). Results are sliced horizontally at the elevation of 1000 m.

**Figure 9.** Bi-objective error landscape of SFSS algorithm during lithotype (**a**) and Au (**b**) estimation.

**Figure 11.** A total of 200 selected features were used for lithotype prediction. The selected features are sliced at an elevation of 1000 m. **Figure 11.** A total of 200 selected features were used for lithotype prediction. The selected features are sliced at an elevation of 1000 m.

The error landscape of the Au grade prediction by SFSS is also shown in Figure 10b. For the last generation, the SFSS algorithm selects 209 out of 227 features for spectral learning (Figure 12). The SFSS can recover a much better Au grade model than conventional SRL estimations. The selected features are also shown in Figure 13. In this case, SFSS detects 18 extremely redundant features in the prediction of Au concentrations that help to facilitate the machine learning predictions. The error landscape of the Au grade prediction by SFSS is also shown in Figure 10b. For the last generation, the SFSS algorithm selects 209 out of 227 features for spectral learning (Figure 12). The SFSS can recover a much better Au grade model than conventional SRL estimations. The selected features are also shown in Figure 13. In this case, SFSS detects 18 extremely redundant features in the prediction of Au concentrations that help to facilitate the machine learning predictions.

**Figure 12.** Au grade predictions: Results of 3D SRL predictions (**a**) versus 3D SFSS predictions (**b**). SFSS resulted in the better reconstruction of Au-grade distributions with a smaller number of input spectral features (*Nf*). Results are sliced horizontally at the elevation of 1000 m. **Figure 12.** Au grade predictions: Results of 3D SRL predictions (**a**) versus 3D SFSS predictions (**b**,**c**). SFSS resulted in the better reconstruction of Au-grade distributions with a smaller number of input spectral features (*N<sup>f</sup>* ). Results are sliced horizontally at the elevation of 1000 m.

**Figure 13.** A total of 209 selected features were used for Au prediction. The selected features are sliced at an elevation of 1000 m. **Figure 13.** A total of 209 selected features were used for Au prediction. The selected features are sliced at an elevation of 1000 m.

#### **4. Conclusions 4. Conclusions**

We propose a spectral feature selection algorithm for supervised learning of geological patterns to perform 3D predictive modeling from 3D inverted physical properties. Our work is based on the synergy between the ICA and CWT feature extraction methods and multi-objective machine learning optimization through GA for feature selection. The spectral feature extraction method provided the inputs of the feature selection algorithm, and we show that our self-proposed SFSS algorithm can pick the relevant spectral features necessary for 3D predictive modeling of the targets. The practical implementation of the SFSS algorithm is also evaluated for an epithermal Au/Ag deposit in British Columbia, Canada. The results show that the spectral learning scheme proposed can efficiently learn geological patterns to make predictions based on 3D physical property inputs. The SFSS also minimizes the number of extracted spectral features and tries to pick the best representative geophysical features for each target learning case. This automated dimensionality reduction strategy also gives interpreters a precise predictive model and an understanding of the relevant and irrelevant selected geological features at the end, which adds value to the interpretability of the machine learning process. Although tested on Newton's epithermal deposit, the proposed feature selection approach should be applicable to similar mineral deposits. Other 3D geophysical imageries acquired with inversion of seismic, gravity, and electromagnetic data sets could be considered as well to enhance the accuracy of the feature extraction and the subsequent feature selection. Future research will be focused on the automatic fine-tuning of SFSS hyperparameters and adapting the 3D SFSS algorithm to provide a practical tool for integrated 3D imaging of mineral deposits. We propose a spectral feature selection algorithm for supervised learning of geological patterns to perform 3D predictive modeling from 3D inverted physical properties. Our work is based on the synergy between the ICA and CWT feature extraction methods and multi-objective machine learning optimization through GA for feature selection. The spectral feature extraction method provided the inputs of the feature selection algorithm, and we show that our self-proposed SFSS algorithm can pick the relevant spectral features necessary for 3D predictive modeling of the targets. The practical implementation of the SFSS algorithm is also evaluated for an epithermal Au/Ag deposit in British Columbia, Canada. The results show that the spectral learning scheme proposed can efficiently learn geological patterns to make predictions based on 3D physical property inputs. The SFSS also minimizes the number of extracted spectral features and tries to pick the best representative geophysical features for each target learning case. This automated dimensionality reduction strategy also gives interpreters a precise predictive model and an understanding of the relevant and irrelevant selected geological features at the end, which adds value to the interpretability of the machine learning process. Although tested on Newton's epithermal deposit, the proposed feature selection approach should be applicable to similar mineral deposits. Other 3D geophysical imageries acquired with inversion of seismic, gravity, and electromagnetic data sets could be considered as well to enhance the accuracy of the feature extraction and the subsequent feature selection. Future research will be focused on the automatic fine-tuning of SFSS hyperparameters and adapting the 3D SFSS algorithm to provide a practical tool for integrated 3D imaging of mineral deposits.

**Author Contributions:** Conceptualization, B.A. and L.-Z.C.; methodology, B.A.; software, B.A.; validation, B.A. and L.-Z.C.; formal analysis, B.A.; investigation, B.A.; resources, B.A.; data curation, B.A.; writing—original draft preparation, B.A.; writing—review and editing, B.A., L.-Z.C., M.J., and

**Author Contributions:** Conceptualization, B.A. and L.-Z.C.; methodology, B.A.; software, B.A.; validation, B.A. and L.-Z.C.; formal analysis, B.A.; investigation, B.A.; resources, B.A.; data curation, B.A.; writing—original draft preparation, B.A.; writing—review and editing, B.A., L.-Z.C., M.J., and D.L.; visualization, B.A.; supervision, L.-Z.C.; project administration, L.-Z.C.; funding acquisition, L.-Z.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by grants from the Natural Sciences and Engineering Research Council of Canada (Grant No. 43505) and Fonds de recherche du Québec—Nature et technologies (Grant No. 133896).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Minerals* Editorial Office E-mail: minerals@mdpi.com www.mdpi.com/journal/minerals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-7437-0