A Modiﬁed Geometrical Optical Model of Row Crops Considering Multiple Scattering Frame

: The canopy reﬂectance model is the physical basis of remote sensing inversion. In canopy reﬂectance modeling, the geometric optical (GO) approach is the most commonly used. However, it ignores the description of a multiple-scattering contribution, which causes an underestimation of the reﬂectance. Although researchers have tried to add a multiple-scattering contribution to the GO approach for forest modeling, di ﬀ erent from forests, row crops have unique geometric characteristics. Therefore, the modeling approach originally applied to forests cannot be directly applied to row crops. In this study, we introduced the adding method and mathematical solution of integral radiative transfer equation into row modeling, and on the basis of improving the overlapping relationship of the gap probabilities involved in the single-scattering contribution, we derived multiple-scattering equations suitable for the GO approach. Based on these modiﬁcations, we established a row model that can accurately describe the single-scattering and multiple-scattering contributions in row crops. We validated the row model using computer simulations and in situ measurements and found that it can be used to simulate crop canopy reﬂectance at di ﬀ erent growth stages. Moreover, the row model can be successfully used to simulate the distribution of reﬂectances (RMSEs < 0.0404). During computer validation, the row model also maintained high accuracy (RMSEs < 0.0062). Our results demonstrate that considering multiple scattering in GO-approach-based modeling can successfully address the underestimation of reﬂectance in the row crops.


Introduction
In recent canopy modeling studies, row modeling has been extensively studied. The canopy modeling of crops based on the row model can more accurately estimate the seasonal changes of biophysical canopy parameters in remote sensing [1]. In the inversion of remote sensing, physical modeling is key, and its approach can be separated into three categories based on different physical mechanisms, name the computer simulation (CS) approach, radiative transfer (RT) approach, and geometric optical (GO) approach [2]. The GO approach can describe the geometric characteristics of an individual canopy. It is most suitable for heterogeneous canopy modeling [3]. The CS approach has the highest accuracy and has been applied primarily for understanding radiation regimes. It has been used as a "truth value" ("surrogate truth" or benchmark in Radiation Transfer Model Intercomparison (RAMI)) [4][5][6] to validate GO and RT approaches. However, compared to the GO approach, the computational time of the CS approach is too slow; hence, it is rarely used in remote sensing inversion. The RT approach is based on volume scattering in the radiative transfer equation. It is very suitable for describing the scattering issue in canopies [7]. However, compared to the GO approach, the RT approach lacks a description of the geometric characteristics of an individual canopy and is mostly used for high-density homogeneous (or continuous) canopy modeling [2]. The row crop the reflectance at the top of the between-row area (r between_row ). The equation of reflectance at the top of the canopy of row crops is r = A 1 A 1 +A 2 r canopy_closure + A 2 A 1 +A 2 r between_row = A 1 r canopy_closure +A 2 r between_row A 1 +A 2 (1)

A A +
Here, A1 is the row width and A2 is the between-row distance. To facilitate a full understanding, r, rcanopy_closure, and rbetween_row are, respectively, called the sum of the reflectance of row crops, the reflectance of the canopy closure, and reflectance of the between-row area for short. In Equation (1), there are differences in the radiation mechanism of the canopy closure and between-row area. Therefore, the sum of the reflectance of row crops (r) is separated into the reflectance of the canopy closure (rcanopy_closure) and between-row area (rbetween_row) for consideration. According to Equation (1), the width of the row (A1) and between-row distance (A2) can be measured as input parameters. Therefore, rcanopy_closure and rbetween_row are key in row modeling. When rcanopy_closure and rbetween_row are calculated, it is implied that r can be calculated. As a result, we can establish a row model. Moreover, the sum of the reflectance of row crops (r) is composed of single-and multiple-scattering contributions [30,31]. Therefore, more detailed modeling examples for the single-and multiple-scattering contributions in the two areas are presented in the next two sections. Here, A1 is the row width, A2 is the between-row distance, and h is the canopy height.

Reflectance of the Canopy Closure
The reflectance of the canopy closure is the sum of the single-and multiple-scattering contributions inside the canopy [30,31]. We therefore have the following equation Here, rcc_1 is the single-scattering contribution of the canopy closure, and rcc_m is the multiplescattering contribution of the canopy closure. In the next two sections, we explain the modeling for rcc_1 and rcc_m.

Single-scattering contribution of the canopy closure
In the GO approach, the single scattering of the canopy closure is a linear combination of the four-component (illuminated leaf, shaded leaf, illuminated soil, and shaded soil [16]) area fraction and the corresponding representative reflectance for each component in the viewing direction [3,24,25]. Therefore _1 cc c c i i z z g g r S r Sr S r S r = + + + Here, rc is the reflectance of the illuminated leaf, ri is the reflectance of the shaded leaf, rz is the reflectance of the illuminated soil, and rg is the reflectance of the shaded soil. Sc is the area fraction of Here, A 1 is the row width, A 2 is the between-row distance, and h is the canopy height.
Here, A 1 is the row width and A 2 is the between-row distance. To facilitate a full understanding, r, r canopy_closure , and r between_row are, respectively, called the sum of the reflectance of row crops, the reflectance of the canopy closure, and reflectance of the between-row area for short. In Equation (1), there are differences in the radiation mechanism of the canopy closure and between-row area. Therefore, the sum of the reflectance of row crops (r) is separated into the reflectance of the canopy closure (r canopy_closure ) and between-row area (r between_row ) for consideration. According to Equation (1), the width of the row (A 1 ) and between-row distance (A 2 ) can be measured as input parameters. Therefore, r canopy_closure and r between_row are key in row modeling. When r canopy_closure and r between_row are calculated, it is implied that r can be calculated. As a result, we can establish a row model. Moreover, the sum of the reflectance of row crops (r) is composed of single-and multiple-scattering contributions [30,31]. Therefore, more detailed modeling examples for the singleand multiple-scattering contributions in the two areas are presented in the next two sections.

Reflectance of the Canopy Closure
The reflectance of the canopy closure is the sum of the single-and multiple-scattering contributions inside the canopy [30,31]. We therefore have the following equation r canopy_closure = r cc_1 + r cc_m (2) Here, r cc_1 is the single-scattering contribution of the canopy closure, and r cc_m is the multiple-scattering contribution of the canopy closure. In the next two sections, we explain the modeling for r cc_1 and r cc_m .

1.
Single-scattering contribution of the canopy closure In the GO approach, the single scattering of the canopy closure is a linear combination of the four-component (illuminated leaf, shaded leaf, illuminated soil, and shaded soil [16]) area fraction and the corresponding representative reflectance for each component in the viewing direction [3,24,25]. Therefore r cc_1 = S c r c + S i r i + S z r z + S g r g Remote Sens. 2020, 12, 3600

of 26
Here, r c is the reflectance of the illuminated leaf, r i is the reflectance of the shaded leaf, r z is the reflectance of the illuminated soil, and r g is the reflectance of the shaded soil. S c is the area fraction of the illuminated vegetation in the canopy closure, S i is the area fraction of the shaded vegetation in the canopy closure, S z is the area fraction of the illuminated soil in the canopy closure, S g is the area fraction of the shaded soil in the canopy closure.
We modified the area fraction equations of the four-component area fraction proposed by Verheof [33] in continuous crops and derived an expression suitable for the four-component area fraction in the canopy closure. The specific derivation can be seen in Supplementary Materials A, and the equations are Here, P o (θ o ,x,h) is the gap probability of the canopy closure in viewing direction, P so (θ s ,θ o ,x,h) is the bidirectional gap probability of the canopy closure, and h 0 P so (θ s , θ o , x, z)dz is the bidirectional vegetation probability of the canopy closure. k is the extinction coefficient of the canopy closure in the viewing direction Here, h is the canopy height, L is the leaf area index, and f (θ l ) is the leaf inclination distribution function (LADF). This study used an elliptic distribution function [34,35]. Combining Equations (4)-(7), P o (θ o ,x,h), P so (θ s ,θ o ,x,h) and h 0 P so (θ s , θ o , x, z)dz are the most important parameters for the calculation of area fraction of the canopy closure (S c , S i , S z , and S g ).
In the calculation of gap and vegetation probabilities, we proposed a new approach to calculate gap probabilities. In this new approach, we used a penetration function [21] to calculate gap probabilities in each path length. In this step, we considered the overlapping relationship between the leaves to calculate an average value of gap probabilities of the canopy closure (Supplementary Materials B-2). Furthermore, the average value of gap probabilities of the canopy closure has been used to represent the whole geometric characteristics of the canopy closure and was substituted into the calculation to analyze the overlapping relationship between the average value of gap probabilities of the canopy closure in the solar direction or the viewing direction. Therefore, the overlapping relationship between individual canopy closures could be considered in the calculation of gap probabilities (Supplementary Materials B-3). According to the above calculation ideas, we can consider both the overlapping relationship between leaves and individual canopy closures. In order to describe the bidirectional gap probability and vegetation probabilities, we modified the hotspot kernel function [25] originally used in forests to be suitable for row crops (Equation (B-25) in Supplementary Materials B-3), which can control the peak value near the hotspot. Based on the above, we attempted to address the computational deviations in the single-scattering contribution near the hotspot. For a detailed mathematical derivation, please refer to Supplementary Materials B.

2.
Multiple-scattering contribution of the canopy closure To calculate the multiple-scattering contribution, we assumed that the canopy closure is an isotropic scattering layer. The principle of adding in the RT approach (adding method) [36][37][38] is introduced (Figure 2). The derived reflectance of the canopy closure is Here, rs is the reflectance of soil and rs = (Szrz+ Sgrg)/(Sz + Sg). τs is the transmittance of the canopy closure in the solar direction, and τo is the transmittance of the canopy closure in the view direction. Equation (10) reflects the interaction (scattering) between light, vegetation, and soil, including the single-scattering contribution of the canopy closure and the multiple-scattering contribution of the canopy closure. Figure 2. Sketch of radiative transfer inside an isotropic scattering layer and the soil. Es is the downward irradiance on the horizontal plane, rcc_1 is the single scattering contribution of the canopy closure, τs is the transmittance of the canopy closure in the solar direction, and τo is the transmittance of the canopy closure in the solar direction. τ is the optical thickness, k is the extinction coefficient, and s is the path length.
We removed the single-scattering contribution (rcc_1) in Equation (10), hence the multiplescattering contribution of the canopy closure is S r S r S r S r S S r S r S r S S r S r S r r S S τ τ τ τ To address the lack of leaf transmittance in the optical input parameters of the GO approach, we introduce the study by Lang [39]. In this study, the transmittance of the canopy is approximately the same as the gap probabilities and replaces τs and τo in Equation (11). Therefore, the multiple scattering of the canopy closure is Previous studies have pointed out that there is a radiation energy exchange between vegetation and soil in the between-row area, which is further affected by multiple scattering between the soil and adjacent canopy closure between the rows [30][31][32]. The between-row area is an area where both the soil (C' in Figure 3a) and the vegetation (A' and B' in Figure 3a) exist, hence the multiple-scattering equation is very complicated. To establish the multiple-scattering equation of the between-row area, we introduced the mathematical solution of integral radiative transfer equation derived by Ma et al. [30]. We parameterized this equation according to the requirements of the GO approach. Then, an equation of reflectance in the between-row suitable for the GO approach was derived. Here, r s is the reflectance of soil and r s = (S z r z + S g r g )/(S z + S g ). τ s is the transmittance of the canopy closure in the solar direction, and τ o is the transmittance of the canopy closure in the view direction. Equation (10) reflects the interaction (scattering) between light, vegetation, and soil, including the single-scattering contribution of the canopy closure and the multiple-scattering contribution of the canopy closure.
We removed the single-scattering contribution (r cc_1 ) in Equation (10), hence the multiple-scattering contribution of the canopy closure is To address the lack of leaf transmittance in the optical input parameters of the GO approach, we introduce the study by Lang [39]. In this study, the transmittance of the canopy is approximately the same as the gap probabilities and replaces τ s and τ o in Equation (11). Therefore, the multiple scattering of the canopy closure is r cc_m = P o (θ o , h)P s (θ s , h) S z r z + S g r g S z + S g − r cc_1 S z r z + S g r g Remote Sens. 2020, 12, 3600 7 of 26

Reflectance of the Between-Row Area
Previous studies have pointed out that there is a radiation energy exchange between vegetation and soil in the between-row area, which is further affected by multiple scattering between the soil and adjacent canopy closure between the rows [30][31][32]. The between-row area is an area where both the soil (C' in Figure 3a) and the vegetation (A' and B' in Figure 3a) exist, hence the multiple-scattering equation is very complicated. To establish the multiple-scattering equation of the between-row area, we introduced the mathematical solution of integral radiative transfer equation derived by Ma et al. [30]. We parameterized this equation according to the requirements of the GO approach. Then, an equation of reflectance in the between-row suitable for the GO approach was derived.  < . Here, αo is the inclined angle projected in the perpendicular plane in the viewing direction, α1 is the openness angle of the between-row area, and α2 is the nonopenness angle of the between-row area. sbr is the path length of the light of the canopy closure for the between-row soil between being observed, r ϕ is the row azimuth angle, A1 is the row width, A2 is the between-row distance, and h is the canopy height.

Multiple-scattering contribution of the between-row area
Multiple scattering of the between-row area (rbr_m) occurs between the soil and adjacent canopy closure between the rows. Ma et al. gave a multiple-scattering equation of the between-row area based on the operation of the differentia integral operator [40], and its solution is Here, Kbr is the transfer probability of the collision, and Kbr = kbrPbr. Here, kbr is the extinction coefficient of the between-row area, and Pbr is the visual probability of each surface in the betweenrow area. The between-row area consists of four surfaces (vegetation surface (A' in Figure 3a To calculate Pbr, we defined the openness angle of the between-row area (α1) and nonopenness angle of the between-row area (α2) (Figure 3a,b). Then, we can calculate the visual probability of the escape surface (Popen) as α1/(α1 + α2), e.g.,  Here, α o is the inclined angle projected in the perpendicular plane in the viewing direction, α 1 is the openness angle of the between-row area, and α 2 is the nonopenness angle of the between-row area. s br is the path length of the light of the canopy closure for the between-row soil between being observed, ϕ r is the row azimuth angle, A 1 is the row width, A 2 is the between-row distance, and h is the canopy height.
The reflectance of the between-row area (r between_row ) is divided into two components: one is the single scattering of the between-row area (r br_1 ) and the other is the multiple scattering of the between-row area (r br_m ), and the expression is 1.

Single-scattering contribution of the between-row area
The single-scattering contribution of the between-row area is the reflectance of the soil that can be observed. According to the projection relationship between canopy closures in the GO approach, the single scattering of the between-row area is Remote Sens. 2020, 12, 3600

of 26
Here, l is the horizontal projected length of the row height on the ground in the solar or viewing directions (for its expression, please refer to Supplementary Materials B-2). When the viewing direction is considered (Figure 3b,c), Equation (14) becomes where P o_br (θ o , x, h) is the average viewing probability of the between-row soil, and Here, s br (θ o ,x,h) is the path length of the light in the canopy closure in the between-row soil between being observed (Figure 3b,c), and where α o is the inclined angle projected in the perpendicular plane in the viewing direction (its specific description can be found in Supplementary Materials B-1), β o is the azimuth of the inclined angle in the viewing direction, and sinβ o = sinϕ ro sin|θ o |/sinα o . x r is a remainder on the x-axis, and Here, mod is the mathematical symbol of modulus, and · is the mathematical notation for rounding down. Their sketches are shown in Figure 3b,c. Similarly, the average viewing probability at z of the between-row area (P o_br (θ o , x, z)) only needs to replace h with z.

2.
Multiple-scattering contribution of the between-row area Multiple scattering of the between-row area (r br_m ) occurs between the soil and adjacent canopy closure between the rows. Ma et al. gave a multiple-scattering equation of the between-row area based on the operation of the differentia integral operator [40], and its solution is Here, K br is the transfer probability of the collision, and K br = k br P br . Here, k br is the extinction coefficient of the between-row area, and P br is the visual probability of each surface in the between-row area. The between-row area consists of four surfaces (vegetation surface (A' in Figure 3a Figure 3a), and escape surface). To calculate P br , we defined the openness angle of the between-row area (α 1 ) and nonopenness angle of the between-row area (α 2 ) (Figure 3a,b). Then, we can calculate the visual probability of the escape surface (P open ) as α 1 /(α 1 + α 2 ), e.g., ∠dac/(π/2), ∠dgc/π, and ∠dbc/(π/2), shown in Figure 3a Here, ϕ r is the row azimuth angle. Equation (19) is the radiation escape probability at the bottom of the between-row area ( Figure 3a). For the escape probability at z of the between-row area, h in Equation (19) needs to be replaced by z (Figure 3b). Since different positions (different coordinate Remote Sens. 2020, 12, 3600 9 of 26 points (x,y) in Figure 3b) in the between-row area have a corresponding radiation escape probability, there are many radiation escape probabilities in the between-row area. In terms of modeling, we focus on the average radiation escape probability value in the between-row area, which is When the viewing direction is considered, Equation (18) can become The average visual factor of the vegetation surface (A' and B' in Figure 3a) and one soil surface (C' in Figure 3a) is defined as P other , and P other = 1 − P open . Combining Equations (18) and (21), the multiple scattering of the between-row area becomes From Equation (22), the extinction coefficient of the between-row area is key in calculating the reflectance of the between-row area. To apply Equation (22) in the GO approach, we assumed that the extinction coefficient of the between-row area is the weight sum of the two mediums, soil and leaf. The weight is determined by the length ratio of the medium in this area, which is Here, k is the extinction coefficient of the canopy closure in the viewing direction (A' and B' in Figure 3a), and its equation is Equation (8). k s is the extinction coefficient of the between-row soil in the viewing direction (C' in Figure 3a). According to [30], the extinction coefficient of the between-row soil in the viewing direction is Here, p(δ) is the scattering phase function of the soil particle, which represents the second-order Legendre polynomial (an approximation of spherical function) [41], p(δ) = 1 + b cos δ + 0.5c 3 cos 2 δ − 1 , and cos δ = cos θ s cos θ o + sin θ s sin θ o cos ϕ so . Here, b and c are the adjustment parameters for the second-order Legendre polynomial in the scattering phase function of the soil particle, which can be determined by [42].

Preparations for Validation-Based Computer-Simulated Data
We used a 3D computer simulation to validate (or compare) the row models. In this study, the 3D computer simulation used an extended 3D Radiosity-Graphics Combined Model (RGM) [43,44]. The RGM model is a radiosity model based on a bilinear equation (a simple nonlinear differential equation) [2]. It uses a numerical calculation method (Gauss-Seidel) to calculate the scattering of polygons in a scene constructed by a computer graphics method with high calculation accuracy [43,44]. To calculate the sum of the reflectance of row crops (r) for the RGM model, we divided it into two steps. In the first step, we used computer graphics to establish a computer scene similar to the assumption of our row model (the turbid medium bound in the periodic box-shaped vegetation material) (abstract scenes in Figure 4). For the abstract scene, we generated four representative scenes of row crops, namely the proportion of between-row dominance (Stage_rc1), proportions of between-row and canopy closure equality (Stage_rc2), the proportion of the canopy closure dominance (Stage_rc3), and continuous crops (Stage_cc). The parameters constructed in the abstract scene are shown in Table 1.
In the second step, based on the established abstract scene, we used the RGM model to calculate polygonal scattering in the abstract scene, and finally calculated the sum of the reflectance of row crops (r). Moreover, the sum of the reflectance of row crops (r) calculated by the RGM model can be used as a "true value" to validate the sum of the reflectance of row crops calculated by our model with the same parameters in Table 1. In setting the angle, the solar zenith angle was 25 • and the solar azimuth angle was 130 • for both models. Moreover, to keep the computer scene consistent with the periodic box-shaped plant materials assumed by our model, we used the infinite canopy of the RGM model in this study (Appendix B in [43]). This set implies that the reflectance calculated by RGM is not the reflectance of the one-or two-row cycle shown in Figure 4, but the reflectance of the scene where the row cycle is infinitely extended. The output results of reflectance are shown in Section 3.1. material) (abstract scenes in Figure 4). For the abstract scene, we generated four representative scenes of row crops, namely the proportion of between-row dominance (Stage_rc1), proportions of betweenrow and canopy closure equality (Stage_rc2), the proportion of the canopy closure dominance (Stage_rc3), and continuous crops (Stage_cc). The parameters constructed in the abstract scene are shown in Table 1. In the second step, based on the established abstract scene, we used the RGM model to calculate polygonal scattering in the abstract scene, and finally calculated the sum of the reflectance of row crops (r). Moreover, the sum of the reflectance of row crops (r) calculated by the RGM model can be used as a "true value" to validate the sum of the reflectance of row crops calculated by our model with the same parameters in Table 1. In setting the angle, the solar zenith angle was 25° and the solar azimuth angle was 130° for both models. Moreover, to keep the computer scene consistent with the periodic box-shaped plant materials assumed by our model, we used the infinite canopy of the RGM model in this study (Appendix B in [43]). This set implies that the reflectance calculated by RGM is not the reflectance of the one-or two-row cycle shown in Figure 4, but the reflectance of the scene where the row cycle is infinitely extended. The output results of reflectance are shown in Section 3.1.

Preparations for Validation-based In Situ Data
The measurements include two sites in arid Northwestern China: Zhangye, Gansu Province, and Zhongwei, Ningxia Hui Autonomous Region. Field and satellite measurements were performed in Zhangye and Zhongwei, respectively.    1 Here, the inclined leaf angle and azimuth leaf angle are set to a random distribution. By counting the leaf inclination angle and width of the leaves for each polygon, we obtained the average leaf inclination angle (θ l ) and the characteristic width of leaves (W p ).

Preparations for Validation-Based In Situ Data
The measurements include two sites in arid Northwestern China: Zhangye, Gansu Province, and Zhongwei, Ningxia Hui Autonomous Region. Field and satellite measurements were performed in Zhangye and Zhongwei, respectively.

1.
Field measurement data in Zhangye The in situ data of Zhangye came from Watershed Allied Telemetry Experiment Research (WATER) [45,46] and were measured from 20 May to 11 July 2008. A plot with an area of 180 × 180 m was selected, and the center coordinates of the plot were 38.857056 • N, 100.410444 • E ( Figure 5). The type of crop in the plot was corn, and four quadrats were randomly selected to measure the required parameters for validation. In the measurement of optical parameters, the reflectances of the illuminated leaf (r c ), shaded leaf (r i ), illuminated soil (r z ), and shaded soil (r g ) were measured by the ASD FieldSpec Pro Spectrometer [47,48]. The sum of the reflectance of row crops (r) (blue line in Figure  9) was measured by the ASD FieldSpec Pro Spectrometer, and the distance from the sensor to the top of the canopy was about 1 m, ensuring that at least one row cycle was observed. The measured spectral curve was from 400 to 2500 nm. The measurement time was 22 May, 25 May, 1 June, 16 June, 22 June, and 1 July (Table 2) [47,48]. Since the row structure was obvious on 22 June (Table 2), we performed a multiangle observation. The sum of the reflectance of row crops (r) (blue triangle scatter and blue line in Figure in Section 3.2.2) was measured by the ASD FieldSpec Pro Spectrometer combined with the multiangle frame. The distance between the instrument and ground was 5 m, the field of view was 25 • , and a complete row cycle was observed. The zenith angle was from −60 • to 60 • with 10 • for an interval [49]. The study considered the anisotropy of reflectance. Measurements were separated by four modes in the azimuth: the principal plane (PP), orthogonal plane (OP, viewing azimuth angle perpendicular to the sun azimuth angle), along-row plane (AR, the plane along the row direction), and orthogonal row plane (OR, the plane perpendicular to the row direction). The other structure parameters used the direct measurement method [50][51][52]. These parameters were the leaf area index (L), average leaf inclination angle (θ l ), row width (A 1 ), between-row distance (A 2 ), canopy height (h), characteristic width of leaves (W p ), and row azimuth angle (ϕ r ). The measurement results are shown in Table 2. The output results of the spectral curve for the growing season are displayed in Figure 9. Corresponding output results of the distribution of reflectances in the multiangle observation (22 June, in Table 1) are displayed in Section 3.2.2. Figure 9) was measured by the ASD FieldSpec Pro Spectrometer, and the distance from the sensor to the top of the canopy was about 1 m, ensuring that at least one row cycle was observed. The measured spectral curve was from 400 to 2500 nm. The measurement time was 22 May, 25 May, 1 June, 16 June, 22 June, and 1 July (Table 2) [47,48]. Since the row structure was obvious on 22 June (Table 2), we performed a multiangle observation. The sum of the reflectance of row crops (r) (blue triangle scatter and blue line in Figure in Section 3.2.2) was measured by the ASD FieldSpec Pro Spectrometer combined with the multiangle frame. The distance between the instrument and ground was 5 m, the field of view was 25°, and a complete row cycle was observed. The zenith angle was from −60° to 60° with 10° for an interval [49]. The study considered the anisotropy of reflectance. Measurements were separated by four modes in the azimuth: the principal plane (PP), orthogonal plane (OP, viewing azimuth angle perpendicular to the sun azimuth angle), along-row plane (AR, the plane along the row direction), and orthogonal row plane (OR, the plane perpendicular to the row direction). The other structure parameters used the direct measurement method [50][51][52]. These parameters were the leaf area index (L), average leaf inclination angle (θl), row width (A1), between-row distance (A2), canopy height (h), characteristic width of leaves (Wp), and row azimuth angle (φr). The measurement results are shown in Table 2. The output results of the spectral curve for the growing season are displayed in Figure 9. Corresponding output results of the distribution of reflectances in the multiangle observation (22 June, in Table 1) are displayed in Section 3.2.2.

Satellite measurement data in Zhongwei
The satellite measurement was performed on 16 July, 2016 at Zhongwei city. In the optical parameter measurements, the sum of the reflectance of row crops (r) was measured via satellite with a high spatial resolution (World-View 3 image in Figure 4b). Therefore, the measured height was about 617 km. The World-View 3 image has 8 bands from 400 to 1040 nm, and the spatial resolution is 1.24 m. To obtain the sum of the reflectance of row crops, we preprocessed the image, involving radiation calibration, geometric correction, and atmospheric correction [53]. The surface objects in the image have three types of row crops: rice, corn, and matrimony vine (strictly speaking, matrimony vine is an economic forest planted in row form). To obtain the optical parameters required for the row model, the seven typical quadrats with the same size as the resolution of the World-View 3 image (Figure 5b) were directly selected to measure the reflectance of the illuminated leaf (r c ), shaded leaf (r i ), illuminated soil (r z ), and shaded soil (r g ). For these optical parameters, we used the ASD FieldSpec Pro Spectrometer on the ground. The other structure parameters used the direct measurement method, and the measurement results are shown in Table 3. The corresponding output results simulated by our model are shown in second part of Section 3.2.1. Table 2. Structural parameters of corn in the growing season.

Phenology
Date (2008)  Table 3. Structural parameters of crops in the quadrats.

Validation of Row Model Using Computer-Simulated Data
In Figure 6, the distribution of the sum of the reflectance of two models in four scenes (Stage_rc1, Stage_rc2, Stage_rc3, and Stage_cc) are shown. With an increase in LAI (leaf area index) and a change in row structure (increase in the row width and height, decrease in the between-row distance) (Table 1), the distribution of the sum of the reflectance in the red band and NIR band simulated by the row model and the distribution of the sum of the reflectance in the red band and NIR band simulated by the RGM model have high consistency in four viewing modes (PP mode, OP mode, AR mode, and OR mode). The correlation coefficients (R) are greater than 0.9281. The root mean square errors (RMSEs) are less than 0.0012 in the red band and less than 0.0095 in the NIR band (Table 4). Moreover, the differences between the reflectance simulated by the two models are less than 5%. The consistency of the red band (mean absolute deviation (MAD) = 0.25 × (1.98% + 3.08% + 4.64% + 4.03%) ≈ 3.43%) is greater than that of the near-infrared band (MAD = 0.25 × (1.87% + 2.99% + 1.68% + 1.82%) ≈ 2.08%) ( Table 4). With an increase in LAI and a change in row structure, the sum of the reflectance in the red band and NIR band shows opposite trends, i.e., the sum of the reflectance in the red band is decreasing while the sum of the reflectance in the NIR band is gradually increasing. Our model is the same as the RGM model in the inverted v-shaped area where the hotspot is located (PP in Figure 6a,b). When Stage_rc1 changes to Stage_cc, the width of hotspots (slope of an inverted v-shaped area) gradually narrows. Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 27    Figure 7 shows the decomposition of the reflectance in the NIR band in Figure 6b,d,f,h, namely the single-scattering contribution in Figure 7a,c,e,g and multiple-scattering contribution in Figure 7b,d,f,h. In Figure 7a,c,e,g, the overall difference between the single-scattering contribution in the NIR band simulated by the row model and the single-scattering contribution in the NIR band simulated by the RGM model is small. The R-value is greater than 0.8974 and RMSEs are less than 0.0057 for single scattering in the NIR band in Table 4. Moreover, the differences between the single scattering in the NIR band simulated by the two models are less than 3% (Table 4). Figure 8 shows the dense points in the single scattering near the hotspot in PP, with our model and RGM model also having high consistency. In addition to the change in the width of the hotspot, the changing trend of the single-scattering hotspot is more obvious than the sum of the reflectance in Figure 6a,b.
In Figure 7b,d,f,h, the multiple-scattering contribution in the NIR band simulated by our model and multiple-scattering contribution in the NIR band using the RGM model have high consistency. The R-value is greater than 0.7971 and RMSEs are less than 0.0062 for multiple scattering in the NIR band in Table 4. Moreover, the differences between the multiple scattering in the NIR band simulated by the two models are less than 10% (Table 4). Comparing Figure 6b,d,f,h and Figure 7b,d,f,h, we can find that the proportion of the multiple scattering in the NIR band in the sum of the reflectance gradually changes from 20% to 50% as LAI increases and changes in row structure (from Stage_rc1 to Stage_cc).
Here, R_NIR_1 is the single scattering in the NIR band, and R_NIR_m is the multiple scattering in the NIR band. R is the correlation coefficient, RMSE is the root mean square error, and MAD is the absolute value of average difference, which represents the difference between Data A and Data B,  Validation of the sum of the reflectance during the growing season using field data In Figure 9, the sum of the reflectance simulated by the row model and the field data have high consistency in the growing season of row crops (it can be described as an increase of LAI and change in row structure in Table 2), except for the spectrum less than 1500 nm measured on 20 May (Figure 9a). As the crop grows over time, the photosynthetically active radiation (400-700 nm) gradually decreases. This phenomenon is especially noticeable at the "red edge" (700-750 nm, the area where the vegetation reflectance from the red band to the near-infrared band increases sharply). On 25 May (Figure 9b), the "red edge" gradually formed, and it was increasingly obvious with the growth of crops (Figure 9b-f). The sum of the reflectance in the near-infrared (NIR) region (700-1100 nm) gradually increases with the growth of crops; after that, it reaches a maximum value (0.51) on 1 July (Figure 9f).

Spectral Curve
1. Validation of the sum of the reflectance during the growing season using field data In Figure 9, the sum of the reflectance simulated by the row model and the field data have high consistency in the growing season of row crops (it can be described as an increase of LAI and change in row structure in Table 2), except for the spectrum less than 1500 nm measured on 20 May ( Figure  9a). As the crop grows over time, the photosynthetically active radiation (400-700 nm) gradually decreases. This phenomenon is especially noticeable at the "red edge" (700-750 nm, the area where the vegetation reflectance from the red band to the near-infrared band increases sharply). On 25 May (Figure 9b), the "red edge" gradually formed, and it was increasingly obvious with the growth of crops (Figure 9b-f). The sum of the reflectance in the near-infrared (NIR) region (700-1100 nm) gradually increases with the growth of crops; after that, it reaches a maximum value (0.51) on 1 July (Figure 9f). 2. Validation of the reflectance for different types of row crops using satellite data We used remote sensing data to validate the row model. The sum of the reflectance simulated by the row model and the sum of the reflectance measured by the World-View 3 satellite have high

2.
Validation of the reflectance for different types of row crops using satellite data We used remote sensing data to validate the row model. The sum of the reflectance simulated by the row model and the sum of the reflectance measured by the World-View 3 satellite have high consistency. In the simulation of each species, the accuracy of the simulation of corn and wheat is the best, while the simulations of wolfberry had a slight deviation (Figure 10). June, and (f) 1 July. Here, the noise in the water vapor absorption was removed.

Validation of the reflectance for different types of row crops using satellite data
We used remote sensing data to validate the row model. The sum of the reflectance simulated by the row model and the sum of the reflectance measured by the World-View 3 satellite have high consistency. In the simulation of each species, the accuracy of the simulation of corn and wheat is the best, while the simulations of wolfberry had a slight deviation (Figure 10).

Distribution of the Sum of the Reflectance on the Multiangle Observation
In the distribution of the sum of the reflectance (Figure 11), the results simulated by the row model and the measured data have a slight systematic deviation from some angles. In the forward direction of the PP mode in the NIR band, the reflectance has the largest computational deviation between the simulation and measurement. However, from the statistical results, we can find that the R-values are greater than 0.8317 (except R = 0.5392 in the PP mode in the NIR band) and RMSE values less than 0.0403 ( Figure 12). The MAD results show that the difference between the sum of the reflectance simulated by the row model and field data is less than 13%, except for the OR mode in the red band (24.11% in Figure 12). Similar to the computer simulation in Figures 6 and 7, the consistency of the red band (MAD = 0.25 × (10.71% + 12.4% + 6.29% + 24.11%) ≈ 13.38%) is greater than that in the NIR band (MAD = 0.25 × (9.65% + 5.29% + 1.61% + 3.76%) ≈ 5.08%) ( Figure 12). Generally speaking, the consistency between the two is still relatively high, simulating multiangle changes in the PP mode, OP mode, AR mode, and OR mode. of the red band (MAD = 0.25 × (10.71% + 12.4% + 6.29% + 24.11%) ≈ 13.38%) is greater than that in the NIR band (MAD = 0.25 × (9.65% + 5.29% + 1.61% + 3.76%) ≈ 5.08%) (Figure 12). Generally speaking, the consistency between the two is still relatively high, simulating multiangle changes in the PP mode, OP mode, AR mode, and OR mode.

Multiple Scattering of Row Crops in the GO Approach
This study substantiates similar findings obtained in earlier studies by showing that it is necessary to consider the multiple-scattering contribution in the NIR band in the GO approach [13,26,27]. Moreover, we considered the multiple-scattering contribution in row crops as an extension of the previous study field (forest field) to the object of row crops. In the row modeling based on the

Multiple Scattering of Row Crops in the GO Approach
This study substantiates similar findings obtained in earlier studies by showing that it is necessary to consider the multiple-scattering contribution in the NIR band in the GO approach [13,26,27]. Moreover, we considered the multiple-scattering contribution in row crops as an extension of the previous study field (forest field) to the object of row crops. In the row modeling based on the GO approach, its initial purpose was to provide a physical basis to calculate temperature (or emissivity). However, due to the different physical properties between temperature (or emissivity) and reflectance, the previous studies rarely considered multiple scattering issues in row modeling [14,28,29]. Although subsequent studies attempted to consider multiple scattering issues in the reflectance modeling of row crops, these models usually consider a row structure to modify existing continuous crop models based on the RT approach [17,31,54,55]. Therefore, their physical basis was still the volume scattering of the radiative transfer equation, and their essence was the RT model. Different from the previous modeling, our model took the multiple-scattering equation constructed by the RT approach in GO modeling further to modify the GO approach. The multiple-scattering equation achieved a good calculation accuracy (Figure 7b,d,f,h). In the multiple-scattering contributions of row crops, the radiation energy of multiple scattering in the NIR band was very strong, accounting for a large proportion of the sum of the reflectance (comparing Figure 6b,d,f,h and Figure 7b,d,f,h). Specifically, as crops gradually grew from row crops to continuous crops, the proportion of multiple scattering in the NIR band in the sum of the reflectance gradually increased from 20% to 50% (comparing Figure 6b,d,f,h and Figure 7b,d,f,h). These results imply that the sum of the reflectance would be significantly underestimated if the multiple-scattering contribution is not considered, and the underestimation becomes more obvious as LAI gets bigger. Therefore, the multiple-scattering equation is considered in GO modeling and can be accurately calculated, which can effectively improve the accuracy of the sum of the reflectance in the NIR band. However, in Table 4, the mean absolute deviation (MAD) of R_NIR_m gradually decreased from 9.82% in the Stage_rc1 to 0.99% in the Stage_cc. This result implies that the row structure influences the calculation accuracy of multiple scattering, whereby, when the row width decreases, the between-row distance increases, and the difference in the multiple-scattering contribution between our model and the RGM model will increase. This conclusion is consistent with the simulation results in [30]. The primary reason is that our model describes Lambertian soil (soil reflectance as a single value, see Equations (13) and (22)), but the RGM model describes non-Lambertian soil (the brown squares shown in Figure 4 are approximated as an anisotropic soil, hence soil reflectance is a multivalue consisting of multiple squares). Therefore, the increase in the between-row distance will cause a difference in the multiple scattering between our model and the RGM model.
In the GO approach, we studied an adding method (Section 2.1.2) and a mathematical solution of integral radiative transfer equation (Section 2.1.3) to derive the multiple-scattering equation. This is a typical method for the multiple-scattering equation, which was established by the RT approach, to be coupled into the GO approach for a modified GO model. This method is based on the RT approach to derive the multiple-scattering equation. However, in a previous study (mainly in the heterogeneous forest), this method often ignored the geometric characteristics between individual canopies [26], which was reflected in the results as a computational deviation in the sum of the reflectance [13,27]. The between-individual canopies mentioned in [13,27] are the between-row area in row crops. Our research focused on this point in row modeling. In our model, the openness angle of the between-row area (α 1 ) and nonopenness angle of the between-row area (α 2 ) were introduced to calculate an average value of the escape probability of radiation in the between-row area (P open ). The above treatment showed the scattering direction (or scattering angle with typical geometric characteristics) in the multiple-scattering contribution of the between-row area, and further reflected the influence of geometric characteristics (A 2 , ϕ r and h in Equation (20)) on the multiple-scattering equation (Equation (22)). The multiple-scattering contribution simulated by our model had higher R-values and lower RMSEs compared to the RGM model ( Table 4). The modified treatments of these equations and results further illustrated that our model considered how the geometric characteristics between individual canopies can improve the calculation accuracy of the multiple-scattering contribution. We further analyzed the sum of the reflectance simulated by our model compared to those of the RGM model was of high consistency ( Figure 6). The simulation results and in situ measurement of our model also have good accuracy overall (Figures 9-11). These results show that geometric characteristics considered in the multiple-scattering equation established by the RT approach are necessary for GO modeling. This consideration not only improves the accuracy of multiple scattering but also the accuracy of the sum of the reflectance.
When we analyzed the multiple-scattering equation in the previous model in the heterogeneity forest (Equation (27) in [27] and Equations (8)- (11) in [13]), we found that the single-scattering albedo of the leaf was added to derive the multiple-scattering equation. The purpose of this treatment in the previous models was to remedy the defects of the GO approach without the scattering phase function in the calculation of scattering. However, the single-scattering albedo of the leaf is difficult to determine in the GO approach without leaf transmittance [7]. Therefore, the single-scattering albedo of the leaf needs to adopt an empirical approach, which leads to the lack of a physical mechanism at the foundation of multiple scattering. In this study, we focused on the above limitation. In the multiple-scattering equation of row crops, we used a single-layer adding method ( Figure 2) and mathematical solution of integral radiative transfer equation (its mathematical properties are also cumulative, Equation (22)) to avoid introducing the single-scattering albedo of the leaf to calculate multiple scattering. We tried to use gap probabilities to approximate the transmittance in canopy closure (P o (θ o ,x,h) and P s (θ s ,x,h) in Equation (12)) and the transmittance in the between-row area (P o_br (θ o , h) and P open in Equation (22)). In this way, the gap probabilities become a connection that can seamlessly connect single-scattering contributions and multiple-scattering contributions together.

Hotspot of Row Crops in the Single-Scattering Contribution
The single-scattering contribution was accurately calculated, which was the basis of the multiple scattering calculations. Therefore, the issues of single scattering cannot be ignored in the calculation of the multiple-scattering contribution. Equations (12) and (22) require a single-scattering contribution as the original function to calculate the multiple-scattering contribution. In this study, we extended the modeling ideas of two-overlapping relationship (the overlapping relationship between leaves and canopy closures [21,22]) previously used for heterogeneity forest modeling to row modeling (Supplementary Materials B). By comparing single-scattering contributions near the hotspot in four canopies (the four growth stages of crops in Figure 4) simulated by our model and the RGM model, we found that the simulation results between the two models were very consistent (Figure 8), which indicates that our study achieved good accuracy. These results substantiate that our model considered the two-overlapping relationship to calculate gap probabilities (especially bidirectional gap probabilities and vegetation probabilities, Equations (B-24) and (B-32) in Supplementary Materials B-3), which effectively improved the calculation accuracy of single scattering near the hotspot in the GO approach. The mechanism for these results in Figure 8 shows that the shadow produced by the overlapping relationship between the canopy closures is gradually replaced by the shadow produced by the overlapping relationship for leaves when the row structure changes from row crops to continuous crops. In addition, the shadow area produced by the overlapping relationship between canopy closures is larger than the shadow area produced by the overlapping relationship for leaves. Therefore, the shadow near the hotspot in the crops will be reduced and single scattering near the hotspot will be enhanced. These results also imply that gap probabilities are calculated and the two-overlapping relationship is necessary.
The cause of the peak value near the hotspot involves bidirectional gap probability and bidirectional vegetation probability issues [2]. The reason for this phenomenon is that the shadow created by the overlap between the medium is not seen when the solar direction and the viewing direction are very close, and the reflectance will have a maximum value (the peak value near the hotspot) [56,57]. In continuous crops, the peak value near the hotspot is mainly influenced by LAI, leaf size, leaf orientation mode, and sun geometry [20]. However, in addition to the above factors, the peak value near the hotspot in row crops is also influenced by the row structure [31]. In our study, we found that as the row structure changes from row crops to continuous crops (from Stage_rc1 to Stage_cc), the width of the hotspot (slope of an inverted v-shaped area) gradually narrows, and the peak value becomes more obvious (Figure 8 for the single-scattering contribution and Figure 6a,b for the sum of the reflectance). Our study accurately simulated the width of the hotspot, which not only considered the two-overlapping relationship but also considered the hotspot kernel function. In order to describe the width of the hotspot, cross-correlation [19] or overlap functions [57] have been commonly used as hotspot kernel functions in the RT approach. However, our study was based on the GO approach, hence the above functions were not applicable. In our study, we used the hotspot kernel function (Equation (B-25) in Supplementary Materials B-3) proposed by Chen et al. [25] in the heterogeneity forest modeling based on the GO approach. To make it suitable for the canopy of row crops, our study analogized the internal relationship between forest parameters and row parameters to the modified C 2 in Equation (B-25), where C 2 is the control factor of the width of the hotspot. From C 2 (Equation (B-25)), it can be found that L, h, θ s , W p , and k (note that k is a function of A 1 , A 2 , h, θ o , θ l , and f (θ l ) in Equation (8)) were considered, which also implies that our study involved equation modeling that takes into account the main influencing factors described above, such as LAI, leaf size, leaf orientation mode, and the sun geometry, row structure. These results and equations in the modeling show that the hot kernel function is very important in simulating the peak value near the hotspot.

Analysis of the Row Model
As a GO model, our model mainly uses the geometric relationship between light and medium to calculate the sum of the reflectance. The RGM model is a computer model based on radiosity. It uses the geometric probability between leaves to calculate the sum of the reflectance. Specifically, a view factor describing the overlapping relationship between polygons (leaves) was previously introduced into the RGM model [43]. If we exclude the numerical calculation of reflectance (Gauss-Seidel algorithm [44]) in the RGM model, the geometric principles involved in our model are very similar to those of the RGM model. Therefore, based on similar physical principles, the consistency between the simulation results of our model and the RGM model is high (Figures 6-8, MAD is less than 10% in Table 4). However, our model still shows calculation deviation in the in situ validation, especially in the distribution of the sum of the reflectance on the multiangle observation ( Figure 11, MAD is less than 24% in Figure 12). These results imply that the calculation deviation may come from other sources. At present, it is difficult for multiangle instruments to accurately measure directional reflectance (specifically, bidirectional reflectance distribution function (BRDF) [58]) and obtain an approximate reflectance distribution. In theory, directional reflectance (BRDF) refers to the measurement value at the same time (the same second), and there is no shadow of the instrument (note: the instrument is nontransparent) during the measurement [58], which is hard to achieve with the current instruments. In the multiangle observation, instruments take at least 10 minutes or more to complete the four modes (PP, OP, AR, and OR). Even if the measurement is performed in strict accordance with the measurement specification, changes in environmental variables over time will have a strong influence on measurement [59]. Therefore, the results in the measurements are only an approximation (abnormal points in the black box in Figure 11a,c,e). The instrument had a field of view (FOV) of 25 • (common FOV of ASD spectrometer), and the observation distance was 5 m (this height is different from the current canopy reflectance model assuming that the sensor is located at infinity) in our study. Therefore, we observed the sum of the reflectance of a limited row cycle (a canopy closure plus a between-row area) in FOV (according to Table 1, we used the triangular relationship to calculate that the observed row cycle changes from 2.2 to 10.4 and from 0 • to ±60 • in the zenith angle), which is different from our model assumptions (periodic box-shaped plant materials). However, the validations of computer simulation have not been influenced by time in the multiangle observation. The infinite canopy is set in the RGM model in this study [43], which is the same as our model assumptions. This implies that the validation of computer simulation has more advantages than the current lack of accurate instruments. Our study reconfirmed the conclusion (or objective) from RAMI [4][5][6]60] and a wide range of mathematical modeling in earth sciences [59], and computer simulation is one of the most effective ways to validate the model of a nonnumerical calculation.
Since the leaves of crops are not randomly distributed in the real world, if the uncertainty in the measurement process is excluded, the most likely cause of this phenomenon is that the state of leaf aggregation on the horizontal plane has not been described, i.e., the clumping index [61]. The calculations of the clumping index and gap probabilities are inseparable. For the calculation of gap probabilities, we used a penetration function (P = e −ks in Supplementary Materials B-2). This function is based on the assumption that leaves are randomly distributed in the canopy. Therefore, there may be deviations between this assumption and the actual canopy, which implies that the equation describing the gap probabilities may need to be further refined considering the clumping index (P = e −ksΩ , where Ω is the clumping index). The clumping index is a hot topic in the field of radiation regimes and agricultural measuring equipment [62][63][64][65]. Compared with coniferous forests, the degree of leaf aggregation of crops is not very obvious. Therefore, considering the computational complexity, and whether the clumping index should be the main focus in improving the model, requires further study.

Conclusions
In this study, we established a row model to accurately estimate the reflectance by considering the multiple scattering framework. We validated the row model using both computer simulations and in situ measurements. The validation results show that the row model can be used to simulate crop canopy reflectance at different growth stages, with an accuracy comparable to the computer simulations. The study modified the calculation accuracy of single scattering near the hotspot by considering the calculation of the two-overlapping relationship and hot kernel function. Moreover, the row model can be successfully used to simulate the multiple-scattering contribution by considering multiple-scattering equations based on the RT approach. Therefore, this can address the underestimation problem in row crops by ignoring multiple scattering calculations in the GO approach. Our results demonstrate that the multiple-scattering contribution is a very important process that needs to be considered in row modeling based on the GO approach. This study provided a potential mechanism for remote sensing inversion based on the physical model.