A Pixel Dichotomy Coupled Linear Kernel-Driven Model for Estimating Fractional Vegetation Cover in Arid Areas From High-Spatial-Resolution Images

With the increased use of high-spatial-resolution (HSR) images for vegetation monitoring in arid areas, more details of the low vegetation coverage and interference from the land “background” are captured in the corresponding images. From computational time and accuracy, the multiangle method (MAM) in the pixel dichotomy model is a potential algorithm to apply in arid areas, but MAM needs the multiangle vegetation index (VI) as the driver parameters. However, most HSR images are obtained in nadir mode, and the multiangle information of reflectance is difficult to obtain, which limits the estimation of multiangle VI from HSR images. To address this issue, this study used a “graphical method” to modify the radiation influence caused by the canopy structure and land “background.” We developed an inversion method of the linear kernel-driven model (KDM) and designed a random sampling method to estimate multiangle VI from HSR images. Then, we proposed a new pixel dichotomy coupled linear KDM (PDKDM), validated using simulated, field-measured, and reference data. The results showed that the FVC in arid areas estimated by PDKDM was highly consistent with “true” data, with root-mean-square error (RMSE) < 0.062, RMSE < 1.125, and RMSE < 0.027 for comparison with simulated, field-measured, and reference data, respectively. PDKDM addressed the issue with the previous MAMs to estimate FVC from HSR images in arid areas. This study provides a useful algorithm with high computational efficiency for producing HSR FVCs in arid areas.

One of adjustment factor. L One of adjustment factor. r RED Black-sky albedo in red band. r NIR Black-sky albedo in near-infrared (NIR) band. θ i Solar zenith angle. θ l Leaf inclination angle. P(θ i ) Gap fraction. E Clumping index. G(θ i ) G function, it denotes the projection function of leaves area in the θ direction. L AI Leaf area index (LAI). r (θi) Black-sky albedo. f iso Parameters of isotropic in the linear KDM. f vol Parameters of RossThick in the linear KDM. f geo Parameters of Lisparse-R in the linear KDM. g * polynomial coefficients in the linear KDM, including g 0iso , g 1iso , g 2iso , g 0vol , g 1vol , g 2vol , g 0geo , g 1geol , and g 2geo .

I. INTRODUCTION
F RACTIONAL vegetation cover (FVC) is a canopy architecture parameter representing the fraction of land surface covered by green foliage in the 2-D plane [1]. In arid areas, FVC is closely related to water distribution [2], [3], [4]. Therefore, FVC monitors water resource distribution, describes desertification, and helps with desertification monitoring in arid areas. The prevalent algorithms include: 1) an empirical method [5], [6]; 2) a physical model-based method [7], [8], [9], [10]; and 3) a spectral mixing analysis (SMA) [7], [8], [11], [12], [13], [14] to estimate FVC from remote sensing images. Considering the integrated efficiency of computational time and accuracy, SMA is the optimal method to invert FVC at a large regional scale. SMA is a semiempirical model which solves a spectrum linear equation to analyze the probability relationship of signals between endmembers (the decomposed "basic component" of a pixel) [14], [15]. However, SMA also involves a complex matrix operation and an analysis of the matrix-vector space [12], [16], which consumes an enormous time. Therefore, to further accelerate the computational speed, the endmember reflectance in SMA has been replaced by the normalized difference vegetation index (NDVI) [17], and a PDM extending from SMA has been proposed [18], [19], [20], [21]. The PDM ignored the clumping index and used an approximate equation to estimate FVC. In PDM, the commonly used model is the model proposed by Gutman and Ignatov [21], i.e., linear PDM, given by FVC = VI − VI min VI max − VI min (1) where FVC denotes the green FVC, with values ranging from 0 to 1; VI min and VI max denote the maximum and minimum NDVI, respectively; and VI denotes the rescaled NDVI value. Equation (1) is a general form of PDM, and VI min and VI max are key parameters to calculate FVC. Currently, there are three methods to determine VI min and VI max : 1) the observation data method (ODM); 2) the direct analysis method from images (DAM); and 3) the MAM. ODM is performed via ground measurement or high-resolution products [6], [22], with high computational accuracy. However, this method requires numerous samples to ensure the accuracy of VI min and VI max . Therefore, it is difficult to promote them in a large-scale application. The proposed DAM addresses this issue [21], [23]. This method has two techniques to analyze VI min and VI max from an image directly: one is to analyze the VI min and VI max by clustering analysis [21] and the other is to analyze changes in the time series of VI combined with the classification frame of land cover to determine VI min and VI max [24]. VI max is the NDVI of a high-density vegetation area in a scene, and hence, can be easily obtained [23], [25]. VI min is the NDVI of land "background" (i.e., soil in the land, sand in the desert, and gravel in the Gobi), which is affected by the component, particle size, and moisture [26] and has large uncertainty. Thus, the determination of VI min is difficult [26], [27]. Therefore, an MAM based on a fixed angle method in photography was developed [28], [29].
In the fixed angle method, when the zenith angle is 57.5 • , the projection function of the area of leaves within the solid angle (i.e., the G function in the radiative transfer of vegetation) is equal to 0.5, and Beer's law-related FVC (A-2) can be solved [28], [29], [30]. In general, MAM uses a zenith angle of 57.5 • to calculate the multiangle VI and then solves Beer's law-related FVC to determine VI min and VI max . At an early stage, an MAM only uses a zenith angle of 57.5 • . Thus, some parameter products involving Beer's law-related FVC need to be considered to solve the Beer's law-related FVC, e.g., LAI and clumping index ( E ) products [31], [32]. A proposed multi-VI algorithm reduces MAM's dependence on the products. In this algorithm, two viewing zenith angles approximate the viewing zenith angle of 57.5 • to establish a multiangle calculation equation to calculate multiangle VI [33]. The multiangle VI is calculated from a medium-coarse spatial resolution product, i.e., the products of bidirectional reflectance distribution function (BRDF)/albedo model parameters (MCD43A1). Therefore, the multiangle VI images are also medium-coarse spatial resolution images. MCD43A1 is a product for calculating black-sky albedo, and its physical properties are determined by solar direction. However, the multi-VI algorithm ignores the difference between the reflectance factors (the black-sky albedo is one of the reflectance factors) on the properties of light (specular and diffuse) and the geometric relationships (incidence and outward) and modeling from the viewing angle. This phenomenon causes an angle singularity between the algorithm and the black-sky albedo calculated from MCD43A1. In addition, the sensitivity of the zenith angle of NDVI between 55 • and 60 • in the forward and backward directions also is ignored in the multi-VI algorithm. Afterward, a multi-HSVI considering the coarse-spatial resolution product of the directional reflectance factor (BRF) of "hotspot" was proposed [34]; however, the design goal of this model is still based on viewing the zenith angle, but the data of physical properties at the "hotspot" are related to solar and viewing zenith angle, which avoids the angle singularity. With the popularity of HSR images in vegetation monitoring, a challenge will arise for MAM. Since most HSR images are obtained in nadir mode, the multiangle reflectance is difficult to obtain, which limits the estimation of multiangle VI from HSR images. To the best of our knowledge, the most critical problem about the absence of an estimation method to obtain the multiangle VI from HSR images has not been addressed.
The previous studies on the PDM were mainly concerned with medium-coarse spatial resolution images and rarely involved estimation of FVC from HSR images. Therefore, the deviation of FVC estimated by the previous multi-VI algorithm was not obvious when the radiation change caused by the land "background" and the canopy structure (i.e., canopy height, leaf inclination angle, and LAI) was ignored [27], [30], [35], [36], [37], [38]. Different from the soil "background" in moist areas, the land "background" in arid areas is distributed with a large area of desert, Gobi, and saline-alkali land, and there is a phenomenon of biological crusts [6]. Moreover, numerous desert shrubs are distributed in interlaced areas between oasis and desert, most of their leaves are needles, and their woody areas are higher than the leaf areas. This special canopy structure causes difficulty in estimating FVC. With the remarkable development in remote sensing with HSR images, radiation information on land "background" and canopy structure are easy to capture. In previous studies [22], [31], [33], [35], the linear PDM ignored these influencing factors to accelerate the computational speed. This process causes mathematical physics defects in the linear PDM, limiting its computational accuracy in estimating FVC from the HSR images, especially for estimating FVC in arid areas. Therefore, it is necessary to modify the linear PDM and further improve the PDM.
To address the two core limitations of the PDM on HSR images mentioned earlier, we proposed a PDKDM to extend the PDM to estimate FVC in arid areas from HSR images. In particular, we modified the linear PDM considering the radiation influence caused by the canopy structure and land "background." We also rederived the multiangle calculation equation. An inversion method of the linear KDM and a random sampling method were introduced to solve the equation. We successfully extended the PDM to estimate FVC from HSR images in arid areas with these two steps.

II. METHODOLOGY
In this study, PDKDM was developed to extend the MAM to estimate FVC in the applications of HSR images. In establishing PDKDM, two steps are considered. First, we modified the linear PDM by considering the radiation influence between the canopy structure and the land "background." Second, we rederived the multiangle calculation equation, and the corresponding solution method was proposed. The implementation flowchart of PDKDM is shown in Fig. 1.

A. Modified Linear PDM
Because the linear PDM ignores the radiation change caused by the land "background" and canopy structure, it is difficult to use the PDM to estimate FVC in the interlaced areas between oasis and desert from HSR images. In this study, a "graphical method" is used to analyze the deformation of FVC images caused by the land "background" and canopy structure (Fig. 12). We introduced the Lagrangian polynomial of three points to fit the nonlinear function of Beer's law-related FVC and modified the mathematical defects of the linear equations involved in the linear PDM. According to the derivation in Appendix A, we developed a modified linear PDM, as shown in the following: where l and L are adjustment factors; they physically represent comprehensive radiation information (multiple-scattering) influenced by the land "background" and canopy structure. Some VIs have adjustment factors, e.g., the soil-adjusted VI (SAVI) and the perpendicular VI (PVI) [39], [40], which would bring some uncertainties in solving (2). In addition, some VIs are not normalized, e.g., the ratio VI (RVI), which would result in outliers (FVC appearing greater than 1 or less than 0) in (2). Therefore, we use the most common index, NDVI, in this study where r NIR and r RED denote the black-sky albedo in the NIR and red bands, respectively. In (2), VI(0 • ) can be obtained directly from an image. Therefore, VI min , VI max , L, and l are key parameters to calculate FVC.

B. Multiangle Calculation Equation
To solve the modified linear PDM, we need to determine the input parameters (VI min , VI max , L, and l). The Beer's lawrelated FVC (A-3) is substituted in (2) to address this issue in the modified linear PDM We introduced a conclusion in the fixed angle method in photography [28], [29], i.e., when θ i = 57.5 • and G(θ i ) = 0.5. Here, θ i is the solar zenith angle. In a previous study [33], the multi-VI used 55 • and 60 • to approximate 57.5 • , and G(57.5 • ) ≈ G(55 • ) ≈ G(60 • ). Therefore, (4) was derived as follows: Because of G(57. where VI(55 • ) and VI(60 • ) denote the images of NDVI when the solar zenith angle is 55 • and 60 • , respectively. Equation (7) is based on the solar zenith angle to establish a multiangle calculation equation and tries to address an angle singularity between the multi-VI algorithm and black-sky albedo calculated from MCD43A1 and expands it to estimate FVC from HSR images.  (3) and (B-1) are combined, and the HSR images of VI(55 • ) and VI(60 • ) can be obtained. When HSR images of VI(55 • ) and VI(60 • ) are calculated, the next step is to calculate the four unknown numbers in (7): VI min , VI max , l, and L. These four unknown numbers are input parameters in the modified linear PDM (2). Therefore, obtaining the values of these unknowns is also a key to solving (2). Considering the amount of data of HSR image is large, we proposed a random sampling method, which used a chessboard block to split an image into subimages, used the Monte Carlo method to randomly obtain n pixels in the subimages to obtain reflectance, and established a multiangle calculation equation [see (B-7)], then find an optimal solution using (B-8)-(B-10) in Appendix B2, i.e., VI max , VI min , L, and l. Finally, FVC can be calculated when four unknown numbers are solved.

III. VALIDATION DATA
To validate PDKDM, simulated, field-measured, and reference data were used. We chose the area of Ebinur Lake [ Fig. 2(b)] as a case study on the interlaced area between the oasis and desert to validate the FVC estimated by PDKDM coupled IA in arid areas, showing the computational accuracy of the canopy of desert shrub in the desert area. In addition, to show the advantages of PDKDM for estimating FVC from

A. Simulated Data
We chose the 4SAIL model (the third generation of the scattering by arbitrary inclined leaves (SAIL) model considers "hotspot" and equation singularities) [41], [42], [43] to generate the datasets comprising FVC, VI(0 • ), VI(55 • ), VI(60 • ), and input parameters. The 4SAIL model is the radiative transfer model with the highest accuracy, which estimates multiangle black-sky albedo and FVC by substituting input parameters and further calculates multiangle VI and FVC. The input parameters in the datasets include radiation, biophysical, and angle parameters. The leaf hemisphere reflectance (ρ), leaf hemisphere transmittance (τ ), and background reflectance (r b ) are radiation parameters, we determined using various measured data of vegetation and land "background" as the standard spectrum in arid areas (Table I). LAI is a biophysical parameter set from 0.01 to 5.01 to show the canopy in different leaf densities. Leaf inclined angle (θ l ) and viewing zenith angle (θ o ) are the angle parameters that determine the interception area of radiation by leaves; Table I also shows their settings. Based on the settings, we obtained 159 600 VI-FVC datasets from the 4SAIL model. Because the VI-FVC pairs in the VI-FVC datasets are relatively large, they are divided into groups to estimate FVC. Each group of VI-FVC pairs is generated by randomly selecting 760 samples within the various ranges of ρ, τ , r b , and θ o , i.e., (4 × 2 × 5 × 19 = 760), which are the corresponding numbers of ρ, τ , r b , and θ o (Table I). According to the arrangement between LAI and θ l , we produced 210 groups of datasets (159 600/760 = 210). In the 210 groups, VI(0 • ), VI(55 • ), and VI(60 • ) can be used as input parameters to establish (B-7) to calculate FVC using the modified linear PDM in PDKDM. The FVC estimated by the 4SAIL model can be used as a "surrogate truth" to validate the FVC estimated by PDKDM. The implementation flowchart was shown in Fig. 3.

B. Field-Measured Data
To validate the computational accuracy of FVC in the interlaced area between the oasis and desert, we designed 69 plots with 10 × 10 m 2 . Based on the resolution of the remote sensing image, i.e., Sentinel-2 (Level-2A) image with a spatial resolution of 10-m [ Fig. 2(b)], five quadrats (approx-imately 2 × 2 m 2 ) were selected in each plot to obtain the average FVC. The FVC measurements were performed from July, 15 to 21, 2021 (note: we ignored the vegetation growth within several days). In the FVC measurement, the SONY-SLT-A37 digital camera (CCD) with a field of view of 46 • was used. We designed two types of measurement methods based on the canopy distribution density. Because the plots are located in arid areas, most vegetation is sparsely distributed. The camera was mounted on an unmanned aerial vehicle (UAV) to capture the sparse canopy. FVC images that covered a large area were acquired, and an average FVC with representativeness was obtained. For the thick canopy, the camera was mounted on a telescopic frame, and nadir images were shot; the measuring height was determined according to the quadrat size (the calculation method is described in [44]). For a forest higher than the telescopic rod, we determined FVC using images captured by shooting upward and downward (refer to (8) in [45]). In processing images captured by digital cover photography (nonfisheye lens), we used the laboratory's spatial conversion technology to convert the captured color images into binary images. In the binary images, a pixel value of 0 represented the background (soil, saline, desert, or Gobi), and a pixel value of 255 represented an area covered by vegetation (Fig. 4). We counted the number of pixels in the vegetation components to obtain FVC in the images in the following: where N v and N s denote the number of vegetation and background pixels, respectively, fvc up and fvc down are FVCs extracted from images captured in the bottom-up and top-down directions.

C. Reference Data
PDKDM is one of the MAMs in the PDM, which is suitable for estimating FVC from HSR images. To further compare the difference between PDKDM and the existing MAM, i.e., multi-VI [33], we chose 30 plots gathered in the subimages in Ebinur Lake, i.e., the oasis and desert-edge areas in Fig. 2(c) and (d). Because the computational accuracy of the physical model-based method is very high, we chose a physical model-based method, PROSAIL coupling backward propagation neural networks (BPNNs) [46]. The LAI estimated from PROSAIL coupling BPNNs was used as a "surrogate truth" to validate the incremental relationship of the exponential function between FVC and LAI in PDKDM. For the dataset involved in PROSAIL coupling BPNNs, see [7], [46].  ]. In addition, the LAI and FVC estimated by 4SAIL have an incremental relationship to the exponential function, whereas the LAI estimated by 4SAIL and FVC estimated by PDKDM have an approximate exponential form with slight fluctuations (Fig. 5). When LAI increases, FVC exhibits a saturated phenomenon, and the saturated phenomenon reduces as the average leaf inclination angle increases. The statistical results also show similar trends; when the average leaf inclination angle is between 10 • and 70 • , root-mean-square error (RMSE) is lower and less than 0.062 with only less than −7.45% difference    Fig. 7(g)-(h) are greener than Fig. 7(i), and Fig. 7(g) and (h) is very close, which implies that the BRDF effect of NDVI in the solar zenith angle is very small. Fig. 7(j) shows the image of the FVC estimated by PDKDM in the area of Ebinur Lake; there is a high consistency between it and the field measurement (RMSE = 0.124). The FVC of the oasis in summer was above 0.7 [green in Fig. 7(j)]. The areas of Ebinur Lake, the saline-alkali land around the lake, desert areas, and the southern slope of the Tianshan Mountains are indicated by white in Fig. 7(j), and the FVC is close to 0. In Fig. 8, the points of desert shrubs (including salt claws, Tamarix ramosissima, shuttle, and other desert shrubs) are closer to the 1:1 line than the point of other vegetation (including farmland, deciduous broadleaf, and aquatic plant). Because of the lack of HSR products to validate the FVC and parameters of the KDM, we used Sentinel-2 with upscaling to compare commonly non-HSR products in Xinjiang, which further demonstrates PDKDM's ability to eliminate the dependence on remote sensing products. The results are shown in the Supplemental Material A-2. Fig. 9 shows the FVCs of the subimages in the area of Ebinur Lake estimated by different algorithms from the Sentinel-2 images. There is a high consistency between the FVCs estimated by PDKDM and multi-VI with a 3 × 3 moving window, but multi-VI with a 3 × 3 moving window has plaque phenomena. These plaque phenomena are particularly Fig. 8. Statistics of FVC in the area of Ebinur Lake estimated by PDKDM from Sentinel-2 against measured FVC. FVC_M denotes the measured FVC. Note: salt claws, Tamarix ramosissima, and shuttle are desert shrubs, and another desert shrub is also a desert shrub without a name. obvious in the oasis area [ Fig. 9(c)], implying that the 3 × 3 moving window is unsuitable for the FVC calculation in HSR images, especially calculating maximum VI and minimum VI in the farmland [ Fig. 9(e) and (f)]. In the desert-edge area, the FVC estimated by multi-VI with a 3 × 3 moving window is greater than 0.3 [ Fig. 9(d)], as well as greater than the FVC estimated by PDKDM, showing a slight overestimation [0.028 > 0.027 in RMSE in Fig. 10(b) and (d)]. In the relationship between FVC and LAI, the FVC estimated from all algorithms and LAI shows the exponential function's  relationship (Fig. 11). The relationship of the exponential function between FVC estimated by PDKDM and LAI is the most obvious. Therefore, this further implies that the "graphical method" considered in the modified linear PDM can simulate the nonlinear function of Beer's law-related FVC in Fig. 12 and has great potential in the applications of HSR images.

A. Advantages of Modified Linear PDM
In the modified linear PDM, we introduced a Lagrangian polynomial of three points (Appendix A) to fit the incremental relationship of the exponential function between FVC and VI(LAI) (the curve of the FVC and LAI in Figs. 5 and 11), and derive two parameters [l and L in (2)] in the "graphical method." Different from modified linear PDM, the linear PDM ignored the multiple scattering caused by the canopy structure and land "background," which caused the different curves of  FVC. In this study, the simple parameters were added in the linear PDM, which obtained high computational accuracy in Figs. 10 and 11. These results indicated that l and L simulated nonlinear curves of FVC and improved the computational accuracy. Therefore, the computational deviation caused by the canopy structure was modified in the modified linear PDM. In the simulated data, L was much larger than l (Fig. A-1 in the Supplemental material A-1). L was located at the denominator, and l was located on the numerator, implying that the FVC simulated by the modified linear PDM should be less than that of the linear PDM, further addressing the overestimated issue of FVC in the linear PDM reported in [27].
In addition, Figs. 5 and 11 showed that the modified linear PDM captures the curve of FVC simulated by Beer's law. These results imply that the simple parameters analyzed by the "graphical method" can well represent its physical properties to simulate the function of the change curve of FVC. This "graphical method" used constants instead of complex equations expression, which made iterative computation involved in the solution of the multiangle calculation equation [see (B-7)-(B-10)] possible, e.g., l and L in Table IV and was very suitable for the design of accelerating algorithms.

B. Advantages of Calculation Methods of Multiangle VI
In HSR images, in addition to ignoring the information on canopy structure in the linear PDM, the maximum obstacle of the multi-VI in improving computational accuracy is the lack of effective methods to calculate multiangle VI from HSR images, i.e., VI(55 • ) and VI(60 • ). At present, the nadir observation mode is a common mode for most HSR satellites. For HSR images captured in a nadir observation mode, we developed an inversion method of the linear KDM to calculate kernel parameters from HSR images in this study, i.e., IA [ Fig. 7(a)-(f)]. IA is used, and the multiangle black sky albedo from HSR images can be calculated as possible.
The multiangle black sky albedo can be used to calculate multiangle VI [ Fig. 7(g)-(h)]. This calculation indicates that FVC can be calculated from HSR images when the multiangle calculation equation [see (7)] is used, see Fig. 7 and RMSE = 0.124 in Fig. 8. In the calculation of kernel parameters, how to calculate each pixel of kernel parameters in the 2 × 2 moving window is difficult. Because the four pixels involved in the 2 × 2 moving window can only calculate a pair of kernel parameters, i.e., f iso , f vol , and f geo [ Fig. 15(b)]. Therefore, the iterative method is used to obtain each average value in the 2 × 2 window [ Fig. 15(d)]. This modeling thought considers upscaling and downscaling in the moving window (see Fig. 15).
In this study, the KDM parameters estimated by IA and MCD43A1 had a high consistency (Figs. A-3 and A-4 in the Supplemental material A-2), implying that IA can calculate the KDM parameters independently. PDKDM eliminated excessive dependence on remote sensing products and further addressed the common limitation issues for existing MAMs. Although IA has more advantages, it is difficult to obtain a correct value when the solar zenith angle equals 0 • , and there is no method for this deficiency before. From the model simulated results (Figs. 13 and 14), we can find that the reflectance is very close when the solar zenith angle is equal to 5 • and 0 • , respectively. The zenith angle of 5 • can be used as approximate 0 • . However, this phenomenon is the deficiencies of the linear KDM itself (B-1), which implies that the linear KDM needs further modification in the future.
In comparison between reference data, there were plaque phenomena in multi-VI, which is similar to the issue of the blocks effect in image fusion [47] caused by a 3 × 3 moving window [ Fig. 9(c)]. The primary reason is that the distance weight difference between the eight adjacent pixels involved in the 3 × 3 window is not considered, i.e., the distance between the diagonal pixel and rescaled pixel is greater than the distance between the adjacent pixel and the rescaled pixel [ Fig. 16(a)]. Therefore, we used a random sampling method (Appendix B2) to address these plaque phenomena [ Fig. 9(a)]. Ensuring certain accuracy, the random sampling method speeded up the computational time. The random sampling method as a subalgorithm in PDKDM based on images captured in the nadir observation mode [r (θ i ) in Section B-I] has strong robustness.

C. Advantages of Calculation Methods of Multiangle VI
In this study, PDKDM estimated FVC based on the solar zenith angle, which addressed the angle singularity between the black-sky albedo determined solar direction [48], [49] and the MAM design from viewing direction [31], [33], [36]. In fact, due to the existence of reciprocity in BRDF, the radiation energy between black-sky albedo (directional-hemispherical reflectance factor, DHR) and hemispherical-directional reflectance factor (HDRF) is very close. The VI calculated based on the reflectance factor also has a very similar value, i.e., VI(θ i = 55) ≈ VI(θ o = 55) and VI(θ i = 60) ≈ VI(θ o = 60) [50], [51]. Therefore, the multiangle calculation equation based on the view zenith angle still maintains a good result when FVC is calculated. However, the principle of reciprocity is just established in theory. In the actual measurement, the principle of reciprocity has a deviation, which is influenced by atmospheric scattering, and adjacent pixels [50]. Fig. 5 shows that there are only some fluctuations in some angles with low LAI for FVC simulated by PDKDM, e.g., average leaf inclination angles of 0 • , 80 • , and 90 • [RMSE > 0.062 in Fig. 6(a), (i), and (j)]. When the distribution of leaf inclination angles is planophile (horizontal leaves are most frequent, i.e., average leaf inclination angles are equal to 80 • and 90 • , respectively) and erectopile (vertical leaves are most frequent, i.e., the average leaf inclination angle is equal to 0 • ), the PDKDM is unstable in these cases when the vegetation is sparse. However, contrary to the simulation experiment's validation (Figs. 8 and 10), PDKDM was applied to estimate HSR images, which has excellent computational accuracy (Fig. 11). The primary reasons are that planophile and erectopile are extreme cases, which happens less in nature. The average leaf inclination angle for vegetation in nature is located between 30 • and 60 • [52], and the fluctuation of radiation variables is not large.

D. Potential of PDKDM in Arid Areas
Estimating FVC in areas with sparse vegetation became one of the key issues in vegetation monitoring. We chose the interlaced area between an oasis and a desert in this study (Fig. 7). In a previous study [6], FVC was overestimated at the early growing stage of vegetation or vegetation in arid areas, and some FVC produce in the interlaced between desert and oasis with desert shrub have no value, i.e., GEOV2 FVC. PDKDM considered the radiation influence on canopy structure and land "background," which had high computational accuracy, especially in desert areas [Figs. 8 and 10 and red points in These results also imply that PDKDM is not only suitable for arid plant communities with low FVC but also for other plant communities with high FVC, which has excellent robustness and has the potential to make up for the lack of value for current products in arid areas.
In addition, because we modified the linear PDM using the "graphical method," stable values of FVC [ Fig. 6(b)-(h)] can also be obtained by PDKDM from VI-FVC pairs comprising different land "backgrounds" (r b in Table I). In the interlaced area between the oasis and the desert, it is the main area is covered by desert shrubs. We further investigated the performance of PDKDM in the interlaced area between the oasis and desert, i.e., the area of Ebinur Lake, and found that PDKDM was suitable for estimating the FVC of desert shrubs, e.g., salt claws (the FVC is less than 0.25 points in Fig. 8). This is because PDKDM considers the canopy structure and land "background" and integrates the linear KDM. The linear KDM involves a certain physical mechanism (B-1), and it has strong universality and can be further extended to HSR images.

VI. CONCLUSION
Considering the information on canopy structure and land "background," we used the "graphical method" to modify the linear PDM and derived a modified linear PDM. By establishing the relationship between linear PDM and Beer's law-related green FVC, we have established a multiangle calculation equation. To solve this equation, a random sampling method as well as proposed in the inversion method of the KDM. Then, we proposed a PDKDM to calculate FVC from HSR images. PDKDM was validated using simulated, field-measured, and reference data. There is high consistency between FVC estimated by PDKDM and simulated data. Compared with the field-measured data, FVC estimated by PDKDM also showed good accuracy. The results showed that the lack of a method to estimate the multiangle VI in the MAM was addressed in this study, and PDKDM could estimate FVC from HSR images in arid areas. In addition, we realized a high computational efficiency for producing HSR FVC.

A. Derivation of Modified Linear PDM
Beer's Law (some references are called a penetration function, Fig. 12) can well describe the relationship between radiation and canopy structure [53]. Therefore, it can be used as a reference to modify the linear PDM. The expression of Beer's Law is where P(θ i ) is the gap fraction and E is a clumping index. G(θ i ) is the G function, it denotes the projection function of the leaves area in the θ direction [30], [54], L AI is LAI, which is a function of VI. Since the gap fraction and FVC together form the scene, the expression of FVC is is Beer's law-related FVC, which shows an incremental relationship of the exponential function between FVC and LAI (Fig. 12). The relationship between LAI and VI also shows an incremental relationship. Therefore, based on the logical derivation for the relationship among FVC, LAI, and VI, the relationship between FVC and VI shows an incremental relationship of nonlinear function least (Fig. 12).
The linear PDM is a typical linear function (Fig. 12), which ignores the multiple scattering, the total radiation received by the satellite sensor is a linear combination between each component, there is In radiative transfer, the interaction between radiation and canopy is divided into single scattering and multiple scattering [30], [38]. Equation (A-3) is established using VI calculated by the single scattering. The multiple scattering is caused by the canopy structure and land "background." Some nonlinear changes of the FVC curve caused by canopy structure and land "background" are not considered in (A-3) [55]. Therefore, (A-3) is difficult to fit the nonlinear relationship between radiation and components, that is, the nonlinear relationship between VI and FVC in the linear PDM.
To address this problem, we use the "graphics method" to analyze graphics differences and derive the parameters with certain physical mechanisms. In this study, we introduce the Lagrangian polynomial to approximate the FVC curve simulated by Beer's law. In mathematics, a Lagrangian polynomial of three different points can fit a nonlinear incremental function [56], and the general equation is Here, P 2 (x) is the mathematical symbol for fitting the FVC sequence using a polynomial, (x 1 , y 1 ), (x 2 , y 2 ), and (x 3 , y 3 ) are three random points for FVC calculated by Beer's law, i.e., (A-3). Let x 1 be VI min , x 2 be VI random , and x 3 be VI max , here VI is NDVI. In general, when NDVI is the minimum value, FVC is approximately equal to 0. Therefore, y 1 is approximately equal to 0. Similarly, y 3 is approximately equal to 1, and FVC random is equal to y 2 . Finally, three points are determined, i.e., (VI min , 0), (VI random , FVC random ), and (VI max , 1). According to these assignments, (A-4) became Here, we used the method of symbolic substitution in mathematics. Let FVC radom ((VI − VI max )(VI − VI min ))/((VI min − VI max )(VI min − VI radom )) = m, VI − VI max = o, and VI max − VI radom = n, then (A-5) becomes Equation (A-6) has the multiplication symbols in the numerator and denominator, and hence, the FVC calculated based on this equation cannot achieve normalization, i.e., the value range of FVC from 0 to 1. Therefore, (A-6) is further derived as Let (nm/o)(VI max − VI min ) = l, and (no − 1) (VI max − VI min ) = L, then (A-7) becomes Here, l and L are adjustment factors; they physically represent comprehensive radiation information (multiple-scattering) influenced by the canopy structure and land "background." These factors are used and they can simplify complex mathematical expressions and ensure the physical characteristics. Meanwhile, these factors are very suitable for the design of accelerating algorithms, especially iterative computation in the solution of the multiangle calculation equation (B-7)-(B-10).

B. Solution of Multiangle Calculation Equation
To solve the multiangle calculation equation [i.e., (7)], we divided components, i.e., acquire VI(55 • ) and VI(60 • ), and calculated four unknowns about input parameters in the modified linear PDM: VI max , VI min , L, and l.
1) Method to Acquire VI(55 • ) and VI(60 • ): The images of VI(55 • ) and VI(60 • ) are calculated using the images of r NIR and r RED at the solar zenith angles of 55 • and 60 • , respectively, (3). The images of r NIR and r RED at the solar zenith angles of 55 • and 60 • are calculated from the images of the parameters of the linear KDM, i.e., f iso , f vol , and f geo [57]. At present, the spatial resolution of the parameter product of linear KDM is more than 500 m, e.g., the BRDF/albedo model parameter product [58]. To further calculate the HSR images of these parameters, the linear KDM was introduced where r (θ i ) denotes the black-sky albedo at θ i , and g 0iso , g 1iso , g 2iso , g 0vol , g 1vol , g 2vol , g 0geo , g 1geol , and g 2geo are polynomial coefficients of the linear KDM; their values are presented in Table II [57]. Generally, an HSR satellite is dominated by nadir observation; we let the initial r (θ i ) be a reflectance image with HSR in the nadir observation. To calculate the HSR images for f iso , f vol , and f geo , we need to perform an inversion of (B-2). Therefore, we proposed an IA to invert (B-2), which introduces the idea of combining upscaling and downscaling to avoid fluctuations caused by high vegetation abundance.
We assumed that the 2 × 2 moving window [vector on the right-hand side of (B-2)] has the same parameters of the linear KDM (i.e., F in Fig. 15). Then, we extracted the reflectance in the 2 × 2 moving window in the HSR image with the nadir observation, i.e., r (x,y) (θ i ), r (x+1,y) (θ i ), r (x,y+1) (θ i ), and r (x+1,y+1) (θ i ) [ Fig. 15(a)]. A 2 × 2 moving window is used to avoid the blocks effect caused by the weight of the distance from the adjacent pixel to the rescaled pixel in the previous 3 × 3 moving window [ Fig. 16(a)] [47]. According to (B-1), we established four equations for three unknowns, i.e., f iso ,  f vol , and f geo in which  where r (x,y) (θ i ), r (x+1,y) (θ i ), r (x,y+1) (θ i ), and r (x+1,y+1) (θ i ) can be obtained from an HSR image of reflectance in the nadir observation, see Fig. 15; θ i denotes the solar zenith angle, which is obtained from the metadata. When the solar zenith angle equals 0, g iso = g vol = g geo = 0, then The reflectance calculated at this time has a problem for (B-2). Therefore, we need to find an approximate solar zenith angle to approximate the solar zenith angle of 0 • , then (B-2) can be solved. According to Table III, the dataset of black-sky albedo simulated by the PROSAIL model is established [41], [42], [43], [59], and the difference in the black-sky albedo is small [root-mean-square error (RMSE) < 0.0046 in Fig. 13(a) and Fig. 14(a)] when the solar zenith angle is less than 5 • performing the sensitivity analysis. Therefore, when the solar zenith angle is less than 5 • , we approximate it to 5 • ; then, (B-2)-(B-5) yield nonzero values and (B-2) can have an approximate solution.
Here, other parameters in the PROSAIL model that have little influence on the black-sky albedo are fixed, such as the hotspot factor (l) = 0.1, sky visual factor (Skly) = 0.1, relative azimuth angle (ϕ) = 0 • , and brown pigment content (C brown ) = 0. Here, C ab is chlorophyll content, C ar is carotenoid content, C w is equivalent water thickness, C m is leaf mass per unit leaf area, and N is structure coefficient.
Because (B-2) is an overdetermined system, its solution is where F denotes a vector of the parameters of the linear KDM, [ f iso , f vol , f geo ] T . According to (B-6), the parameters of the linear KDM in the 2 × 2 moving window are calculated, which is similar to upscaling [ Fig. 15(c)]. With the application of the 2 × 2 moving windows, the step size of the movement is set to a pixel. In the 2 × 2 moving windows, the repeated pixels are assigned an average value [ Fig. 15(b) and (d)]. This step is similar to downscaling. According to these calculation steps, the reflectance of HSR images in the nadir observation is used, and the HSR images of the parameters of the KDM are acquired in (B-6). Finally, (B-1) are combined, the HSR images of r (55 • ) and r (60 • ) are calculated, then the HSR images of VI(55 • ) and VI(60 • ) can be obtained according to (3).
2) Calculation of Parameters Involved in the Modified Linear PDM: Once the HSR images of VI(55 • ) and VI(60 • ) are calculated, the next step is the calculation of the four unknown numbers in (2): VI min , VI max , l, and L. To determine the value ranges of VI min , VI max , l, and L, a previous study used a cyclic method to establish a dataset [26], [33], [36] (Table IV).  Fig. 16(a)] is used to acquire the pixels in the images of VI (55 • ) and VI (60 • ). A multiangle calculation equation [see (7)] with the unknowns is established to solve VI min , VI max , l, and L. However, similar to the 2 × 2 moving windows [Fig. B-4(b)] in the parameters of the linear KDM, the 3 × 3 moving window involves numerous iterative calculations, and the computational efficiency is low. With the improvement of spatial resolution, the number of pixels in the HSR images increases geometrically. Therefore, the computational time of the 3 × 3 moving window with the numerous iterative calculations is very long. Therefore, we used a chessboard block to split an image into subimages. Then, the Monte Carlo method is used to randomly obtain n pixels in the subimages to obtain reflectance and established multiangle calculation equations. This is called the random sampling method, i.e., chessboard block + Monte Carlo method [ Fig. 16(c)]. This method replaced the 3 × 3 moving window in the previous study to solve the multiangle calculation equation [33], and the equations are VI min , VI max , l, and L in (B-7) have certain physical meanings with certain value ranges. Therefore, we referred to previous studies [21], [26], [27], [33] and employed prior knowledge to determine the value ranges of VI min , VI max , l, and L, as shown in Table IV. He is currently an Associate Professor with the Faculty of Geo-Information Science and Earth Observation, University of Twente, Enschede. His research interest includes the development and use of advanced remote sensing and machine learning technologies for understanding, monitoring, and predicting biodiversity changes. His research helps inform environmental and wildlife management, supporting sound decision-making on the conservation and sustainable use of biodiversity in a rapidly changing world, which, in turn, contributes directly and indirectly to two United Nations Sustainable Development Goals (SDGs), namely, SDG 13 Climate Action and SDG 15 Life on Land.
Lei Lu received the Ph.D. degree from the Xinjiang Institute of Ecology and Geography, Chinese Academy of Science, Urumqi, China, in 2011.
She is currently an Assistant Professor with the College of Earth and Environmental Sciences, Lanzhou University, Lanzhou, China. Her research interests include land surface temperature and inversion of vegetation parameters.
Hui Sun received the bachelor's degree from the Hebei University of Technology, Tianjin, China, in 2021. She is currently pursuing the M.S. degree from Xingjiang University, Urumqi, China.
Her research interest includes quantitative remote sensing of biological soil crusts.
Fei Zhang is currently a Professor with the Xinjiang Key Laboratory of Oasis Ecology, Xinjiang University, Urumqi, China. His research interest includes remote sensing in the ecological environment of arid areas.
Xiao Cheng is currently a Professor with the Office of Scientific Research and Development, School of Geospatial Engineering and Science, Polar Science Center, Sun Yat-sen University, Zhuhai, China. His research interests include polar remote sensing, including regional climate system interaction between Antarctic ice and sea ice, inversion of key parameters in the sea ice, and monitoring technology of UAVs.
Ilyas Nurmemet is currently an Associate Professor with the Xinjiang Key Laboratory of Oasis Ecology, College of Geography and Remote Sensing Sciences, Xingjiang University, Urumqi, China. His research interests include microwave remote sensing and the applications of remote sensing technology in arid areas.