A Generic Method to Derive Coastal Bathymetry From Satellite Photogrammetry for Tsunami Hazard Assessment

In this study, we use Pleiades satellite stereo imagery to map shallow‐water bathymetry specifically for the purpose of assessing tsunami hazard in coastal regions. In order to calculate the refraction correction factor, we independently measure the refractive indices of seawater mixed with a variety of phytoplankton, which was previously unquantified. We find that the refractive index of seawater increases with cell density of phytoplankton, ranging from 1.3410 to 1.3425. Through rigorous error analysis, we demonstrate that the change in refraction due to heterogeneous distribution of phytoplankton is negligible, suggesting that the technique is applicable to the global ocean. Using the Penghu Islands as a case example, we show that photogrammetric bathymetry is able to reveal tsunami impacts more accurately and reasonably than the commonly used General Bathymetric Chart of the Oceans data, highlighting the importance of high‐resolution and high‐accuracy bathymetry for tsunami hazard assessment.

situ depth data. Empirical method does not usually require atmospheric correction, but does need a large amount of previously existing in situ data. The SDB method is efficient but less accurate than acoustic and LiDAR surveys as it depends largely on water quality and atmospheric calibration of multispectral images. Recently, photogrammetry has been used to map bathymetry of shallow waters (Casella et al., 2017;Dietrich, 2017;Hodúl et al., 2018;Maas, 2015;Mandlburger, 2019), reaching an accuracy of ∼ 1 m down to 10 m depths (Hodúl et al., 2018), providing a new means of obtaining high-resolution and high-accuracy bathymetric data in coastal regions.
Stimulated by Hodúl et al. (2018)'s study, and other published studies (Dietrich, 2017;Maas, 2015;Mandlburger, 2019) that used photogrammetry to derive bathymetry, we are interested to see to what extent this new data set, in comparison to the widely used General Bathymetric Chart of the Oceans (GEBCO, Weatherall et al., 2015) data, can contribute to improving our ability of tsunami hazard assessment. In this study, we use Pleiades stereo imagery to derive coastal bathymetry for such purpose. In order to measure actual water depths, through-water refraction correction has to be taken into account, in which the refractive index of the media is a key variable. Although using satellite photogrammetry to retrieve bathymetry has been well established in recent studies (e.g., Hodúl et al., 2018;Mandlburger, 2019), none of them quantified the potential effect of phytoplankton that are abundant throughout the ocean (Leblanc et al., 2012;O'Brien et al., 2013) on refraction. Here, we measure for the first time the refractive indices of seawater with different species and densities of phytoplankton and quantify the associated errors. We then use the newly derived bathymetry to calculate local tsunami wavefield generated by a synthetic earthquake on the Manila trench (Megawati et al., 2009). We select the Penghu Islands (Figure 1), located in the northeastern corner of the South China Sea, as a test area because it is densely populated and threatened by tsunamigenic earthquakes that may occur along the Manila subduction zone. A recent study by Lu et al. (2019) has revealed geological records of three possible tsunami events since 400 AD on the islands with a recurrence interval of 400-500 years, highlighting the importance of evaluating the tsunami risks of the region. We present the first case study of tsunami hazard assessment using photogrammetric bathymetry, showing the improvement in hazard assessment. Meanwhile, this paper provides the fundamental parameter-the refractive index of natural seawater-for the derivation of photogrammetric bathymetry, and demonstrates through rigorous error analysis that heterogeneous distribution of phytoplankton in the oceans is negligible, indicating that the technique is generic can be applied to the global ocean.

Pleiades Imagery and Processing Methods
We purchased one pair of Pleiades bistereo imagery that was acquired on November 17, 2017 (listed in Table S1 in Supporting Information S1). The product consists of one 50-cm panchromatic band and four 2-m multispectral (red, green, blue, and near-infrared) bands. The multispectral bands are used to derive photogrammetric bathymetry.

Photogrammetric Process
The quality of photogrammetric process is determined by the goodness of image features. Usually, panchromatic images are used to extract topography given their high spatial resolution, but since phytoplankton-rich water has a very different spectral reflectance to soil (Han, 1997), panchromatic images may not be suitable for extracting bathymetry. It is worth analyzing the characteristics of each Pleiades band to explore which band can well represent seafloor features (Legleiter et al., 2009). Note that sun glint is a common phenomenon in remotely sensed data in shallow-water environments under conditions of clear skies and nonflat water surfaces, which needs to be removed from the images beforehand. This can be done using a simple linear correction method proposed by Hedley et al. (2005). Figure 2 presents an example of seafloor in different bands, including red (600-720 nm), green (490-610 nm), blue (430-550 nm), near-infrared (750-950 nm), and panchromatic (480-830 nm). Red and near-infrared bands have a very low reflectance beneath the water, resulting from a combined effect of pigment and water absorption. Given its extremely low reflectance of water, the near-infrared band can be used to separate land and water areas. Green and blue bands show seafloor features well; this is because phytoplankton have a low absorption rate of light around 550 nm (Gitelson, 1992), leaving a broad reflectance maximum in this wavelength range. Because the panchromatic band covers a much wider spectrum (480-830 nm), the features are less clear in the image, caused by the loss of reflectance between 400-500 nm, 650-690 nm, and 750-900 nm. To further analyze the applicability of each band, we used all of them to extract bathymetry. The Leica Photogrammetry Suite (LPS) built in ERDAS was used for the photogrammetric extraction (Zhou et al., 2015). A standard bundle block adjustment was performed using the rational function models (RFMs) (Fraser et al., 2006;Tao & Hu, 2002). We checked the tie points that were automatically generated, and removed those that were generated on waves or floating objects. A dense point cloud was produced by a pairwise pixel-by-pixel matching procedure with a search window size of 7 pixels × 7 pixels in the new automatic terrain extraction module (eATE) using a Dense Image Matching algorithm (Hirschmuller, 2007;Rothermel et al., 2012). As shown in Figure 2, the green band yields the densest point cloud owing to the peak reflectance of seawater around 550 nm. We took this point cloud for postprocessing to generate bathymetric data. The near-infrared image was used to mask out land areas before implementing refraction correction.

Refraction Correction
When the ray of light passes the water/air interface, it bends, that is, refraction occurs, as shown in Figure 3. Traditional photogrammetric process does not take refraction into account, so the depth extracted directly from photogrammetry is apparent water depth ( * ) that needs to be corrected (Mandlburger, 2019). Note that the topography derived from satellite images is normally referenced to the WGS-84 ellipsoid (i.e., ellipsoidal heights). Before refraction correction, we need to convert ellipsoidal heights to water depths by adding a constant to make sea surface (or seaside for approximation) elevations zero.
In the geometry as shown in Figure 3, the distance between the two intersections (A and B) of light and water surface can be calculated by (Maas, 2015): Diagram illustrating the geometry of through-water photogrammetry. Q is a point on the seafloor, and P is its apparent location. A and B are the intersections of light and water surface. * is the apparent water depth, is the actual water depth, and are the incidence angles of the left-view and right-view images, respectively, and are the corresponding refraction angles, and are the view azimuth angles, and Ω and Ω cross-track and in-track view angles.
where * is the apparent water depth, is the actual water depth, and are the incidence angles of the left-view and right-view images, respectively, and are the corresponding refraction angles, and and are the view azimuth angles.
Therefore, can be derived from * via The view azimuth angle ( ) of each image can be calculated by its cross-track ( Ω ) and in-track ( Ω ) view angles: and are related by the refractive index of water (Snell's Law): Combining Equations 3, 6, and 7, the actual water depth can be calculated by Equation 8 serves as the basis for deriving photogrammetric bathymetry, in which the refractive index of seawater is a key parameter. may vary in the global ocean due to heterogeneous distribution of phytoplankton, which was poorly understood. In this study, we put a firm constraint on through laboratory measurements.

Measurements of Refractive Index of Seawater With Phytoplankton
We measured the refractive indices of seawater mixed with a range of different species of phytoplankton that represent the major groups in the global ocean, including cyanobacteria Synechococcus, chlorophyte Chlamydomonas concordia, dinoflagellate Thoracosphaera heimii, diatom Phaeodactylum tricornutum, and coccolithophore Emiliania huxleyi (see Figure S1 in Supporting Information S1). The reason for doing this is that phytoplankton are unevenly distributed in the global ocean, both horizontally and vertically, and we must evaluate their effects on the refractive index of natural seawater in order to make sure that this technique is generic and applicable worldwide.
The phytoplankton were cultured in f/2 medium (Guillard & Ryther, 1962). They were grown at a constant temperature of 20 ± 1 • C and illuminated with 100 mmol photons m −2 S −1 on a light to dark cycle of 12:12 hr (Zhang & Rickaby, 2020). The cells were harvested at exponential phase by centrifugation at 10,000 g for 10 min. They were then resuspended in filtered natural seawater (0.22m) and diluted to different cell densities. The number of cells in the samples were determined by using an electronic particle counter Coulter Counter (W. Wilson et al., 2019). Refractive indices of seawater with phytoplankton at different cell densities were then measured by using an Abbe refractometer ( Figure S1 in Supporting Information S1).

Converting Photogrammetric Elevations to Bathymetry
The point cloud extracted from the green band was used in the derivation of bathymetry. The elevations calculated based solely on satellite orbits are in the WGS-84 ellipsoid reference system, so the first step is to convert ellipsoidal heights to water depths. To do so, we manually selected 19,536 waterline points (their distribution is shown in Figure 1) by visual looking at the true-color image to compute the conversion factor (i.e., the mean shoreline height). We subtracted the mean value (22.14 m, as indicated in Figure S2a in Supporting Information S1) from the elevations of the original point cloud to get the apparent water depths. After waterline correction, we separated land and water points using the mask generated from the near-infrared image, and applied Equation 8 to water points for refraction correction. The parameters in Equation 8 are listed in Table S1 in Supporting Information S1. After refraction correction, the bathymetric data are referenced to the local sea level at the time of image acquisition. In order to make the data compatible with other data sets, we converted this "arbitrary" vertical datum to the same datum as the global topographic data. The ALOS Global Digital Surface Model (referred to as the "ALOS DEM" hereafter) was used for the conversion. We calculated the differences between the Pleiades and ALOS elevations of the waterline points, and obtained a mean value of 2.37 m ( Figure S2b in Supporting Information S1). We applied this value to bring the Pleiades elevation points into the same vertical datum as the ALOS DEM and the EGM96 Geopotential Model. We then gridded the point cloud using a Triangulated Irregular Network (TIN) (Van Kreveld, 1996). We set 500 m as the maximum length of the triangles to delineate the TINs in order to minimize interpolation errors (Figure 4).

Refractive Index of Seawater
For calibration of the refractometer, we first measured the refractive index of Milli-Q water ( = 1.3345 ) and seawater ( = 1.3405 ). The results are listed in Table S2 in Supporting Information S1. The refractive index of seawater increases slightly with the density of phytoplankton, from 1.3410 (low cell density) to 1.3425 (ultra-high cell density in natural ocean). The changes in the refractive index of natural seawater (on the order of 0.001) due to heterogeneous distribution of phytoplankton can be ignored without causing significant errors in the derivation of bathymetry (discussed in Section 4), so a constant refractive index ( = 1.34 ) is used in our calculations. This value can also be applied to the global ocean.

Bathymetry Derived From Pleiades Stereo Images
The spatial resolution of the gridded DEM (Figure 4c) was set to be 20 m given the average point density of the water area, ∼ 0.0025 points/ m 2 , that is, 1 point/400 m 2 . Although the point cloud for the land area is much denser, we used 20 m to obtain a uniform DEM. For comparison, we plotted the GEBCO data (500 m resolution) in Figure 4d. GEBCO data have been widely used in tsunami-wave field calculation, for example, Fujii and Satake (2008); Heidarzadeh et al. (2016); Watada et al. (2014), but it may be too coarse to represent the near-shore bathymetry, as shown in our case, as well as other recent studies (Heidarzadeh et al., 2019;Kulikov et al., 2016). Photogrammetry provides a new means of mapping high-resolution coastal bathymetry, and it works well for clear waters with depths of 0-20 m. With increasing depth and turbidity, the reflectance of seawater decreases, making photogrammetric process challenging.

Error Analysis of Photogrammetric Bathymetry
Errors in the incidence and view angles, and the refractive index will result in errors in the measured bathymetry. We evaluate the contribution of each component to the final error (for simplicity, we assume that the error in the apparent water depth is independent of the incidence and view angles), with the focus on the influence of the observed changes in the refraction index of seawater on the calculation of bathymetry.
Since we do not have in situ measurements for calibration, we obtained an independent set of bathymetry using NASA's Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) space-based altimeter for comparison with our photogrammetry-derived results. It has been tested that the ICESat-2 derived bathymetry can reach an accuracy of 0.4-0.6 m (Parrish et al., 2019). Following the routine processing procedure of ICESat-2 data as outlined in many studies such as Babbel et al. (2021); Thomas et al. (2021);Hsu et al. (2021);and Parrish et al. (2019), we computed ICESat-2 bathymetry for our study area. The difference between the photogrammetric and ICESat-2 bathymetry is −0.4 ± 2.1 m (see Figure S4 in Supporting Information S1), suggesting that the accuracy of the photogrammetric-derived bathymetry is 1-2 m, consistent with the accuracy ( ∼ 1.2 m) estimated by Hodúl et al. (2018) where they have in situ measurements for verification.

Application of High-Resolution Coastal Bathymetry in Tsunami Modeling
We use a simple example to show to what extent the new coastal bathymetry can improve our ability of assessing tsunami hazard. We simulate a tsunami scenario generated by a synthetic earthquake with 9.3 occurring along the Manila Trench (Megawati et al., 2009). The setup of the computational grids is shown in Figure 5a. The details of tsunami modeling can be found in the Supporting Information S1. The inundation map calculated from high-resolution bathymetry ( Figure 5) shows that the low-elevation area along the southern coastline of Penghu would experience the worst tsunami impact. Some key infrastructures, including desalination plant and harbors could expect maximum tsunami wave heights of 3-5 m. Coastal regions in the east and northwest would suffer minor tsunami impact with only limited areas being inundated. In contrast, the inundation map calculated based on the GEBCO data ( Figure 5c) is not able to reveal the inundation area in a reasonable way. Large chunks of land mass that appear in the region are supposed to be ocean. Shallow waters, especially in near-shore regions with water depths of shallower than 20 m, cannot be clearly delineated by GEBCO data, creating challenges for tsunami hazard evaluation as features of tsunami waves, including their height, speed, and direction are strongly affected by the location and orientation of coastline, and water depth, especially in near-shore regions (e.g., Schambach et al., 2019).
The importance of near-shore bathymetry has also been highlighted in previous studies. By comparing simulations based on multibeam and GEBCO data, Xie et al. (2019) pointed out that the tsunami wave height is incredibly sensitive to the resolution of bathymetry at water depths of shallower than 50 m. Even small variations in the near-shore water depths could affect the behavior of tsunamis. For example, the existence of fringing coral reefs has a debatable role in increasing or decreasing tsunami wave heights (Dilmen et al., 2018). Such detailed information can only be resolved by high-resolution and high-accuracy bathymetry. Accurate water depth is a critical prerequisite for accurate tsunami hazard assessment, effective evacuation planning, and appropriate response during a real event. Unfortunately, high-resolution bathymetric data rarely exist in most of the coastal regions that are threatened by tsunami hazard. Satellite photogrammetry provides an effective tool to close the data gap.

Conclusions
In this study, we map coastal bathymetry using high-resolution satellite stereo imagery. The green band (490-610 nm) of Pleiades data outperforms the other bands in generating seafloor topography, owing to the maximum reflectance of seawater in this wavelength range. In order to calculate the refraction correction factor, we measure for the first time the refractive indices of seawater mixed with various phytoplankton that represent the major groups in the ocean. Changes in the refractive index of seawater as a result of heterogeneous distribution of phytoplankton are negligible in the derivation of bathymetry, confirming that the technique is generic and applicable to the global ocean. This method is limited to shallow waters ( ≤ 20 m). With increasing water depth, reflectance decreases, making image matching less reliable. The accuracy of photogrammetric bathymetry is estimated to be ∼ 1-2 m. The newly derived bathymetric data, 20 m spatial resolution, was applied to tsunami hazard assessment. Using the Penghu Islands as a case example, we demonstrate that high-resolution and high-accuracy photogrammetric bathymetry captures near-shore tsunami impacts more accurately than the commonly used GEBCO data, providing an effective tool for hazard evaluation. Such high-resolution and high-accuracy bathymetric data could potentially be used in various applications such as hazard assessment, coastal zone management, and navigation worldwide.