Integrative 3D Geological Modeling Derived from SWIR Hyperspectral Imaging Techniques and UAV-Based 3D Model for Carbonate Rocks

: This paper demonstrates an integrative 3D model of short-wave infrared (SWIR) hyperspectral mapping and unmanned aerial vehicle (UAV)-based digital elevation model (DEM) for a carbonate rock outcrop including limestone and dolostone in a ﬁeld condition. The spectral characteristics in the target outcrop showed the limestone well coincided with the reference spectra, while the dolostone did not show clear absorption features compared to the reference spectra, indicating a mixture of clay minerals. The spectral indices based on SWIR hyperspectral images were derived for limestone and dolostone using aluminum hydroxide (AlOH), hydroxide (OH), iron hydroxide (FeOH), magnesium hydroxide (MgOH) and carbonate ion (CO 32 − ) absorption features based on random forest and logistic regression models with an accuracy over 87%. Given that the indices were derived from ﬁeld data with consideration of commonly occurring geological units, the indices have better applicability for real world cases. The integrative 3D geological model developed by co-registration between hyperspectral map and UAV-based DEM using best matching SIFT descriptor pairs showed the 3D rock formations between limestone and dolostone. Moreover, additional geological information of the outcrop was extracted including thickness, slope, rock classiﬁcation, strike, and dip.


Introduction
Sedimentary rocks, one of the three major rock types, are one of the most widely distributed rocks on Earth's surface covering 66% [1]. They are mainly produced in a water environment such as an ocean basin as a sedimentation originated from weathering, erosion, and transportation of continental rocks. The main components describing the sedimentary rocks are lithology, thickness, strike, and dip of bedding planes. Lithology and thickness describe the characteristics of environment and duration of sedimentation. Strike and dip depict how a sedimentary layer is positioned and its inference of the history of Earth's dynamics caused by tectonic events. To figure out those properties, geological field surveys were conducted including visual inspection of rock types, acid reaction test, field tremolite/dolomite map on a UAV 3D model for an open pit mine. Teodoro et al. [15] used a thermal infrared sensor on a UAV for environmental monitoring such as temperature and topography in a coal mine waste pile and the surrounding areas. Kirsch et al. [11] used a UAV SWIR hyperspectral imager to identify hydrothermal alteration minerals (smectite-goethite-jarosite) and developed a 3D mineral distribution model for a quarry. These previous studies have classified specific minerals using hyperspectral mapping for fresh open pit mines without derivation of characteristics of rock layers such as thickness, strike, and dip. The outcrops in an open pit mine and the field outcrops are different where relatively impure lithology occurs in the field for most cases, and the geological models without strike and dip of rock layers have limited information. Moreover, the extendibility of the previous classification models is questionable as the models were specified and trained on their sites. Therefore, the more generalized models for detection of limestone and dolostone are in a field environment with specified equations.
This study introduces a 3D integrative geological model of carbonate rock outcrops using ground-based SWIR hyperspectral images and UAV-based 3D model. An outcrop with alternative limestone and dolostone layers was scanned by a Hyspex SWIR-384 hyperspectral imaging camera. Importance of spectral variables for detection of limestone and dolostone were identified by a spectroscopic analysis and a random forest classification model from the images. Univariable logistic regression models were tested using the spectral indices for classification of lithological layers from the images. The derived lithological image was draped on UAV-based 3D model for a perspective visualization of geological distribution. The geological information including lithology, thickness, strike, and dip as well as topographic characteristics were extracted from the 3D integrative model.

Study Area
The study area, Gumunso is in Taebaek city, Gangwon province (37 • 5 38"N and 129 • 2 29"E), 195 km east from Seoul in South Korea ( Figure 1). The target outcrop has an area of 1881 square meters and exposed as a rocky outcrop along a stream. The outcrop shows an alternative rock formation between limestone and dolostone with a dimension of 22 by 42 m. The outcrop belongs to Makgol formation consisting of dominant lithofacies of dark gray massive dolostone, sometimes intercalated with bioturbated wackestone/grainstone and limestone conglomerate [16][17][18][19]. The Makgol formation is characterized by cyclic sedimentation of meter-scale shallowing-upward units in a peritidal environment in the latter part of the Paleozoic era [20]. Due to the stratigraphic distribution with topographic variation, the study site is ideal for developing a 3D integrative geologic model using hyperspectral imaging analysis and a UAV-based 3D surface model.

Hyperspectral Image Acquisition and Preprocessing
The study outcrop was scanned by a HySpex SWIR-384 SWIR hyperspectral camera in 288 spectral channels, ranging from 930 to 2500 nm (Table 1). For field operations, to ensure stable and reliable acquisition of images on a corrugated surface, a rugged yet portable tripod was used. Barium sulfate white and black reference panels were stationed in the field of view for radiometric calibration. The images were acquired from 14:00 to 17:30 on 28 May 2020, under a natural illumination condition with sufficient incident energy to ensure a continuous spectrum captured by the camera. Then, the integration time and the camera focus were adjusted for preferable image acquisition. After adjusting these presets, the hyperspectral image was referenced to panels followed by outcrop scanning [21]. The target outcrop was scanned 6 times for 2 sections at 60 m distance from the target for maximum image coverage. After radiometric calibration, we selected the two best hyperspectral images (1616 × 620 pixels) for further analysis (Figure 2).

Hyperspectral Image Acquisition and Preprocessing
The study outcrop was scanned by a HySpex SWIR-384 SWIR hyperspectral camera in 288 spectral channels, ranging from 930 to 2500 nm (Table 1). For field operations, to ensure stable and reliable acquisition of images on a corrugated surface, a rugged yet portable tripod was used. Barium sulfate white and black reference panels were stationed in the field of view for radiometric calibration. The images were acquired from 14:00 to 17:30 on 28 May 2020, under a natural illumination condition with sufficient incident energy to ensure a continuous spectrum captured by the camera. Then, the integration time and the camera focus were adjusted for preferable image acquisition. After adjusting these presets, the hyperspectral image was referenced to panels followed by outcrop scanning [21]. The target outcrop was scanned 6 times for 2 sections at 60 m distance from the target for maximum image coverage. After radiometric calibration, we selected the two best hyperspectral images (1616 × 620 pixels) for further analysis ( Figure 2).
The SWIR hyperspectral images were smoothed by the Savitsky-Golay smoothing filter to remove noisy signals [22]. After that, dark subtraction and minimum noise fraction (MNF) transformation were applied to remove random noise and correct atmospheric noise in the images. The MNF transformation is composed of two cascaded principal component transformations (PCT). The first PCT uses the sampled noise covariance to transform noise to white noise (unit variance and no band-to-band correlation). The second PCT transforms the noise-whitened data to two parts: large eigenvalues corresponding to coherent images and near-unity eigenvalues corresponding to noise dominated images. The noise dominated images were treated by low-pass filters to suppress the noise. A backward transform was performed to restore the spectral bands with minimal noise content [23,24]. The noise-reduced hyperspectral images were used for further analyses [25]. The image pro-processing and further imaging analysis were conducted by ENVI (Harris Geospatial, Boulder, CO, USA) and Rstudio (Integrated Development for R. RStudio, PBC, Boston, MA, US).  The SWIR hyperspectral images were smoothed by the Savitsky-Golay smoothing filter to remove noisy signals [22]. After that, dark subtraction and minimum noise fraction (MNF) transformation were applied to remove random noise and correct atmospheric noise in the images. The MNF transformation is composed of two cascaded principal component transformations (PCT). The first PCT uses the sampled noise covariance to transform noise to white noise (unit variance and no band-to-band correlation). The second PCT transforms the noise-whitened data to two parts: large eigenvalues corresponding to coherent images and near-unity eigenvalues corresponding to noise dominated images. The noise dominated images were treated by low-pass filters to suppress the noise. A backward transform was performed to restore the spectral bands with minimal noise content [23,24]. The noise-reduced hyperspectral images were used for further analyses [25]. The image pro-processing and further imaging analysis were conducted by ENVI (Harris Geospatial, Boulder, CO, USA) and Rstudio (Integrated Development for R. RStudio, PBC, Boston, MA, USA).

UAV Survey and DEM Processing
This study selected seven ground control points (GCP hereafter) over the UAV survey site for geometric calibration and accuracy assessment of the output orthorectified image and DEM. To avoid bias in geometric correction, the positions of GCPs were uniformly spread over the study area ( Figure 3). The location of GCPs were measured by Trimble R2 GNSS Receiver based on the virtual reference station (VRS) system with max horizontal accuracy of 10 mm and vertical accuracy of 20 mm.
This study used DJI's Phantom 4 Pro+ UAV to obtain an orthorectified image and DEM. The UAV data was acquired simultaneously with hyperspectral scanning in a sunny weather condition with no cloud cover. Over 100 UAV images were acquired at an altitude of 40 m with 65% overlap at image resolution of 4000 × 3000 pixels. Point clouds were extracted from the UAV images by the SfM algorithm using Agisoft PhotoScan software [26,27]. The UAV data processing for construction of orthorectified image and digital elevation model (DEM) was carried out by the procedure in Figure 4. This study selected seven ground control points (GCP hereafter) over the UAV survey site for geometric calibration and accuracy assessment of the output orthorectified image and DEM. To avoid bias in geometric correction, the positions of GCPs were uniformly spread over the study area ( Figure 3). The location of GCPs were measured by Trimble R2 GNSS Receiver based on the virtual reference station (VRS) system with max horizontal accuracy of 10 mm and vertical accuracy of 20 mm. This study used DJI's Phantom 4 Pro+ UAV to obtain an orthorectified image and DEM. The UAV data was acquired simultaneously with hyperspectral scanning in a sunny weather condition with no cloud cover. Over 100 UAV images were acquired at an altitude of 40 m with 65% overlap at image resolution of 4000 × 3000 pixels. Point clouds were extracted from the UAV images by the SfM algorithm using Agisoft PhotoScan software [26,27]. The UAV data processing for construction of orthorectified image and digital elevation model (DEM) was carried out by the procedure in Figure 4.

Ground Truthing of the Outcrop
Ground truth data were obtained by field inspections and air-photo delineation in a GIS environment. The field survey was conducted at the same date of the hyperspectral and UAV data acquisition. Each rock layer was inspected by a field geologist about lithology, and the geometric properties including thickness, strike and dip of each rock layer were measured. Microscopic observation of rock layers was additionally carried out for

Ground Truthing of the Outcrop
Ground truth data were obtained by field inspections and air-photo delineation in a GIS environment. The field survey was conducted at the same date of the hyperspectral and UAV data acquisition. Each rock layer was inspected by a field geologist about lithology, and the geometric properties including thickness, strike and dip of each rock layer were measured. Microscopic observation of rock layers was additionally carried out for definition of rock layers composing the out crop. The field survey map was, then, reinterpreted to the UAV air-photo to be used as ground truth data to train and validate the hyperspectral image classification (

UAV-Based Orthorectified Image and Digital Elevation Model of the Outcrop
The UAV survey on the target outcrop constructed an orthorectified i tial resolution of 1.7 cm with 7360 × 8448 pixels. The accuracy of orthore sessed based on the GPS survey results of 7 GCPs. The overall RMSE of the image was 21.1 cm. The orthorectified image revealed that the target outcro ric properties of 228 m and 2175 m 2 for perimeter and area, respectively (F tably, the outcrop showed an alternative rock formation between two dif lithology with colors ranging from white gray to dark gray. The DEM derived from the UAV survey had a spatial resolution of 2990 pixels. In general, the overall horizontal and vertical RMSEs of the DE 21.1 and 37.4 cm. The derived DEM revealed that the target outcrop (ma square) is topographically positioned at a minimum elevation of 526.3 m

Spectral Analysis
In the shortwave infrared wavelength region, many mineral groups have distinct absorption features caused by their chemical component. For instance, the spectral characteristics of carbonate rocks such as limestone and dolostone have distinctive absorption features at 2340-2345 and 2320-2328 nm due to CO 3 2− ion [28]. These spectral features can be used to distinguish calcite, dolomite, and other non-carbonate minerals [29]. The absorption feature at 2300-2350 nm is the most noticeable for carbonate rocks when using ground-based Hyperspectral Camera under natural illumination [30].
This study analyzed the spectral characteristics of rock formations in the target outcrop based on pre-processed reflectance spectra and hull quotient corrected spectra extracted from the hyperspectral image pixels. The hull quotient correction technique maximizes and characterizes absorption features, and, thus, is useful for detecting the position and absorption depth characteristics [31].

Band Selection and Spectral Index Derivation
The random forest (RF) classification algorithm was applied to the preprocessed hyperspectral images to select the most effective bands to distinguish different carbonate rock types. Although RF can provide accurate classification results, the extendibility of such model is questionable: the model trained in one study site cannot be used in others. Therefore, the purpose of using random forest is to use the Gini index provided by the RF model as a guidance to build a more straightforward logistic regression model with better extendibility. The training and validation pixels were randomly selected from the hyperspectral image based on the field survey data. In total, 70% of the extracted image pixels (24,198 pixels) were used as training data for RF model development, and the remaining 30% of pixels (10,370 pixels) were used for validation of the classification model. Moreover, 482 mineral and 76 rock spectra from USGS and JPL spectral libraries were included for RF and spectral indices development to identify spectral variables effective to distinguish limestone and dolostone from commonly found rocks and minerals. The mineral spectra include rock forming minerals, a total of 482 mineral spectra and 76 rock spectra. The library mineral spectra include hydrothermal alteration minerals and ore minerals, which are considered the most possible mineral occurrences including alunite, beryl, biotite, calcite, dolomite, chlorite, epidote, fluorite, gypsum, kaolinite, muscovite, pyrophyllite, quartz, tourmaline, and so on. The library rock spectra also include common rock occurrences in the field such as pyroxenite, gneiss, gabbro, sandstone, marble, basalt, limestone, dolostone, shale, granite, diorite, pegmatite, diabase and so on.
RF is a classification algorithm based on decision trees and improved bagging and bootstrapping techniques [32]. The bagging and bootstrapping selected 2/3 samples for training a tree each time and reserved 1/3 samples as out of bag (OOB) validation. Trees were generated by repeating randomized formation of sample and variable sets using the bagging and bootstrapping strategy [33][34][35]. Classification decisions are made by major votes from the trees. The two parameters to be specified in RF classification algorithm are ntree (number of trees to grow) and mtry (number of variables to divide at each node). The RF algorithm selects important bands based on the variable importance by the mean decrease Gini (data purity). The lower the Gini index, the more important the variable in the classification model [36]. We analyzed the spectral bands most useful for the rock classification and used them for further model development. The RF classification was assessed with the overall accuracy, kappa coefficient, commission and omission errors compared to the field survey map [37].
To derive the indices for detection of carbonate rocks including limestone and dolostone for field application, this study used binary logistic regression (BLR). The BLR is composed of one or multiple independent variables that describe a relationship to a dependent response variable [38,39]. In this study, the limestone or dolostone occurrence (variable Y) is regulated by the spectral value of specific bands (variable X) and used to predict an event occurring with only 2 numbers indicating limestone (1) and non-limestone (0) and dolostone (1) and non-dolostone (0). We derived a binary logistic regression equation with probability value of 0 or 1 based on the standard cutoff of 0.5. The equation can be written as a function of P as follows (1) [40][41][42].
where P limestone or dolostone is the probability of limestone or dolostone occurrence, C is the constant (or intercept) of the equation, β is the coefficient (or beta weight) of the predictor variables and x is the reflectance value of the selected band n . Evaluation of the indices was based on the results of Hosmer and Lemeshow test, the Wald statistic, and pseudo-R 2 values of "Nagelkerke" and "Cox and Snell". The Hosmer-Lemeshow test is a Goodness-of-Fit test for logistic regression, especially for models of risk prediction. This test divides the samples to deciles (10-fold) and assesses the model performance in each subgroup by the probability of matching between observed and Furthermore, Cox-Snell and Nagelkerke R 2 statistics were used to evaluate the BLR models, following Equations (2) and (3) below [43][44][45].
where L(M Intercrpt ) is the log-likelihood of the null model and L(M Full ) is the log-likelihood of the full model. The ratio of the pseudo-R 2 closer to 1 indicates the model has more potency [46]. Finally, the best equation for carbonate rock prediction was selected with the highest prediction accuracy.

Fusion of Hyperspectral Imaging and UAV-Based 3D Model
The integrative 3D geologic model for carbonate rocks was constructed using Autodesk 3ds Max (Autodesk, Inc., San Rafael, CA, USA) software based on the scale-invariant feature transform (SIFT hereafter) algorithm to find key-points between two images [47].
The key-points are conjugate points easily recognizable between the classification map and the orthorectified image where geometric properties including coordinate, scale, and orientation were extracted for each key-point. The SIFT algorithm analyzed geometric properties of each key-point based on a 128-dimensional eigenvector (SIFT descriptor) and extracted SIFT descriptors for orthorectified image and the rock type map. Finally, the SIFT algorithm compared between SIFT descriptors of the orthorectified mosaic image and the rock type map to identify the best matching coefficients for image co-registration [48,49]. Once the coordinates are matched, the classified rock type image is draped on the 3D surface of the rock to compute the geometric measures of the stratifications such as thickness, strike, dip, slope, and aspect.

UAV-Based Orthorectified Image and Digital Elevation Model of the Outcrop
The UAV survey on the target outcrop constructed an orthorectified image at a spatial resolution of 1.7 cm with 7360 × 8448 pixels. The accuracy of orthorectified was assessed based on the GPS survey results of 7 GCPs. The overall RMSE of the orthorectified image was 21.1 cm. The orthorectified image revealed that the target outcrop has geometric properties of 228 m and 2175 m 2 for perimeter and area, respectively (Figure 5a). Notably, the outcrop showed an alternative rock formation between two different types of lithology with colors ranging from white gray to dark gray.
The DEM derived from the UAV survey had a spatial resolution of 9 cm at 2828 × 2990 pixels. In general, the overall horizontal and vertical RMSEs of the DEM image were 21.

Spectral Characteristics of Limestone and Dolostone
The spectral characteristics of carbonate rocks forming the target outcrop were analyzed based on reflectance and hull quotient spectra. From the previous studies, the average surface reflectance of carbonate rocks in the outcrop when using ground-based Hyperspectral Camera under natural illumination can observe from 2000 to 2400 nm wavelength range, because other spectral regions are blurred and degraded by strong atmospheric scattering [5,9,28]. The occurrence of the carbonate rocks in the study area could be recognized from the weak absorption features at 2140 nm for the dolostone and 2150 nm for the limestone manifested by CO3 2− ion ( Figure 6) [50]. The limestone shows additional absorption features around 1990, 2150, 2295, 2340 and 2480 nm, where the strongest absorption features at 2340 nm and minor absorptions at 1990, 2150, and 2480 nm are manifested by CO3 2− ion of calcite and that at 2295 nm is associated with FeOH ( Figure 6). The spectral characteristics of the target limestone coincides well with the reference spectra of limestone and calcite ( Figure 6) except the FeOH absorption. Given that FeOH absorption is mainly associated with clay minerals, it indicates that the limestone of the target outcrop is mainly composed of calcium carbonate mineral (calcite) with minor content of clay minerals. The hull quotient corrected spectra of limestone in the target outcrop shows strong absorption features at 2340 nm followed by 2480 nm where the absorption features were identical to the library references of limestone (Figure 6b) confirming the calcite as major mineral component.
The absorption features of dolostone in the target outcrop are at 1990, 2140, 2260, 2320 and 2470 nm, where absorptions at 1990, 2140, 2320 and 2470 nm are manifested by Mg and CO3 2− component, and that at 2260 nm is associated with FeOH ( Figure 6). The distinctive differences between the limestone and dolostone in the target outcrop can be detected at their absorption features. Limestone in the target outcrop shows strong absorption at 2340 nm while dolostone shows weaker absorption at 2320 nm which shifted toward shorter wavelengths. This phenomenon is a characteristic spectral behavior between calcite and limestone reported by previous studies [5,[9][10][11]50,51]. The spectral characteristics of the target dolostone coincides well with the reference spectra of dolomite and dolostone (Figure 6), while the absorption depth of target dolostone is relatively weaker. It indicates that the target dolostone is not a pure dolostone but contains clay minerals as accessary minerals [52].

Spectral Characteristics of Limestone and Dolostone
The spectral characteristics of carbonate rocks forming the target outcrop were analyzed based on reflectance and hull quotient spectra. From the previous studies, the average surface reflectance of carbonate rocks in the outcrop when using ground-based Hyperspectral Camera under natural illumination can observe from 2000 to 2400 nm wavelength range, because other spectral regions are blurred and degraded by strong atmospheric scattering [5,9,28]. The occurrence of the carbonate rocks in the study area could be recognized from the weak absorption features at 2140 nm for the dolostone and 2150 nm for the limestone manifested by CO 3 2− ion ( Figure 6) [50]. The limestone shows additional absorption features around 1990, 2150, 2295, 2340 and 2480 nm, where the strongest absorption features at 2340 nm and minor absorptions at 1990, 2150, and 2480 nm are manifested by CO 3 2− ion of calcite and that at 2295 nm is associated with FeOH ( Figure 6). The spectral characteristics of the target limestone coincides well with the reference spectra of limestone and calcite ( Figure 6) except the FeOH absorption. Given that FeOH absorption is mainly associated with clay minerals, it indicates that the limestone of the target outcrop is mainly composed of calcium carbonate mineral (calcite) with minor content of clay minerals. The hull quotient corrected spectra of limestone in the target outcrop shows strong absorption features at 2340 nm followed by 2480 nm where the absorption features were identical to the library references of limestone (Figure 6b) confirming the calcite as major mineral component.
The absorption features of dolostone in the target outcrop are at 1990, 2140, 2260, 2320 and 2470 nm, where absorptions at 1990, 2140, 2320 and 2470 nm are manifested by Mg and CO 3 2− component, and that at 2260 nm is associated with FeOH ( Figure 6). The distinctive differences between the limestone and dolostone in the target outcrop can be detected at their absorption features. Limestone in the target outcrop shows strong absorption at 2340 nm while dolostone shows weaker absorption at 2320 nm which shifted toward shorter wavelengths. This phenomenon is a characteristic spectral behavior between calcite and limestone reported by previous studies [5,[9][10][11]50,51]. The spectral characteristics of the target dolostone coincides well with the reference spectra of dolomite and dolostone (Figure 6), while the absorption depth of target dolostone is relatively weaker. It indicates that the target dolostone is not a pure dolostone but contains clay minerals as accessary minerals [52].

Carbonate Rock Classification by Random Forest Classification
A total of 34,569 pixels of the hyperspectral images were extracted from limestone and dolostone in the target outcrop for training and validation data of the RF and BLR algorithms. For each classification algorithm, the dataset was randomly partitioned into 70% to train the model and the remaining 30% for validation. Based on the training data,

Carbonate Rock Classification by Random Forest Classification
A total of 34,569 pixels of the hyperspectral images were extracted from limestone and dolostone in the target outcrop for training and validation data of the RF and BLR algorithms. For each classification algorithm, the dataset was randomly partitioned into 70% to train the model and the remaining 30% for validation. Based on the training data, the classification map of limestone and dolostone was derived by RF classification from the SWIR hyperspectral images. In this study, we set the ntree value to 250 trees and mtry to 400 based on the OOB error analysis (Figure 7).
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 18 the SWIR hyperspectral images. In this study, we set the ntree value to 250 trees and mtry to 400 based on the OOB error analysis (Figure 7). The accuracy assessment of the classification on limestone and dolostone showed high overall accuracy and kappa coefficient of 91.71% and 0.89, respectively ( Table 2). The RF algorithm showed good performance for carbonate rock classification in the target outcrop with high overall accuracy, kappa coefficient, and low commission and omission errors for both training and validation data ( Table 2).  The accuracy assessment of the classification on limestone and dolostone showed high overall accuracy and kappa coefficient of 91.71% and 0.89, respectively ( Table 2). The RF algorithm showed good performance for carbonate rock classification in the target outcrop with high overall accuracy, kappa coefficient, and low commission and omission errors for both training and validation data ( Table 2).

Band Selection and Derivation of Carbonate Rock Indices from Binary Logistic Regression
The variable importance plot shows the most important variables marked at the wavelength from 2300 to 2362 nm followed by 2450-2495 nm, 970-1026 nm and 1548-1656 nm (Figure 8). The first tier is the major absorption features of limestone and dolostone from 2320 to 2345 nm manifested by MgOH and/or CO 3 2− . In particular, the peak importance occurs at 2336 nm where absorption features of 70% pure calcite are found [9], followed by a second peak at 2341 nm manifested of CO 3 2− ion. The third peak is at 2331 nm related to MgOH and/or CO 3 2− . The second tier is the absorption features of limestone and dolostone from 2450 to 2495 nm. The third and the last tiers are from 970 to 1026 nm and 1548 to 1656 nm related to absorption features of clay minerals. Through the RF's variable importance index, 30 bands were selected because they also gave the lowest OOB errors [11]. These 30 bands were used to construct the BLR models.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 18 from 2320 to 2345 nm manifested by MgOH and/or CO3 2− . In particular, the peak importance occurs at 2336 nm where absorption features of 70% pure calcite are found [9], followed by a second peak at 2341 nm manifested of CO3 2− ion. The third peak is at 2331 nm related to MgOH and/or CO3 2− . The second tier is the absorption features of limestone and dolostone from 2450 to 2495 nm. The third and the last tiers are from 970 to 1026 nm and 1548 to 1656 nm related to absorption features of clay minerals. Through the RF's variable importance index, 30 bands were selected because they also gave the lowest OOB errors [11]. These 30 bands were used to construct the BLR models. The spectral indices for detection of carbonate rocks derived by BLR model showed that the best detection model of limestone was the one using 15 spectral bands associated with OH, AlOH, FeOH, MgOH and/or CO3 2− absorption features, among which the MgOH and/or CO3 2− absorptions are associated with eight bands (Table 3). Different from limestone, the dolostone model was developed based on nine spectral bands associated with three bands at MgOH and/or CO3 2− absorptions and six bands of clay minerals (Table  3).  The spectral indices for detection of carbonate rocks derived by BLR model showed that the best detection model of limestone was the one using 15 spectral bands associated with OH, AlOH, FeOH, MgOH and/or CO 3 2− absorption features, among which the MgOH and/or CO 3 2− absorptions are associated with eight bands (Table 3). Different from limestone, the dolostone model was developed based on nine spectral bands associated with three bands at MgOH and/or CO 3 2− absorptions and six bands of clay minerals (Table 3).
Both models are statistically significant in the Pseudo-R 2 goodness of fit test 0.734 and 0.979 for limestone, and 0.712 and 0.956 for dolostone. Generally, Cox and Snell pseudo-R 2 values higher than 0.2 are a meaningful fit [53]. The results of the Hosmer and Lemeshow test showed that p values of the BLR model were 0.662 for limestone and 0.07 for dolostone. In the Hosmer and Lemeshow test, the p value higher than 0.05 cannot reject the null hypothesis, which indicates the model is acceptable [54] (Table 4).
Accuracy assessment of validation data showed an overall accuracy of 87.91% and a kappa coefficient of 0.77 (Table 5). It indicates that these models can distinguish limestone and dolostone with high accuracy and acceptable statistical significance. Compared to the RF model, although the accuracy of the BLR models is lower than the RF model, they are more generalizable because there are only several coefficients to be adjusted rather than the complex decision making structure of the RF model. Furthermore, the detection model in this study also included various rocks and mineral spectra commonly found in geological units. Therefore, the detection model derived by BLR is easier to interpret for geological survey.   , 2362 and 2480 nm and nine spectral bands for dolostone detection at  1968, 2009, 2024, 2055, 2091, 2132, 2300, 2336 and 2341 nm based on MgOH, CO 3 2− and clay mineral absorption features (Table 3 and Figure 6). The model accuracy using field data is as good as laboratory samples, which is a generalization example of laboratory experiments to real-world cases.

Fusion of Classification Map and UAV-Based 3D Model
The co-registration between HIS and UAV images ensures the match up of the BLR result with the rasterized DEM from the UAV image point cloud. The integration between BLR image and the DEM was completed by draping the BLR image on the DEM. The geometric properties were then computed from the DEM. As lithology is extracted from hyperspectral images, the 3D model shows the position of the extracted sedimentary rock layers as well as bedding planes (Figure 9). The 3D coordinates (x, y, z) of the elements in the rock layer map can be used to calculate 3D geometric properties of the layers such as strike and dip, where strike is the azimuthal direction of a bedding plane, and dip is the slope angle of it ( Figure 9). The rock layer thickness ranges from 0.4 to 4.7 m, which matches with field survey results. Furthermore, our 3D model revealed that the strike and dip of the rock layers are N85 • E and 52 • NW ( Figure 9). The geological and topographic information extracted from the model was verified by the field survey data. The results suggest that our approach using HIS and UAV data can obtain reliable measurements of rock layer thickness, strike, and dip for understanding geological positioning of different rock layers, geological profile, and geological relationship between rock layers without contact surveys.  (Table 3 and Figure 6). The model accuracy using field data is as good as laboratory samples, which is a generalization example of laboratory experiments to real-world cases.

Fusion of Classification Map and UAV-Based 3D Model
The co-registration between HIS and UAV images ensures the match up of the BLR result with the rasterized DEM from the UAV image point cloud. The integration between BLR image and the DEM was completed by draping the BLR image on the DEM. The geometric properties were then computed from the DEM. As lithology is extracted from hyperspectral images, the 3D model shows the position of the extracted sedimentary rock layers as well as bedding planes (Figure 9). The 3D coordinates (x, y, z) of the elements in the rock layer map can be used to calculate 3D geometric properties of the layers such as strike and dip, where strike is the azimuthal direction of a bedding plane, and dip is the slope angle of it ( Figure 9). The rock layer thickness ranges from 0.4 to 4.7 m, which matches with field survey results. Furthermore, our 3D model revealed that the strike and dip of the rock layers are N85° E and 52° NW (Figure 9). The geological and topographic information extracted from the model was verified by the field survey data. The results suggest that our approach using HIS and UAV data can obtain reliable measurements of rock layer thickness, strike, and dip for understanding geological positioning of different rock layers, geological profile, and geological relationship between rock layers without contact surveys.  One limitation found by this case study is that because of the roughness of the rock, it was difficult to find an ideal illumination condition to cover the entire rock body during one time of survey ( Figure 9). There were always shadows on the rock surface, which Remote Sens. 2021, 13, 3037 16 of 18 must be removed from the classification. The 3D point cloud was also degraded by the shadows because the SfM algorithm relies on image matching between features excluding shadows. Therefore, multiple trips are needed to obtain a shadow-free, full coverage of the rock outcrop, and there would be another problem of dealing with the changing illumination conditions.

Conclusions
This work introduced an integrative 3D model between hyperspectral mapping and a UAV-based 3D model for a carbonate rock outcrop. A carbonate rock outcrop composed of alternative layers of limestone and dolostone was scanned with a SWIR hyperspectral imaging system and classified with an RF algorithm to derive spectral bands sensitive to mapping of carbonate rocks for constructing logistic regression models. A DEM and ortho-rectified images were generated by the SfM algorithm from UAV images. The hyperspectral classification map was co-registered with the ortho-rectified image using the best-matching SIFT descriptor pairs as control points, and the integrative 3D geologic model was developed.
The spectral characteristics of limestone and dolostone showed variations in MgOH and CO 3 2− absorption features. The spectral analysis revealed that the target limestone is relatively purer than dolostone as dolostone has weaker absorption of dolomite with weak clay absorptions. The RF classification achieved an overall accuracy of 91.71% and a kappa coefficient of 0.89. As an outcome, 30 bands with the highest Gini index from the RF model were extracted from OH, AlOH, FeOH, MgOH and CO 3 2− absorption features. The logistic regressions using the spectral indices had an overall accuracy of 87.73% for limestone and 88.02% for dolostone. Given that the spectral indices were derived from field data with consideration of commonly found geological materials, the indices have better applicability in the real-world cases than laboratory samples.
Using the 3D mineral map produced by draping the mineral classification result on the UAV digital surface model, detailed measurement of lithologic units such as rock type, thickness, strike, and dip were extracted from the model. Therefore, the method introduced in this study opens a fast, reliable, and convenient option for geological surveys, and provides a fast, convenient solution to expediting geological surveys. Furthermore, the model can be easily extended for a subsurface geological model and mineral resources assessment.