The spatial distributions of water depth and bottom type are important information for coastal zone management. In shallow and undulating areas such as coral reefs, exhaustive in-situ surveys of these variables are time consuming, costly, and sometimes even hazardous. As a low-cost complementary technique for mapping water depth and bottom type in these areas, various passive remote sensing methods applicable to multi-spectral satellite imagery in the visible region have been proposed.18.104.22.168.22.214.171.124.10.11.12.13.–14
Among these methods, that documented by Lyzenga et al.1 is one of the most widely used methods for water depth mapping and that proposed by Lyzenga2 is one of the most widely used methods for bottom type mapping. Both of these methods, as well as many others,34.5.–6 are based on the following shallow-water reflectance model (for each visible band) proposed by Lyzenga3 [Eq. (1) in his paper with the modification given in Eq. (6) in his paper] and others:
Here, is the remote-sensing reflectance just above the water surface, where and are the upwelling radiance (including the surface-reflected radiance) and downwelling irradiance just above the surface, respectively. is defined as , depends on the bottom reflectance and surface transmittance, is the effective attenuation coefficient, and is the water depth.
The beauty of this model is that it can be transformed to the linear form
Photons from the sun go through various underwater pathways (Fig. 1) before contributing to (the upwelling radiance just above the surface). Model (1) was derived from radiative transfer theory by ignoring light internally reflected at the surface and by simplifying multiple scattering in water.3 It is recognized that ignoring internal reflection may make the model inaccurate for very shallow water and high bottom reflectances.3 However, the inaccuracy has not been quantified: despite the wide application of the model, no study has quantitatively validated the model itself on the basis of radiative transfer theory. Although many years have passed since the model was first proposed, an understanding of its accuracy and limitations are necessary for appropriate application of the still-popular remote sensing methods based on it.
In this study, we examined the accuracy of the model and its application under various idealized optical conditions using Monte Carlo radiative transfer simulations. Specifically, we evaluated how accurately the model describes the relationship between and calculated by the simulation for each condition. Then, we investigated the error caused by using the model in remote sensing of water depth. Finally, we provided a formula that enables readers to reproduce the simulation results obtained in this study for arbitrary conditions.
Note that we do not intend to develop a new shallow-water reflectance model considering internal reflection. Such a model already exists.15 However, it is complex and cannot be linearized like Eq. (2), which is the basis of the popular linear methods: water depth mapping using a linear predictor1 and bottom type mapping using a linear index.2 In addition, its application requires numerical optimization and is vulnerable to statistical overfitting.16 Because complexity has these disadvantages, it is important to know the inaccuracy of the simple model (1) well before making it more complex; that is our strategy.
Radiative Transfer Simulation
The derivation of model (1) by Lyzenga3 presumes a level surface, homogeneous optical properties of the water, and a level Lambertian bottom with homogeneous reflectance. It also assumes that the effects of polarization, fluorescence, and inelastic scattering are negligible. We simulated radiative transfer in order to estimate the “true” value of in light of the radiative transfer theory under these assumptions: our simulations are based on the same assumptions as model (1) except for the ignorance of internal reflection and the simplification of multiple scattering. This enables us to examine only the inaccuracy caused by the ignorance or simplification.
Under these assumptions, depends only on the optical conditions listed in Table 1. Specifically, these conditions are the incident angle of the incident beam at the surface; the zenith angle of ; the single scattering albedo, defined as (: the scattering coefficient; : the beam attenuation coefficient, defined as , where is the absorption coefficient); the scattering phase function ; the optical depth of water; the irradiance reflectance of the bottom; and the relative refractive index of the air–water interface .
Optical conditions and values set in radiative transfer simulation.
|Optical condition||Values set in simulation|
|Incident angle of incident beam at the surface (θi)||0, 30, 60 (deg) and uniforma|
|Zenith angle of upwelling radiance L (θL)||0 (deg)|
|Single scattering albedo (b/c)||0.2, 0.5, and 0.8|
|Scattering phase function (β˜)||• Mob: The function described in Refs. 17 and 18 obtained by averaging the measurements of Ref. 19. The probability of backscattering B is 0.0178.|
|• M01 and P02: The functions presented in Ref. 20 and named in Ref. 21, characterized by small (0.00453) and large (0.0445) B, respectively.|
|Optical depth of water (ch)||100.2n(n=−10,−9,⋯,3),∞b|
|Irradiance reflectance of the bottom (rb)||0.1 to 0.6 (0.05 intervals)|
|Relative refractive index of the air–water interface (naw)||1.333|
The “uniform” case approximates skylight: the incident radiance is uniformly distributed over the sky.
Simulated R or R′ for c·h=∞ was used as R∞ or R∞′.
We calculated on the basis of the radiative transfer simulation for all the combinations of optical conditions shown in Table 1 ( combinations). values of 0, 30, and 60 deg represent direct sunlight with different solar elevations, and (the incident radiance is uniformly distributed over the sky) approximates diffuse sky light. was fixed at zero to simulate observation by quasi-zenith satellites. Three values of (0.2, 0.5, and 0.8) and three types of (“Mob,” “M01,” and “P02”) were used to represent various in-situ measurements.1920.–21 Here, “Mob” indicates the scattering phase function described in Refs. 17 and 18 obtained by averaging the measurements of Ref. 19. The probability of backscattering (the ratio of the backward scattering coefficient to ) was 0.0178 in our implementation. “M01” and “P02” are the scattering phase functions presented in Ref. 20 and named in Ref. 21, characterized by small (0.00453) and large (0.0445) , respectively. We set in geometric progression from to with a ratio of . This is because the scale of of interest in shallow-water remote sensing using multispectral satellite imagery is diverse: sometimes very shallow ranges such as 0 to 0.3 m13 are discussed in a centimeter scale; in other cases, wide ranges of about 1 to 20 m1,22 are targeted. The range handled in this paper corresponds to 0.05 to 19.91 m when (a possible value for clear seawater at blue and green wavelengths).18,19,23 We set to cover a wide range of 0.1 to 0.6 with a small interval of 0.05. An value of 0.6 is possible for carbonate sand at green and red wavelengths.24 Although an of , such as 0.05, is common for algae and corals,24 this condition was not considered because it is not suitable for the use of model (1) in the form of Eq. (2): becomes nonpositive when , and we cannot calculate the logarithm in Eq. (2). Because is rather stable in natural waters,18 it was fixed at a typical value (1.333).
Our simulation was based on the forward Monte Carlo method described in Ref. 25. For each combination of , , , and , we injected photons downward to the water surface and traced their coordinates until they disappeared by absorption or flew into the air. In this process, we used MT1993726 pseudo-random numbers with the optical conditions to determine the behaviors of each photon when it interacts with the water surface, water body, and bottom: whether the photon is reflected or refracted when it hits the water surface, the distance it travels before interacting with water, the type of interaction (scattering or absorption), and the direction it travels after being scattered in water or reflected at the bottom. Then, we tallied the photons just above the surface to estimate . Here, because we fixed at zero, only the upwelling photons with zenith angles smaller than 5 deg were counted to calculate .
The simulation code was developed in the Fortran 90 language and validated using the approaches described in Refs. 27 and 28. The former approach is based on a rigorous relationship among the vector and scalar irradiances and (the first formula of Ref. 27). The latter approach is based on an analytical solution (shown in Table 1 of Ref. 28) for a collimated beam normally incident on a homogeneous single-layer slab with isotropic scattering. For further validation, we also developed a backward Monte Carlo code and confirmed that its output was consistent with that of the forward Monte Carlo code.
Gordon and Brown25 showed that any radiometric quantity for an arbitrary bottom reflectance () can be exactly computed from the contribution to from the photons that strike the bottom zero, one, and two times when . Therefore, we performed the simulation with for each combination of , , , and , and the values for various listed in Table 1 were then calculated using this time-saving technique. In addition to , the remote-sensing reflectance without the contribution of photons internally reflected at least once () was also tallied.
The upward reflection at the water surface was not simulated directly except for the loss of incident photons in the case of . Instead, the effect was corrected afterward based on Fresnel equation.
First, we evaluated how accurately model (1) can describe the relationship between and calculated by the radiative transfer simulation for each combination of , , , and . Model (1) can be transformed to the following form via Eq. (2):
For each combination, we prepared a dataset consisting of the values listed in Table 1 (except for ; 14 in total) and the corresponding values of the simulated . Here, the simulated for was used as . Next, Eq. (3) was fitted to the dataset by least squares fitting, treating and as free parameters. Then, the determination coefficient () and root mean square residual (RMSR) of the fitting were evaluated. We also performed the same evaluation for (by replacing and with and in the above procedure) in order to examine the effect of internal reflection on the model’s accuracy.
In most of the applications of model (1) where is assumed to be a constant, the relationship between (or ) and is essentially equivalent to that between (or ) and . Therefore, the above procedure is equivalent to an accuracy evaluation of model (1) in terms of the relationship between (or ) and for arbitrary .
Second, we investigated the error caused by using model (1) in remote sensing of water depth when the depth range of the calibration data is limited. We considered the simplest case in which the bottom type is uniform and thus only one visible band is used. First, Eq. (3) was calibrated with (fitted by least squares to) the data for limited ranges selected from the dataset described above. Then, the calibrated Eq. (3) was used to estimate from the simulated values for the entire range, and the estimation error was evaluated.
Results and Discussion
Accuracy of Model (1)
Figure 2 demonstrates the relationships between or and calculated by the radiative transfer simulation for several combinations of , , , and . Obviously, or increases with decreasing and increasing . Naturally, this is due to the increase in the bottom reflection component of or .
Figure 2 also shows the curves of model (1) as expressed in Eq. (3) fitted to the plotted data as described above. Overall, the fitting is fairly good for both and . The model appears to underestimate small values when is large, but the error (residual of the fitting) is actually small: the error in the region is exaggerated because the axis is in log scale.
Table 2 lists the statistics of the and RMSR values of the fitting for all the combinations of , , , and . According to this table, the overall mean is 0.9991 for and 0.9997 for . The overall mean RMSR is 0.02861 for and 0.01274 for . A RMSR of 0.02861 is just 0.72% of the target range (). The overall minimum for is 0.9935, whereas the overall maximum RMSR for is 0.09158. Even this maximum RMSR is just 2.3% of the target range. These results indicate that when model (1) is fitted to data that cover the entire target range of , it is reasonably accurate in describing the relationship between and although the internal reflection slightly degrades the accuracy on average. Figure 3 shows (the contribution of internally reflected photons to ) as a function of for the same combinations of optical conditions as in Fig. 2. We can observe that increases with decreasing and increasing . In fact, increased monotonically with decreasing for all the combinations of , , , and ( combinations) except for two combinations with , in which for was slightly larger than that for . increased monotonically with increasing for all the combinations of , , , and ( combinations) without exception.
Statistics of R2 and RMSR for fitting model (1) [in the form of Eq. (3)] to the dataset of c·h and simulated log[R−R∞] for all the combinations of θi, b/c, β˜, and rb. The values for R′ (results obtained when log[R′−R∞′] is used instead of log[R−R∞]) are also shown.
Figure 4 is a contour plot of the averaged ratio of to versus and . We can see from this figure that internal reflection significantly increases even in common conditions. For example, the averaged ratio of to exceeds 0.05 (5%) when and , as indicated in the figure by the symbol “.” Here, a value of 0.63 corresponds to an value of 2.1 m when , and both of these values are not too small to be common: has been observed in various locations in oceans and the Mediterranean for green wavelengths19,23 and in a coral reef even at 660 nm.29 An value of 0.4 is not a large value for carbonate sand in coral reefs: the average of 670 measurements by Hochberg24 falls within 0.4 to 0.6 at wavelengths of 475 to 700 nm. The maximum value of the ratio of to for all the combinations of optical conditions listed in Table 1 was 0.28 (28%).
Figure 3 also shows the curve of the following exponential function fitted to the plotted data:
Error Caused by Model (1) When Calibration Data are Limited
Figure 5 shows the relationship between and calculated by the radiative transfer simulation for the same combinations of optical conditions as in Fig. 2. However, the curves in Fig. 5 are not those of model (1) calibrated using all the data plotted. The red curve shows the model as expressed in Eq. (3) fitted to only the five data points with the smallest (), and the blue curve shows model (1) fitted to only the five data points with the largest (). This figure shows what happens when the depth range of the calibration data does not fully cover the target depth range for depth estimation using model (1). For medium and large (, 0.5), the red curve underestimates in the region with large , whereas the blue curve underestimates in the region with small . This shows that when model (1) is used to estimate from , it significantly underestimates if it is calibrated using data for a different depth range. In short, model (1) cannot accurately extrapolate using .
Figure 5 also shows the maximum error (“m.e.” in the figure) of estimating the plotted from using each curve. In other words, the maximum error is the maximum difference in between the curve and the plotted data points. Naturally, for the red curve, which was calibrated using the data points with the smallest , the maximum error is the error for the data point with the largest . Similarly, for the blue curve, which was calibrated using the data points with the largest , the maximum error is the error for the data point with the smallest . For example, when , , =“Mob,” and , the maximum error for the red curve is 1.382. This means that the red curve underestimated the maximum by 1.382. When , this error is as large as , which is 35% of the target depth range (13.24 m).
Table 3 lists the top three maximum errors for all the combinations of , , , and . The largest error, 1.705, was observed when , , =“M01,” and . This error corresponds to of depth estimation error when and accounts for as much as 43% of the target depth range (13.24 m).
The top three maximum errors in c·h estimation using model (1) [in the form of Eq. (3)] fitted to only five data points with smallest or largest c·h. The corresponding optical conditions and the maximum error for R′ (the maximum error obtained when log[R′−R∞′] is used instead of log[R−R∞]) are also shown.
|Data used for fitting||Rank||Maximum error for R||θi||b/c||β˜||rb||Maximum error for R′|
|Five data points with smallest c·h||1||1.705||0||0.5||M01||0.6||0.062|
|Five data points with largest c·h||1||0.345||0||0.8||M01||0.6||0.023|
Table 3 also shows the maximum error for (the maximum error obtained when is used instead of ). Considering that the maximum errors for are of those for , the large errors for can mostly be attributed to the fact that model (1) ignores internal reflection. Although both and decrease exponentially with as discussed above, their sum is only approximate and not strictly an exponential function of , and hence the exponential model (1) cannot be used for extrapolation.
Mathematical Presentation of Simulation Results
To enable readers to reproduce the simulation results obtained in this study for arbitrary conditions, we empirically modeled and as follows:
Here, is the contribution of specular reflection of incident light at the surface, which is nonzero only for and . Based on Fresnel equation, is for and 0.02037 for . is the contribution of in-water scattering for infinitely deep water, which was calculated on the basis of Ref. 30 as7) is the Fresnel reflectance of upwelling radiance at the water surface.
The decimal coefficients in models (5) and (6) were determined by least squares fitting of the models to the simulated and for all the combinations of optical conditions shown in Table 1. The values of the fitting of models (5) and (6) were as large as 0.999854 and 0.999152, respectively. The RMSR values were as small as 0.001132 and 0.000895, respectively. These results indicate that the models can accurately reproduce the simulated and .
We examined the accuracy of the widely used shallow-water reflectance model (1) using Monte Carlo simulations and found that the internal reflection at the water surface significantly increases (the remote sensing reflectance just above the surface) at small depths and large bottom reflectances. However, this does not make the shallow-water reflectance model (1) significantly inaccurate. The exponential model, if fitted to (calibrated with) data covering the entire target depth range, describes the relationship between and depth reasonably accurately. This is because the internally reflected component of , as well as the other component, decreases exponentially with depth.
However, because the sum of two exponentially decreasing functions is not strictly exponential, the model cannot be applied accurately when the calibration data do not cover the entire depth range of interest. For example, when used to estimate depth from , the model significantly underestimates the depth (on the order of meters in some cases) if it is calibrated using data for a different depth range.
Ariyo Kanno received his PhD of environmental studies from The University of Tokyo in 2010 after working as a research fellow of the Japan Society for the Promotion of Science. Since then, he is working as an assistant professor at Graduate School of Science and Engineering, Yamaguchi University. He is interested in the statistical analysis of environmental data including satellite images and meteorological data. He has expertise in multispectral remote sensing, radiative-transfer simulation, numerical flow simulation, and field observation in coastal waters. He was awarded two incentive prizes for his papers in multispectral bathymetry: one in 2009 and another in 2011.
Yoji Tanaka received his PhD of environmental studies from The University of Tokyo in 2007. He began his professional career as a postdoctoral fellowship in 2007 at the Port and Airport Research Institute (PARI). He has been a research associate at Yokohama National University since 2012. He has expertise in computational fluid dynamics, coastal environment, and statistics. He is currently working on the water quality problems in semi-enclosed seas. He was awarded the incentive prize of Annual Journal of Hydraulic Engineering from the Japan Society of Civil Engineers in 2013 for his paper regarding the effect of increasing solar radiation on the environment of a semi-enclosed bay.
Masahiko Sekine received his PhD of engineering from Kyoto University in 1991. He is currently a professor at Graduate School of Science and Engineering, Yamaguchi University. He has worked in the field of environmental engineering using his wide range of expertise in sanitary, ecological, and river engineering. His current research topics include evaluation and design of river fish habitat and quick testing of environmental water toxicity. He has been a member of various committees in Japan Society of Civil Engineers, Japan Society of Water Environment, and International Society for Ecological Modeling. He has also contributed to a number of local and international, governmental and nongovernmental projects to improve aquatic environments.