Deriving Bathymetry from Multispectral Remote Sensing Data

Abstract

The use of passive satellite sensor data in shallow waters is complicated by the combined atmospheric, water, and bottom signals. Accurate determination of water depth is important for monitoring underwater topography and detection of moved sediments and in support of navigation. A Worldview 2 (WV2) image was used to develop high-resolution bathymetric maps (four meters) that were validated using bathymetry from an active sensor Light Detection and Ranging (LiDAR). The influence of atmospheric corrections in depth retrievals was evaluated using the Dark Substract, Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) and the Cloud Shadow Approach (CSA) atmospheric corrections. The CSA combined with a simple band ratio (Band2/Band3) provided the best performance, where it explained 82% of model values. The WV2 depth model was validated at another site within the image, where it successfully retrieved depth values with a coefficient of determination (r²) of 0.90 for all the depth values sampled, and an r² of 0.70, for a depth range to 20 m. The WV2 bands in the visible region were useful for testing different band combinations to derive bathymetry that, when combined with a robust atmospheric correction, provided depth retrievals even in areas with variable bottom composition and near the limits of detection.

1. Introduction

Accurate determination of water depth is important for monitoring underwater topography and detection of moved sediments, and for producing nautical charts in support of navigation [1]. Water depth information is fundamental in discriminating and characterizing coral reef habitats, and allows estimation of bottom albedo, which can improve benthic habitat mapping [2]. The use of satellite and airborne sensors in shallow waters is complicated by the combined atmospheric, water, and bottom signals. This includes variability in the signal from water column scattering and absorption due to dissolved and suspended materials (i.e., sediments, chlorophyll, and colored dissolved organic matter). Some of these limitations have been overcome by the introduction of high-resolution multispectral sensors that use light reflected from the seafloor to extract benthic information [3,4].

Previous researchers have mapped water depth but, in many cases, the methods that were used required the input of known depth values [5], or assumptions that a pair of spectral bands can be identified such that the ratio of the reflectance in these two bands was the same for all bottom types [6]. The limitations and validity of obtaining water depth from passive sensors was demonstrated by Maritorena et al. [7]. Lyzenga [8] further enhanced the methodology using multiple linear regressions, which required knowing the optical properties of the water at the time of image acquisition. Stumpf [9] expanded on the strategy developed by Lyzenga [5] of using the blue/green spectral bands for depth estimates by reducing the number of parameters that needed to be estimated. The ratio transform proposed by Stumpf [9] assumed that depth-driven change is significantly larger than the corresponding benthic albedo-driven change. These authors used the ratio transform with two Ikonos satellite sensor wavebands, characterized by differential water attenuation.

Mishra [10,11] estimated water depth for each pixel based on a site-specific polynomial model using high-resolution multispectral IKONOS data in a site near Roatan Island, Honduras. A ratio of wavebands (blue and green) were identified that were constant for all bottom types, and it was found that the correlation coefficient between actual depth and estimated depth was 0.94 m, with a Root Mean Square Error (RMSE) of 2.71 m. Based on this approach, the model overestimated depths beyond 21 m.

Collin and Hench [12] presented another approach in which water depth was retrieved from the WV2 imagery based on different band combinations from all bands (including 5 visible bands) provided by the sensor. This built on the Stumpf [10] method by enhancing the digital depth models and increasing the range of depth estimation by testing different atmospheric corrections and spatial resolutions.

Hyperspectral sensors have been used to derive properties of the water column and bottom, which includes bottom albedo and water depths. The Airborne Visible Infrared Image Spectrometer (AVIRIS) hyperspectral sensor was used to obtain water depth as well as other properties based on model-derived image optimization techniques [13]. The results suggested that the model and methods work well for extracting subsurface water properties and demonstrated that the model-derived depths agree with depths measured from 0 to 4.6 m. The AVIRIS sensor was also used in deriving bottom depths in a shallow water marine environment utilizing a remote sensing reflectance model and comparing the model-derived depths to high-resolution LiDAR bathymetry data [4]. The method is based on the assumption that the use of two bands allows separation of variations in depth from variations in bottom albedo and compensates implicitly for variable bottom [11].

Bathymetric maps for our study area, La Parguera Reserve in southwestern Puerto Rico, have relied on active sensors such as airborne LiDAR and ship acoustic surveys. The use of active sensors were mainly focused in coral reef mapping using side scan sonar [14] comparison of active sensors (LiDAR and multibeam sonar) for coral reef mapping [15], or to develop morphometrics from LiDAR to predict the diversity and abundance of fish and corals [16]. Torres-Madronero [17] performed a fusion of Airborne Imaging Spectrometer for Applications AISA hyperspectral data and LiDAR SHOALS using Lee’s bio-optical model [13] to derive depths and other water and bottom optical properties for a coral reef within our study area. This was the only study that used passive sensors to derive bottom depth but only focused on a small offshore cay within the La Parguera Reserve.

The present work has similarities as well as differences in the methods utilized from previous remote sensing-derived bathymetry approaches. This study presents the first high-resolution, large-scale water depth map for La Parguera Reserve derived from passive sensors. This bathymetric map was developed to a significant high spatial resolution (four meters) and was validated using bathymetric data from an active sensor (LiDAR). The project includes an analysis of methods and sensors to derive optimum bathymetry values from passive sensors and the influence of atmospheric corrections in depth retrievals that can be implemented throughout the Caribbean Region and elsewhere.

2. Methods

2.1. Study Area and Materials

La Parguera Natural Reserve is a marine protected area located in the southwest coast of Puerto Rico and it extends from the coast to the shelf edge with an area of approximate 168 square kilometers (Figure 1). The area of La Parguera is recognized for the exceptional value of its marine resources, which includes extensive coral reef ecosystems, seagrass beds, coastal mangrove fringe and mangrove islands, and two bioluminescent bays [18]. The shelf is also characterized by an irregular and complex physiography, with submerged patch reefs extending from 1 to 9 m off the bottom, an emergent reef with a characteristic reef crest, and a well-developed coral reef formation at the shelf edge [19]. The lack of significant rainfall and large rivers in the southwest coast, combined with the oligotrophic waters of the northern Caribbean contribute to the high levels of water transparency of the insular shelf and results in an optimization of growth rates and depth extent of coral reefs and seagrass beds [19]. The average annual rainfall at the Isla Magueyes Station in La Parguera is 74.52 cm. The average water depth is 18–20 m in the outer and middle shelf, while the shallow inner reef lagoon is less than 6 m deep.

Figure 1. WV2 image of study area of La Parguera Reserve, Southwestern Puerto Rico.

Figure 1. WV2 image of study area of La Parguera Reserve, Southwestern Puerto Rico.

2.2. Image Resources and Pre-Processing

2.2.1. WV2 Image

WV2 is a high-resolution 8-band multispectral commercial satellite operating at an altitude of 770 km with a spatial resolution at nadir of 46 cm in panchromatic mode and 1.84 m in multispectral mode at nadir. This sensor can collect a 16.4-km swath with 11-bit radiometric resolution. The advantage of this sensor for coastal areas is that it includes five bands in the visible region with a purple “coastal” band that provides better water penetration capabilities. The image was collected on 4 December, 2011 with a viewing angle of 3.27° off-nadir.

The image used in this study was acquired over La Parguera on 4 December, 2011 without apparent clouds or “sun glint” present (Figure 1).

The WV2 image was radiometrically corrected before any additional processing was performed. This radiometric correction used the WorldView Radiance calibration routine in ENVI 5.0, which converts relative radiance into absolute radiance in units of (μW·sr^-1·cm^-2·nm^-1) based on the calibration factor for each band. The image was evaluated for data gaps and a subset of the image was made based on the La Parguera Reserve polygon. The original spatial resolution of the WV2 image was 1.84 m, but the image was resampled to 4 m to match the spatial resolution of the LiDAR image, because there was no higher resolution bathymetry available. The WV2 image was warped for the co-registration with the LiDAR image and both were set to the coordinate system World Geodetic Survey 84 (WGS-84). A total of 40 points were used as ground control points for the registration and the total RMSE for the co-registration was 0.5 m.

2.2.2. LiDAR SHOALS Bathymetry

LiDAR is an active remote sensing instrument that transmits high-power pulses of laser light and in the case of LiDAR for bathymetry, the instrument uses a green laser that penetrates the water column to measure the distance from the surface of the water to the bottom [20].

LiDAR topographic and bathymetric data were acquired over La Parguera Reserve region between April and May 2006 (Figure 2) for coastal elevations of 50 m above sea level to 70 m water depths using a Laser Airborne Depth Sounder (LADS) Mk II Airborne System [21]. This airborne system uses a 900 Hz Nd: YAG (neodymium-doped yttrium aluminum garnet) laser, which is split by an optical coupler into an infrared (1064 nm) beam and a green (532 nm) beam. The final product (16-bit Geotiff image) produced a 4 × 4 m bathymetry surface where the soundings were positioned relative to the NAD83 UTM 19 N horizontal coordinate system and to the Mean Lower Low Water (MLLW) vertical tidal coordinate system [22]. According to Costa [15], the horizontal accuracy of the dataset is better than ±5 m, and the vertical accuracy or maximum total vertical uncertainty (TVU) of the dataset is better than ±0.82 m.

Figure 2. LiDAR bathymetry image, subset of La Parguera Reserve.

Figure 2. LiDAR bathymetry image, subset of La Parguera Reserve.

2.2.3. Ancillary Data

The images were corrected for fluctuations in tide readings measured at the Magueyes Island Tide Station (Station ID 9759110) in La Parguera, Puerto Rico [22]. The WV2 image was acquired on 4 December 2011 (15:25 GMT) with the tide reading for this station at 0.249 m above mean lower low water (MLLW) (15:24 GMT) . The LiDAR SHOALS bathymetry data was processed and corrected for errors in position and tidal changes [21].

3. Methods

3.1. Atmospheric Corrections

An important step in multispectral and hyperspectral remote sensing of ocean targets is to correct for atmospheric effects. To analytically derive water and/or bottom properties from any satellite ocean color data, the first step is to get high-quality spectral remote sensing reflectance (Rrs) that contains water and/or bottom information [23]. There are several methods for atmospheric correction currently available for both multispectral and hyperspectral imagery. We applied three different atmospheric corrections to our imagery including Dark Pixel Subtraction, Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH), and Cloud Shadow Approach (CSA) to evaluate the influence of the atmospheric correction in bathymetry retrievals.

3.1.1. Dark Pixel Subtraction

The Dark Pixel Subtraction or Dark Subtraction Routine [24] is based on the assumption that dark objects reflect no light and any values greater than zero must result from atmospheric scattering. The scattering is removed by subtracting these values from every pixel in the band. An area of deep oceanic pixels was selected and a dark subtraction was performed on the WV2 image.

3.1.2. FLAASH

FLAASH is a first-principle atmospheric correction tool that it is based on the MODTRAN4 radiation transfer code [25]. It is a method used to correct sensor radiance for atmospheric effects by mathematically modeling the physical behavior of radiation as it passes through the atmosphere. It can be used in multispectral and hyperspectral imagery, and basically consist of defining the tropical atmosphere over a maritime area and solve the radiative transfer equation. One of the limitations of FLAASH is that requires that the image contains bands in appropriate wavelength positions for water vapor and aerosol retrieval corrections. For the WV2 image, these bands were not available for this correction to be applied with the longest wavelength available centered at 950 nm. The FLAASH atmospheric correction was evaluated for the WV2 image.

3.1.3. Cloud-Shadow Approach (CSA)

CSA is a practical image-driven method for correcting the effects of atmosphere and obtains remote sensing reflectance from the image. This technique was used by Lee [23] in which the atmospheric radiance (La) is calculated from a pair of sun and shadow pixels, where the product of (t) transmittance and down-welling irradiance (Ed) is estimated using the reflected radiance from the top of clouds [26]. One advantage of this image-driven approach is that all radiance measurements come from the same sensor and are collected simultaneously.

For this application of the CSA we used a modified approach from Lee [23]. Since the image for the selected study area had no clouds, we used the shadow of a tall building in our study area to obtain the sun-shadow pixels. To calculate the reflectance we collected in situ spectra of a homogenous area on the rooftop of the building using a field spectroradiometer (GER 1500, SpectraVista Corp., Poughkeepsie, NY, USA) within the sun-shadow pixels selected. For the down-welling irradiance (Ed) sky (diffuse)/Ed total ratio we used the EKO Tracker STR-22 Multiple Radiometers instrument at the Bio-Optical Oceanography Laboratory in Isla Magueyes, which collects irradiance from direct and diffuse solar radiation and had data for the exact image day/time acquisition. The final product was an atmospheric corrected WV2 image with values of remote sensing reflectance (Rrs, sr⁻¹).

3.2. Bathymetry Retrieval

In the water column, the spectral signal undergoes selective exponential attenuation in the visible region as mediated by absorption and scattering processes [27]. There have been different approaches to map water depth based on the use of wavebands or combinations of wavebands. Our approach was focused on deriving depth using band combinations (band ratios) to maximize the use of visible bands available from WV2. This approach used a simple blue/green band ratio for the WV2 image. The technique was also used to evaluate the influence of the atmospheric correction in bathymetry retrievals. The selected band combinations were calibrated using LiDAR SHOALS to derive bathymetry.

3.2.1. Depth Estimated from Band Ratios

The bands selected were visible bands in the blue/green region that provide higher water penetration for improved bathymetry retrieval. This approach assumes that a pair of wave bands can be identified and that the ratio of the reflectance in these bands is the same for all bottom types within a given scene [6]. All five visible WV2 bands were evaluated but only combinations of bands 1, 2, and 3 were selected for the estimation of bathymetry.

Once we established the best band pairs from the imagery, a simple band ratio was used to estimate bathymetry. According to Mishra [11], wavebands (such as blue and green) have different water absorptions, one band will have arithmetically lesser values than the other while changes in bottom albedo affects both bands similarly. In addition, as depth increases, radiance of the band with higher absorption (green) decreases proportionally faster than the band with lower absorption (blue); and the radiance ratio of the blue to the green increased [28].

3.2.2. Ground-Truth from LiDAR SHOALS

To establish a relationship between bottom depth from WV2 band ratios and depth from LiDAR SHOALS bathymetry, a random sample was used to extract the values from both images. This random sample consisted of 5000 points and accounted for all bottom types and depths across the images. These random points were constrained to the La Parguera Reserve Boundary and used as ground-truth data for the calibration (Figure 3).

Figure 3. LIDAR SHOALS random points selected for La Parguera Reserve.

Figure 3. LIDAR SHOALS random points selected for La Parguera Reserve.

4. Results

4.1. Atmospheric Corrections

To evaluate the atmospheric influence in the remote sensed signal, radiance values of water over different bottom types were analyzed before and after the atmospheric correction. These bottom types included seagrass, sand, corals, mud and deep-water areas (Figure 4). The results after the atmospheric correction show that sand had the highest values when compared to other substrates (e.g., coral, seagrass, mud). The atmospheric contribution to the radiance signal for all substrates was estimated at 83% and the “coastal blue”, blue and green bands were found to have the highest values after atmospheric correction, while the near infrared band values were all reduced to near zero.

Figure 4. Comparison of top of atmosphere (TOA) and water leaving radiances (Luw) derived from WV2 image over different substrates that include seagrass, sand, coral reef, mud and deep water after the CSA atmospheric correction. Notice that the scale in the y-axis of the sand graph is higher.

Figure 4. Comparison of top of atmosphere (TOA) and water leaving radiances (Luw) derived from WV2 image over different substrates that include seagrass, sand, coral reef, mud and deep water after the CSA atmospheric correction. Notice that the scale in the y-axis of the sand graph is higher.

4.2. LiDAR Depths

LiDAR bathymetry data for the study area ranged from 0 to 46 m. The inner shelf area can be distinguished by its complex physiography of submerged patch and emergent reefs, inshore cays, and mud and seagrass plains and depths ranging from 0 to 20 m. The mid and outer shelf have a wider range of depths ranging from 10 to 47 m, with the deeper areas located in the southwest area of the reserve (>30 m) near the shelf edge. This area is characterized by extensive coral reef development, especially near the shelf edge, where depth of this habitat can range from approximately 12 to 35 m. The average depth is approximately 16 m for the selected areas (Figure 5).

Figure 5. LiDAR SHOALS depth frecuency histogram for the selected points in La Parguera Reserve.

Figure 5. LiDAR SHOALS depth frecuency histogram for the selected points in La Parguera Reserve.

4.3. Simple Band Ratios

Different band ratios were developed to determine the best ratio to retrieve bathymetry from the WV2 image. The band selection was based in selecting bands in the blue/green region that provided the best water penetration for successful bathymetry retrieval. Based on this approach 3, 2 band ratios (B1/B2, B1/B3, B2/B3) from WV2 visible bands were selected for six depth ranges. The correlation coefficient was calculated to test the performance of band ratios for different atmospheric corrections and was also divided in six depth ranges to evaluate the influence of depth. (Table 1)

Table 1. The coefficient of determination (r²) for six depth ranges for the different atmospheric corrections. The performance of each band ratio on (B1/B2, B1/B3, B2/B3) was evaluated when compared with LiDAR bathymetry per pixel values.

		No Atmospheric Correction	Dark Subtract	FLAASH	Cloud Shadow
Depth	Band Ratio	r2	r2	r2	r2
25–35 m	B1/B2	0.061	0.146	0.000	0.073
	B1/B3	0.154	0.229	0.174	0.175
	B2/B3	0.037	0.121	0.063	0.054
20–30 m	B1/B2	0.015	0.096	0.008	0.020
	B1/B3	0.277	0.315	0.103	0.316
	B2/B3	0.239	0.318	0.145	0.265
15–25 m	B1/B2	0.031	0.099	0.002	0.033
	B1/B3	0.169	0.248	0.047	0.222
	B2/B3	0.212	0.272	0.076	0.276
10–20 m	B1/B2	0.042	0.052	0.008	0.044
	B1/B3	0.086	0.134	0.024	0.121
	B2/B3	0.100	0.141	0.026	0.156
5–15 m	B1/B2	0.132	0.188	0.044	0.161
	B1/B3	0.254	0.465	0.067	0.408
	B2/B3	0.289	0.425	0.063	0.524
0–10 m	B1/B2	0.098	0.097	0.035	0.182
	B1/B3	0.283	0.596	0.090	0.576
	B2/B3	0.350	0.359	0.095	0.723

Table 1. The coefficient of determination (r²) for six depth ranges for the different atmospheric corrections. The performance of each band ratio on (B1/B2, B1/B3, B2/B3) was evaluated when compared with LiDAR bathymetry per pixel values.

		No Atmospheric Correction	Dark Subtract	FLAASH	Cloud Shadow
Depth	Band Ratio	r2	r2	r2	r2
25–35 m	B1/B2	0.061	0.146	0.000	0.073
	B1/B3	0.154	0.229	0.174	0.175
	B2/B3	0.037	0.121	0.063	0.054
20–30 m	B1/B2	0.015	0.096	0.008	0.020
	B1/B3	0.277	0.315	0.103	0.316
	B2/B3	0.239	0.318	0.145	0.265
15–25 m	B1/B2	0.031	0.099	0.002	0.033
	B1/B3	0.169	0.248	0.047	0.222
	B2/B3	0.212	0.272	0.076	0.276
10–20 m	B1/B2	0.042	0.052	0.008	0.044
	B1/B3	0.086	0.134	0.024	0.121
	B2/B3	0.100	0.141	0.026	0.156
5–15 m	B1/B2	0.132	0.188	0.044	0.161
	B1/B3	0.254	0.465	0.067	0.408
	B2/B3	0.289	0.425	0.063	0.524
0–10 m	B1/B2	0.098	0.097	0.035	0.182
	B1/B3	0.283	0.596	0.090	0.576
	B2/B3	0.350	0.359	0.095	0.723

The best overall correlation was by using B2/B3 with the CSA atmospheric correction for the depth range of 0–10 m (r² = 0.72). The best performance for the image with no atmospheric correction was the band ratio B2/B3 for the 0–10 m depth (r² = 0.35). For the dark subtract atmospheric correction the band ratio B1/B3 provided the best performance and was for the depth range of 0–10 m (r² = 0.60). For the atmospheric correction FLAASH, the best performance for WV2 was the band ratio B1/B3 for the 25–35 m depth band (r² = 0.17).

The band ratios were also evaluated including all sampled points and depth ranges and the performance for the no atmospheric correction was B1/B2 = 0.04, B1/B3 = 0.48, B2/B3 = 0.60; for dark subtract correction B1/B2 = 0.69, B1/B3 = 0.86, B2/B3 = 0.85; for the FLAASH atmospheric correction were B1/B2 = 0.29, B1/B3 = 0.45, B2/B3 = 0.46, and for the CSA atmospheric correction were B1/B2 = 0.64, B1/B3 = 0.86, B2/B3 = 0.88.

Based on the Pearson-product moment relationships of band ratios for different atmospheric corrections, we selected the best models to determine a best-fit equation for the relationship. For the band ratio with no atmospheric correction the equation of best fit was a second order polynomial (r² = 0.65), for the dark subtract atmospheric correction was a lineal equation (r² = 0.74), for the FLAASH atmospheric correction the best fit was an exponential equation (r² = 0.48), and for the CSA atmospheric correction the best fit was a second order polynomial (r² = 0.82) (Figure 6).

Figure 6. Coefficient of determination and best models to determine an equation of best fit for the relationship between the band ratios of WV2 image and LiDAR SHOALS data. The correlations presented are for no atmospheric correction (Top Left), the dark subtract atmospheric correction (Top Right), the FLAASH atmospheric correction (Bottom Left), and the CSA atmospheric correction (Bottom Right).

Figure 6. Coefficient of determination and best models to determine an equation of best fit for the relationship between the band ratios of WV2 image and LiDAR SHOALS data. The correlations presented are for no atmospheric correction (Top Left), the dark subtract atmospheric correction (Top Right), the FLAASH atmospheric correction (Bottom Left), and the CSA atmospheric correction (Bottom Right).

After the selection of the best band ratio and model for the relationship, a site validation was used to evaluate the model performance in retrieving depth for other locations within the WV2 image.

This model was applied to the image and analyzed for accuracy in depth estimation when compared with LiDAR depths. We also evaluated the influence of the atmospheric correction in the bathymetry retrievals.

4.4. WV2 Model Validation

The model from the B2/B3 ratio from CSA with a second order polynomial equation (y = −0.0006x² + 0.0359x + 1.1498) was applied to another location within the WV2 image (Figure 7) This model was applied to the image and analyzed for accuracy of the model in depth estimation when compared with LiDAR depths. The validation area shows an average depth of 15.7 m with a depth range from 1 to 43 m (Figure 8).

Figure 7. Validation site within the WV2 image off Southwestern Puerto Rico.

Figure 7. Validation site within the WV2 image off Southwestern Puerto Rico.

Figure 8. LiDAR SHOALS depth frecuency histogram for the selected points in the site validation area in Southwestern Puerto Rico.

Figure 8. LiDAR SHOALS depth frecuency histogram for the selected points in the site validation area in Southwestern Puerto Rico.

The values collected from WV2 Depth Model and LiDAR bathymetry were log transformed and then plotted to validate the model performance in estimating depth (Figure 9). The coefficient of determination (r²) was 0.90 for all the depth values sampled. The values were also analyzed to determine the model performance by depth range (Table 2).

Figure 9. Linearized (Ln) Relationship between values of the WV2 Depth Model and the LiDAR SHOALS depth (meters).

Figure 9. Linearized (Ln) Relationship between values of the WV2 Depth Model and the LiDAR SHOALS depth (meters).

Table 2. A coefficient of determination (r²), standard error, standard deviation and the Root Mean Square Error (RMSE) for six depth ranges for the WV2 depth model for validation. WV2 Depth Model vs. LiDAR SHOALS.

Depth Range	r2	St. Dev. (m)	N	Standard Error	RMSE
1–10 m	0.72	3.54	1921	1.58	1.26
5–15 m	0.70	6.42	3591	2.16	1.47
10–20 m	0.69	9.54	3672	2.47	1.57
15–25 m	0.44	14.36	2760	2.80	1.67
20–30 m	0.18	17.36	2180	3.00	1.73
25–35 m	0.07	19.67	1101	3.16	1.78
All depths	0.90	10.54	8110	2.43	1.56

Table 2. A coefficient of determination (r²), standard error, standard deviation and the Root Mean Square Error (RMSE) for six depth ranges for the WV2 depth model for validation. WV2 Depth Model vs. LiDAR SHOALS.

Depth Range	r2	St. Dev. (m)	N	Standard Error	RMSE
1–10 m	0.72	3.54	1921	1.58	1.26
5–15 m	0.70	6.42	3591	2.16	1.47
10–20 m	0.69	9.54	3672	2.47	1.57
15–25 m	0.44	14.36	2760	2.80	1.67
20–30 m	0.18	17.36	2180	3.00	1.73
25–35 m	0.07	19.67	1101	3.16	1.78
All depths	0.90	10.54	8110	2.43	1.56

The WV2 depth model successfully retrieved depth values with an r² = 0.70 to a depth range of 20 m. Additionally, of the total values considered (n = 8110), 69% of these were between the depth range of 1–20 m (n = 5600). The aggregated RMSE for the validation values of the WV2 depth model was 1.56 m. These values of RMSE increase from 1.26 m for the 0–10 m depth range to 1.78 m for the 25–35 m depth range. A derived-bathymetry image was generated from the model used (Figure 10).

Figure 10. (Top) Depth estimation based on WV2 depth model for validation site within the WV2 image off southwestern Puerto Rico. (Bottom) Depth values (in meters) for WV2 depth model and LiDAR depths with 95% confidence and prediction bands. Note values dispersion after 20-m depth.

Figure 10. (Top) Depth estimation based on WV2 depth model for validation site within the WV2 image off southwestern Puerto Rico. (Bottom) Depth values (in meters) for WV2 depth model and LiDAR depths with 95% confidence and prediction bands. Note values dispersion after 20-m depth.

5. Discussion

The water depth models tested here showed the importance of the atmospheric correction used in high-resolution multispectral image analysis. The data followed the results of Collin and Hench [12] where the relation between modeled and actual depths varied as a function of the atmospheric correction used. This empirical model using regression-based modeling requires to take into account the local set of conditions of the study area and the atmospheric effects on the signal path length [29]. Since the sensible depth is a function of the wavelength used and water clarity [1], different band ratios were evaluated. The waters of La Parguera Reserve although considered oligotrophic, are susceptible to sediment resuspension, CDOM input from streams, and higher chlorophyll concentrations especially in the inner-shelf area [30].

5.1. WV2 Depth Model

For the WV2, the cloud shadow atmospheric correction provided the best results (r = 0.875) of all models tested including all the depth values and depth ranges that were analyzed. Both the cloud shadow and the dark subtract corrected values (r = 0.772) show an improvement from the values using no atmospheric correction (r = 0.592), confirming that the atmospheric correction is an essential pre-processing step in the image analysis process [8]. As expected with bathymetry models from passive sensors [31], the shallower values from the depth band of 0–10 m provide the best correlation between modeled and actual values (Table 1). This correlation presented a gradual decrease with increasing depth, and for the depth range 25–35 m, these values were significantly reduced. However, for this maximum depth range the correlation was improved from 0.397 from the no atmospheric correction, to 0.479 and 0.418 for the dark subtract correction and cloud shadow corrections, respectively.

For the dark subtract the best model used the combination of the “coastal” band (425 nm) with the green band (B1/B3 = 0.863) for all depth values. The use of the “coastal” band in combination with the green band outperformed the traditional band combination that other high spatial resolution sensors provide (i.e., QuickBird, IKONOS) where the blue band is centered at 480 nm. For the cloud shadow atmospheric correction the best model used the combination of the blue and the green bands (B2/B3 = 0.823) for all depth values.

The coefficient of determination was obtained and an equation of “best fit” and applied to our depth models. The model with the cloud shadow atmospheric correction provided the best performance (r² = 0.823). For each of the depth ranges the r² decreased with depth from the highest r² = 0.732 at 0–10 m, 0.532 at 5–15 m, to 0.162 at 10–20 m, suggesting that the model performance is reduced for values deeper than 15 m deep, as expected.

For the band ratio with no atmospheric correction the equation of best fit was a linear (r² = 0.638), and for the FLAASH atmospheric correction the best fit was an exponential equation (r² = 0.481) (Figure 5). The scatterplot for the FLAASH atmospheric correction show the values dispersed from the exponential model for all depths up to the 23 m depth band. A present limitation of the FLAASH model for atmospheric correction is that WV2 imagery does not include bands for determining water vapor and aerosols [12]. Our results showed that the analytical atmospheric correction (FLAASH) provided poorer results than no atmospheric correction and the dark subtract correction.

5.2. WV2 Depth Model Validation

The WV2 depth model successfully retrieved depth values with an r² = 0.90 for all depth ranges, and of 0.70 to a depth range of 20 m at the validation site. The derived bathymetry map is important since it provides good depth estimates to 20 m, which represents 75% of the total depth range of the study area.

Other researchers have reported various depth ranges and values using band ratios but these were limited to shallow areas with very high water clarity [4,8], which contrast significantly with our study area that consist of variable bathymetry, heterogeneous bottom composition and geographic extent. Additionally, the random point’s selection for the model was not limited to a specific depth range (i.e., 15 m) or a geographic location (i.e., near-shore). By including this near-shore areas with higher turbidity and deeper areas the model represents more realistic conditions and the limitations of the image derived bathymetry using the WV2 image. Similar to our results, other researchers have used a site-specific model and obtained good results to 20 m [10,11], and to 30 m [12].

The WV2 depth model was adjusted using LiDAR SHOALS bathymetry data an assumed that the LiDAR values were obtained with minimal error. According to Costa [15], the pixels with the largest depth differences estimated by the LiDAR bathymetry for our site (that were greater than the maximum allowable vertical error) occurred primarily along the shelf edge in water deeper than 35 m. Although approximately only 25% of our study area exceeds the 20-m depth range, this limitation in the accuracy of the bathymetry from the LiDAR image in deeper areas, combined with the significant water attenuation effects of the passive signal, limits the retrieval of more accurate depths in these areas.

6. Conclusions

This study presents the first broad scale (168 km²) bathymetry map for La Parguera Reserve derived from a passive sensor over a broad depth range and to the limits of detection of WV-2 data in these waters, which is also related to the bottom albedo signal. It also provided a comprehensive evaluation of the different atmospheric corrections and the influence in the water depth retrievals. This study also confirms that atmospheric correction methods should be applied to WV2 imagery prior to the development of bathymetric maps in order to increase the accuracy of the final product [1,12,32]. Additionally, the selection of the atmospheric correction method is extremely important in the performance of the depth model, even for simple blue/green ratios, and can be determinant based on the final product that will be developed from the imagery (i.e., benthic habitat map, depth retrievals).

In our case the best model was the one that use the Cloud Shadow Approach with a Band 2/Band 3 combination for the study area of La Parguera Reserve. This atmospheric correction proved successful even when it was modified due to cloud absence in the scene. Also, the use of simple band ratios provided the first high-resolution bathymetry map for La Parguera Reserve from passive sensors. The dark subtract correction provided the second best results and proved to be an excellent first approach for the correction of the atmospheric influence. The temporal difference in bathymetry of five years between the LiDAR bathymetry (2006) and the WV2 image (2011) were assumed minimal for the study area at this spatial scale. Additionally, the water quality conditions allowed for successful retrieval of bathymetry using passive sensor, which may not be the case in very high turbid waters. Using these methods, high-resolution bathymetry can be obtained from passive multispectral sensors with minimum processing while utilizing the information contained in the imagery.