Analytical Relationship between Two-Band Spectral Vegetation Indices Measured at Multiple Sensors on a Parametric Representation of Soil Isoline Equations

Abstract

Differences between the wavelength band specifications of distinct sensors introduce systematic differences into the values of a spectral vegetation index (VI). Such relative errors must be minimized algorithmically after data acquisition, based on a relationship between the measurements. This study introduces a technique for deriving the analytical relationship between the VIs from two sensors. The derivation proceeds using a parametric form of the soil isoline equations, which relate the reflectances of two different wavelengths. First, the derivation steps are explained conceptually. Next, the conceptual steps are cast in a practical derivation by assuming a general form of the two-band VI. Finally, the derived expressions are demonstrated numerically using a coupled leaf and canopy radiative transfer model. The results confirm that the derived expression reduced the original differences between the VI values obtained from the two sensors, indicating the validity of the derived expressions. The derived expressions and numerical results suggested that the relationship between the VIs measured at different wavelengths varied with the soil reflectance spectrum beneath the vegetation canopy. These results indicate that caution is required when retrieving intersensor VI relationships over regions consisting of soil surfaces having distinctive spectra.

1. Introduction

Environmental studies often require data records spanning several decades of time intervals to validate a research hypothesis. Because the designed lifespan of an Earth observation satellite is generally about five years, long-term data records usually comprise multiple datasets acquired using several Earth observation satellites. Therefore, intersensor calibrations among the past, current, and future satellite sensors play an important role [1] in determining the quality of such data.

A primary purpose of intercalibration processes is to minimize biases caused by differences in the sensor specifications, such as, but not limited to, the spectral response functions [2,3,4,5], the spatial resolution [6,7,8,9], and the illumination geometries [10,11], in addition to differences in the absolute calibration of the systems [12], algorithms used for parameter retrievals, and atmospheric correction schemes [13]. Significant efforts have been devoted to understanding the mechanisms underlying such complex processes [14,15].

The intercalibration of radiances and/or reflectance spectra obtained from two similar instruments has been conducted [16,17,18], and several techniques have been proposed to minimize certain bias variations among the sensors [2,4,18]. The differences among reflectance spectra obtained from different sensors have been examined using images acquired nearly simultaneously [19,20,21]. The outputs from one sensor have been plotted against the corresponding outputs from other sensors to characterize the biases to some extent [22,23]. These studies often reach similar conclusions, for example, that the relationship (or, rigorously, the coefficients that describe the relationship) between two spectral vegetation indices (VIs) vary from site to site [17,21,24].

Recently, a series of intensive studies [25,26,27,28] and comprehensive review [29] regarding intersensor calibration methods for VIs have been reported. These studies addressed clearly the importance of resolving land cover dependencies upon intersensor calibration from a theoretical point of view. Most remarkable was Fan et al. (2014) introducing and employing a theory of radiative transfer and showing its potential as a feasible approach for intercalibration among multiple sensors [25,28,30,31]. Their work inspired us to further investigate the issue of land cover dependency over the intersensor VI-to-VI relationships. This study mainly took an analytical approach to better understand the dependency mechanism.

Several intercalibration studies in the land characterization discipline have selected VIs as an example of satellite data products [28,32,33,34,35,36,37]. Because the VI has been used as a proximity measure of biophysical parameters, the impacts of systematic errors due to differences in the band configuration may be further propagated into the downstream data records. In fact, widely used indices, such as the Normalized Difference Vegetation Index (NDVI) [38], can contribute to such biases [39,40,41,42]. Thus, the selection of a VI, as an example, provides an initial and meaningful step toward tackling bias reduction in data using complex algorithms.

Limiting our discussion to the investigation of the relationship between the VI values of different sensors, one difficulty in this investigation involves the selection (and often identification) of a model that provides a set of reflectance spectra under any desirable conditions. In any investigation based on a numerical model, this selection depends simply on the parameter range covered by a model. On the other hand, the availability of analytical and convenient models is very limited. One possibility is to use the analytical relationship between two reflectances of different wavelengths, known as the isoline equations [43,44,45,46]. The “soil isoline equatio”’ introduced in a previous study [47] is one such choice [27,29].

Soil isolines introduce several advantages into investigations of the VI relationships, as discussed elsewhere [48,49,50]. One advantage to this approach is that it can clarify the influences of the pure soil reflectance spectrum on the intersensor VI relationship. Because the soil isoline is a set of reflectance spectra obtained from canopies with a constant soil content, any relationships along a soil isoline can be purely attributed to a specific soil spectrum. Another advantage is that the use of the isoline enables the derivation of the relationship between two VI values obtained from different sensors. Two VIs may be related directly by eliminating the reflectance variables from a system of equations. One such example has been introduced elsewhere [51]. Soil isolines are good candidate measures for examining sensor-to-sensor variability. The advantages of the soil isoline equations may provide a new perspective on the intersensor calibration of VIs, to which this study contributes.

The purpose of this study is to clarify the influences of the soil reflectance spectrum on the biases that exist when two sensors capture the same target spectra at different wavelengths. To this end, we set three objectives: First, we derived a relationship between two VI values measured from different sensors based on the recently introduced soil isoline equations. Second, we validated the derived results using a radiative transfer model to describe the coupling between the vegetation canopy and leaf reflectance. The third objective was to identify the dependencies of the derived relationships on the soil reflectance spectra obtained using both analytical and numerical approaches.

The use of the soil isoline for the intercalibration of VI data has been examined briefly in pilot studies [52,53,54,55], which indicated that the soil isoline equation may provide good insights into the intersensor spectral differences. This study attempts to further advance the derivation and numerical validation studies. The theoretical background will be explained in Section 2. Then, we introduce the derivation steps conceptually without using any numerical models in Section 3. Some practical considerations are provided in Section 4. The results of numerical experiments are presented in Section 5. Finally, the discussion and conclusions summarize our findings in Section 6 and Section 7, respectively.

2. Background

A soil isoline is defined as a set of reflectance spectra, primarily in the red and near-infra red (NIR) reflectance subspaces, obtained by assuming fixed soil surface conditions beneath a vegetation canopy. Such isolines may be simulated using a radiative transfer model of the vegetation canopy, e.g., PROSAIL [56,57,58]. The conditions assumed during the simulation comprise a fixed soil spectrum that creates a set of spectra considered to be an isoline. Soil isoline equations comprise a system of equations that describe the relationship between the red and NIR reflectances, along with a soil isoline [47,59]. A formal derivation of the soil isoline has been described in previous studies [47]. Below, we simply summarize the derived isoline equation.

2.1. Soil Line Concept in the Red-NIR Reflectance Subspace

The isoline derivation requires a soil line [60], which is a well-established concept. A soil line is essentially a zero-vegetation isoline [61] and is often represented using a linear form of the relationship between the red ($R_{sr}$) and NIR reflectances ($R_{sn}$) of the soil surface,

$$\begin{array}{c}\hfill R_{sn}=s_0+s_1{R}_{sr},\end{array}$$

(1)

where $s_1$ and $s_0$ are the slope and offset, respectively. The slope of a soil line is used as information for calculating the parametric representation of soil isoline equations, which is mentioned below. Here, the angle is defined.

$$\begin{array}{ccc}\hfill \theta & =& \arctan \left(s_1\right).\hfill \end{array}$$

(2)

2.2. Soil Isoline Equation in the Red-NIR Reflectance Subspace

According to the derivation steps described in Taniguchi et al. [47], the form of the soil isoline is expressed as a parametric representation between two reflectances. Through this study, reflectance vector $\mathit{\rho }$ stores red (${\rho }_r$) and NIR (${\rho }_n$) reflectance spectra. The derivation process of parametric representation for soil isolines is summarized in Appendix A. The final form is expressed as a product of a matrix and vector:

$$\begin{array}{c}\hfill \mathit{\rho }=\mathrm{S}{\mathit{\rho }}_n^{\prime },\end{array}$$

(3)

where $\mathrm{S}$ represents a matrix whose elements are functions of soil brightness, and ${\mathit{\rho }}_n^{\prime }$ represents a vector whose elements are series of a parameter ${\rho }_n^{\prime }$, which is a transformed NIR reflectance. The above equation is equivalent to following form:

$$\begin{array}{ccc}\hfill \mathit{\rho }& =& \left(\begin{array}{ccccc}{\mathit{c}}_0,& {\mathit{c}}_1,& \cdots & {\mathit{c}}_i,& \cdots \end{array}\right)\left(\begin{array}{c}1\\ {\rho }_n^{\prime }\\ ⋮\\ {\rho }_n^{\prime i}\\ ⋮\end{array},\right)\hfill \end{array}$$

(4)

where the vector c for the ith order terms, ${\mathit{c}}_i$, in matrix $\mathrm{S}$ is constructed as two coefficients, namely ${\alpha }_i$ and ${\beta }_i$ (${\mathit{c}}_i={({\alpha }_i,{\beta }_i)}^t$), which are defined as

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\alpha }_i\\ {\beta }_i\end{array}\right)=\left(\begin{array}{cc}\cos \left(\theta \right)& -\sin \left(\theta \right)\\ \sin \left(\theta \right)& \cos \left(\theta \right)\end{array}\right)\left(\begin{array}{c}{p}_i\\ {\delta }_{1i}\end{array}\right)+\left(\begin{array}{c}0\\ s_0{\delta }_{0i}\end{array}\right),\end{array}$$

(5)

with the Kronecker delta $\delta$. The coefficients $p_i$s are obtained from a series of numerical simulations carried out using canopy radiative transfer models. The $p_i$s depend only on the soil surface reflectance. Because the red and NIR reflectances of the soil surface assumed during the derivation and numerical simulation can be described by the soil line represented in Equation (1), the only parameter that characterizes the $p_i$ values is either $R_{sr}$ or $R_{sn}$ in Equation (1) in this study. The details of the derivation are summarized in Appendix A.

The derivation of the intersensor VI relationships introduced in this study depends on the parameter ${\rho }_n^{\prime }$. Although ${\rho }_n^{\prime }$ is simply a parameter in the soil isoline equations (Equation (4)), its meaning in the derivation of the soil isoline equations was explained in our previous study [47]. The definition of ${\rho }_n^{\prime }$ is equivalent to the difference VI (DVI) [62] or, more rigorously, to the weighted DVI (WDVI), fully accounting for the soil line slope and offset [63].

2.3. A General Form of the Assumed VI Model

A general form of the VI equation was employed in this study to cover a variety of two-band VI models. We employed a form expressed as the ratio of two linear sums of the red and NIR bands with a constant term. This form may be expressed by

$$\begin{array}{ccc}\hfill v& =& F(\mathit{Q},\mathit{\rho })\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill & =& q_0\frac{{\mathit{q}}^{\mathrm{U}}\mathit{\rho }+q_3^{\mathrm{U}}}{{\mathit{q}}^{\mathrm{D}}\mathit{\rho }+q_3^{\mathrm{D}}},\hfill \end{array}$$

(7)

where v represents a VI value, and a vector Q is composed of seven coefficients specific to each VI model. Note that the row vector ${\mathit{q}}^z=(q_1^z,q_2^z)$ with index z meaning $\mathrm{U}$ or $\mathrm{D}$ is defined to simplify its representation. The coefficients are determined by the specific choice of VI, such as the NDVI; hence, these values characterize the performances of the VIs. The derivation of the intersensor VI relationships will be introduced based on this general form. Although this form cannot cover all two-band VIs, some of the used two-band VIs may be represented by choosing the coefficients of Q. Several examples of well-known VIs are summarized in

Table 1. The seven coefficients of $\mathit{Q}$ provide a generalized form of the ratio-based two-band vegetation indices (VIs). Note that the vector ${\mathit{q}}^z=(q_1^z,q_2^z)$ of Equation (7) (where z represents $\mathrm{U}$ and $\mathrm{D}$) is defined to simplify the VI representation.

v	${\mathit{q}}_0$	${\mathit{q}}_1^{\mathbf{U}}$	${\mathit{q}}_2^{\mathbf{U}}$	${\mathit{q}}_3^{\mathbf{U}}$	${\mathit{q}}_1^{\mathbf{D}}$	${\mathit{q}}_2^{\mathbf{D}}$	${\mathit{q}}_3^{\mathbf{D}}$
NDVI	1	$-1$	1	0	1	1	0
DVI	1	$-1$	1	0	0	0	1
SAVI	$1.5$	$-1$	1	0	1	1	$0.5$
EVI2	$1.5$	$-1$	1	0	$2.4$	1	1

.

3. Derivation Steps for Obtaining the Intersensor VI Relationships

This section symbolically describes the derivation steps applied to obtain the intersensor VI relationships based on the soil isoline equations. The derivation involves the non-unique choice of variables and relationships and the consideration of several terms in the derivations of individual cases. A variety of choices during the derivation can decrease the clarity of the model. These choices may be avoided by implementing the symbolic form introduced in this section.

This study assumed that the sensor differences were attributed only to the wavelengths during the derivation and numerical simulations. Our intention was to eliminate the influences of the spatial resolution, the viewing and illumination geometries, and so on, from the influences of the sampling wavelength. Similar approaches have been explored in our pilot studies [52,53,54,55].

The first step of the derivation involves the elimination of one of the two reflectances. This step assumes that an intersensor VI relationship may be obtained from each soil isoline. In other words, we are attempting to characterize the relationship along a soil isoline. This approach enables us to use the relationship between the red and NIR reflectances expressed by Equation (4). This first step is explained in the following subsection.

3.1. Relationship between the VI and ${\rho }_n^{\prime }$

In this subsection, we derive the VI along with a single soil isoline as a function of the parameter ${\rho }_n^{\prime }$. This derivation may be performed simply by substituting the soil isoline equations, Equation (4), into the red and NIR reflectances of the general form of the VI model equation, Equation (6). As a result, Equation (6) becomes

$$\begin{array}{c}\hfill v=F_v(\mathit{Q},{\rho }_n^{\prime }).\end{array}$$

(8)

3.2. Symbolic Form of the Intersensor VI Relationship

In this subsection, we explicitly differentiate two sensors, denoted sensor A and sensor B. The VI model equation, Equation (8), is applied to two sensors, indicated by the subscripts A and B throughout this study. Two definitions of the VI model equations are used here:

$$\begin{array}{ccc}\hfill v_A& =& F_{vA}\left({\rho }_{nA}^{\prime }\right),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill v_B& =& F_{vB}\left({\rho }_{nB}^{\prime }\right),\hfill \end{array}$$

(10)

where the vector Q is omitted for brevity.

Our goal is to directly relate $v_A$ and $v_B$, defined by Equations (9) and (10). To this end, we will require one condition in addition to Equations (9) and (10). From a physical point of view, the required relationship is one that relates the variables of the two sensors. Our choice of condition is the relationship between the parameters of sensor A (${\rho }_{nA}^{\prime }$) and sensor B (${\rho }_{nB}^{\prime }$), represented by

$$\begin{array}{c}\hfill {\rho }_{nB}^{\prime }=f_{AB}\left({\rho }_{nA}^{\prime }\right).\end{array}$$

(11)

In summary, three relationships, Equations (9)–(11), were the key to the derivation of the relationship between $v_A$ and $v_B$.

Equation (9) was solved for $v_A$ symbolically,

$$\begin{array}{c}\hfill {\rho }_{nA}^{\prime }=F_{vA}^{-1}\left(v_A\right).\end{array}$$

(12)

The relationship between $v_A$ and $v_B$ was given by

$$\begin{array}{c}\hfill v_B=F_{vB}\circ f_{AB}\circ F_{vA}^{-1}\left(v_A\right),\end{array}$$

(13)

where the circle “∘” represents the composition of functions. Equation (13) describes the relationship between $v_A$ and $v_B$, which is one objective of this study. This form may be further modified to relate $v_B$ to the spectrum of the original sensor ${\mathit{\rho }}_A$ by inserting Equation (6) into Equation (13),

$$\begin{array}{c}\hfill v_B=F_{vB}\circ f_{AB}\circ F_{vA}^{-1}\circ F\left({\mathit{\rho }}_A\right).\end{array}$$

(14)

Note that Equation (14) directly relates the reflectance spectrum of the original sensor (${\mathit{\rho }}_A$) to a VI value of the destination sensor ($v_B$), which enables us to understand the derivation steps. These derivation steps are summarized in the flowchart shown in Figure 1.

4. Practical Consideration of the Intersensor VI Relationship

The derivation explained symbolically in the previous section involves the inversion of one of the functions, $F_{vA}$. Thus, the analytical approach is restricted by the choice of function in practical applications. Specifically, although $F_{vA}$ may be approximated using a geometric series, the order of the polynomial is limited to a small number. Practical considerations must be addressed to ensure model accuracy while limiting the order of the polynomial to relatively low values, preferably, to first-order polynomials. This section addresses several points that should be treated carefully for practical applications.

4.1. Treatment of Higher-Order Terms of Soil Isoline Equations

We used a simple analytical approach to achieve the derivation of the intersensor VI relationship. The approach is based on an operation that separates linear and non-linear terms of soil isoline equations. Let us redefine soil isoline equations as

$$\begin{array}{c}\hfill \mathit{\rho }=\left(\begin{array}{cc}{\mathit{c}}_0,& {\mathit{c}}_1\end{array}\right){\mathit{\rho }}_{n,L}^{\prime }+\left(\begin{array}{ccc}{\mathit{c}}_2,& {\mathit{c}}_3,& \cdots \end{array}\right){\mathit{\rho }}_{n,H}^{\prime },\end{array}$$

(15)

where ${\mathit{\rho }}_{n,L}^{\prime }={}^t\left(\begin{array}{cc}1,& {\rho }_n^{\prime }\end{array}\right)$ and ${\mathit{\rho }}_{n,H}^{\prime }={}^t\left(\begin{array}{ccc}{\rho }_n^{\prime 2},& {\rho }_n^{\prime 3},& \cdots \end{array}\right)$. One can concatenate ${}^t{\mathit{\rho }}_{n,L}^{\prime }$ with ${}^t{\mathit{\rho }}_{n,H}^{\prime }$ to yield ${}^t{\mathit{\rho }}_n^{\prime }$.

In order to derive the soil isoline-based VI relationship, the present study deals with the higher-order terms as a residual term.

$$\begin{array}{c}\hfill \mathit{\epsilon }=\left(\begin{array}{ccc}{\mathit{c}}_2,& {\mathit{c}}_3,& \cdots \end{array}\right){\mathit{\rho }}_{n,H}^{\prime }.\end{array}$$

(16)

Because it is considered as a residual vector, the above higher-order term is included in the zeroth-order term of the isoline equations.

$$\begin{array}{c}\hfill \mathit{\rho }=\left(\begin{array}{cc}{\mathit{c}}_0+\mathit{\epsilon },& {\mathit{c}}_1\end{array}\right){\mathit{\rho }}_{n,L}^{\prime }.\end{array}$$

(17)

Equation (17) can be interpreted as an implicit form of Equation (3), which is a matrix-vector representation of soil isoline equations in the red-NIR reflectance subspace.

Subsequently, one can yield an equation corresponding to Equation (8) from Equations (6) and (17):

$$\begin{array}{c}\hfill v=q_0\frac{{\gamma }_0^{\mathrm{U}\ast }+{\gamma }_1^{\mathrm{U}}{\rho }_n^{\prime }}{{\gamma }_0^{\mathrm{D}\ast }+{\gamma }_1^{\mathrm{D}}{\rho }_n^{\prime }},\end{array}$$

(18)

where

$$\begin{array}{ccc}\hfill {\gamma }_i^z& =& {\mathit{q}}^{\mathrm{z}}{\mathit{c}}_i+q_3^z{\delta }_{0i},\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\gamma }_0^{z\ast }& =& {\mathit{q}}^z\mathit{\epsilon }+{\gamma }_0^z.\hfill \end{array}$$

(20)

Recall that the vectors $q_0$ and ${\mathit{q}}^z=(q_1^z,q_2^z)$ are both vegetation index model parameters. The former variable is a scaling factor of the vegetation index, while the latter vector contains two elements, with index z meaning $\mathrm{U}$ or $\mathrm{D}$ as its numerator and denominator part, respectively. The coefficients ${\gamma }_0^{z\ast }$ are not generally constant due to the non-linear effect of ${\rho }_n^{\prime }$ in the residual term against the variability of the pixel condition.

Note that the second- and higher-order terms included in the zeroth-order term could be computed numerically for the original sensor (sensor A) from the physical definition of the parameter ${\rho }_n^{\prime }$, which is discussed in our previous work [47],

$$\begin{array}{c}\hfill {\rho }_n^{\prime }=-\sin \left(\theta \right){\rho }_r+\cos \left(\theta \right)({\rho }_n-s_0).\end{array}$$

(21)

Furthermore, the value of ${\rho }_n^{\prime }$ obtained from the destination sensor, ${\rho }_{nB}^{\prime }$, could be estimated from the known variable ${\rho }_{nA}^{\prime }$ based on the relationship between ${\rho }_{nA}^{\prime }$ and ${\rho }_{nB}^{\prime }$, as represented by the function $f_{AB}$ of Equation (11).

Equation (18) is applicable to both sensor A and sensor B; however, the coefficients represented by the $\gamma$s are unique to each sensor because both coefficients are actually functions of the coefficients $p_i$ and, hence, the values of ${\alpha }_i$ and ${\beta }_i$ in the soil isoline equations, Equation (4). Therefore, Equation (18) should be distinct from the equation obtained from the other sensors, as achieved using the subscripts A and B for sensor A and sensor B, respectively. We now have

$$\begin{array}{ccc}\hfill v_A=F_{vA}\left({\rho }_{nA}^{\prime }\right)& =& q_0\frac{{\gamma }_{0A}^{\mathrm{U}\ast }+{\gamma }_{1A}^{\mathrm{U}}{\rho }_{nA}^{\prime }}{{\gamma }_{0A}^{\mathrm{D}\ast }+{\gamma }_{1A}^{\mathrm{D}}{\rho }_{nA}^{\prime }},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill v_B=F_{vB}\left({\rho }_{nB}^{\prime }\right)& =& q_0\frac{{\gamma }_{0B}^{\mathrm{U}\ast }+{\gamma }_{1B}^{\mathrm{U}}{\rho }_{nB}^{\prime }}{{\gamma }_{0B}^{\mathrm{D}\ast }+{\gamma }_{1B}^{\mathrm{D}}{\rho }_{nB}^{\prime }}.\hfill \end{array}$$

(23)

The inverse of $F_{vA}$ may be simply defined by solving Equation (22) for ${\rho }_{nA}^{\prime }$,

$$\begin{array}{ccc}\hfill {\rho }_{nA}^{\prime }& =& F_{vA}^{-1}\left(v_A\right)=\frac{-{\gamma }_{0A}^{\mathrm{D}\ast }{v}_A+q_0{\gamma }_{0A}^{\mathrm{U}\ast }}{{\gamma }_{1A}^{\mathrm{D}}{v}_A-q_0{\gamma }_{1A}^{\mathrm{U}}}.\hfill \end{array}$$

(24)

This function will be used later in this section.

4.2. Intersensor Relationship between the Parameter ${\rho }_n^{\prime }$

Our next step is to prepare the relationship between ${\rho }_{nA}^{\prime }$ and ${\rho }_{nB}^{\prime }$, represented by the function $f_{AB}$. In this study, the relationship is represented as a polynomial with s separation of linear and residual components. These assumptions on linkage between parameter ${\rho }_n^{\prime }$ in soil isoline equations in red-NIR space permits the following equations:

$$\begin{array}{ccc}\hfill f_{AB}\left({\rho }_{nA}^{\prime }\right)& =& \sum _{i=0}^{\infty }{u}_i{\rho }_{nA}^{\prime i}\hfill \end{array}$$

(25)

$$\begin{array}{cc}=& u_0+u_1{\rho }_{nA}^{\prime }+{\epsilon }_u\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill & =& u_0^{\ast }+u_1{\rho }_{nA}^{\prime },\hfill \end{array}$$

(27)

where ${\epsilon }_u$ represents the second and higher order term, and $u_0^{\ast }$ is defined by

$$\begin{array}{c}\hfill u_0^{\ast }=u_0+{\epsilon }_u.\end{array}$$

(28)

Recall that the relationship between the red and NIR reflectance spectra measured by a specific sensor may be approximated by a soil isoline equation for a soil spectrum that remains constant throughout the numerical simulations. The coefficients $u_i$ of Equation (25) depend on the soil spectrum because ${\rho }_n^{\prime }$ is defined along with a soil isoline whose coefficients show variation with the changes of the soil spectrum. This analogy suggests that the accuracy of Equation (25) may be improved by expressing the coefficients as functions of the soil spectrum.

4.3. Intersensor VI Relationship

According to the above two subsections, one can yield an intersensor VI relationship via the composition of three functions shown in Figure 1. Before the description of the practical derivation step, we introduce and define a list $G_{\epsilon }$ to represent three residual terms, namely, vectors ${\mathit{\epsilon }}_A$ and ${\mathit{\epsilon }}_B$ and a scalar ${\epsilon }_u$. The vectors ${\mathit{\epsilon }}_A$ and ${\mathit{\epsilon }}_B$ represent residual terms of $\mathit{\epsilon }$ in Equation (16) for sensor A and B, respectively. The scalar ${\epsilon }_u$ is the residue in Equation (28). The definition of the list is represented by $G_{\epsilon }=({\mathit{\epsilon }}_A,{\epsilon }_u,{\mathit{\epsilon }}_B)$ throughout this study. Note that $G_{\epsilon }$ is not a vector with a fixed length or having a certain order of components. Rather, $G_{\epsilon }$ simply represents a combination of the three residual terms.

From Equations (23), (24) and (27), the intersensor VI form is derived as a fractional form with four functions of the residual list $G_{\epsilon }$:

$$\begin{array}{c}\hfill v_B=q_0\frac{{\psi }^{\mathrm{UD}\ast }\left(G_{\epsilon }\right)v_A-q_0{\psi }^{\mathrm{UU}\ast }\left(G_{\epsilon }\right)}{{\psi }^{\mathrm{DD}\ast }\left(G_{\epsilon }\right)v_A-q_0{\psi }^{\mathrm{DU}\ast }\left(G_{\epsilon }\right)},\end{array}$$

(29)

where four ${\psi }^{\mathrm{xy}\ast }$ functions with $\mathrm{x}$ and $\mathrm{y}$ for indices of $\mathrm{U}$ or $\mathrm{D}$ are written as the following expression.

$$\begin{array}{c}\hfill {\psi }^{\mathrm{xy}\ast }\left(G_{\epsilon }\right)={\gamma }_{1A}^{\mathrm{y}}({\gamma }_{0B}^{\mathrm{x}\ast }+{\gamma }_{1B}^{\mathrm{x}}{u}_0^{\ast })-{\gamma }_{1B}^{\mathrm{x}}{u}_1{\gamma }_{0A}^{\mathrm{y}\ast }\end{array}$$

(30)

The residual terms in list $G_{\epsilon }$ are implicitly included in ${\gamma }_{0A}^{\mathrm{z}\ast }$, ${\gamma }_{0B}^{\mathrm{z}\ast }$, and $u_0^{\ast }$. One can discriminate the resultant function ${\psi }^{\mathrm{xy}\ast }\left(G_{\epsilon }\right)$ into two types of variables. One type includes pure soil-dependent variables, while another corresponds to reflectance-dependent variables. In other words, the latter is the effect of higher-order terms in soil isoline equations in the intersensor VI relation model, which are mathematically defined as elements of $G_{\epsilon }$ in Equation (30). The effect in the VI relation model is unavoidable factors propagated from soil isoline non-linearity under a separating operation for soil isoline equations proposed in this study.

There are several characteristics in the derived VI relationship. The residual-based VI relationship keeps its rational function under the condition that all the residual terms in $G_{\epsilon }$ are zero. The finding shows that the nonlinearity of the intersensor VI relationship obtained from soil isoline equations without higher-order terms holds from a mathematical point of view. In contrast, such nonlinearity doe not occur for VI models whose coefficients is ${\mathit{q}}^{\mathrm{D}}=\mathbf{0}$. The VI group is the so-called distance-based index, such as DVI, which meets ${\psi }^{\mathrm{DD}\ast }(G_{\epsilon }=\mathbf{0})=0$ from its definition.

It is the case that ${\psi }^{\mathrm{DD}\ast }\left(G_{\epsilon }\right)$ is constant zero for the DVI; nevertheless, the other three coefficients retain their unique isoline parameters. Therefore, it is worth noting that the potential for analysis of relationships between distance-based indices derived from distinct sensors also remains, as well as ratio-based indices such as NDVI.

4.4. Intersensor VI Relationship after Approximation of $G_{\epsilon }$

As the goal of the present study, we finally examine the effect on the approximation of higher-order terms in soil isoline equations. The effect can be easily investigated by introducing an operation that replaces residues, $G_{\epsilon }$ with a function of ${\rho }_{nA}^{\prime }$. The operation may lead to a practical form of the intersensor VI relationship. During the replacement of residues for sensor B (${\mathit{\epsilon }}_B$), it is necessary to further replace ${\rho }_{nB}^{\prime }$ as ${\rho }_{nA}^{\prime }$ because we cannot use any product derived from reflectances of the destination sensor B in operational cases.

When translation from observed data from sensor A to that from sensor B is performed, residual sets need to be approximated as functions of ${\rho }_{nA}^{\prime }$. In this study, we deal with three types of residual terms in ${\psi }^{\mathrm{xy}\ast }\left(G_{\epsilon }\right)$ as introducing truncation orders $N_1$ for soil isoline equations of Equation (4) and $N_2$ for a relationship between ${\rho }_n^{\prime }$ in Equation (25). Appendix B is allocated for details on their treatment.

The intersensor VI relationship after replacement of $G_{\epsilon }$ with functions of ${\rho }_{nA}^{\prime }$ is obtained as the following:

$$\begin{array}{c}\hfill {\widehat{v}}_B=q_0\frac{{\psi }_{(N_1,N_2)}^{\mathrm{UD}}\left({\rho }_{nA}^{\prime }\right)v_A-q_0{\psi }_{(N_1,N_2)}^{\mathrm{UU}}\left({\rho }_{nA}^{\prime }\right)}{{\psi }_{(N_1,N_2)}^{\mathrm{DD}}\left({\rho }_{nA}^{\prime }\right)v_A-q_0{\psi }_{(N_1,N_2)}^{\mathrm{DU}}\left({\rho }_{nA}^{\prime }\right)},\end{array}$$

(31)

where

$$\begin{array}{c}\hfill {\psi }_{(N_1,N_2)}^{\mathrm{xy}}\left({\rho }_{nA}^{\prime }\right)={\gamma }_{1A}^{\mathrm{y}}(({\gamma }_{0B}^{\mathrm{x}}+f_{{\mathit{\epsilon }}_B(N_1,N_2)}^{\mathrm{x}}\left({\rho }_{nA}^{\prime }\right))+{\gamma }_{1B}^{\mathrm{x}}(u_0+f_{{\epsilon }_u\left(N_2\right)}\left({\rho }_{nA}^{\prime }\right)))-{\gamma }_{1B}^{\mathrm{x}}{u}_1({\gamma }_{0A}^{\mathrm{y}}+f_{{\mathit{\epsilon }}_A\left(N_1\right)}^{\mathrm{y}}\left({\rho }_{nA}^{\prime }\right)).\end{array}$$

The derivation steps for three residual terms approximated by ${\rho }_{n,A}^{\prime }$, that is, $f_{{\mathit{\epsilon }}_B(N_1,N_2)}^{\mathrm{x}}$, $f_{{\epsilon }_u\left(N_2\right)}$, and $f_{{\mathit{\epsilon }}_A\left(N_1\right)}^{\mathrm{y}}$ are described in Appendix B.

A special case is the relationship for $N_1=N_2=1$. The scenario means that ${\psi }_{(1,1)}^{\mathrm{xy}}\left({\rho }_{nA}^{\prime }\right)$ is equivalent to ${\psi }^{\mathrm{xy}\ast }(G_{\epsilon }=\mathbf{0})$, where all four coefficients relating $v_A$ with $v_B$ are fully independent of ${\rho }_{nA}^{\prime }$. In other words, the relationship is constructed by coefficients as a function of soil brightness, $R_s$. Therefore, this drastic case has a clue to explaining the dependency of soil on the intersensor VI relationship, which is an issue that the present study attempts to address.

5. Results of Numerical Simulations

This section demonstrates the validity of the relationships derived and introduced in this study. All uncertainties associated with the spectral acquisition were eliminated by applying the framework to a data set calculated numerically using radiative transfer (RT) models of the vegetation canopy and leaf layers. The validity of the derived relationships was then investigated by comparing the VI values of the two sensors to the VI values translated from the other sensor. The translation from the VI value of one sensor to the value of the other sensor was performed by applying the relationship derived using the coefficients computed from the coefficients of the soil isolines. The following subsections describe the numerical simulation conditions (parameter settings) and obtained results.

5.1. Numerical Simulations of the Inter-VI Relationships

A combined leaf and canopy RT model, PROSAIL, was employed to simulate the top-of-the-canopy (TOC) reflectance spectra under a variety of conditions. The numerically obtained TOC reflectance spectra were then used to simulate the band reflectances of four hypothetical sensors by selecting different pairs of wavelengths as the red and NIR bands (shown below). One of the four sensors was assigned as the “original” sensor (sensor A), the VI values of which were translated into the values of the other three sensors (sensor B), denoted as the “destination” sensors during the derivation.

The numerical simulations were conducted by selecting three different parameters to characterize the VI differences. The first parameter was the leaf area index (LAI), a representative biophysical parameter. The value of the LAI varied from 0 to 4.0 at 0.5 intervals. The second parameter was the soil reflectance spectrum, which can disrupt the VI values. The soil reflectance spectra were obtained by linearly combining the wet and dry soil spectra stored in the PROSAIL model. The third parameter was the wavelength, which introduced differences in the band reflectances measured by the different sensors, thereby producing differences in the VI values as well. The other parameters in the leaf model PROSPECT and the canopy model SAIL were held constant during the experiments.

The parameters that characterized the leaf chemical content in the PROSPECT model were fixed at the standard (average) values reported in [47]. The leaf angle distribution, an input of the SAIL model, was fixed at a spherical distribution. Under these assumptions, the LAI was the only parameter of biophysical properties in this study.

The band center wavelengths of the four hypothetical sensors were set to values equal to the band center wavelengths of the GOSAT-CAI [64], LANDSAT8-OLI [65], Suomi NPP-VIIRS [66], and Terra-MODIS [67,68]. The wavelength pairs of the red ${\rho }_r$ and NIR bands ${\rho }_n$ of the four hypothetical sensors were (${\rho }_r$, ${\rho }_n$) = (674, 870) for sensor A, (655, 865) for sensor B1, (672, 865) for sensor B2, and (645, 869) for sensor B3, in [nm]. Figure 2 shows the reflectance spectra and soil isolines retrieved from sensor A in the red and NIR reflectance space. The figure shows good agreement among the retrieved soil isolines, which will be used to translate the VI values. The detailed steps used to obtain the soil isoline retrievals have been described previously [47].

The translation of the VI values from one sensor to another sensor was performed by regarding the CAI-like sensor (sensor A) as the original sensor, and the other three sensors (sensor B1, B2, and B3) were the destination sensors. We considered three translation cases from sensor A (original) to sensors B1, B2, and B3 (destinations). We labeled these three pairs of sensors as translation case 1, case 2, and case 3, respectively. Again, the translation process involved translating a set of reflectance spectra (of the two bands) from sensor A (CAI-like sensor); the VI values were then translated into sensor B1 (case 1), sensor B2 (case 2), and sensor B3 (case 3). The VI translation was performed in all cases simply by using the derived VI relationships. The validity of the derived relationship was then evaluated based on a comparison between the simulated VIs obtained from sensor B and the translated VIs obtained from sensor A.

5.2. Dependence of the Coefficient $\psi$ on Soil Reflectance

We firstly investigated the effect of soil reflectances on the intersensor VI relationship. In the context of derivation results introduced in the present work, we distinguished two types of translation coefficients. One is residual-based translation with ${\psi }^{\mathrm{xy}\ast }\left(G_{\epsilon }\right)$, and another is its modification/approximation by a usable parameter ${\rho }_n^{\prime }$ with ${\psi }_{(N_1,N_2)}^{\mathrm{xy}}\left({\rho }_{nA}^{\prime }\right)$. The difference between two coefficients is investigated to evaluate the effectiveness of the derived VI relationship after modification.

$$\begin{array}{c}\hfill {\delta }_{(N_1,N_2)}^{\mathrm{xy}}={\psi }^{\mathrm{xy}\ast }\left(G_{\epsilon }\right)-{\psi }_{(N_1,N_2)}^{\mathrm{xy}}\left({\rho }_{nA}^{\prime }\right)\end{array}$$

(32)

Figure 3 shows translation coefficients as a function of soil reflectance. In the figure, results of the four coefficients represented by ${\psi }^{\mathrm{xy}\ast }\left(G_{\epsilon }\right)$ for NDVI relations under three sensor pairs are plotted. Solid lines and gray x marks depict trends of linear and nonlinear components of ${\psi }^{\mathrm{xy}\ast }\left(G_{\epsilon }\right)$, respectively. As shown in the figure, all four coefficients varied mostly linearly with the soil brightness variations. In addition, one can see small residual-induced variations in individual soil brightness relative to those of changing soil brightness conditions, except for ${\psi }^{\mathrm{DD}\ast }$ in case 2.

Solid lines for the y-axis on the left depict coefficients under the special condition of a residual $G_{\epsilon }$ as zero. It means there exists the effect of only pure soil variation in the function ${\psi }^{\mathrm{xy}\ast }$. In all panels in the figure, a single line is observed because no vegetation variable is contained in the component. As can be seen, differences in Equation (32) at the y-axis on the right side are relatively small compared with pure the soil-oriented effect in intersensor VI relationship at that y-axis on the left side. The comparison result of variation ranges implies the importance of soil consideration for intersensor VI calibration.

Various truncation scenarios are also examined (Figure 4). Figure 4 shows a plot of differences ${\delta }_{(N_1,N_2)}^{\mathrm{xy}}$ in Equation (32) for a sensor pair for case 1. The figure is constructed with four panels. The panels located at the upper-left, upper-right, lower-left, and lower-right positions correspond to ${\delta }^{\mathrm{UD}}$, ${\delta }^{\mathrm{UU}}$, ${\delta }^{\mathrm{DD}}$, and ${\delta }^{\mathrm{DU}}$, respectively. The individual panels show results for ${\delta }^{\mathrm{xy}}$ with truncation combinations ($N_1,N_2$). The gray, red, and green plots correspond to ($N_1,N_2$) that equals (1,2), (2,1), and (2,2), and solid and dashed lines are for LAIs of 1.0 and 3.0, respectively. In either case, we observe higher accuracy in the case of symmetric (2,2) compared to asymmetric cases such as (1,2) and (2,1). This is expected from the fact that the accuracy of soil isolines depends on the order of polynomials.

These results indicate that the intersensor VI relationship depended on the soil reflectance beneath the canopy layer, a representative parameter of the land surface conditions. The implication of these results are serious (and important): Whn determining the regression coefficients without accounting for the variations in the soil spectrum, the intersensor VI model could suffer from the variations in the soil spectrum.

5.3. Accuracy of the Derived Translation Function

The validity of the derived relationships with truncation orders, $N_1$ and $N_2$, were evaluated by comparing the VI values of the two sensors before and after the translation. The VI value translated from sensor A to sensor B was denoted by ${\overline{v}}_B$, which was computed according to Equation (29). The aim of this section is to validate the derived expressions by confirming the accuracy improvement of $v_B$ after the translation from $v_A$ (hence ${\overline{v}}_B$). The comparison between ($v_A-v_B$) and ($v_B-{\overline{v}}_B$) was facilitated by representing the differences as ${\epsilon }_v$ and ${\overline{\epsilon }}_v$, respectively, such as

$$\begin{array}{ccc}\hfill {\epsilon }_v& =& v_B-v_A,\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill {\overline{\epsilon }}_v& =& v_B-{\overline{v}}_B,\hfill \end{array}$$

(34)

where $v_A$ is the VI value of the CAI-like sensor A, and $v_B$ represents the VI value obtained from the destination sensors, namely, OLI-like sensor B1, VIIRS-like sensor B2, or MODIS-like sensor-B3. The values of ${\epsilon }_v$ and ${\overline{\epsilon }}_v$ were negative if the VI values measured using the original sensor ($v_A$) and the values obtained from the translated VI values ($\overline{v_B}$) exceeded $v_B$. The latter case also indicated that translation by Equation (29) overestimated the value of the destination sensor.

Figure 5 summarizes the plots of ${\epsilon }_v$ and ${\overline{\epsilon }}_v$ for the NDVI as a function of a LAI characterized by one of two different soil spectra (wet or dry conditions). The calculated values of ${\epsilon }_v$ and ${\overline{\epsilon }}_v$ are indicated by the gray area and lines, respectively. The differences between the soil reflectance spectra of the wet and dry soils are represented by the solid and dashed lines, respectively. The upper three plots correspond to case 1 (CAI vs. OLI) and retain terms $N_1$ up to the first-, second-, or third-order terms of a polynomial in the soil isoline equations, while red and green lines distinguish the other terms $N_2$ up to the first and third orders. The middle and bottom rows of the figure show the results obtained from case 2 (CAI vs. VIIRS) and case 3 (CAI vs. MODIS), respectively.

This figure clearly shows that the translation error ${\overline{\epsilon }}_v$ (lines) was smaller than the original differences ${\epsilon }_v$ (gray colored zones) across nearly the entire LAI range in all cases. This fact indicates the validity of the expression derived in this study.

A comparison of overall accuracy for sensor-to-sensor NDVI translations by intersensor relationships obtained from various orders ($N_1$,$N_2$) is shown in Figure 6. The bar plot visualizes trends of the root mean square error (RMSE) for overall LAI and soil mixture conditions. From the figure, it is revealed that the translation performance enhances with increments of higher-order terms taken into account in the translation. Furthermore, the translation performance becomes more stable with both orders ($N_1$ and $N_2$) considered as second-order terms.

The trends in the other VI models were investigated by creating similar plots for the SAVI, EVI2, and DVI, as shown in Figure 7. These trends were similar to those obtained using the NDVI, as observed in the SAVI for soil isolines approximated using a first-order polynomial (left column) or a third-order polynomial (second-left column). The third column and the left column plot the results obtained from the use of the EVI2 and DVI, approximated using third-order soil isolines, respectively. Regardless of VI types, accuracy improvement is observed in the numerical simulation environment.

Table 2. Normalized root mean square error in the derived translation functions (in percent).

VI		Combination of Orders (${\mathit{N}}_1,{\mathit{N}}_2$)
NDVI		(1,1)	(1,3)	(3,1)	(2,2)	(3,3)
	Case 1	20.6	9.6	15.5	1.9	0.6
	Case 2	10.9	18.3	11.7	2.2	1.2
	Case 3	17.3	5.3	12.6	1.3	0.3
SAVI		(1,1)	(1,3)	(3,1)	(2,2)	(3,3)
	Case 1	29.5	11.1	17.9	3.3	0.9
	Case 2	12.1	6.2	11.7	2.9	1.1
	Case 3	22.0	7.5	13.3	1.8	0.4
EVI2		(1,1)	(1,3)	(3,1)	(2,2)	(3,3)
	Case 1	34.0	13.0	20.2	4.1	1.1
	Case 2	13.6	6.3	13.3	3.5	1.3
	Case 3	24.6	8.7	14.3	2.2	0.5
DVI		(1,1)	(1,3)	(3,1)	(2,2)	(3,3)
	Case 1	30.9	14.8	18.0	5.3	1.1
	Case 2	12.4	3.6	9.9	3.2	0.7
	Case 3	22.0	10.1	13.5	3.0	0.6

summarizes the RMSE for the four VI values over the entire range of the LAI and soil brightness. The table shows normalized values by the original difference for the five truncation scenarios, along with the three cases individually. These results show that the original differences ($\epsilon$) decreased significantly after the translation ($\overline{\epsilon }$), especially increasing the retained terms of ${\rho }_{nA}^{\prime }$, thereby validating the derived expressions.

Finally, we examined the land cover dependence on the derived intersensor VI relationships. Figure 8 shows the normalized RMSE of NDVI under various combinations of chlorophyll content (20, 40, and 80 $\mathsf{\mu }{\mathrm{g}/\mathrm{cm}}^2$) and leaf angle distribution (planophile, erectophile, extremophile, and uniform distributions). Although variation of errors under various conditions inner a canopy layer, such as for chlorophyll and leaf angle distribution (LAD), was observed here, the error decreased as the truncation order became higher. Thus, this implies some extent of land cover dependence on translation accuracy, but the validity of the derived relationship is consistently shown in all results. Further comprehensive investigations will be needed to further understand the land cover dependence on the derived relationships.

6. Discussion

This study explained that the intersensor VI relationship varies with the soil reflectance spectrum beneath the canopy layer. The mechanism by which this relationship varies may be inferred from the fact that the four coefficients of the intersensor VI relationship may all be expressed as functions of the soil reflectance spectrum, as a parameter. This analytical result shows the following important facts: (1) a model of the intersensor VI relationship that does not address soil surface classification may retrieve a relationship that suffers from soil surface spectrum variations within target regions; (2) the retrieved intersensor VI relationship across a region characterized by a specific soil spectrum cannot be directly applicable to other regions characterized by a different soil spectrum; (3) soil spectrum estimates using any means must be accomplished prior to intersensor calibration when the calibration is based on formulations similar to the expressions derived in this study.

A different approach to performing the intersensor calibration may rely on cross-calibrating similar bands, one by one, measured in two sensors. This approach is by far the simplest, although it has several disadvantages. First, the number of bands that must be cross-calibrated equals the number of bands required by the algorithm employed. Doing so can unfortunately be difficult for large numbers of bands. These cross-sensor band-to-band relationships vary with the land surface conditions, depending on the canopy and soil spectrum. These dependencies present the major obstacles to this type of cross-calibration. Recall that the intersensor calibration is not really the same as the absolute calibration due to differences in the band configurations (band center and width characterized by the spectral response functions). By contrast, the approach introduced in this study requires only the relationships between ${\rho }_{nA}^{\prime }$ and ${\rho }_{nB}^{\prime }$, instead of requiring all the relationships of similar bands measured in the two sensors. The simpler approach described here is advantageous if a product algorithm requires a greater number of bands.

One drawback to the present approach is that the accuracy of the derived relationship may depend on the accuracy of the model of the relationship between ${\rho }_{nA}^{\prime }$ and ${\rho }_{nB}^{\prime }$. This relationship is expressed by a polynomial of ${\rho }_{nA}^{\prime }$ and must be solved for ${\rho }_{nB}^{\prime }$ afterward. The use of more suitable functional forms might improve the accuracy of the intersensor VI relationship; however, such steps would require additional approximations at some point during the derivation. From a practical perspective, the accompanying loss of model simplicity reduces the practicality of the model. Therefore, there is a trade-off between practicality and accuracy in modeling a relationship between ${\rho }_{nA}^{\prime }$ and ${\rho }_{nB}^{\prime }$. The choice of relationship must be optimized for ${\rho }_n^{\prime }$, which would require a significant effort beyond the scope of the present discussion. This matter is worth investigating in a separate study.

This study focused on the influence of the soil spectrum as a source of variation. The derived VI relationships will also be influenced by the biophysical parameters. This issue should be investigated thoroughly in a future study.

7. Conclusions

This study introduced an analytical technique for relating the VIs measured by two sensors. The relationship was derived using the soil isoline equation, which consists of a set of reflectance spectra obtained under a constant soil spectrum. First, the relationship between the VIs measured by two sensors was explained conceptually. The conceptual relationship was applied to several realistic cases by truncating the order of the polynomials used to describe the soil isoline equations. Finally, the derived relationships were numerically examined using common radiative transfer models. The results validated the derived expressions and the applicability of the truncations. The derived VI relationships and the numerical results also indicated an important fact: that the intersensor relationship among the measured VIs can be influenced by the soil reflectance spectrum. Although this implication may be obtained from numerical simulations, the derivations introduced in this study confirmed this fact analytically.

Further studies are needed to enhance the practicality of the proposed method by combining estimation processing of translation parameters (soil reflectance in this study) with isoline-based translation functions. Such efforts should be devoted mainly to solving the inverse problem of soil isoline equations. This investigation is left to future studies.