2.2.1. Photometric Models for Planetary Surfaces
Shepard [
39] gives an overview of photometric models in planetary remote sensing. The simplest but crudest approximation is to assume an isotropic surface and apply Lambert’s reflectance model, as done by Dorrer and Zhou [
13] and Hartt and Carlotto [
12].
The empirical Lunar-Lambert model is a weighted superposition of Lambert’s law and the Lommel–Seeliger law with an additional phase function to address anisotropy, and was originally tailored to the Lunar surface [
40]. Wu et al. [
9] and Alexandrov and Beyer [
10] use the Lunar-Lambert model for the Moon and Lohse et al. [
41] and Gehrke [
15] adopt it to Mars. The models of Shkuratov et al. [
42] and Hapke (summarized in Hapke [
43]) rely on a more rigorous derivation. They are advantageous since they take physical properties into account and are proven to work well for particulate materials. The Hapke model, especially, has sound quantitative support and can be regarded as the standard model for planetary photometry. Nevertheless, some shortcomings are discussed by Shepard and Helfenstein [
44] and Shkuratov et al. [
45].
In the earth remote sensing community, the semi-empirical kernel method is commonly employed to approximate the bi-directional distribution function (BRDF) of a surface by a weighted superposition of kernels [
46]. These kernels are chosen from a catalog according to the scene type under investigation. Jiang et al. [
18] combine Ross-Thick and Li-Sparse (RTLS) kernels for Mars as described by Ceamanos et al. [
28]. The performance of this method was compared to the Hapke model, showing a good agreement in common cases, but stronger deviations, especially under large incidence angles [
28].
This approach relies on the Hapke model as the de-facto standard of photometric planetary analysis, which has successfully been employed for SfS [
1,
8]. Several formulations of the Hapke model exist. The Anisotropic Multiple Scattering Approximation (AMSA) of the Hapke model is given by:
The model is a function of a variety of photometric and geometric parameters. The main parameter is the single-scattering albedo
w. Additionally, the opposition effect is modeled with the terms
denoting the Shadow Hiding Opposition Effect (SHOE), and
, for the Coherent Backscatter Opposition Effect (CBOE). Both terms depend on two parameters each: one parameter for the strength of the opposition surge (
/
) and another parameter that describes the shape or width of the effect (
/
). In the model, multiple scattering within the medium is accounted for by the approximation
, which is given by:
The geometry is defined by the cosines of the incidence and the emission angle
and
, respectively. The angle
g is termed the phase angle and is measured between the incoming light direction and the direction of the emitted light. The material dependency is addressed by the phase function
and the single-scattering albedo
w. The phase function can be expressed in terms of Legendre polynomials with the material specific Legendre coefficients
:
Alternatively, the double-Henyey–Greenstein function can be used to approximate the phase function with only two material-specific parameters
b and
c:
The Ambartsumian–Chandrasekhar
H-function is defined as
It is not analytically solvable, but viable approximations [
43] exist:
with
The remaining terms for
and
are defined by:
The coefficient
is calculated by evaluating the following definition
2.2.2. Atmospheric Model
Mars is encircled by a thin, but non-negligible atmosphere, which alters the radiance emerging from the planet’s surface. The mean surface pressure on Mars is approximately 6.3 hPa, but is highly variable compared to Earth [
47]. The atmosphere consists almost entirely of CO
2 (≈95%) [
47], but the largest effect on the optical opacity of the atmosphere is caused by finely grained dust (effective radius approximately 1.4–1.7
m [
48]). Especially during the latter part of the Martian year, winds faster than 30 m/s [
49] can raise large amounts of dust, leading to a higher optical depth. The effects that arise from the presence of an atmosphere are addressed by embedding the reflectance model into an atmospheric model. The atmospheric model, which serves the purpose of accurate surface estimation, must be physically plausible and accurate, while staying at a manageable level of complexity. These requirements can be balanced by making a set of assumptions, which can roughly be grouped into either arising from the geometry of the atmospheric model or from the physical effects that govern the interaction between light and the atmosphere.
The Martian atmosphere can be modeled as a horizontally stratified layer on top of a flat planetary surface. The geometrical constraints for our atmospheric model are accompanied with three simplifying assumptions. (1) Because the region of interest is rather small compared to its distance to the orbiter and to the sun, we omit small angular changes and assume constant orbiter zenith angles and constant solar zenith angles. (2) Further, the topographical changes within the region of interest are small compared to the thickness of the atmosphere, which allows us to neglect them as well. (3) Unlike Earth, Mars usually does not exhibit strong weather changes on local scales and hence we can assume constant atmospheric parameters for our region of interest. If cases arise in which local weather events are present or very strong topographic changes occur, they may be addressed individually by dividing the region of interest into smaller segments. A sketch of the atmospheric slab model can be found in
Figure 1.
The total radiation that arrives at the image sensor is, in principle, a superposition of reflected and scattered radiation and thermally emitted radiation. Because the spectral ranges of CTX (500–700 nm [
37]) and HiRISE (400–1000 nm, more specifically, 550–850 nm for the red channel [
36], exclusively used in this work) imagery are in the visible region, thermal contributions in the infrared range are not relevant. The reflected and scattered light of a planet
which arrives at the sensor is commonly modeled [
50] as the superposition of three contributing processes, i.e., directly reflected light
, reflected skylight
and path-scattered light
The processes that govern the contributions are sketched in
Figure 1 (right) and are outlined in the following paragraphs.
The Lambert–Beer law [
51,
52] states that the incoming light
, which travels along a path
S through the atmosphere with extinction
, is attenuated exponentially. If the optical depth
of the atmosphere is defined as the line integral of the extinction E evaluated along a line
N normal to the surface, we need to divide by the cosine of the solar zenith angle
to project it to the path of oblique illumination
S. The top-of-atmosphere irradiance is denoted by
and the irradiance that reaches the Martian surface
is given by
The radiance
is reflected by the Martian surface and travels back to the orbiter’s image sensor along the line
V under the orbiter zenith angle
. The cosine of this angle is denoted as
. En route, the light is attenuated a second time, such that
becomes
Parts of the incoming radiance
are scattered by the atmosphere and act as a diffuse source of illumination, termed skylight. The skylight
illuminates the Martian surface diffusely and is reflected into the direction of the sensor. This process can be directly modeled by employing the hemispherical-directional reflectance
, which links diffuse illumination to a collimated detector:
The factor
describes the fraction of light that is actually scattered by the atmosphere to serve as diffuse skylight. The factor is obtained empirically in the parameter estimation step of our method. Following Nicodemus [
53] and Hapke [
43], the hemispherical-directional reflectance can be obtained by integrating the bi-directional reflectance
r over the upper half sphere:
The bi-directional reflectance
r can be implemented by Hapke’s model (Equation (
1)). It is feasible to assume that the narrow peaks of the opposition effects have a neglectable contribution to the whole integral over the upper half sphere. Therefore, we can set
and
and arrive at a simpler formulation of Hapke’s model. The integral then evaluates to
Finally, some parts of the light are scattered by the atmosphere into the direction of the sensor due to
“Rayleigh scattering and aerosol and particulate Mie scattering and [it] is reasonably assumed to be constant over the entire scene” [
50]. Thus, we state
2.2.3. Full Model and Parameter Estimation
Summarizing the components (Equations (
16), (
17) and (
20)) yields the full atmospheric model:
The full reflectance model
depends on 16 parameters summarized in
Table 1, which can be grouped into three distinct sets, i.e., geometric, atmospheric and material parameters. To accurately perform SfS, it is crucial to correctly obtain the reflectance parameters for a given scene and avoid interdependencies between parameters. The parameter determination is done in three ways. One set of parameters is given by a priori knowledge and, therefore does not need to be estimated. The second set of parameters is determined using a one-shot estimation based on the initial DEM and the image (see
Section 2.2.4 for details). These are the atmospheric parameters in particular. The remaining geometric and material parameters are iteratively updated throughout the SfS method.
The geometric parameters
and
describe the cosine of the sun zenith angle and the orbiter’s zenith angle. Given the small angle approximation, they can be assumed to be constant over the entire image and are inferred from the ancillary data of the CTX or HiRISE imagery. The geometric parameters
,
and
g are derived from the illumination vector, the observation vector and the normal vector of a pixel of the DEM. All three parameters require a DEM, but the parameters are also the input for the reflectance model, which is actually used to estimate the DEM. Due to this coupling, the geometric parameters
,
and
g, as well as the albedo
w, are estimated during the iterative shape and albedo from shading procedure (
Section 2.3).
The single-scattering albedo
w is the most important material parameter and describes the ratio of reflected and incoming light that interacts with a single particle. It is generally assumed that the albedo varies locally, and hence it is estimated pixel-wise within the iterative shape and albedo from shading procedure. The initial albedo
is assumed to be constant over the entire image and estimated once before the SfS procedure. Johnson et al. [
54] retrieved Hapke’s single-scattering albedo
w for different surface materials from pancam measurements of the Opportunity rover. In the red wavelength range for the CTX camera (500–700 nm) and the HiRISE camera (550–850 nm), the albedo is in the interval
. Laboratory measurements of Martian dust analogs (HWMK919) even yield albedos of up to
in the crucial wavelength range [
55]. Thus, physically plausible albedo values on Mars are somewhere in the broad interval of
. The phase function
describes the phase-angle dependent scattering behavior of particles and can be approximated by the double-lobed Henyey–Greenstein function. This function depends on the backscattering parameter
b and the asymmetry parameter
c, which are related by the so-called hockeystick relation [
56]. The parameters
b and
c are located on this curve depending on the material. However, as done in previous work (e.g., [
1]), there is a working assumption to equip the phase function with constant parameters. For the Moon, we used globally averaged
and
from Warell [
57]. For Mars, Fernando et al. [
58] set
and
for various materials observed by the pancam of Opportunity and Spirit. Johnson et al. [
54] found slightly smaller values:
. A globally averaged phase function based on CRISM and OMEGA observations was determined by Vincendon [
59], who found
and
. Given these data, we set
and adopt
. Fernando et al. [
58] also determined the roughness for various Martian materials, as seen by Spirit’s pancam. The unweighted average roughness of these materials is
. Vincendon [
59] derived a global roughness of
. For this study, we omit the roughness to be able to compare the model with and without atmospheric compensation, because the roughness is not considered in the hemispherical-directional reflectance.
The optical depth
and the weights
for the diffuse illumination and
for the path scattered radiance are assumed to be constant over the entire image and the estimation is performed once before the SfS algorithm is carried out. The optical depth of the Martian atmosphere undergoes seasonal variations and is usually in the range of
[
60]. The pancam of Opportunity rover reported values in the range of
at its exploration site [
60] and the 880 nm channel of the Spirit rover camera encountered values around
[
61] for roughly a quarter of the Martian year. During the dust storm season later in the year, the optical depth rises well above
and may reach levels up to
[
61], such that the atmosphere becomes opaque. The region of interest investigated in this work is at a higher latitude than the rovers, and the optical depth is usually lower compared to the equator [
33]. For a reasonable reconstruction, we set
as boundaries. The weight
is approximately the relative importance of the hemispherical reflected light compared to the directly reflected radiance. We do not expect the hemispherical reflected light to be stronger than 20% of the directly reflected light. Consequently, we set
to be between 0 and
. In contrast to
,
is an additive term. In the case of a global dust storm, almost all received radiance can be attributed to path scattered radiance—the surface is invisible. In most cases, however, it should only be a minor factor compared to the total radiance, which is in the order of
to
. Therefore, we expect
to be smaller than
.
2.2.4. One-Shot Parameter Estimation
The atmospheric parameters and the initial mean albedo are assumed to be constant over the scene and are estimated once before the SfS method and are not updated in the SfS scheme. Then, the parameters are adjusted such that the rendered image matches the measured data from the CTX or HiRISE instrument. These atmospheric parameters that best reconstruct the image data are used throughout the SfS. The resulting albedo serves as a global initialization, which is then refined locally within the SfS algorithm. Some challenges were encountered when devising a scheme for parameter retrieval and, as follows, the solutions to these challenges form the final scheme:
Firstly, pixel-locking is a stereo artifact that yields stair-like structures at slopes, e.g., at crater walls. Consequently, the gradient of these areas oscillates between a value that is either too steep or too flat compared to the real slope. The surface reflectance model depends on the surface inclination towards the sun, which is derived from the gradients of the DEM. If the reflectance takes oscillating gradients, it predicts an oscillating reflectance, which is not observed on the planet. To mitigate this effect, a Gaussian filter with a standard deviation of 2 is applied to the initial DEM in MATLAB® such that the stairs are smoothed out and the gradient becomes oscillation-free.
Secondly, the problem is ill-posed. Changes in the observed value can originate from atmospheric effects, albedo changes or changes in the surface slope. For example, choosing a lower albedo or higher optical depth can darken the overall scene. The model relies on the assumption that the scene is described by constant atmospheric parameters and that the albedo changes are of lower spatial frequency. Most of the variations in brightness are explained by the shape of the surface. To disentangle these effects, it is best to use many different data points to cover the variations in the model. This is achieved by sampling areas with significant slope variations. Nonetheless, albedo changes can account for a change in . To improve the parameter retrieval, it is best to sample areas with nearly constant albedo.
This also means that the parameters cannot be obtained entirely automatically. Experience and prior knowledge about the formation and the atmospheric conditions of the scene at time of acquisition are necessary. Physically reasonable boundaries for the parameters need to be chosen to constrain the optimization procedure. To address the given challenges, the following scheme for parameter retrieval is devised:
- (1)
The initial stereo DEM is filtered with a Gaussian filter to suppress pixel-locking;
- (2)
An image is rendered from the smoothed, downsampled DEM using the combined atmospheric and reflectance model;
- (3)
Training profiles are sampled, which should:
- (a)
not cross severe stereo artefacts such as spikes and holes,
- (b)
cover enough variations of slopes (e.g., shadows, as well as bright facets oriented towards the sun),
- (c)
have a constant albedo. This is seen if the image has brightness variations which are not correlated with slope variations;
- (4)
A Bayesian optimization algorithm in MATLAB
® (bayesopt) is used to impose constraints on the parameters and to avoid terminating in a local minimum (e.g., [
62,
63,
64]).