*Article* **The Use of Machine Learning and Satellite Imagery to Detect Roman Fortified Sites: The Case Study of Blad Talh (Tunisia Section)**

**Nabil Bachagha 1,\*, Abdelrazek Elnashar 2, Moussa Tababi 3, Fatma Souei <sup>4</sup> and Wenbin Xu <sup>1</sup>**


**Abstract:** This study focuses on an ad hoc machine-learning method for locating archaeological sites in arid environments. Pleiades (P1B) were uploaded to the cloud asset of the Google Earth Engine (GEE) environment because they are not yet available on the platform. The average of the SAR data was combined with the P1B image in the selected study area called Blad Talh at Gafsa, which is located in southern Tunisia. This pre-desert region has long been investigated as an important area of Roman civilization (106 BCE). The results show an accurate probability map with an overall accuracy and Kappa coefficient of 0.93 and 0.91, respectively, when validated with field survey data. The results of this research demonstrate, from the perspective of archaeologists, the capability of satellite data and machine learning to discover buried archaeological sites. This work shows that the area presents more archaeological sites, which has major implications for understanding the archaeological significance of the region. Remote sensing combined with machine learning algorithms provides an effective way to augment archaeological surveys and detect new cultural deposits.

**Keywords:** archaeological; machine learning; Google Earth Engine; remote sensing; fortified sites

#### **1. Introduction**

Similar to most research-driven fields, archaeology fundamentally depends on the following two precious and scarce resources: time and money [1]. Archaeologists often travel long distances to reach areas of interest and devote a considerable amount of time to excavations and surveys. Additionally, the discovery of archaeological sites is among the most time-consuming and labour-intensive activities. Archaeologists often use advanced technologies to search for less expensive and faster methodologies for archaeological research. Southern Tunisia is a vast region with difficult access to the investigation area and limited opportunities for in-person data collection. Consequently, land managers spend precious and dwindling resources conducting expensive surveys that result in very few representative samples. One way to address this problem is to develop survey strategies that focus on archaeological potential. Satellite image analysis is a relatively low-cost method with great potential for addressing these needs. The application of remote sensing technology is an effective method for producing relatively complete records of archaeological settlement patterns and human activity at the landscape scale. The literature concerning archaeological remote sensing (RS) has shown that multispectral satellite sensors and airborne LiDAR are particularly useful for addressing survey coverage and positioning issues [2–7]. In many regions worldwide, RS has been used to detect and map archaeological proxy indicators [8–12] by leveraging multisource imagery and spatial pattern analyses. RS technologies enable the discovery of new sites and maps

**Citation:** Bachagha, N.; Elnashar, A.; Tababi, M.; Souei, F.; Xu, W. The Use of Machine Learning and Satellite Imagery to Detect Roman Fortified Sites: The Case Study of Blad Talh (Tunisia Section). *Appl. Sci.* **2023**, *13*, 2613. https://doi.org/10.3390/ app13042613

Academic Editor: Tung-Ching Su

Received: 20 January 2023 Revised: 10 February 2023 Accepted: 12 February 2023 Published: 17 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of ancient remnants and can be used to assist with the digital reconstruction of ancient monuments and their environmental backgrounds. For example, in Mediterranean regions, such as Greece and Croatia, several studies have demonstrated excellent performance using a variety of passive and active remote sensing data (including multispectral imagery, thermal images, and UAV-based LiDAR data) to identify ancient Roman remains [13–17]. Fonte et al. [18] evaluated and debated the use of digital modelling tools in northwestern Iberia to identify and define Roman roads through GIS-based spatial analyses.

In recent years, RS-based progressive archaeological studies have steadily integrated machine-learning algorithms that enable automated feature finding [19]. Most applications of these methods have focused on small-scale feature detection using LiDAR [20–22] or WorldView satellite data [8,23,24]. In the east, Menze and Ur [25] applied a random forest (RF) classifier to a multitemporal collection of ASTER imagery to classify probable anthrosols. In the Cholistan Desert, Hector et al. [19] applied an RF classifier that integrates a multisensor, multitemporal machine learning approach to detecting archaeological mounds, and Abdulmalik et al. [26,27] used an artificial neural network model combined with a geographic information system (GIS) and remote sensing for surface water detection. Other studies have used multiple types of information to identify archaeological remains or threats to archaeological heritage, such as urban sprawl and archaeological looting [28–30]. Combined with machine learning algorithms, remote sensing provides an effective way to expand investigation areas and detect new cultural deposits [31–34]. In particular, the application of this technology allows researchers to (1) investigate remote areas that are not easy to access due to harsh geography and (2) survey vast and complex environments (characterised, for example, by otherwise challenging topography). Platforms, such as Google Earth Engine (GEE) [35], provide a scalable, cloud-based, geospatial retrieval and processing platform, particularly in remote areas with little or almost no ground information. Google Earth Engine is a cloud-based platform that makes it simple to access high-performance computing resources for covering geospatial datasets. Once an algorithm is developed on GEE, users can produce regular data products or deploy interactive applications powered by Google Earth Engine resources. Archaeologists have only recently begun to search for more complex models, drawing on innovative analysis of spatial and statistical computing and machine learning methods [36,37]. To date, these studies are increasing when dealing with the detection of archaeological sites based on high-resolution satellite or drone imagery and in pottery classification [38–40]. Several studies have applied different methods for the identification of buried structures [5,41]. For instance, in arid or desert environments, satellite remote data provides very reliable results that can be used to analyse archaeological sites and their environments. According to previous research, there are more Roman civilisation sites in this area [8], but only a few have been uncovered by remote sensing methods and identified the major distinguishing features of Roman settlements, such as villa and tower walls. In this context, the discovery of archaeological features using machine learning classification algorithms has emerged as one of the most important topics regarding remote sensing; it is an important method to better understand and reconstruct ancient landscapes. Thus far, these topics have mainly been studied by experts, sometimes with the help of automated tools, which results in a very costly and time-consuming process. Therefore, in recent years, many automation methods combining satellite image processing with machine learning models have been developed [19,42–45].

Due to the difficult geography, the investigation of ancient sites in this region is not easy; therefore, machine learning via satellite can play an important role in the detection and documentation of archaeological sites. The aim of this study is to evaluate the potential offered by the combined use of remote sensing data and machine learning to discover fortified sites in arid areas. An innovative aspect of our approach is the combination of very high spatial resolution satellite (PB1) technology with SAR data for better detection and to help obtain enhanced results.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area is positioned in Blad Talh (see Figure 1). The weather is hot, and the average temperature is 19.6 ◦C; the location receives less than 100 mm of precipitation annually [46,47]. The climate is characterised by extremely irregular rainfall that varies across seasons and within the same season. The alternation of rainy and dry years occurs without a defined pattern, and periods of prolonged drought often occur. The location of Blad Talh (between 34◦26 8.18 N and 34◦21 36.52 N and between 9◦20 3.93 E and 9◦30 56.12 E) is shown in Figure 1. It is formed by a series of high plains that are separated by mountains and gradually descend towards the south. Blad Talh has a maximum height of 1165 metres and is characterised by a unique morphology. Blad Talh has been the subject of several studies, mostly geological in nature, that assessed the tectonic and structural formation of two mountain ranges, the Orbita-Bouhedma and Chemsi-Belkhir. Despite its geographical position, located less than 40 km from the Mediterranean and more than 100 km north of the Sahara, Blad Talh is part of pre-Saharan Tunisia, bordered to the north by the Orbita-Bouhedma link, and is classified in the arid bioclimatic stage. This region played an important historical role in the Roman era and was part of the sensitive defence system (Limes) located along the desert border, as attested by the several forts that provided protection [8,46–49]. Many sites that appear to be military structures—forts or towers of various sizes—have been identified at these locations based on their layouts and the methods used to construct them. The study area, Blad Talh, was a region of great importance in Roman times. In a report by Toussaint on the campaign of 1902–1903 [50], we found that the fortification of KSAR GROUECH was especially important. The general assessment of this report leads to the conclusion that the Roman ruins in the region of Blad Talh are very numerous in the plains and much rarer and even less important in the mountainous regions.

**Figure 1.** (**A**) Study area location in Tunisia (North Africa); (**B**) Archaeological sites in the Roman Fort area; (**C**) Fort of KSAR GRAOUECH surrounding wall and remains.

#### *2.2. Methodology*

The study describes the process of archaeological site detection and the development of an integrated approach that includes remote sensing, a geographic information system, and random forest to achieve these objectives. First, we selected Blad Talh for this research because of the authenticated archaeological remains and the opportunity to gather survey information for this project in November 2021. In total, 101 known archaeological sites are located in Sebkhet Nouail (Figure 1B). Based on a field survey of the entire study area, the archaeological sites in this area were classified. After the field sample collection was carried out by GNSS, satellite imagery was used to locate the archaeological sites at each selected location using ArcGIS based on the actual locations identified by the GNSS instrument. The second step is to use very high resolution (VHR) data to delineate and map study area sites, as well as Sentinel-1 data, which is used as input to help achieve better results. For geometric corrections, these images were imported into the ArcGIS 10.4 package using WGS 84 and UTM zone 32. In the third step, a random forest was applied using a combination of VHR and SAR data in this study, while employing Google Earth Engine as a platform to decrease the computational time and avoid any other barriers. These sensors were a crucial factor in ensuring that the features that would generate values analogous to those of spots in one source were distinguished from the other features. The final step was survey validation and determining archaeological site significance using remote sensing-based archaeological surveys. Remote sensing, machine learning, and GNSS allow us to identify unknown sites and provide new directions for future archaeological research in our study area in the ancient Roman period. The overall methodology is presented in Figure 2, while the specifics of each major section are given in the subsequent sections.

#### **Figure 2.** Framework approach.

#### *2.3. Field Investigation Sample Collection*

We selected Blad Talh for this research because of the authenticated archaeological remains and the opportunity to gather survey information for this project. The database recaptured from the historic record of the region allowed an overview of the literature review of the Roman period. The GIS integration of the results of these different steps was critical to planning a systematic field survey focused on specific areas of interest, for which the primary results are presented. In November 2021, the search team arrived at the study area of interest. A number of ceramic African sigils were found at these sites. Along these roads, we have classified these sites of interest as being Roman fortified structures based on the ancient remains observed at these sites. By coordinating our findings with GIS, remote sensing, and machine learning, the importance of these new fortifications allowed

speculation regarding their ancient conditions. By confirming our results with remote sensing images and machine learning, we aimed to study the spatial distributions of the archaeological remains and reconstruct the Roman settlement in the region.

#### *2.4. Data Preprocessing*

#### 2.4.1. Pleiades

These images were obtained from the Pleiades sensor (P1B) in a strip of the Blad Talh region with the simultaneous acquisition of panchromatic and multispectral information in the Universal Transverse Mercator (UTM) planar projection, Zone 32 N, and WGS 84 datum on 16 July 2021. Of course, the two sets of bands had two different geometric resolutions (Table 1). While the multispectral bands had a ground resolution of 2.14 m, we used 2.14 \* 2.14 m/px. The acquisition in the multispectral mode was subdivided into the following four bands of the spectrum: NIR, red, green, and blue. In the first stage of the analysis, P1B images were used to study the archaeological landscape [51] and highlight the most important areas where the cartographic features were determined (boundaries, forts, sites, etc.).

**Table 1.** Pleiades data from Blad Talh acquired on 16 July 2021.


#### 2.4.2. Sentinel-1

Sentinel-1 (S1) in the band C imagery contains a 12- or 6-day revisit cycle depending on the availability of images 1A and 1B with 10 m spatial resolution. The analysis presented here selected the interferometric wide swath mode, which is the main operational mode and produces data in single (HH or VV) or double (HH + HV or VV + VH) polarisation in both ascending and descending modes. Each scene available at GEE was preprocessed using the Sentinel-1 Toolbox to (1) remove low-intensity noise and invalid data from scene edges and (2) remove thermal noise [19]. The edge of S1 was omitted using the 'connected components' method in GEE, which hides sets of contiguous pixels processed with values less than −25 dB in the VV polarisation [52,53]. The VV and VH noises were filtered using the refined Lee filter [54]. We collected S1 data on 13 and 18 July 2021. In this study, the mean of both images was combined with the Pleiades image. In addition, Sentinel-1 data helped obtain better results. SAR is able to reflect the soil's roughness, texture, and the different physical conditions of the ground. Another benefit of using SAR is that it has a certain amount of soil penetration in dry sand, which makes it particularly suitable for this specific area.

#### 2.4.3. Processing Environment

The Google Earth Engine (GEE) platform is currently a free answer to the matter of restricted access to processing power [19,35]. GEE is a cloud-based platform that enables a planetary-scale analysis of petabytes of freely obtainable satellite imagery. Combining this large amount of information with the parallel computing strength of Google's infrastructure facilitates the rapid and easy evaluation of satellite images. GEE parallelises and executes code developed in JavaScript or Python, exploiting Google's cloud computing infrastructure and allowing work with intensive processes at new scales. This study used GEE for the following three reasons: (1) it is a free platform that has many built-in functions for geospatial analysis and provides quick results without downloading datasets to a local storage system; (2) it has been widely used in several environmental studies [55–57]; and (3) it uses one processing environment and may reduce the uncertainties introduced by various techniques, such as temporal aggregation, resampling (e.g., downscaling and upscaling), and data reprojection. Notably, the final results from the GEE could be used to create maps (raster.tif) and tables (table.csv) in Google Drive for further processing in the Google Collaboratory (Colab) cloud service [58]. Pleiades bands (blue, green, red, and NIR) were uploaded to the cloud asset of the GEE environment because they are not yet available in the Earth Engine. Thus, the overall framework of the research can be applied to other regions to achieve related aims.

#### 2.4.4. Random Forest Classifier and Accuracy Report

Several artificial intelligence models can be used to process geospatial data when searching for sites in different environments [34,59–62]. The present study aims to detect and map fortified archaeological sites using an integrated method composed of RS data, GIS, field surveys, and a random forest algorithm. The steps of this model followed to classify archaeological sites included categorising the image composite, collecting training data, training the classifier model, and validating the classifier with a separated dataset. We use a selection of 101 archaeological sites linked in Blad Talh's investigation on predesert zones as our training (*n* = 76) and confirmation (*n* = 25) sets. Notwithstanding the existence of quality problems in Blad Talh's data, we sort out those sites that might be evidently specified and precisely situated in the high-resolution imagery obtainable in GEE. These matched up with large and well-preserved sites. Polygons were drawn in GEE to outline the areas of the designated sites from which pixel values in the image composite could be extracted for algorithm training. The random forest classification algorithm was applied in this study [56,63,64]. An RF classifier was chosen for the GEE machine-learning implementation. The RF establishes a number of 128 trees on preliminary training samples, but it takes every split in a tree into consideration. For instance, a new subset of predictors is considered randomly in every split when splitting nodes. In addition, the average of the resulting trees helps to evade overfitting and is therefore more stable and reliable than the other decision tree-based classifiers [63,65]. Moreover, RF classifiers can handle a small number of training samples, and it is possible to get the number of votes for each class. These are two advantages that are particularly useful for RS-based archaeology with limited land-use and cover information. In the current GEE algorithm, the RF is combined with trees, which is considered a sufficient value to gain perfect results without unnecessarily increasing the computational cost. The probability mode of the RF was set to evaluate, filter, and refine the results to ameliorate the algorithm's detection potential.

Usually, 75% of the samples are used for the training tree, while the remaining 25% are used for confirmation [63,65]. Herein, the machine-learning process went through three iterations. The original iteration of the algorithm generated pleasant results in that it had the ability to explicitly categorise known sites used as training data and many potential sites by high RF probability values. However, two additional iterations were necessary to adjust the pixels to higher percentages that did not correspond to archaeological sites. RF classifier analysis is a probability domain in bitmap format, where each value records the probability that a given pixel is a "vulnerable site." The RF probability limit for the site values was determined after careful examination of the test data on high-resolution images, which produced a raster map of the clusters ("fortified sites"). A higher threshold resulted in better identification of large, clear sites, but many small pixel clusters corresponding to the small sites covered were lost. Confusion metrics were used to calculate RF accuracy, including the overall accuracy and kappa coefficient. Several researchers demonstrated the perfect application of the GEE platform for training and applied an RF classifier due to its high performance [19,66].

#### **3. Results**

Using the proposed methodology, we succeeded in detecting fortified sites that were previously unknown in Blad Talh. The automated detection of fortified sites exhibited high accuracy, equal to that of human detection. The results show that machine learning performed with high accuracy (see Table 2).

**Table 2.** Classification results of the study area.


The image combines SAR bands (a single VV and VH and a dual VH–VV band) and P1B bands (b1–b4) (as shown in Figure 3). The high-quality P1B and Sentinel-1 bands are not affected by specific environmental or vision conditions.

**Figure 3.** Feature importance in the classification model.

To the best of our understanding, as of mid-2021, the P1B image received on 16 July 2021 was the most current VHR optical satellite image covering Blad Talh that was accessible online through the GEE. The GEE and machine learning have proven to be applicable instruments for automatically detecting archaeological sites in satellite images. The GEE, machine learning, and remote sensing were highly effective at the site scale and identified pixels of archaeological remains. Furthermore, pixels of fortified structures were successfully detected in Blad Talh, confirming the ability of the GEE to quickly run an analysis on a substantial scale. The discovery and distribution of archaeological sites extended towards the southern part of the Blad Talh area and the surrounding pre-desert towards the east, which was a region that was relatively unreachable in past investigations. Predictably, the allocation of the fortified sites detected by the algorithm was concentrated in the Blad Talh area. Sentinel-1 texture differences between layers were accentuated in the area (see Figure 3).

The algorithm's ability to detect fortified archaeological sites is likely due to polarisation bands and the capacity of the SAR C-band to penetrate soil. Because of vegetation and sand cover, the RF potential produced only some well-illustrated square shapes. The majority of the newly proposed fortified sites presented segmented square shapes or amorphous groups of pixels. It is quite possible that in addition to these five detection sites, there are also remains of sand-covered archaeological sites with low RF probability values. Notably, when high-resolution visible images were used, the algorithm still assisted in identifying small pixels or backscattering. SAR alone showed no significant change in land cover. As a result, SAR alone cannot give enough information to identify archaeological sites. Also, the Pleiades images are not capable of segregating the spectral structures of the sites on their own using an RF classifier. Furthermore, our algorithm was designed to integrate Sentinel-1 and PB1 images for better archaeological discovery, an approach that showed accurate results (see Figure 4). Comparing the results between high-resolution imagery PB1 without SAR images and the virtual absence of sites, their combination for archaeological detection of the known fortified site was used as a validation set (see Figure 4). When possible, the new detection was matched to previous research on past archaeological sites [8,47,67]. The

archaeological sites relative to the existence of sand dunes suggested that many more sites can be covered or located underground in the pre-desert, and thus, the aeolian sediments in southern Tunisia may have played an important role in both the field records of the area and the discernibility of the archaeological ruins. Regions with a large number of sites were located in open mudflats that were distributed throughout the area, while many sites were still relatively covered underground. Five newly detected fortified sites were identified (see Figure 5). The distribution of the archaeological sites extended towards the southern part of the Blad Talh area and the surrounding pre-desert towards the east, which was a region that was relatively unreachable in past investigations. The resulting cluster of high-RF probability pixels showing fortified sites was used to rebuild the regions. Image interpretation was performed using P1B, SAR data, and machine learning. The blend of these data secured enough environmental variability and time to estimate and determine the potential extent of the fortified sites that the algorithm specified. Because of the sand cover, the RF probability yielded only a few well-defined square shapes. Most of the newly proposed sites present fragmentary square shapes and elongated bands of pixels. It is quite possible that in addition to these discovered sites, there are also the remains of archaeological sites covered by sand or shrub vegetation with low RF probability values. It is worth emphasizing that the algorithm helped to identify small groups of site pixels.

**Figure 4.** (**A**) The study area, showing the distribution of sites. (**B**) A Google Earth base map locating a well-preserved fort. (**C**) Dual Sentinel-1 band [VV, VH] in ascending mode. (**D**) Single Sentinel-1 band [VV] in ascending mode. Results of the RF classifier: (**E**) Visible high-resolution PB1 imagery without SAR images compared to the virtual absence of sites. (**F**) An example of SAR images combined with PB1 that detected a known fortified site used as a validation set; note that the white dots scattered through the region indicate high-RF probability sites.

**Figure 5.** RF probability results. (**A**) Illustration of fortified detection sites. (**B**) Clustering of highprobability pixels in the area showing the detection structures of fortified sites. (**C**) Pan sharpened Pleiades imagery. Image was interpreted as signifying the existence of a fort hidden underground.

#### **4. Discussion**

The RF method, which has been used in the classification of remote sensing data since the early 2000s, was chosen as the classification algorithmic rule for this study due to its simplicity of implementation, robustness, and possible predictability [63,68,69]. The challenge we faced was the discovery of fortified sites. Our objective was to develop a straightforward yet efficient methodology that combined remote sensing techniques and automated archaeological remains discovery, as well as to evaluate the potential offered by the combined use of remote sensing data and machine learning to discover fortified sites in arid areas. Using a substantive approach, we confirmed an advantage similar to that of other applications of machine learning techniques, and we show that the RF potential field-detected sites may be considered ancient settlements. In this research, given the high likelihood of improvement in the classification accuracy, we assessed the RF's performance in a wider context (different archaeological features or sites). As shown in the results section, the performance of the RF classifier is highly satisfactory, with a kappa coefficient of 0.91 and a precision and overall accuracy of 0.93 in the "fortified site group." The RF potential area highlights five unknown locations shown on historical maps. The sites are partly hidden by sand, and the RF potential only shows different parts of the structures' fortified sites, which was verified by a field survey. The discovery of these sites is very important because it documents the existence of ancient Roman fortifications. Based on previous studies, there are more sites of Roman civilisation in this area [8,47], and only some small sites have been uncovered by remote sensing methods. Currently, new sites could be included with those previously identified and dated by Bachagha et al. [8,47,67]. To try to understand the changes in the settlement distribution over time, research in southern Tunisia documented all sites encountered, including those of the Roman era. Considering that most of the previously reported sites were settlements attributed to periods of Roman civilisation, there is a high possibility that most of the newly discovered sites were also occupied during the same period. Fortified structures such as fortresses and defensive works are fairly common features of the archaeological landscape in southern Tunisia. Such features have previously been uncovered in the city of Gafsa. In this study, we were able to locate five fortified sites situated south of the Blad Talh area. The fortified sites were built on approximately level ground in a nearly square form that is covered by sand dunes. The layout showed visible signs of erosion on all sides, and loot pits were visible, especially on the south side. Based on our findings in the literature, we propose a suitable reconstruction of the ancient settlement in Gafsa, indicating that the fortifications in the locale were a defence project planned to defend, protect, and control the Limes [8,70–74]. We found fortified sites located in the same area and close to each other, in contrast to previous research [8]. The fortified sites were separated by a distance of approximately 2.5 km. From the machine learning and P1B data of the study area, a red-dashed line rectangle can be seen in Figure 5. Using machine learning and P1B data surveying, we demonstrate the detection of hidden fortified sites in Gafsa, which is well-situated across the frontier of the Roman or Byzantine era in southern Tunisia.

#### *Surveys: Validation and Archaeological Site Significance*

Building on recent advances in digital image analysis and feature detection, we developed an RF algorithm that allows the automated detection of these types of enclosures in Pleiades images with high resolution and Sentinel-1 for better archaeological detection. We also used SAR images to increase the accuracy of classification, which provided good results. In addition, we only used the Google Earth Engine for the view, so we cannot use it for building a machine learning detection model. We successfully used the RF model to detect new sites in satellite imagery. The model detected these sites, which we verified and confirmed in the field. We built the training dataset using both Pleiades and Sentinel-1 images in the Google Earth Engine. Of course, the identification of any object or location as archaeological remains based on images requires confirmation from the ground. The discovered sites described above were confirmed in the area of Blad Talh in southern

Tunisia by machine learning, remote sensing data, and field investigation. Other than the surrounding wall and the ruins of the fortified sites, no other archaeological ruins were seen above ground. Almost all the visual features in the analysed RF results were determined to be micro-topographical results of archaeological importance. As a result, our approach reveals a different layout of archaeological sites, and the algorithms identified structures in the area very well. In November 2021, the exploration team arrived at sites approximately 10 km north of the suspected location of the ancient Roman fortification, KSAR GRAOUECH (a square-shaped mark). Each side of the suspected location was measured to be around 50 m long. These anomalies have no intact walls, but the remains of some wall piers were investigated to identify anomalies. Shards of ancient ceramics were scattered among the walls, and their environs date back to the Roman period, as shown in Figure 6G,I,N,L. Buried walls at these five fortified sites were also unearthed, representing the appearance of concealed remains, as shown in Figure 6C,F,H,J,M. The field investigation conducted in this region reveals the existence of a farm dating back between the Roman Empire and Byzantine times based on the discovery of tiles and shards of ancient ceramics from the Roman period. These fortified sites are situated between Capsa (Gafsa) and the coastal towns of Iunci and Lariscus. Controlling one of the passes to the north would have been of obvious strategic importance during the Roman era, a consideration that strongly supports their inclusion in the list of official fortifications of that period. In addition, these sites have common topographical and archaeological features, and their dimensions are relatively similar to each other. According to Trousset [74], these fortified sites are widely spread in the rear Limes area and are located on the immediate edge of a Roman road or in a control location of a passage; some have been identified as military works that are a part of the defensive system. In the report by Rebuffat et al. [75], these were probably fortified farms that adopted a plan inspired by the military project. The parallels between military architecture and these fortified sites, namely Mattingly Qsur, are very evident in the adoption of the rectangular form with towers projecting from two or four facades [76]. The results of studies investigating this type of site in Fazzan and Cyrenaica in eastern and western Libya, respectively, confirm that the peak occupation of fortified sites occurred after the end of the third century and continued in some regions in later times. More specifically, the majority were occupied between the 4th and 6th centuries (350–540 AD) [77].

Indeed, the general landscape of the Blad Talh was a highly important part of the ancient Roman period. As a result, based on our fortified sites and the literature, we firmly conclude that Blad Talh has played a critical historical role in migration between ancient cities located in current-day north and south Tunisia, as attested by the caravan roads. It was also a part of the sensitive defence system (Limes) located along the desert border, as attested by the several forts that protected the area (Figure 7).

**Figure 6.** Magnified views of anomalies 1 to 5. (**A**,**B**) Landscape of the study area. (**G**,**I**,**N**,**L**) Shards of ancient ceramics. (**E**) Ancient Architectute element (**K**) Lintel photographed at the entrance to a fortress. (**C**,**D**,**F**,**H**,**J**,**M**) Field photos of wall remains.

**Figure 7.** Comprehensive archaeological map of the Blad Talh landscape in the ancient Roman period: localities and ancient road networks. Modified from [71].

#### **5. Conclusions**

In summary, we present an integrated workflow that combines SAR images and optical Pleiades imagery with basic spatial analysis in the Google Earth Engine for automatically detecting fortified sites. The random forest-based approach demonstrates that it is a useful tool for overcoming problems with data size and structure. Our research illustrates the applicability of this detection technique to fortifications and demonstrates that it is possible to automate fortification discoveries for the first time in Tunisia. It is possible to summarise insights and knowledge directly from the data, and the algorithm is able to highlight relationships between archaeological sites. The machine-learning algorithm we use is capable of identifying all the currently known sites in the study area and the accurate locations of new fortified sites. The approach provides outputs that are significantly superior to those obtained with a single-source RS technique. RS-based applications in arid and semi-arid regions elsewhere can benefit from the integration of globally available Pleiades and Sentinel data into an accessible, flexible, and repeatable GEE environment to perform and evaluate machine learning workflows. Even though the approach shows accurate results, the limited availability of VHR data and the expense to obtain the data are considered a limitation of this study that can be improved in the future. Archaeological data, along with landscape analysis and mapping items affecting site visibility, indicates that these are only a few of the fortified sites in the region, many of which may be covered by sand in the surveyed area. We discovered fortified sites hidden by thick vegetation and confirmed these sites using field investigation, remote sensing, and machine learning. The outcome of the work presented here provides advances in the technology and archaeology domains. Technologically, using these methods, various derived models are used to evaluate their abilities in investigating the micro-topography of a site of archaeological importance under a sand dune. The results viewed from an archaeological perspective are interesting. The machine learning-based method enabled us to detect five fortified sites for which only some data was obtainable from past records and shed new light on the fortified KSAR GROUECH structure dating back to the Byzantine period (6th century), for which only the presence of a fort was well-known. The findings are validated by on-site processing, demonstrating that the automatic procedure can identify and extract the major distinctive characteristics of Roman settlements, such as villas and tower walls. Notably, the exploration elaborated in this research is promising in terms of technological invention. By providing for previous quantitative examinations in this region, our exploration raises awareness of the need to use quantitative styles to address more pressing questions, similar to those about the protection and preservation of threatened archaeological sites. As a result, one of the benefits of this research is that it demonstrates, from the perspective of archaeologists, the capability of satellite data and machine learning to discover buried archaeological sites. While this study focuses on archaeological sites in southern Tunisia, the proposed approach is effective in terms of time and cost, particularly in locations where data availability is scarce. Validation of this approach was proved using the nature and characteristics of the study area, and we expect that the classification accuracy of our approach would be further improved by using a convolutional neural network (CNN) and images obtained by the Google Earth Engine.

**Author Contributions:** Conceptualisation: N.B., W.X. and A.E.; methodology, A.E. and N.B.; software, N.B. and A.E.; validation, N.B., W.X and A.E.; formal analysis, N.B.,W.X. and A.E.; investigation, N.B., M.T. and A.E.; resources, W.X., N.B., and M.T.; data curation, N.B., M.T., A.E., F.S.; writing—original draft preparation, N.B., W.X., A.E., M.T., W.X. and F.S.; writing and supervision, W.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (No. 42174023).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Sentinel-1 data are available via the European Space Agency's Copernicus Open Access Hub https://scihub.copernicus.eu/ (accessed on 13 July 2021).

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **A Simulation Method of Two-Dimensional Sea-Surface Current Field for Trajectory Crossing Spaceborne SAR**

**Yan Li 1,2,3, Jinsong Chong 1,2,3,\* and Zongze Li 1,2,3**


**Abstract:** The trajectory-crossing method is an important and effective method for spaceborne SAR (synthetic aperture radar) to detect a two-dimensional current field. However, in practical applications, there are few spaceborne SAR data that meet the requirements of close acquisition time and overlapping trajectories, which makes it difficult to comprehensively analyze the impact of different systems and environmental parameters on current measurement results from real data processing. With the proposal of a SAR constellation plan in the future, the data resources will be constantly enriched, and the trajectory-crossing method will be widely used. It is, thus, necessary to lay some theoretical foundation at present. A two-dimensional sea-surface-current-field simulation method for trajectory crossing spaceborne SAR is proposed in this paper. Based on the principles of the trajectory-crossing method, the proposed method can realize the simulation of two-dimensional sea-surface-current field under given spaceborne SAR system and environmental parameters. In this paper, the simulation process of this method is given, and the simulation experiment is performed. Compared with the current measurement results of SAR data, the simulation experiment shows that the current velocity error is less than 0.03 m/s and the direction error is less than 10 degrees, which proves the reliability of the proposed simulation method. The proposed method lays a foundation for analyzing the influence of various parameters in the application of the trajectory-crossing method.

**Keywords:** synthetic aperture radar; trajectory-crossing method; sea-surface current; simulation analysis

#### **1. Introduction**

As an important ocean dynamic parameter, the sea-surface-current field has a profound impact on climate change, ocean engineering, fishery resources, energy, offshore target detection and so on [1,2]. Therefore, monitoring of ocean currents is becoming crucial. Spaceborne synthetic aperture radar (SAR) has the ability to observe the earth all day, in all weather and at high-resolution, which makes it an important remote-sensing tool that can retrieve the sea-surface-current field.

Currently, along-track interferometry (ATI) [3] and Doppler centroid analysis (DCA) [4] are the two main techniques for spaceborne SAR to detect ocean current fields; however, they can only retrieve the range velocity component of the current, i.e., they cannot directly obtain the current vector—that is, the two-dimensional current field. Generally, two sets of SAR data with crossed trajectories can be used to synthesize the velocity components in two different directions on the overlapping sea-surface area to thus obtain the current vector. This is an important two-dimensional current field detection method, which is denoted as a trajectory-crossing method in this paper.

However, due to the requirement that the overlapping area needs to cover the research scope and as the data acquisition time is close, there are fewer real spaceborne SAR data to meet this feature. The current measurement ability of the trajectory-crossing method is only

**Citation:** Li, Y.; Chong, J.; Li, Z. A Simulation Method of Two-Dimensional Sea-Surface Current Field for Trajectory Crossing Spaceborne SAR. *Appl. Sci.* **2022**, *12*, 5900. https://doi.org/ 10.3390/app12125900

Academic Editor: Tung-Ching Su

Received: 17 February 2022 Accepted: 7 June 2022 Published: 9 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

verified in some airborne [5] and spaceborne [6] SAR data processing, and there is a lack of further research for the impact of different system and environmental parameters on the current measurement effect. Therefore, simulation is an effective research method. At present, to the best of knowledge, there is no specific simulation method for the trajectory crossing current measurement process.

Due to the complexity of sea-surface electromagnetic scattering, the simulation of seasurface SAR data usually does not start from the original echo but mainly considers the modeling using wave spectrum and modulation transfer function. In 2000, Romeiser and Thompson [7] proposed a sea-surface Doppler spectrum model based on the improved combined sea-surface model, in which the wave–current interaction and the modulation transmission process of sea-surface backscattering were considered [8,9]. Then, the model was integrated into SAR ocean imaging model and widely used in sea-surface current measurement and simulation [10–12]. Therefore, this model provides a possibility for the simulation of two-dimensional sea-surface-current field for trajectory crossing spaceborne SAR.

Based on the principle of the sea-surface SAR imaging process and trajectory-crossing method, a two-dimensional current field simulation method is proposed in this paper. In the second part, the specific process of the proposed simulation method is given. In the third part, the simulation experiment and verification are conducted. Finally, the above is summarized.

#### **2. Simulation Method of Two-Dimensional Sea-Surface Current Field for Trajectory Crossing Spaceborne SAR**

The trajectory-crossing method needs to use two sets of spaceborne SAR data with overlapping areas, usually from ascending and descending passes, respectively, and their acquisition time should be close. Geometric relationships of different velocity components of sea-surface motion over overlapping area are shown in Figure 1, where the blue and red thick arrows indicate the flight direction of the ascending and descending satellite, respectively. The sea-surface imaging range is *S* and *G*, bounded by blue and red dashed lines, respectively. The middle is the overlapping area of *S* and *G*, where *Vtot* and *Dtot* represent the magnitude and direction of current velocity, respectively. Sea-surface movement produces projection components in four directions, where *VS* and *VG* represent horizontal radial velocity components along beams far from the ascending and descending satellites. *VV* and *VU* represent surface velocity components along the north–south and east–west directions, respectively. In addition, *α* represents the intersection angle between the direction of the beam center line and the east–west direction on the horizontal plane.

**Figure 1.** Schematic diagram of the geometric relationship between the sea-surface velocity and its components in the trajectory-crossing method.

The simulation method proposed in this paper is mainly divided into three parts: radial velocity simulation, wave velocity simulation and current vector simulation. The flow chart is shown in Figure 2. The first step is to calculate the radial velocity using the SAR ocean imaging model, the second step is to calculate the wave motion component, and the third step is to calculate the two-dimensional current field combined with the trajectory-crossing method. The radar parameters include the signal frequency, look direction of the beam and incident angle. Platform parameters include the satellite orbit height and flight speed.

**Figure 2.** Flow chart of two-dimensional sea-surface-current-field simulation for trajectory crossing spaceborne SAR.

Based on the SAR ocean imaging model [7–9] proposed by Romeiser, University of Hamburg, Germany, the sea-surface horizontal radial velocity components along the spaceborne SAR beam look in directions with the ascending and descending passes simulated. The model supports simulation in the frequency range of 0.4 to 35 GHz. The main processes include the calculation of the wave–current interaction model, sea-surface scattering model and sea-surface SAR Doppler model.

Wave–current interaction is a physical process that must be considered in the real sea motion environment. We need to first give the sea-surface-current field and wind field parameters and then solve the action spectrum balance equation and calculate the modulated wave spectrum according to the weak hydrodynamic interaction theory. The action spectrum balance equation to be solved is [9,13,14]

$$\frac{\partial N(\mathbf{x}, \mathbf{k}, t)}{\partial t} + \left[c\_{\mathcal{S}}(\mathbf{k}) + \mathcal{U}(\mathbf{x}, t)\right] \frac{\partial N(\mathbf{x}, \mathbf{k}, t)}{\partial \mathbf{x}} - k \frac{\partial \mathcal{U}(\mathbf{x}, t)}{\partial \mathbf{x}} \frac{\partial N(\mathbf{x}, \mathbf{k}, t)}{\partial \mathbf{k}} = Q(\mathbf{x}, \mathbf{k}, t) \tag{1}$$

where *t* represents time, *x* = (*x*, *y*) represents the sea-surface spatial position, *k* = - *kx*, *ky* represents the wave number, and *kx* and *ky* represent the wave numbers in the X and Y directions, respectively. *N* represents the action spectral density of sea-surface micro scale wave, *U* represents the surface current field, and *Cg* represents the modulated wave group velocity. These variables about ocean waves are related to the conditions of sea-surface wind field. *Q* represents the source function:

$$Q(\mathbf{x}, \mathbf{k}, t) = \mu(\mathbf{k}) N(\mathbf{x}, \mathbf{k}, t) \left( 1 - \frac{N(\mathbf{x}, \mathbf{k}, t)}{N\_0(\mathbf{k})} \right) \tag{2}$$

where *N*<sup>0</sup> represents the action spectral density in the sea-surface equilibrium state, and *μ* represents the relaxation rate.

The relationship between the wave spectrum and action spectral density is [15]

$$N(\mathbf{x},k,t) = \frac{\rho \omega\_0(k)}{k} \psi(\mathbf{x},k,t) \tag{3}$$

where *ψ* represents the modulated wave spectrum. The natural angular frequency *ω*0(*k*) = *gk* + (*τs*/*ρ*)*k*3, *g* is the gravitational acceleration, *ρ* is the seawater density, and *τ<sup>s</sup>* is the seawater surface tension.

The sea-surface scattering model is based on the improved combined sea-surface model [8,9] proposed by Romeiser and Alpers, which considers the modulation effect of large-scale wave on small-scale Bragg waves. The normalized scattering coefficient can be calculated by taking the modulated wave spectrum and the given radar parameters as the input of the scattering model [8,16]

 *<sup>σ</sup>* <sup>=</sup> *<sup>σ</sup>*(0) <sup>+</sup> *σ*(2) = *σ*|*s*=<sup>0</sup> + <sup>⎛</sup> ⎝ *∂*<sup>2</sup> ∧∨ *σ ∂* <sup>∧</sup> *sp ∂* <sup>∨</sup> *sp* ∧ *s*=0 + *∂*<sup>2</sup> ∨∧ *σ ∂* <sup>∨</sup> *sp ∂* <sup>∧</sup> *sp* ∧ *s*=0 ⎞ <sup>⎠</sup> · *<sup>k</sup>*<sup>2</sup> *<sup>p</sup>ψ*(*k*)d2*<sup>k</sup>* + <sup>⎛</sup> ⎝ *∂*<sup>2</sup> ∧∨ *σ ∂* <sup>∧</sup> *sn ∂* <sup>∨</sup> *sn* ∧ *s*=0 + *∂*<sup>2</sup> ∨∧ *σ ∂* <sup>∨</sup> *sn ∂* <sup>∧</sup> *sn* ∧ *s*=0 ⎞ <sup>⎠</sup> · *<sup>k</sup>*<sup>2</sup> *<sup>n</sup>ψ*(*k*)d2*<sup>k</sup>* + <sup>⎛</sup> ⎝ *∂*<sup>2</sup> ∧∨ *σ ∂* <sup>∧</sup> *sp ∂* <sup>∨</sup> *sn* ∧ *s*=0 + *∂*<sup>2</sup> ∧∨ *σ ∂* <sup>∨</sup> *sn ∂* <sup>∧</sup> *sp* ∧ *s*=0 *∂*<sup>2</sup> ∨∧ *σ ∂* <sup>∨</sup> *sp ∂* <sup>∧</sup> *sn* ∧ *s*=0 + *∂*<sup>2</sup> ∨∧ *σ ∂* <sup>∨</sup> *sn ∂* <sup>∧</sup> *sp* ∧ *s*=0 ⎞ <sup>⎠</sup> · *kpknψ*(*k*)d2*<sup>k</sup>* (4)

where the symbol · represents the statistical average; *<sup>σ</sup>*(0) represents the normalized backscattering coefficient when the sea-surface slope is zero; *σ*(2) is the sum of secondorder Bragg scattering caused by sea-surface slope; *s* = - *sp*,*sn* represents the sea-surface slope vector; *kp* and *kn*, respectively, represent the Bragg wave number components parallel and perpendicular to the radar look direction; *ψ*(*k*) represents the wave number spectrum; ∧ and ∨ denote the Fourier transform and its conjugate of *σ* to wave number *k*; and ∧∨ and ∨∧ denote the Fourier transform and its conjugate of *σ* to combined wave number *k*<sup>1</sup> − *k*2.

The Doppler spectrum and variance can be calculated using the double Gaussian Doppler spectrum model [7] proposed by Romeiser and Thompson, which divides the sea-surface echo into the superposition of the Bragg wave Doppler spectrum far from and towards the radar. Each Doppler spectrum component is defined as a Gaussian shape and expressed as

$$\begin{split} S(f\_D) &= \frac{\langle \sigma\_+ \rangle}{\sqrt{2\pi \gamma\_{D+}^2}} e^{-(f\_D - \langle f\_{D+} \rangle\_{\mathcal{S}'})^2} \Big/ \gamma\_{D+}^2 \\ &+ \frac{\langle \sigma\_- \rangle}{\sqrt{2\pi \gamma\_{D-}^2}} e^{-(f\_D - \langle f\_{D-} \rangle\_{\mathcal{S}'})^2} \Big/ \gamma\_{D-}^2 \end{split} \tag{5}$$

where ± represents the Bragg wave components away from and towards the radar, respectively; *fD*± *<sup>σ</sup>* is the Doppler center weighted by the normalized backscattering coefficient; and *γD*<sup>±</sup> is the variance of Doppler spectrum. The specific calculation formulas of *fD*± *<sup>σ</sup>* and *γD*<sup>±</sup> can be found in [7]. Then, the horizontal Doppler velocity *Vdop* along the beam radial direction can be expressed as

$$V\_{dop} = \frac{\pi \cdot \langle f\_{D\pm} \rangle\_{\sigma}}{k\_{\varepsilon} \sin \theta} \tag{6}$$

where *ke* is the wave number of the electromagnetic waves, and *θ* is the incident angle. Therefore, the radial horizontal sea-surface Doppler velocity along the ascending and descending satellite beam represented by *VS* and *VG* in Figure 2 is calculated successively by the above process.

The radial velocity is actually the Doppler velocity that reflects the average motion state of scatterers in the sea-surface resolution unit. In the range of medium incident angles, the radial Doppler velocity *ur* includes not only the current velocity but also the component of the sea-surface wave motion [17].

$$
\mu\_r = \mu\_\varepsilon + \mu\_\sigma + \mathfrak{u}\_\mathbb{b} \tag{7}
$$

where *uc* is the current velocity, *uo* is the large-scale wave orbital velocity, and *ub* is the net Bragg wave phase velocity.

At present, the related wave motion model is not perfect, and thus it is difficult to extract the current velocity component by directly calculating the accurate wave velocity [10]. The authors in [18] stated that when there is no current field, the sea-surface velocity obtained by using the SAR ocean imaging model comes from the wave motion. Therefore, in the second step, the input current field velocity is set to zero, combined with the same wind field, and through the calculation process similar to the radial velocity, the wave motion velocities in the ascending and descending beam direction are obtained, which are *V<sup>w</sup> <sup>S</sup>* and *<sup>V</sup><sup>w</sup> <sup>G</sup>* , respectively. Since the influence of large-scale waves is considered in the previously used wave spectrum and scattering model, the simulated wave velocities include the sum of the last two velocity components of Equation (7).

By making a difference between the radial Doppler velocity and the wave velocity, the current velocity components *V<sup>c</sup> <sup>S</sup>* and *<sup>V</sup><sup>c</sup> <sup>G</sup>* in the ascending and descending beam direction can be obtained:

$$V\_S^c = V\_S - V\_S^w \tag{8}$$

$$V\_G^c = V\_G - V\_G^w \tag{9}$$

Then, combined with the trajectory crossing current measurement geometric model in Figure 1, the east–west velocity component *VU* and north–south velocity component *VV* of the sea-surface current can be calculated:

$$V\_{ll} = \frac{V\_S^c - V\_G^c}{2\cos\alpha} \tag{10}$$

$$V\_V = \frac{V\_S^c + V\_G^c}{2\sin\alpha} \tag{11}$$

Finally, the magnitude and direction of the current vector are, respectively, expressed as:

$$V\_{tot} = \sqrt{V\_{ll}^2 + V\_{V}^2} = \frac{\sqrt{V\_S^2 + V\_G^2 - 2V\_S V\_G \cos(2a)}}{\sin(2a)}\tag{12}$$

$$D\_{\rm tot} = \arctan\left(\frac{V\_V}{V\_{\rm II}}\right) = \arctan\left(\frac{V\_S + V\_G}{(V\_S - V\_G)\tan\alpha}\right) \tag{13}$$

#### **3. Simulation Experiment and Verification**

In this section, a simulation experiment was conducted based on the spaceborne SAR system and sea environment parameters. The effectiveness of the simulation method was verified by comparing the simulation results with the measurement results using SAR data. First, the spaceborne SAR data, system parameters and input sea environment data involved in the experiment are introduced. The geographical location of the spaceborne SAR data used in this paper is shown in Figure 3. This is the Mozambique strait between southeast Africa and Madagascar, where the sea surface is affected by the warm current of Mozambique, and the motion state is relatively stable. The red box indicates the stripmap SAR data range of ESA sentinel-1 satellite's ascending pass, the acquisition time is 8 October 2019, 15:35 UTC (orbit number: 018385), the black box indicates the stripmap SAR data range of China Gaofen-3 (GF-3) satellite's descending pass, and the acquisition time is 9 October 2019, 02:43 UTC (orbit number: 016651). Although these two SAR data were not acquired at the same time, the time interval is short, the current motion status is basically unchanged here, and these two SAR images overlap partially in this sea area, which meets the requirements of the trajectory-crossing method.

**Figure 3.** Spatial geographic location of the actual SAR data.

Figure 4 shows the SAR amplitude image of Sentinel-1 satellite and the horizontal radial Doppler velocity image of its secondary product. The dimmer brightness in the lower left corner of Figure 4a is due to the lower wind speed and the weaker echo energy due to the oil spill on the sea, which also results in errors in the Doppler velocity, as shown in the corresponding position in Figure 4b. The absolute values are small in the lower right part of the Doppler velocity map, and the overall velocity includes the influence of sea-surface wind and wave motion, which needs to be further removed.

**Figure 4.** Sentinel-1 SAR (**a**) amplitude image and (**b**) radial Doppler velocity.

Figure 5 shows the SAR amplitude image of the GF-3 satellite and the Doppler velocity estimated from the single look complex (SLC) data. The dark part of the right corner of Figure 5a is obviously affected by the sea oil film, which reduces the signal-to-noise ratio of the image. The Doppler velocity of Figure 5b is estimated from the Doppler processing in reference [19]. It can be seen that the absolute value in the lower left corner of the image is

small, in addition to the anomaly in the oil film position, and the overall velocity overlaps with the influence of the wind and wave motion.

**Figure 5.** GF-3 SAR (**a**) amplitude image and (**b**) radial Doppler velocity.

Although the overlap range of the above SAR data is small, it has clear velocity characteristics, which is conducive to the validation of simulation method. Therefore, a simulation experiment was performed in this sea area. The input system parameters are consistent with the real SAR data parameters, representing the ascending Sentinel-1 satellite and the descending GF-3 satellite, respectively, as shown in Table 1. It should be noted that, at present, there are few ascending and descending data from a single satellite that can meet the requirements of irradiating the same sea area at a short time, and it is difficult to obtain them directly. In future applications, two or more satellites must cooperate to obtain the ascending and descending passes data of the same space-time sea area for current measurement. Therefore, the following dual satellite SAR data processing and system simulation are of practical significance.


**Table 1.** Simulation parameters.

The hybrid coordinate ocean model (HYCOM) reanalysis current field data and the European centre for medium-range weather forecasts (ECMWF) wind field data are obtained for the simulation input, respectively, as shown in Figure 6. The HYCOM model reanalysis data, which combines satellite altimeter data with temperature and salt data obtained by a buoy, is widely used in ocean current research [20–22]. The date selected for the HYCOM data is the same as Sentinel-1 SAR data acquisition date with spatial resolution of 1/12 degree. Wind field data at the same time and place were obtained from the ECMWF website with spatial resolution of 1/4 degree, and interpolated to ensure the same number of data points as the input current field. From the input data, it is known that the velocity in the overlapping area is small, and the sea-surface wind speed is mainly in the range of 2–7 m/s, which belongs to the moderate wind speed.

**Figure 6.** Simulation input sea-surface environment data. (**a**) current field and (**b**) wind field.

The given current and wind data are input into the SAR ocean imaging model to calculate the radial Doppler velocity *VS* under the parameters of ascending Sentinel-1 satellite and the radial Doppler velocity *VG* under the parameters of descending GF-3 satellite, as shown in Figure 7a,b, respectively. It can be seen that, due to the same orbital inclination of these two satellites, the distribution of the relative value of the velocity in the look direction of the two beams is approximately symmetrical; however, the absolute value is different due to the influence of wind and wave motion. Figure 7c,d shows the wave velocity when the current field is set to zero under the same wind field conditions, in which the absolute value of the velocity changes due to the change of the relationship between wind direction and beam look direction. When the wind direction is perpendicular to the beam look direction, the wave velocity is zero, which is consistent with the research conclusion in the literature [11].

**Figure 7.** (**a**) Sea-surface radial velocity along the ascending beam. (**b**) Sea-surface radial velocity along the descending beam. (**c**) Wave velocity along the ascending beam. (**d**) Wave velocity along the descending beam.

According to the flow chart in Figure 2, the corresponding sea-surface radial velocity and wave velocity of ascending and descending passes in Figure 8 are subtracted, and the current vectors are calculated using the trajectory-crossing method. The simulated two-dimensional current field is shown in Figure 8. Comparing with the input current field given in Figure 6a, we found that they were in good agreement. Note that, as for the computational demand, the calculation process of this method takes only a few minutes, the calculation amount of wave spectrum is the largest. The calculation times will be different according to the amount of input current and wind data.

**Figure 8.** Simulation results of two-dimensional current field.

To quantitatively validate the experimental results, we selected eight geographic locations in the overlap area (Figure 9) and compare the simulation results with the SAR data measurement results, as shown in Table 2. As the size of the simulated current field resolution unit is about 8 km and that estimated from Sentinel-1 and GF-3 SAR data is about 1 km, the 8 × 8 Doppler velocities around the corresponding eight latitude and longitude positions in Figures 4b and 5b are averaged. After subtracting the simulated wave velocity in Figure 7 from the average data, the two-dimensional current field results of the SAR data are calculated in combination with Equations (8)–(13). It can be seen that the corresponding current velocity and direction are close. However, because the current and wind data input by the simulation will be different from the instantaneous environment in the real SAR data, there will also be discrepancies between the data in Table 2.

**Figure 9.** Location points in the current field as verification.


**Table 2.** Comparison between the simulated and measured current data.

The statistical analysis of the experimental data is shown in Figure 10. Figure 10a shows a scatterplot of the velocity of the simulated current field compared with that of the measured current field from SAR data, with a correlation coefficient greater than 0.9 and a root mean square error (RMSE) of 0.028 m/s. Figure 10b shows the scatterplot comparing the direction of the simulated current compared with that of the measurement data. The correlation coefficient is greater than 0.98 and the root mean square error (RMSE) is 6.05 degrees. It can be seen that the simulated current field is in good agreement with the measured current of SAR data.

We also notice that there is a visible deviation between these two data, which will not only come from the difference of input data but also be limited by the model. For example, the wave spectrum is different from the real sea-surface state. Nevertheless, we believe that the simulation results are consistent with the measurement results. Therefore, the proposed method in this paper can effectively simulate the two-dimensional sea-surface-current field in the given spaceborne SAR system and marine environment parameters. The results also show that the proposed simulation method is conducive to the parameter analysis process of the trajectory-crossing method in the application of current measurement in the future.

**Figure 10.** Statistical comparison between the simulated current field and SAR data measured current field. (**a**) Velocity and (**b**) direction.

#### **4. Conclusions**

With the continuous enrichment of spaceborne SAR data in the future, the trajectorycrossing method to measure the two-dimensional current field on the sea surface will be widely used. However, at present, it is difficult to comprehensively analyze the influence of system and environmental parameters on the current measurement effect from the real SAR data. In this paper, a simulation method of two-dimensional sea-surface-current field for trajectory crossing spaceborne SAR is proposed, which can realize two-dimensional current field simulation under the given spaceborne SAR system and environmental parameters. In the simulation process, the wave–current interaction and the influence of seasurface backscattering on the SAR Doppler centroid frequency are considered, and the influence of the wave motion on the current measurement results is removed.

The simulation experiment selected the real spaceborne SAR system and sea-surface environment parameters, and the results are in good agreement with the SAR data processing results, which shows the effectiveness of the simulation method in this paper. Therefore, the proposed method is conducive to the parameter analysis process of the trajectory-crossing method in the application of current measurement. However, because the influence of the wave-breaking scattering mechanism under high wind speed and sea conditions was not considered, the simulation method in this paper is not suitable for twodimensional current field simulation in extreme marine environments. In future work, more scattering mechanisms will be considered to improve the scope of application and promote the development of application research for two-dimensional sea-surface-current field measurements.

**Author Contributions:** Conceptualization, Y.L. and J.C.; Methodology, Y.L.; Software, Y.L. and Z.L.; Validation, Y.L. and Z.L.; Formal analysis, Y.L., J.C. and Z.L.; Investigation, Y.L.; Writing—original draft preparation, Y.L.; Writing—review and editing, J.C. and Y.L.; Supervision, J.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


**Jiarui Shi 1,2, Qian Shen 1,2,\*, Yue Yao 1,2, Fangfang Zhang 1,2, Junsheng Li 1,2 and Libing Wang 1,2**


**Abstract:** Remote sensing reflectance (R*rs*), which is currently measured mainly using the abovewater approach, is the most crucial parameter in the remote sensing inversion of plateau inland water colors. It is very difficult to measure the R*rs* of plateau inland unmanned areas; thus, we provide a measurement solution using a micro-spectrometer. Currently, commercial micro-spectrometers are not factory calibrated for radiation, and thus, a radiometric calibration of the micro-spectrometer is an essential step. This article uses an Ocean Optics micro-spectrometer (STS-VIS) and a traditional water spectrometer (Trios) to simultaneously measure the irradiance and radiance of diffuse reflectance plates with different reflectance values for field calibration. The results show the following: (1) different fiber types have different calibration coefficients, and the integration time is determined according to the diameter of the fiber and the type of fiber, and (2) by comparing the simultaneous measurement results of STS-VIS with Trios, the mean absolute percentage difference (MAPD) of both reached 18.64% and 5.11% for Qinghai Lake and Golmud River, respectively, which are accurate R*rs* measurements of water bodies. The R*rs* of the Hoh Xil and Qarhan Salt Lake water bodies in unmanned areas of China was measured, and this was the first collection of in situ spectral information with a micro-spectrometer. This article shows that the micro-spectrometer can perform the in situ measurement of water R*rs* in unmanned inland areas. With this breakthrough in the radiometric performance of the micro-spectrometer, we are able to obtain more accurate remote sensing reflectance results of unmanned water bodies.

**Keywords:** field radiometric calibration; micro-spectrometer; remote sensing reflectance; plateau inland waters

#### **1. Introduction**

Remote sensing reflectance (sr−1) of water bodies is an important apparent optical quantity that can be detected in water color remote sensing and is one of the important parameters to quantify the spectral information of water bodies [1]. The field measurement of water remote sensing reflectance is a key step in the remote sensing of water color, and is widely used in the research and commercialization of satellite authenticity verification, water quality parameter inversion [2], cyanobacterial bloom [3] and black-odor water body identification [4]. The remote sensing of water color can be divided into the remote sensing of marine water bodies and the remote sensing of inland water bodies, with the remote sensing of marine water bodies being utilized to mainly observe the surface chlorophyll-a (Chl-a) concentration [5], suspended particular matter (SPM) concentration [6] and colored dissolved organic matter (CDOM) [7]. The remote sensing of inland water bodies has been used to estimate, in addition to the parameters mentioned above, total phosphorus (TP), total nitrogen (TN), the trophic state index (TSI) and, more recently, dissolved organic carbon in an eutrophic lake [8] and methane emissions in a lake [9].

**Citation:** Shi, J.; Shen, Q.; Yao, Y.; Zhang, F.; Li, J.; Wang, L. Field Radiometric Calibration of a Micro-Spectrometer Based on Remote Sensing of Plateau Inland Water Colors. *Appl. Sci.* **2023**, *13*, 2117. https://doi.org/10.3390/ app13042117

Academic Editor: Tung-Ching Su

Received: 15 November 2022 Revised: 3 February 2023 Accepted: 4 February 2023 Published: 7 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Current spectrometers dedicated to R*rs* based on water bodies include the Analytical Spectral Devices (ASD) series from Malvern Panalytical in the UK [10], the Profiling Reflectance Radiometer (PRR) series from Biospherical in the USA [11], the PAMSES series from Trios in Germany [12] and the Satlantic series from a subsidiary of SeaBird in the USA [13]. Measuring the R*rs* is an essential step in both marine remote sensing and inland water remote sensing. In oceanic remote sensing, the spectrometer is mainly carried on board a research vessel for walk-around measurements, and the same method is applied in inland lakes. There is a great deal of preparation before the measurement is carried out: firstly, the spectrometer is charged; secondly, the instruments are inventoried for transport; and finally, the spectrometer is deployed on board the research vessel. All of the above spectrometers can be utilized for the measurement of parameters in water bodies. However, they all are of a large size and weight and require a power supply, and thus, we cannot bring the spectrometer to lakes or rivers in unmanned areas such as the Tibetan Plateau. Therefore, there has been a gap in our knowledge on the water spectra of lakes or rivers in unmanned areas of the Tibetan Plateau. Furthermore, in inland water environments there are often unpopulated waters and waters that cannot be measured by using large vessels, and the above-mentioned water body R*rs* cannot be measured. Micro-spectrometers have the advantage of being small in size [14], and instead of spectrometers, they have been carried on board unmanned aerial vehicles or unmanned ships for R*rs* measurements to solve these problems. Nowadays, scholars often use imaging spectrometers for spectrometric measurements via unmanned aerial vehicles, and there are empirical linear-based [15] and look-up table methods [16] for radiometric calibrations of imaging spectrometers. The main applications are in vegetation and crop studies [17], and scholars have also used them to perform techniques for estimating water quality parameters [18,19]. However, imaging spectrometers tend to have a low signal-to-noise ratio for measuring the remote sensing spectra of water bodies. However, none of the current commercial micro-spectrometers are radiometrically calibrated; thus, we need to carry out a radiometric calibration of the micro-spectrometer.

The radiometric calibration of the micro-spectrometer can be performed by using a high-precision spectrometer and a micro-spectrometer to simultaneously measure highprecision stable irradiance and irradiance sources and to solve for the conversion relationship between the digital number (DN) recorded by the micro-spectrometer and the radiance and irradiance values received by the micro-spectrometer. The main radiometric calibration methods used thus far for remote sensing sensors are laboratory integrating sphere calibration [20–22], on-station calibration [23], diffuse reflectance plate calibration [24,25], etc. Laboratory integrating sphere calibration is mainly a high-precision test of parameters, such as the remote sensor electronic gain and bias, central wavelength and bandwidth of each channel, signal-to-noise ratio, spatial resolution and spectral response function using relevant laboratory equipment, which is traceable with a low-temperature absolute radiometer or a standard blackbody uniform radiation quantification standard. On-station calibrations use a calibration system mounted on a satellite platform to periodically monitor changes in the radiation response of the satellite during its orbital operation. Diffuse reflectance plate calibration uses a diffuse reflectance plate and a spectrometer probe for simultaneous fixed-point continuous observation. Concerning the radiometric calibration methods for remote sensing sensors, only laboratory integrating sphere calibration and field diffuse reflectance plate calibration methods can meet the requirements for the radiometric calibration of micro-spectrometers. The field diffuse reflector calibration method is simple and easy to operate, whereas the laboratory integrating sphere calibration requires a lot of time and money, and our group does not meet the conditions for conducting laboratory integrating sphere calibration. Therefore, the field diffuse reflector calibration method is chosen for the radiation calibration of the micro-spectrometer.

The main purpose of this article is a preliminary evaluation of the feasibility and accuracy of micro-spectrometer measurements of the R*rs* of plateau inland water colors. Specifically: (1) we calculated the micro-spectrometer radiometric calibration coefficients

using different fiber types; (2) we quantified the accuracy of measuring the R*rs* of clear and turbid waters in plateau inland water bodies and (3) we provided a case study of the application of the micro-spectrometer to R*rs* of plateau inland water bodies in unmanned areas.

#### **2. Data and Methods**

#### *2.1. Equipment and Methods*

The equipment used in the micro-spectrometer radiometric calibration is shown in Figure 1 and includes the following: the micro-spectrometer (STS-VIS, Ocean Optics company, Orlando, FL, USA), an optical fiber, a cosine receiver, a spectrometer (Trios, TriOS Mess, Rastede, Germany), a diffuse reflectance plate with different reflectance values and a target cloth, a micro-spectrometer linked to the optical fiber to measure radiance and a micro-spectrometer linked to the optical fiber and cosine receiver to measure irradiance. The diffuse reflectance plate and the target cloth provide stable and varying values of irradiance and radiance. The measured spectral ranges and spectral resolution of Trios and STS-VIS are shown in Table 1. The Trios can measure spectra from 200 nm up to 1100 nm with a spectral resolution of 3.3 nm, whereas the STS-VIS can measure spectra from 350 to 810 nm with a spectral resolution of 2.2 nm. There are three common methods of field measurements of remotely-sensed reflectance ratios in water bodies: the in-water approach [26], the above-water approach [27] and the skylight-blocked approach [28]. The in-water approach generally uses more expensive instruments, which are complex to operate and deploy, and are subject to a certain amount of self-shadowing and uncertainty in the data results. This method can generally only be used in water depths greater than 10 m, and thus, this method is very widely used in pelagic Class I waters, but rarely in inland waters. The skylight-blocked approach uses a mask to directly block the sky light from the field of view of the observation sensor, thus enabling the direct measurement of the off-water irradiance of a water body. This method avoids the geometrical errors in observation caused by the complexity of the field and the uncertainties associated with the skylight rejection method. However, it requires a continuous power supply, and therefore, cannot measure water spectra in the unmanned areas of the Tibetan Plateau. The above-water method has the advantages of a simple field operation and the low cost of field experiments, and is currently the most commonly used measurement method in the study of the spectral properties of water bodies; however, it is affected by water and weather conditions. The details of this method of measurement can be found in Figure A1 in Appendix A. The main zenith angle conditions met by the instrument are that the downward irradiance sensor should point vertically towards the sky, the radiance sensor for measuring sky light should be equal to 50◦ and the radiance sensor for measuring sky light should be equal to 140◦. The main azimuth angle conditions met by the instrument are that the angle between the left and right sunlight should be equal to 45◦. This helps avoid sun glint. This method is the easiest water body measurement method to implement in unmanned lakes on a plateau; therefore, the above-water method was selected for this paper to carry out spectral measurements of water bodies in unmanned areas on the Tibetan Plateau, and the formula for measuring R*rs* is as follows (Equation (1)):

$$\mathcal{R}\_{rs}(\lambda) = \frac{L\_{\text{w}}(\lambda)}{E\_{\text{d}}(\lambda)} = \frac{L\_{\text{u}}(\lambda) - r\_{\text{sky}}(\lambda)L\_{\text{sky}}(\lambda)}{E\_{\text{d}}(\lambda)}\tag{1}$$

R*rs* is equal to the ratio of the off-water radiance to the downward irradiance, and in the above formula, *Lu*(*λ*) is the total off-water radiance, *rsky*(*λ*) can be derived from the look-up table [29], *Lsky*(*λ*) is the skylight radiance and *E*d(*λ*) is the downward irradiance. They are measured simultaneously using a micro-spectrometer, using (a)-1 for sky radiance, (a)-2 for total off-water radiance and (a)-3 for downward irradiance, so as to avoid the effects of weather variations.

**Figure 1.** (**a**) Ocean Optics micro-spectrometers (STS-VIS); (a)-1 sensor for measuring sky radiance; (a)-2 sensor for measuring water radiance; (a)-3 sensor for measuring downward irradiance. (**b**) Ocean Optics and Hygirel fiber. (**c**) Cosine Correctors. (**d**) Trios. (**e**) Diffuse reflectance plates and target cloth (99%, 95%, 30%, 20%, 5%, and 2%).

**Table 1.** Parameter table of Trios and STS-VIS.


The procedure for field calibration is to use the Trios sensors and micro-spectrometer to simultaneously measure diffuse reflectors with different reflectance values against a target cloth, and to obtain the simultaneously measured Trios radiance, irradiance and microspectra with their corresponding DN values. The formula for field-calibrated radiance is presented in Equation (2) and for field-calibrated irradiance in Equation (3).

$$L\_{T\text{rise}} = \text{gain} \times DN\_{STS-VIS} + offset \tag{2}$$

$$Ed\_{Tiro} = gain \times DN\_{STS-VIS} + offset \tag{3}$$

In Formula (2) above, *LTrios* is the radiance of the different diffuse reflection plates measured by Trios, and *DNSTS-VIS* is the DN of the different diffuse reflection plates measured by the micro-spectrometer. *EdTrios* in Formula (3) is the irradiance of the different diffuse reflectors measured by Trios and *DNSTS*−*VIS* is the the DN of the different diffuse reflection plates measured by the micro-spectrometer.

We used the micro-spectrometer to obtain the R*rs* of plateau lakes as follows: firstly, the DN values obtained with the micro-spectrometers (a)-1, (a)-2 and (a)-3 are converted into radiance and irradiance according to Formulas (2) and (3), and then the R*rs* of the micro-spectrometer is calculated according to Formula (1).

The micro-spectrometer needs to be physically connected to the micro-computer in the unmanned ship using a USB cable; additionally, we also need to use professional software to identify the sensor, and only after the identification is passed can the data acquisition parameters be set. The main parameters to be set are integration time, multiple scan averaging, sliding average width and other parameters. It is also necessary to set the address where the spectrum file is to be saved. Once the settings are complete, the micro-spectrometer can be used for water spectrum acquisition.

#### *2.2. Data Processing and Evaluation Metrics*

We used the Oceanview 6.0 software for hardware control and we needed to set the integration time before collection. The integration times of the micro-spectrometer for measuring radiance and measuring irradiance are described in detail in Section 3.1. The micro-spectrometer and Trios both measure different floating-point wavelength data; therefore, we needed to resample the Trios and STS-VIS spectrometer data to the same integer wavelength for field radiation calibration. The cubic spline interpolation method has a small error (we have also compared the results of the other three resampling methods plotted in Figure A2 Appendix A). Thus, we used this method and obtained a result that is very small, as shown in Figure 2.

**Figure 2.** Comparison of raw and resampled data.

The following statistics were used to evaluate the spectral results, including the correlation coefficient (r), bias and mean absolute percentage difference (MAPD).

$$\mathbf{r} = \frac{\sum\_{i=1}^{n} \left(\mathbf{S}\_{i,1} - \overline{\mathbf{S}\_{1}}\right) \left(\mathbf{S}\_{i,2} - \overline{\mathbf{S}\_{2}}\right)}{\sqrt{\sum\_{i=1}^{n} \left(\mathbf{S}\_{i,1} - \overline{\mathbf{S}\_{1}}\right)^{2} + \sum\_{i=1}^{n} \left(\mathbf{S}\_{i,2} - \overline{\mathbf{S}\_{2}}\right)^{2}}} \tag{4}$$

$$\text{bias} = (\text{S}\_{\text{i},1} - \text{S}\_{\text{i},2}) / \text{S}\_{\text{i},2} \times 100\text{\%} \tag{5}$$

$$\text{MAPD} = \frac{1}{N} \sum\_{i=1}^{N} |(\text{S}\_{i,1} - \text{S}\_{i,2}) / \text{S}\_{i,2}| \times 100\% \tag{6}$$

where Si,1 and Si,2 denote the DN values corresponding to the measured diffuse reflectance plates, the gain and offset observed under different conditions and the R*rs* measure with the Trios and micro-spectrometer.

In this article, we compare the effects of different fibers on the gain and offset, which is followed by a quantitative analysis to measure the spectra of different diffuse reflectors with clear water and turbid water, and finally, we perform a first application demonstration in an unmanned area.

#### **3. Results**

We first determined the important acquisition parameter integration times for the micro-spectrometer under sunny conditions, then compared the gain and offset of two fiber types using a Hygirel fiber-optic micro-spectrometer (detailed parameters of the two types of micro-spectrometers are given in Section 3.2) and Trios under sunny conditions, and finally, we evaluated the simultaneous measurements of 95% and 20% for the diffuse reflectance plates, and 5% and 2% for the target cloth.

#### *3.1. Effect of Integration Time on Spectral Acquisition of Water*

We all know that when a spectrometer performs electro-optical conversion it is mainly influenced by three main parameters [30]: the integration time, aperture and temperature of the micro-spectrometer. The longer the integration time, the more energy is collected by the micro-spectrometer; the larger the aperture, the more energy is collected by the micro-spectrometer; the temperature of the micro-spectrometer creates thermal noise; and different temperatures bring about different thermal noises. When measuring in the field, we pre-warmed the micro-spectrometer for 60 min to keep it below normal operating temperature to avoid the effects of thermal noise caused by temperature. For dealing with the effects of the aperture, micro-spectrometers are often used with a fixed fiber; thus, we use a fixed fiber diameter to control the effect of the aperture on the spectrum.

The signal-to-noise ratio (SNR) of a spectrometer is a very important metric when using micro-spectrometers for the spectral acquisition of water bodies. The current international water color satellites MODIS Terra/Aqua (spatial resolution: 250 m and 500 m), MERIS Envisat (spatial resolution: 300 m), VIIRS (spatial resolution: 750 m) and Sentinel-3A/B (spatial resolution: 300 m) all have SNRs greater than 1000 for measurements that occur in the ocean [31]; measurements in inland lakes mainly use Landsat OLI (spatial resolution: 30 m) HJ-1 (spatial resolution: 30 m), etc., all of which have SNRs greater than 100 [32]. The use of Landsat as a data source for lake water observations on the Qinghai–Tibet Plateau requires that the signal-to-noise ratio of the micro-spectrometer be greater than 100 for application in field measurements.

We believe that in the electro-optical conversion of a micro-spectrometer, the noise of the micro-spectrometer is caused by the dark current; thus, we calculated the SNRs as the total signal ratio over the dark current signal (Equation (7)). The SNRs are the ratio of the useful signal to the noisy signal in the total signal. The higher the signal-to-noise ratio, the better.

$$\text{SNRs} = \frac{L\_{\text{total}}}{L\_{\text{darkcurrent}}} \tag{7}$$

We first measured the dark current at different integration times after warming up the micro-spectrometer for 60 min (Figure 3a) and found that the dark current increased in a logarithmic fashion with integration time. The micro-spectrometer is considered to operate at a stable maximum temperature after 60 min of operation. As we used a reference plate with a maximum reflectance of 95% for field calibrations, but the DN range of the micro-spectrometer is 0–65,535, we first needed to ensure that we did not exceed the DN limit when collecting 95% of the reference plat and that our integration time was not too long, as the state of the water body changes rapidly. We should ideally be able to measure a water body spectrum within 2 s, but in a full water body spectrum measurement, the effects of waves, boat wake, etc., need to be taken into account. We generally took several measurements to average the water body spectra for field measurements.

**Figure 3.** Dark current diagram for different integration times (**a**), SNR diagram for radiance measurement of DN values (**b**), SNR diagram for irradiance measurement of DN values (**c**).

Two schemes were used to measure 95% of the reference plate, the first using a microspectrometer and optical fiber to measure the DN values corresponding to radiance, and the second using a micro-spectrometer, optical fiber and cosine receiver to measure the DN values corresponding to irradiance. The SNR plot of the DN values corresponding to radiance (Figure 3b) shows that the integration time exceeds the SNRs ~100:1 at 500 ms, but that there are higher SNRs at 1000 ms and 1500 ms; therefore, we chose 1000 ms as the integration time for the DN values of radiance. The SNR diagram for the irradiance

corresponding to the DN value (Figure 3c) shows that the integration time exceeds the SNRs ~100:1 at 500 ms, and has similar SNRs at 500 ms, 1000 ms and 1500 ms; thus, we chose 500 ms as the integration time for the DN value of the irradiance.

#### *3.2. Comparison of Gain and Offset Via Two Fiber Types*

We performed field radiometric calibrations using different fiber types (Table 2). The two different fibers are from Ocean Optics (OC) and Hygirel (HY). The OC fiber is manufactured by Ocean Optics and has a spectral range of 200–1000 nm, and its core diameter is 600 μm and its length is 2 m; HY is manufactured by HAIJILEKEJI and has a spectral range of 200–1100 nm, and its core diameter is 1000 μm and its length is 2 m. Figure 4 shows the results of the field radiometric calibration coefficients gain and offset, which we evaluated using *r* for least squares regression coefficient accuracy. We found that (a)-1, (a)-2 and (a)-3 have a large variation in gain coefficients across fiber types, but there is a linear relationship, with r reaching 0.99 for Ocean Optics' gain and Hygirel's gain. The gain data are again noticeably jittery at 760 nm, which we attribute to the different integration times measured by Trios and STS-VIS at the time of measurement and to the variation in atmospheric gases. We post-processed the gain data by means of smoothing in subsequent processing. We also found for (a)-1, (a)-2 and (a)-3 that the offset coefficients vary considerably across fiber types and do not have any linear relationship, with *r* reaching 0.99 for offset obtained on sunny days. We found that the overall offset coefficient of HY is much larger than that of OC, which we believe is due to the different core diameters of the two fibers and the different amounts of sunlight energy entering the micro-spectrometer at the same integration time for HY and OC.

**Table 2.** Parameter table of different fiber types.


**Figure 4.** Graphs of gain and offset of micro-spectrometer connection with different fiber types.

We found that for a fixed integration time and different fiber types, there are different field radiometric calibration coefficients for gain and offset. We have made the gain and offset for sunny days for Ocean Optics and Hygirel fiber types publicly available for use so that a field radiometric calibration can be performed before use if a more accurate R*rs* is required. Different fiber types mainly affect the calibration coefficients of gain and offset, and the field radiometric calibration coefficients of gain and offset for different fiber types can be found at https://github.com/765302995/FRC\_Micro-spec (accessed on 15 November 2022).

#### *3.3. Comparison of Radiance/Irradiance Measured in Four Diffuse Reflector/Target Cloths*

When measuring four different reference plates using Trios and STS-VIS, we used the vertical measurement method because the reference plates are diffuse light sources. The results of the simultaneous measurements of the diffuse reflectance plate/target cloth using the micro-spectrometer and Trios under sunny conditions are shown in Figure 5. It can be seen that the micro-spectrometer and Trios have a large error at 760 nm in the oxygen absorption. In 95% of the diffuse reflectance plate results, the bias error of R*rs* measured by the micro-spectrometer and Trios ranged from 0.5 to 1.6%; in 20% of the diffuse reflectance plate results, the bias error of R*rs* measured by the micro-spectrometer and Trios ranged from −4.15 to −7.44%; in 5% of the target cloth results, the bias error of R*rs* measured by the micro-spectrometer and Trios ranged from −1.87 to −12.63%; and in 2% of the target cloth results, the bias error of R*rs* measured by the micro-spectrometer and Trios ranged from −22.9 to −47.41%. The lower the reflectance, the poorer the signal-to-noise ratio of the micro-spectrometer and, therefore, the greater the error in the measured reflectance.

**Figure 5.** Diffuse reflectance plate results with bias plots for 95%, 20%, 5% and 2% measurements using Trios and micro-spectrometer.

#### **4. Discussion**

The R*rs* accuracy and availability of micro-spectrometer measurements are evaluated here. We firstly evaluate the simultaneous measurements from Qinghai Lake and Golmud River waters using the micro-spectrometer and Trios, and finally, we carry the calibrated micro-spectrometer in the unmanned boat to perform the first measurement of the water body in plateau unmanned lakes.

The water body has a clear bidirectional reflectance distribution, so when we measure the water body, the measurement angle needs to meet the zenith angle of the sun at 40◦ and the solar azimuth angle at 135◦ or 45◦ to eliminate the influence of sun glint. We use the above-water approach to measure water bodies from the same angles as described above.

#### *4.1. Comparison of Rrs Measured in Two Plateau Inland Water Body Types*

We carried out and compared R*rs* measurements by using the micro-spectrometer and Trios on Qinghai Lake and in the Golmud city River, and the comparison results are shown in Figure 6: there is an overestimation of R*rs* when measured using the micro-spectrometer in Qinghai Lake (Figure 6a), where the r of R*rs* measured by the micro-spectrometer and Trios is 0.99; the MAPD reaches 18.64%; and the bias is at 400–450 nm, 590–610 nm and 700–750 nm, all which exceeded 20%. We can see from the in situ photo (a) that this inland lake of the Qinghai–Tibet Plateau is a clean type. The cleaner the water body, the lower the reflectivity of the water body and the lower the SNR of the micro-spectrometer, meaning that the MPAD is larger. However, we can see by the shape of the water body spectrum that we obtained a complete measurement of the spectral shape of the clean water body. There is an underestimation of R*rs* in the Golmud River when measured with the microspectrometer (Figure 6b); the micro-spectrometer and Trios reached an r of 0.99 for R*rs*, an MAPD of 5.11% and a bias of over 10% at 700–750 nm. We can see from the in situ photo (b) that this inland river of the Qinghai–Tibet Plateau has a turbid type. The more turbid the water body, the higher the reflectance of the water body, and the higher the SNR of the micro-spectrometer. We can see by the shape of the water body spectrum that we obtained a complete measurement of the spectral shape of the clean water body. For the turbid water body, not only is the water body shape the same, but the value of the water body spectrum is also very similar; thus, the water body spectrum measured by the micro-spectrometer can be used for quantitative water body parameter research.

**Figure 6.** Micro-spectrometer measurements of clear water and turbid water R*rs* and bias, and in situ photos of (**a**) spectrum of Qinghai Lake, (**b**) spectrum of Golmud city River.

#### *4.2. Unmanned Area Applications*

In field experiments that are carried out in plateau inland waters, many lakes and rivers are in unmanned areas and, therefore, cannot be measured with spectrometers such as Trios. Instead, we can use micro-spectrometers to measure the spectra of unmanned water bodies. The calibrated micro-spectrometer is carried on unmanned boats for the measurements. We designed the route according to the GPS of the unmanned boat and in accordance with the experiment time to ensure that the three probes of the micro-spectrometer met the solar zenith angle of 40◦ and the solar azimuth angle of 135◦, or met the 45◦ water measurement angle requirements. When the unmanned boat arrives at the predetermined point, the unmanned boat will stop and sway according to the water waves; thus, we arrived at the point moving slowly in order to obtain the measurement of the point walking water

body spectrum. Here, the water spectra of the Qarhan Salt Lake (DBX,NHBX) and the Hoh Xil Salt Lake (KKXL) were measured, as shown in Figure 7. From the in situ photos, it can be seen that DBX3 is the most turbid, with the highest water reflectance at 590 nm at 0.055 sr<sup>−</sup>1, and at 500–700 nm exceeding 0.03 sr−1, resulting in a heavy greyish-brown color. DBX8 is also more turbid, with the highest water reflectance at 600 nm at 0.039 sr<sup>−</sup>1, and 550–720 nm. The reflectance of DBX6 reached the highest value at 600 nm at 0.02 sr<sup>−</sup>1, and had a grey-green color; the reflectance of KKXL1 and KKXL2 reached the highest value at 580 nm, and the reflectance was more stable at 400–580 nm, and thus, the color of the water body was blue.

**Figure 7.** Micro-spectrometer measurements of unmanned areas of Qarhan and Hoh Xil Salt Lake R*rs* and in situ photos.

#### **5. Conclusions**

Radiation calibration is needed to carry out irradiance and radiance measurements with micro-spectrometers. This work introduces the method of field radiation calibration for micro-spectrometers, analyzes the effect of different types of connected fibers on the radiation calibration, and analyzes the performance of micro-spectrometers in measuring R*rs* in inland waters. The main findings are as follows: (1) Different fiber types mainly affect the calibration coefficients of gain and offset, and the field radiometric calibration coefficients of gain and offset for different fiber types can be found at https://github.com/ 765302995/FRC\_Micro-spec (accessed on 15 November 2022). (2) The MAPD of the microspectrometer reached 18.64% and 5.11% for clear water and turbid water, respectively, and the water body R*rs* values of unmanned plateau lakes were obtained for the first time using the micro-spectrometer. This article shows that the micro-spectrometer can meet the requirements for field measurements of R*rs* of water bodies in inland unmanned areas, and with this breakthrough in the radiation performance of the micro-spectrometer, we can obtain more accurate R*rs* measurements of water bodies in unmanned areas.

**Author Contributions:** Conceptualization J.S. and Q.S.; project administration, Y.Y. and J.L.; writing original draft preparation, J.S.; writing—review and editing, Q.S.; validation, L.W. and F.Z.; supervision, Q.S.; methodology, J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the National Key Research and Development Program of China under grant number 2021YFB3901101.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank Qian Shen from the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences for help with writing. We are also thankful to all anonymous reviewers for their constructive comments provided on the study. Thanks to Ge Yongquan and Gao Hangyu from Shenyang Jianzhu University for their help in obtaining the data.

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

#### **Appendix A**

**Figure A1.** Conceptual view of the unmanned ship carrying a micro-spectrometer to collect water body spectrum; the side view details the angle of measurement of the zenith angle of the three sensors followed for spectrum collection and the top view details the angle of measurement of the azimuth of the three sensors for spectrum collection.

**Figure A2.** Comparison figure of different spectral resampling methods, (**a**) zero spline interpolation; (**b**) linear spline interpolation; (**c**) quadratic spline interpolation; (**d**) cubic spline interpolation.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
