Wavelength Extension of the Optimized Asymmetric-Order Vegetation Isoline Equation to Cover the Range from Visible to Near-Infrared

: Vegetation isoline equations describe analytical relationships between two reﬂectances of different wavelengths. Their applications range from retrievals of biophysical parameters to the derivation of the inter-sensor relationships of spectral vegetation indexes. Among the three variants of vegetation isoline equations introduced thus far, the optimized asymmetric-order vegetation isoline equation is the newest and is known to be the most accurate. This accuracy assessment, however, has been performed only for the wavelength pair of red and near-infrared (NIR) bands ﬁxed at ∼ 655 nm and ∼ 865 nm, respectively. The objective of this study is to extend this wavelength limitation. An accuracy assessment was therefore performed over a wider range of wavelengths, from 400 to 1200 nm. The optimized asymmetric-order vegetation isoline equation was conﬁrmed to demonstrate the highest accuracy among the three isolines for all the investigated wavelength pairs. The second-best equation, the asymmetric-order isoline equation, which does not include an optimization factor, was not superior to the least-accurate equation (i.e., the ﬁrst-order isoline equation) in some cases. This tendency was prominent when the reﬂectances of the two wavelengths were similar. By contrast, the optimized asymmetric-order vegetation isoline showed stable performance throughout this study. A single factor introduced into the optimized asymmetric-order isoline equation was concluded to effectively reduce errors in the isoline for all the wavelength combinations examined in this study. The results show the superiority of the optimized asymmetric-order vegetation isoline over the other two isolines.


Introduction
Remotely sensed reflectance spectra have been used to estimate biophysical parameters for terrestrial vegetation [1,2]. The simplest and most widely accepted approach for this purpose is to use a spectral vegetation index (VI). Because VIs are often defined as algebraic manipulations of reflectance spectra, VI model equations are also considered a relationship that the reflectance spectrum should satisfy under a fixed biophysical parameter. In particular, a fundamental relationship between the reflectances of two bands (e.g., red and near-infrared (NIR) bands) has been the basis of numerous VIs [3][4][5][6][7][8][9]. This relationship is known as a vegetation isoline equation. The vegetation isoline concept has been investigated intensively over the past two decades [10][11][12][13][14][15] and is considered a useful tool [11,16,17] to explore various applications [18,19].
The concept of vegetation isolines originated from investigations examining the influence of soil brightness on VIs [3,20]. In general, soil brightness varies considerably with changes in soil moisture, surface roughness, shadow amount, and soil characteristics such as organic matter content [21][22][23][24]. These variations affect the soil background component of reflectance spectra, which leads to variations in VI values [16,17]. This effect, known as soil noise, has been a major issue in the field of land remote sensing. Studies seeking to minimize soil noise have paid particular attention analytically to a spectral trajectory observed in the red-NIR reflectance subspace under a constant canopy [3,7]. The derived concept, the vegetation isoline, has been applied to solve the issue of the influence of soil brightness, which has affected various VIs to some extent [3][4][5][6]9,25].
The applications of the vegetation isoline concept are not limited to the analyses of VIs. It can also be used to directly retrieve biophysical parameters such as the fraction of vegetation cover (FVC), leaf area index (LAI), and the chlorophyll-a and -b contents [19,[26][27][28]. In addition, vegetation isoline equations that represent an analytical relationship of two reflectances have been used in investigations of inter-sensor VI translation [18,[29][30][31][32][33] in the context of satellite constellations. The derivations of such isoline equations require several simplifications of the model, which somewhat restricts their range of applications.
A system of vegetation canopy that is assumed for deriving vegetation isoline equations consists of a canopy layer, which is considered a turbid medium, and a soil layer of a Lambertian surface underneath the canopy. For this system of layers, a top-of-thecanopy (TOC) reflectance is represented by an analytical model known as the successive order of scattering, in which the TOC reflectance is modeled by the sum of all the photons categorized by the number of interactions between the soil surface and the canopy layers. For example, the photons scattered only once by the soil surface are labeled as the first-order interaction term. Using this analytical representation of TOC reflectance, the initial attempt to derive a vegetation isoline equation was conducted by truncating the second-order and the higher-order interaction terms between the vegetation canopy and soil surfaces [3,10]. The obtained isoline equation is referred to as the first-order vegetation isoline because the isoline equation accounts only for the first-order interaction term. Although the derived isoline equation is simple and suitable for the analyses of VIs, the simplification associated with the truncation adversely affects its accuracy, especially when the soil brightness is high. Researchers would later develop methods to overcome this weakness.
In our previous work, we improved the accuracy of the first-order vegetation isoline by accounting for the higher-order interaction term during the derivation. This improvement was performed for the widely used red and NIR wavelengths. The derived form is referred to as the asymmetric-order vegetation isoline equation because the second-order interaction term is retained only in the NIR wavelength [13]. In addition, the accuracy of the asymmetric-order vegetation isoline equation was improved dramatically by the introduction of a single factor. The derived new isoline equation was named the "optimized" asymmetric-order vegetation isoline equation [14,15]. Errors in the optimized asymmetricorder isoline equation were reduced to 4% relative to the errors in the first-order isoline equation [13] and were also reduced to 20% relative to the errors in the asymmetric-order isoline equation [14]. After this improvement, the errors in the optimized asymmetricorder isoline equations were smaller in magnitude than the magnitude equivalent to the signal-to-noise ratio (SNR) of some existing sensors.
The accuracy improvements and evaluations of the vegetation isoline equations reported in the previous works were limited to the pair of red and NIR bands-specifically, those at ∼655 nm and ∼865 nm, respectively. This limitation should be extended to make the derived isoline equations more useful for a broader range of applications. One such application is hyperspectral sensors. In recent years, numerous hyperspectral sensors have been available for Earth observation, and new missions are underway [34][35][36][37][38][39][40]. Such hyperspectral sensors enable us to process more than one hundred bands, which could improve the accuracy in biophysical parameter retrievals. For these motivations, the wavelength range should be extended to a broader range of visible-to-NIR wavelengths.
In the present work, the accuracies of the three variants of the vegetation isoline equations are evaluated in the wavelength range from 400 to 1200 nm on the basis of numerical simulations by a radiative transfer model. The explanation that follows starts with the introduction of the three variants of the vegetation isoline equations. The numerical experiments and their results are then described to evaluate their accuracy in a wide range of wavelengths. Finally, the results are further discussed and the study is concluded with descriptions of some unresolved issues.

First-Order Vegetation Isoline Equation
The first derivation of a vegetation isoline equation began with the approximation of a linear relationship between reflectances measured at two distinct wavelengths, λ 1 and λ 2 . During the derivation, the well-known concept of the soil line [41] was assumed. The derived relationship is known as the first-order approximated vegetation isoline equation. During the derivation of the first-order vegetation isoline, the second-and higher-order interaction terms were truncated. The isoline equation is expressed by with γ 1 and D 1 defined as where ω is the FVC defined as a vertically projected area of vegetation canopy [42], ρ λ 1 and ρ λ 2 are the TOC reflectances of the different wavelengths, ρ vλ 1 and ρ vλ 2 are the "pure" canopy reflectances (canopy path radiance), which are independent of the soil surface beneath the canopy layer, and T 2 λ 1 and T 2 λ 2 represent the two-way transmittances of the vegetation canopy, T 2 λ , averaged over a target area with the FVC, ω, as its weight. The contribution of higher-order interaction terms is represented by 1 . The two constants, a and b, are the soil line slope and offset, respectively. Although the derived expression is simple and sufficiently accurate to analyze the performance of VIs, the truncation error represented by 1 tends to be large when the soil reflectance becomes high [10,11]. In recent years, this shortcoming has been overcome through a unique approach described in the next subsection.

Asymmetric-Order Vegetation Isoline Equation
The asymmetric-order form was derived by including the interaction terms up to the second-order in λ 2 ; by contrast, in λ 1 , only the terms up to the first-order interaction were considered [13]. The derived asymmetric-order vegetation isoline equation becomes with the following definitions: where R vλ 2 represents the bi-hemispherical reflectance of the canopy layer at the bottom surface, which appears only in the band λ 2 . The contribution of higher-order interaction terms in this case is represented by 2 , which is considered a truncation error during numerical simulation. The error of the derived isoline equation is decreased dramatically from that of the first-order form [14]. Most notably, this increase in accuracy was achieved without losing its simplicity in the analytical expression. To improve the accuracy one step further, a single factor was introduced into the derived expression, which is considered the "optimized" asymmetric-order vegetation isoline equation.

Optimized Asymmetric-Order Vegetation Isoline Equation
The optimized asymmetric-order vegetation isoline equation was obtained by introducing a single factor k into the asymmetric-order form of the equation. After neglecting the higher-order interaction term 2 in Equation (5) [14], the factor k is introduced in the third term of the right-hand side, The newly introduced factor k can be derived by solving Equation (11) for k: This representation of k implies that the factor depends on biophysical parameters, as well as soil line parameters. That is, the factor is not a constant in this context. At the same time, if we replace k by a constant, there would be an optimum value. Even though the actual k is not a constant, such an "optimum" value for k can be expected to minimize the error to some extent. Again, k in Equation (12) depends on the biophysical and soil line parameters. It is not realistic, however, to find both biophysical-and soil-dependent k in actual satellite data. To solve this problem, we attempted to find a constant factor as k for each "wavelength pair" that minimizes the error of the isoline equation on the basis of numerical simulations. In summary, the optimum constant factor for the replacement of the variable k would depend heavily on the choice of the wavelength pairs (λ 1 and λ 2 ), which is the central theme of the present study.

Demonstration of the Errors in the Vegetation Isoline Equations
In this subsection, we demonstrate the accuracy improvement of the isoline equations by comparing the three forms of the previously introduced isolines. The errors of the three isoline equations were simulated numerically by a radiative transfer (RT) model. In this demonstration, the wavelength pair λ 1 and λ 2 was fixed at 655 nm (red) and 865 nm (NIR), respectively. The PROSAIL model [43] was used to simulate both the spectral reflectance as a function of the LAI and the soil reflectance spectra (from dark to bright soil). The parameter settings of the RT model and the retrieval algorithms of the isoline parameters [13,14] will be further explained in the following section. Figure 1 shows the errors in the three forms of isolines (the top, middle, and bottom plots correspond to the first-order, asymmetric-order, and optimized asymmetric-order forms, respectively) as a function of the LAI and soil reflectance of the red wavelength represented by R sR . The left and right columns in Figure 1 correspond to the cases of FVC = 1.0 and 0.5, respectively. The error observed in the asymmetric-order form is smaller than that in the first-order form, but it increases with increasing soil brightness. On the contrary, the errors in the optimized asymmetric-order form are much smaller than those in the other two forms even in cases where the soil brightness increased. This accuracy improvement is a result of our optimization of the single factor introduced into the asymmetric-order vegetation isoline equation [14]. , asymmetric-order isoline equation (middle panels, (b,e)), and optimized asymmetric-order isoline equation (bottom panels, (c,f)) for red and NIR wavelength pairs as a function of the LAI and soil reflectance. The FVC is 1.0 for the left panels (a-c) and 0.5 for the right panels (d-f).

Parameter Settings of Simulations to Determine Vegetation Isolines
Numerical simulations of vegetation isoline equations (Equations (1), (5), and (11)) were conducted with various settings of the LAI, FVC, and soil brightness for a wide range of wavelength (band) pairs. TOC reflectance spectra were computed using a well-known canopy radiative transfer code, PROSAIL [43], which is a combination of the leaf optical properties model PROSPECT [44] and the canopy bidirectional reflectance model SAIL [45]. The parameter settings in the simulations are summarized in Table 1. The LAI was varied from 0.0 to 4.0 in increments of 0.8 (six variants). The soil brightness is represented by the factor "psoil" in the model, which was varied from 0.0 to 1.0 in increments of 0.2 (six variants). The factor psoil represents the mixture ratio of the two representative soil spectra (i.e., a wet soil spectrum and a dry soil spectrum) implemented in the RT model. The canopy reflectance spectra obtained using PROSAIL were linearly mixed with the soil spectra using the FVC, ω, as the weight, which was varied from 0.0 to 1.0 in increments of 0.2 (six variants).
For the leaf angle distribution (LAD), the spherical distribution model was assumed. The other input parameters for PROSAIL were also fixed throughout the study, as summarized in Table 1. Note that a single mesophyll size was assumed throughout the simulation. Furthermore, note that the hotspot size parameter was assumed to be 0.01 to avoid the strong contribution of this effect during the simulation. The total number of parameter combinations was therefore 216 (6 × 6 × 6). The isoline equations for these conditions were simulated for each pair of wavelengths, λ 1 and λ 2 , where λ 2 was varied from 410 to 1200 nm in 10 nm intervals (total of 80 discrete bands) and λ 1 was varied from 400 to λ 2 minus 10 nm at 10 nm intervals, resulting in a total of 80 C 2 (=3160) combinations. In summary, for each pair of wavelength (band) combinations (3160 pairs), a total of 216 TOC reflectance spectra were obtained by the numerical simulations.  The algorithm used to determine the parameters in the vegetation isoline equations is identical to that used in numerous previous studies [10][11][12][13][14]. The canopy reflectance parameter ρ vλ was computed using the spectrally flat zero reflectance of the soil surface. The two-way transmittances for vegetation canopy, T 2 vλ , were approximated on the basis of the simulated TOC reflectance with a soil reflectance of medium brightness, where ρ vλ was computed in the previous step [14]. In the simulation, the TOC reflectances in the λ 2 were approximated from the first-and the second-order interaction terms between the canopy layer and the soil surface: where R sλ 2 represents the bi-hemispherical reflectance of the soil surface for the λ 2 band. The spherical albedo of the canopy layer, R vλ 2 , was obtained by solving Equation (13) for R vλ 2 , in which the soil spectrum was brighter than that in the simulation used to compute T 2 λ . The soil spectrum was also spectrally flat in this case [13]. The slope and offset of the soil line equation over the λ 1 and λ 2 reflectance space were obtained by linear regression of the reflectance spectra for the wet and dry soils ( Table 1). The TOC reflectance spectra used to depict true isolines were simply the output of PROSAIL under the assumption of zero-brightness soil.

Definition of Errors in Isoline Equations and Determination of the Optimum Value of k opt
Errors of the vegetation isoline equations were measured using the distance between the true and the approximated isolines ( (k)): whereρ(k) represents the spectrum derived by asymmetric-order vegetation isoline equations as a function of k and ρ denotes the spectrum on the vegetation isoline including all the higher-order interaction terms (i.e., the true spectrum). Note that the first-order and the asymmetric-order isoline equations correspond to the cases of k = 0 and k = 1 in Equation (11); thus, the errors in the three isoline equations can be simply measured using Equation (14). We again note that the RT simulations were conducted with 216 different parameter settings, which resulted in 216 different values of (k) for each pair of λ 1 and λ 2 . For the sake of total evaluation, these (k) values were averaged over the 216 simulation cases to obtain a mean error represented by (k). In the following discussion, we use the term "mean error" and "mean epsilon, (k)", interchangeably because, as previously described, Equation (14) (i.e., (k opt )) can be applicable to all the cases when k is set appropriately. The optimum value of k in Equation (11) for each wavelength pair was explored on the basis of (k). First, the values of k were computed from Equation (12) for all 216 parameter combinations. The maximum and minimum were then found among the 216 different k values. Finally, within the range of the maximum and minimum values, the optimum value of k, referred to as k opt , was searched, which resulted in a minimum of the mean (k) (i.e., the minimum of (k)) [14,15]. Apparently, this minimum of (k) is identical to the mean of (k opt ) (i.e., (k opt )), which is defined as the mean error (or mean epsilon) of the optimized asymmetric-order vegetation isoline. Figure 2 depicts the typical variation patterns of k opt as a function of λ 1 and λ 2 . The value of k opt is greater than 1.0 for the combinations visible as λ 1 and NIR as λ 2 . Specifically, these wavelength combinations are located in the upper-left region of Figure 2; λ 1 and λ 2 range from 400 nm to 710 nm and from 720 nm to 1200 nm, respectively. On the contrary, the value of k opt is negative, with a relatively larger magnitude within a limited region indicated in blue at the lower left. This small region of wavelengths represents the combinations of green as λ 1 and red as λ 2 . For these two distinct regions of the subspace, the factor k strongly affects the error reduction. Meanwhile, for the upper-right region and near the lower-left corner of the subspace, the value of k opt is scattered around zero. This result implies that the effect of factor k on the error reduction is relatively less significant than in the two aforementioned regions. These trends shown in Figure 2 are common to all the cases considered in the present study. To observe the details of the variations, we plot the cross-section of the figure. Figure 3a,b show k opt as a function of λ 1 for four different values of λ 2 (810, 860, 910, and 940 nm). In Figure 3a, the variations of k opt are shown within the λ 1 range from 400 to 930 nm. To further investigate the details of the differences in the visible region, we plot the results for λ 1 shorter than 700 nm in Figure 3b. The values of k opt vary from 1.2 to 1.4 when λ 1 is shorter than 700 nm; it then decreases sharply at 700 nm and approaches zero as λ 1 increases. As the next step of the analysis, we referred to a band configuration of an existing satellite sensor. For this purpose, we focused on a new-generation geostationary satellite, Himawari-8 [46]. We applied the band configuration of the Advanced Himawari Imager (AHI) to investigate the behavior of k opt in greater detail to meet the needs of this emerging field of study [42,[47][48][49][50][51][52]. In general, the variation patterns of k opt as a function of λ 2 show a local minimum and/or maximum in the four cases. Regarding the range, the value of k opt falls between −0.5 and 1.4 when λ 1 is in the visible wavelength region. On the contrary, k opt varies in a relatively smaller range (from 0.0 to 0.35) when λ 1 corresponds to the NIR band ( Figure 4). A local maximum and minimum were found in the results in Figure 4 for the case of λ 1 = 470 nm. The local maximum is 0.92 at a wavelength near 550 nm, and the local minimum is 0.36 near 670 nm. Similarly, the local maximum and minimum for the case of λ 1 = 510 nm were found at approximately the same wavelengths as in the case of λ 1 = 470 nm. Interestingly, no local maximum was observed in the case of λ 1 = 640 nm. Although both a local maximum and minimum were observed in the case of λ 1 = 860 nm, their locations differ substantially from those in the other three cases. From these results, special caution is needed when the wavelengths (or bands) from both NIR regions are combined. Several aspects of the variation of k opt can be explained from the spectral shape of the TOC reflectance. One example is the value of k opt when the TOC reflectance of λ 2 is greater than that of λ 1 . Such a case occurs when λ 1 and λ 2 are chosen from the visible and NIR bands, respectively. For example, Figure 5a,b show the TOC reflectance spectra for various LAI values with the bright and the dark soils, respectively. In the figure, the reflectance spectra for the visible bands are much smaller than those for the NIR bands, as expected. In this case, the value of k opt becomes greater than 1.0. On the contrary, k opt is expected to be smaller than 1.0 or even a negative value when the TOC reflectance of λ 2 is smaller than that of λ 1 . Furthermore, the value of k opt tends to be nearly zero when the TOC reflectances of the two wavelengths are chosen from a similar wavelength region. When we chose the two wavelengths from a similar region, the vegetation isoline tended to be linear. Because the higher-order interaction term does not need to be included to obtain a linear relationship, the value of k opt results in a small-magnitude value to minimize the correction term of the isoline equation. Note also that the fundamentals of the error reduction mechanism of the asymmetric-order vegetation isoline equation were described in [13].

Accuracy of the Three Vegetation Isolines
In this subsection, we focus on evaluating the accuracy of the three vegetation isoline equations. Recall that the three forms of isolines are denoted by the first-order, the asymmetric-order, and the optimized asymmetric-order vegetation isoline equations. Comparisons of the accuracy were based on the arithmetic mean of (k) computed from the 216 simulations corresponding to different conditions for each wavelength pair. Figure 6 shows the mean epsilon as a function of λ 1 and λ 2 for the first-order, the asymmetric-order, and the optimized asymmetric-order isolines (Figure 6a-c). The mean epsilon of the asymmetric-order vegetation isoline is smaller than the first-order vegetation isoline when λ 1 is chosen from the range between 400 and 710 nm and λ 2 is chosen from the range between 710 and 1200 nm. By contrast, the mean epsilon of the asymmetric-order vegetation isoline is greater than that of the first-order vegetation isoline in the case where both λ 1 and λ 2 are chosen from the NIR range from 720 to 1200 nm. The mean epsilon of the asymmetric-order isoline is similar to that of the first-order isoline when both λ 1 and λ 2 are chosen from the visible range (400-700 nm). The optimized asymmetric-order vegetation isoline resulted in a much smaller value of the mean epsilon, which is mostly less than 0.001 for all the regions. The results show the superiority of the optimized asymmetric-order vegetation isoline over the other two isolines.  (Figure 7d). The mean epsilon of the first-order isoline is mostly greater than those of the other two isolines when λ 1 is less than ∼700 nm for all four cases of λ 2 . The error in the first-order vegetation isoline, however, decreases with the decreasing difference between λ 1 and λ 2 . The mean epsilon of the first-order isoline becomes smaller than that of the asymmetric-order isoline when λ 1 is greater than ∼750 nm. The mean epsilon of the asymmetric-order vegetation isoline becomes (2.0-3.0) × 10 −4 for λ 1 < 700 nm, but the mean epsilon increases rapidly for λ 1 > 700 nm. This result indicates that the correction term introduced in the asymmetric-order isoline makes the error even larger than that of the simplest isoline equation. This error deterioration in the asymmetric-order isoline is improved by introducing the factor k. The optimized asymmetric-order vegetation isoline results in the smallest values of the mean epsilon for all the investigated wavelength pair combinations. These results indicate that the accuracy of the three types of isolines heavily depends on the selected wavelength combination, which implies the importance of the optimization factor k. Importantly, however, the value of k should be set appropriately for each wavelength pair to reduce the error.   (Figure 8d) in the four cases. The mean epsilon of all three isolines is close to zero in the λ 2 range less than 700 nm, except for the result shown in Figure 8d, which is the case of λ 1 = 860 nm. By contrast, the mean epsilon shows increasing trends when λ 2 is greater than 700 nm. Overall, the mean epsilon of the first-order isoline is the greatest among the mean epsilons of the three isolines. The mean epsilon of the asymmetric-order isoline equation, however, becomes the largest in the case of λ 1 = 860 nm (Figure 8d). This result also indicates the need for special caution when choosing the wavelength pair. The results in Figures 7 and 8 imply that the first-order isoline equation could be more accurate than the asymmetric-order isoline equation in the case of some specific wavelength pairs. The reason for this somewhat unexpected result is carefully investigated in the next subsection.  Figure 9 shows the true vegetation isoline and the spectra predicted by the three vegetation isoline equations plotted in the reflectance subspace. The wavelength for λ 2 was fixed at 860 nm, whereas λ 1 was varied among four discrete values: (a) 690 nm, (b) 710 nm, (c) 730 nm, and (d) 850 nm, covering from the red edge to the NIR region. Note that the last selected λ 1 in Figure 9d is close to λ 2 ; the difference between the two wavelengths is only 10 nm. This choice was made on the basis of the results in Figure 7b, which shows an interesting behavior of the first-order and the asymmetric-order isolines. The values for the LAI and FVC were fixed during the simulation at 1.6 and 1.0, respectively.

Influence of the Higher-Order Interaction Terms Demonstrated in the Reflectance Subspace
In Figure 9, the first-order, the asymmetric-order, and the optimized asymmetric-order isolines are indicated by blue circles, a red dashed line, and green crosses, respectively. In Figure 9a, the asymmetric-order isoline (red dashed line) is closer to the true isoline (black solid line), clearly indicating that the accuracy of the asymmetric-order vegetation isoline is greater than that of the first-order vegetation isoline. The asymmetric-order isoline in Figure 9b is even closer to the true one compared with that in Figure 9a. Interestingly, the superiority of the asymmetric-order isoline over the first-order isoline is inconsistent. For example, in Figure 9c, the spectra approximated by the asymmetric-order isoline (red dashed line) exceed the true isoline (black solid line), and the distance between the true isoline and the asymmetric-order isoline is almost identical to the distance between the true isoline and the first-order isoline (blue circles). The distance between the true isoline and the asymmetric-order isoline becomes even larger in Figure 9d. The first-order isoline obviously accurately describes the true isoline in the case of Figure 9d, indicating that the influence of the last term of the asymmetric-order isoline adversely affects the accuracy in this case. This problem in the asymmetric-order isoline has been dramatically improved by the introduction of the optimization factor k (the optimized asymmetric-order isoline represented by the green crosses) in all the investigated cases. To confirm the influence of the higher-order term and its optimization performance represented by the last term in Equation (11), we conducted further error comparisons among the three isolines. Recall that Equation (11) is identical to the first-order vegetation isoline when k = 0. When k = 1 and k = k opt , Equation (11) results in the asymmetricorder and the optimized asymmetric-order isoline, respectively. Thus, the expression of Equation (11) covers the three types of isoline equations simply because of the selected k factor. We refer to the last term as the "over-correction term" because it is intentionally introduced to overly increase ρ λ 2 to compensate for the higher-order truncations. Figure 10a-d show this overcorrection term computed from the last term for the asymmetric-order isoline, k = 1 (red line), and the optimized higher-order isoline, k = k opt (green crosses). Furthermore, the error of the first-order isoline, which corresponds to the case of k = 0, is shown by the blue line. The combination of the wavelengths, λ 1 and λ 2 , is the same as in the case shown in Figure 9. In the cases of Figure 10a,b, the two wavelengths differ to a certain extent; thus, the asymmetric-order isoline shows high accuracy. As evident in the figures, both the green crosses (the optimized asymmetric-order isoline) and the red line (the asymmetric-order isoline) are close to the blue line (the error of the first-order isoline), as we expected. However, when the two wavelengths are similar, the first-order isoline results in better accuracy significantly than the asymmetric-order isoline. This case corresponds to Figure 10c,d, in which the influence of the correction term of the optimized asymmetric-order isoline (green crosses) is much closer to the error of the first-order isoline (blue line) than the correction term of the asymmetric-order isoline (red line). This result implies that the optimized asymmetric-order isoline includes an appropriate value for k and that, because of this appropriate choice of k, the overcorrection term becomes close to the error of the first-order isoline (the truncation error of the isoline) based on the wavelength pair. In summary, the optimized higher-order term nicely minimized the truncation error of the isoline for all combinations of wavelengths, retaining high accuracy irrespective of the wavelength combination. Figure 10. Plot of the error in the first-order vegetation isoline equation (blue line) and the overcorrection terms of the asymmetric-order vegetation isoline equation (red line) and the adjusted overcorrection term with k opt (green crosses). The wavelength pair of λ 1 and λ 2 is the same as in Figure 9. The wavelength λ 2 was fixed at 860 nm, whereas λ 1 was selected from (a) 690 nm, (b) 710 nm, (c) 730 nm, and (d) 850 nm (λ 2 minus 10).

Discussion
The original motivation for developing the isoline formula (the first-order isoline equation) was to analyze the behaviors of the vegetation isoline itself and the spectral VI [10,11]. The derived first-order isoline equation is the simplest among the three isoline equations; hence, it still has an advantage in various analyses. For example, several studies have demonstrated a method that uses the isoline concept for retrieving biophysical parameters such as the FVC and LAI [19,27,[53][54][55]. If the first-order isoline equation is used to retrieve these biophysical parameters, the errors caused by the truncation of the higher-order interaction terms will be a major source of error in the retrieved results. For instance, the errors of the first-order isoline are not negligible when it is used for high-accuracy parameter retrievals. This problem should be addressed by retaining its simplicity in the isoline formulation, which was the major challenge of the present research. The initial attempt was to include the higher-order interaction term in a unique manner [13]. The new isoline formula (i.e., the asymmetric-order isoline) can achieve better accuracy than the firstorder isoline. However, because the formulation was somewhat unusual (asymmetrical inclusions of the higher-order term), we encountered unexpected behavior in the present study. In particular, the behavior of the errors caused by the choice of the wavelength pair needed to be addressed thoroughly.
In our previous work, we discussed the accuracy of the vegetation isolines for pairs of red and NIR wavelengths by referring to the bands of some existing sensors [14]. For example, when a wavelength pair of 655 and 865 nm is chosen for the red and NIR bands, respectively, we found that the best accuracy was achieved using the optimized asymmetricorder isoline among the three isoline formulas. This study shows that the findings of the previous study remain the same for all the investigated wavelength combinations; the wavelength pair is not limited to the combination of red and NIR wavelengths to support the findings in the previous work [14]. The advantage of the optimized asymmetric-order isoline can also be seen by considering the specifications of the existing hyperspectral sensors. For example, the signal-to-noise ratio (SNR) of EnMAP is 400 in the optical region [38]. The corresponding value of Hyperion is even lower: 144 at 650 nm [56]. Since the reflectance level assumed for the SNR is 0.3 [38], the noise equivalent reflectance becomes 7.5 × 10 −4 for EnMAP and 2.1 × 10 −3 for Hyperion. From Figure 7a-c, the mean epsilon of the optimized asymmetric-order isoline (black line) is well below the noise equivalent reflectance of those sensors.
On the contrary, the magnitude of the errors in the asymmetric-order isoline is not always smaller than the errors in the first-order isoline, which is a somewhat unexpected result from our previous work [14]. For example, in cases such as those shown in Figure 9c,d, the asymmetric-order isoline showed larger errors than the first-order isoline despite the higher-order interaction term being considered during the derivation. This result implies that special caution is needed to include the higher-order interaction term asymmetrically. Interestingly, the optimization factor k can overcome this difficulty in a wide range of wavelength pairs when the optimum value for k is determined algorithmically. From these results, by using the optimized asymmetric-order isoline, the wavelength region can be extended without exacerbating the errors.
The advantages of the optimized asymmetric-order isoline are improvements in both the accuracy and the flexibility of the wavelength choice. These advantages will lead to better accuracy in parameter retrievals if the first-order isoline is replaced by the optimized asymmetric-order isoline. Suppose that one attempts to retrieve a biophysical parameter using the first-order vegetation isoline equation. Because the error of a retrieved biophysical parameter is expected to be proportional to the error of the isoline equation, an error of nearly 8% [10,11] relative to the retrieved value would be observed when the soil of the canopy background is bright. (This error is due to the relatively larger contribution of the higher-order interaction terms.) This error could be reduced to 1/25 by replacing the first-order isoline by the optimized asymmetric-order isoline. Although the magnitude of the error reduction depends on the wavelength pair, replacement of the isoline may provide a great advantage in various fields of application in which the isoline concept plays an important role [18,33,[57][58][59]. This study investigated the errors in the isoline equations within a limited set of parameter settings. This limitation should be noted to consider the results of this study.

Conclusions
In our previous work, accuracy assessments were limited to a common pair of red and NIR wavelengths. This limitation was problematic when the isoline concept was used in a broader wavelength region. This study addressed this problem by extending the wavelength range from 400 to 1200 nm. We conducted a set of numerical simulations to characterize the behavior of the three isoline equations. Specifically, the investigation focused on the optimization factor k of the optimized asymmetric-order vegetation isoline.
The factor k introduced in the asymmetric-order isoline equation was simulated under various conditions for each wavelength pair, and its optimum value, k opt , was determined algorithmically for each wavelength pair. The value of k opt was found to strongly depend on the choice of the wavelength pair. In general, three trends were observed in the relationships between the k opt and the wavelength choice. The trends were identified on the basis of the relative differences between the two reflectances. When the reflectance of λ 2 was greater than that of λ 1 , the value of k opt tended to be greater than 1. The value of k opt was mostly between 0 and 1 when the reflectances of λ 1 and λ 2 were similar. Finally, the value of k opt tended to be less than 0 when the reflectance of λ 2 was less than that of λ 1 .
With respect to the accuracy of the three variants of the vegetation isoline, the optimized asymmetric-order isoline equation was the most accurate for all the wavelength pairs investigated in the present study. The asymmetric-order isoline equation, in general, showed better accuracy than the first-order isoline equation; however, this trend did not hold for some wavelength pairs. Specifically, when the reflectances of the two wavelengths were similar, the accuracy of the asymmetric-order isoline equation was worse than that of the first-order isoline equation. This result could have been caused by the effects of the asymmetrical inclusion of the higher-order interaction. This negative influence of the asymmetric term was clarified in the present study. We expect the developed isoline formula to be useful in a wide range of applications involving analytical and numerical tools. The application of the optimized asymmetric-order isoline equation to actual data processing of existing hyperspectral sensors remains undone, which should be considered as a future effort. Further efforts will be needed to facilitate the isoline formula as a practical tool for processing satellite data.