Spatial Estimation of Daily Growth Biomass in Paddy Rice Field Using Canopy Photosynthesis Model Based on Ground and UAV Observations

1. Introduction

To ensure stable food production in response to global population growth, it is necessary to increase unit yield in agriculture adapted to limited farmland resources and local environments such as climate and meteorology [1,2]. For this purpose, monitoring of crop growth and adapting appropriate cultivation management are required. Monitoring the growth of crops and properly managing their cultivation are skills that farmers have long developed empirically. However, despite rising global food consumption and production, the number of farmers is declining [3]. As a result, worldwide efforts are being made to overcome this problem through efficient and effective precision farming.

Efficient crop phenology observation and growth monitoring at the field level are critical in precision farming. This requires not only labor-saving operations but also extremely accurate yield estimation and prediction. In cultivation management, it is critical to understand the differences in growth stages due to differences in environmental factors (e.g., weather and field conditions in the year of planting), and different production management (e.g., varieties and fertilizer amounts). Monitoring the growth condition of the entire field is also a useful way for determining the effects of abnormal and extreme weather on crops, which are expected to become more frequent as a result of global warming. Precision farming, which replaces traditional farming methods, (e.g., use of fixed-point cameras and ITC) makes it possible to quantitatively grasp the growth status of fields at any time, either by visual identification or by automatic digital processing [4,5]. Furthermore, regular UAV observation allows for detailed field monitoring [6]. UAVs have enabled accurate identification of agricultural growth conditions at the individual level by taking photographs from low altitudes, however, the finer spatial scale requires equally detailed information on the time-series changes in growth. Despite the advancements in UAV technology, the realistic observation frequency is still once a week in consideration of flight planning, operation, data processing, and other aspects.

On the other hand, crop growth, especially biomass, is affected by temporal changes in sunlight, which is most closely related to the rate of photosynthesis. Precision farming, which requires detailed information on both spatial and temporal scales, has the potential to estimate and predict a field's growth conditions spatiotemporally by combining fixed-point observations and UAV observations with meteorological data that can be input into crop growth models.

There are many cases related to the estimation of leaf area index (LAI) in applied research on remote sensing for precision farming [7,8,9,10,11,12,13,14]. This is because LAI is an essential parameter for input into crop models for biomass and yield estimation and prediction. Traditionally, obtaining LAI verification data has not been easy because it requires actual measurement by cutting samples, however there is a method to indirectly estimate LAI based on the light environment (light transmission by leaves) under the canopy [15]. LAI can be estimated non-destructively and effectively by estimating the light transmittance within a canopy using ground-based or UAV-based remote sensing observation data. Furthermore, the canopy photosynthesis model [15,16,17,18] used in several crop models allows photosynthetic photon flux density (PPFD) to correspond to diurnal changes in sun altitude. It allows the geometric calculation of light absorption and transmission received by leaves, as well as the simulation of the amount of growing biomass per day by inputting diurnal changes in PPFD and time-series LAI [19,20,21]. Therefore, the integrated use of observed data in the field and numerical models can be expected to provide continuous growth monitoring both spatially and temporally.

The purpose of this study was to investigate observational methods for efficient indirect estimation of LAI, which is essential for biomass estimation and yield prediction. Furthermore, daily biomass spatial estimation of daily biomass through application of canopy photosynthesis model. Data from ground- and UAV-based observations at two paddy rice field sites and rice cultivated in different years were used to assess a method for measuring the light environment under rice canopy (relative light intensity at the top and bottom of the canopy). In addition, the time-series LAI estimated indirectly from the relative light intensity and the observed light intensity used in photosynthesis (PPFD) were applied to a canopy photosynthesis model for calculating the amount of growth biomass per day. The validity of this method was further examined by comparing the accumulative biomass with the above-ground biomass.

2. Materials and Methods

2.1. Model description

2.1.1. Light distribution under canopy

The photosynthesis rate is determined by how much light is absorbed by the leaves under the canopy layer. Lambert-Beer’s law can describe the light attenuation absorbed by leaves from the top to under the canopy using this equation [15]:

$$I_i=I_0\mathrm{e}\mathrm{x}\mathrm{p}(-K{F}_i)$$

(1)

Whereas I_i is the horizontal light intensity at layer i in the canopy, I₀ is the horizontal light intensity, K is the extinction coefficient, and F_i is the leaf area index (LAI, m²/m²). I₀ denotes the photon flux density per unit leaf area per unit time on the top of the canopy. The F_i reaches a maximum at the bottom layer of the canopy and shows LAI under all layers. The relationship of the relative light intensity logarithm (I_i/I₀) in layers i and F_i is linear, with K as the slope. The extinction coefficient K is closer to 1 for horizontal leaves.

2.1.2. Canopy photosynthesis model

The canopy photosynthesis model based on Monsi and Saeki (1953) [15] has been widely used and modified by considering the light environment and the leaf morphology e.g.,[16,17,18,22]. It is based on a mathematical model that analyzes the light response to changes in the sun's elevation during the day and the absorption and transmission processes of light received by leaves for each variety of plant. Daily productivity (biomass) can be estimated in time series by inputting the incident light intensity and LAI e.g., [20,21].

In this study, the canopy photosynthesis model modified by Anten (1997) [18,19] is used, in which the total incident light on the top of the canopy is divided into direct and diffuse lights. This model is based on the calculation of the photosynthetic rate in each layer by separating the sunlit and the shaded leaves.

Light intensity received by sunlit leaves and shaded leaves: The canopy is divided vertically into multiple layers, and the intensity of light received by the leaves in each layer is calculated. The sunlit leaves receive both direct and diffuse light, while shaded leaves receive only diffuse light. The absorbed light intensity of the sun leaves (I_sl,i, mmol/m²/s) in layer i is expressed by:

$$I_{sl,i}=I_{sh,i}+\frac{O{I}_{0b}}{\sin {\beta }_s}$$

(2)

Where I_sh,i (mmol/m²/s) is the absorbed light intensity received by shade leaves in layer i, O is the projected area of leaves from the view of the sun, I_ob is the direct light (mmol/m²/s) received at the horizontal plane above the canopy, and b_s is the sun elevation angle. O is a parameter that varies depending on the leaf slope and the sun elevation, O is calculated by dividing the slopes into three classes [17]:

$$O=f_{15}{O}_{15}+f_{45}{O}_{45}+f_{75}{O}_{75}$$

(3)

Whereas O₁₅, O₄₅, O₇₅ are the projected area of slope leaves from 0 to 30, 30 to 60 and 60 to 90 degrees, , respectively and f₁₅, f₄₅, f₇₅ are the fraction of three slope classes. In this study, f₁₅, f₄₅, f₇₅ were set based on the growing conditions of rice leaves to 0.6: 0.3: 0.1, respectively.

The projected area of the leaves, O₁₅, O₄₅, O₇₅ of the 3 classes can be calculated using Eq. 4 when the sun elevation angle B_s is higher than the leaf slope while Eq. 5 when it is lower than the leaf slope. O₄₅ and O₇₅ are similarly calculated.

$$O_{15}=\sin {\beta }_s\cos \left(15\right)$$

(4)

$$O_{15}=\frac{\pi }{2}[\sin {\beta }_s\cos \left(15\right)\left.{\sin }^{-1}\left(\frac{\tan {\beta }_s}{\tan \left(15\right)}\right)+{\left({\sin }^2{\beta }_s+{\mathrm{s}\mathrm{i}\mathrm{n}}^2\left(15\right)\right)}^{0.5}\right]$$

(5)

Furthermore in Equation 1, I_sh,i is expressed by:

$$I_{sh,i}=\frac{K_d}{{\left(1-\sigma \right)}^{0.5}}(I_{d,i}+I_{bd,i})$$

(6)

Where K_d is the extinction coefficient for diffuse light. s is the leaf scattering coefficient and was set as 0.3 by considering that the reflectance and transmittance of rice leaf is relatively high. I_d,i is the light intensity of the diffuse light received by the horizontal plane in layer i:

$$I_{d,i}=I_{0d}\mathrm{e}\mathrm{x}\mathrm{p}(-K_d{F}_i)$$

(7)

Whereas I_0d (mmol/m²/s) is the light intensity of the diffuse light received by the horizontal plane above the canopy. I_bd,i, which is the diffuse light derived from the direct light at layer i, is calculated by:

$$I_{bd,i}=I_{d,i}-I_{bb,i}$$

(8)

I_bb,i indicates the light intensity of non-diffusive direct light in layer i. I_b,i and I_bb,i are indicated by Eq. 9, and 10:

$$I_{b,i}=I_{0b}\mathrm{e}\mathrm{x}\mathrm{p}(-K_b{{\left(1-\sigma \right)}^{0.5}F}_i)$$

(9)

$$I_{bb,i}=I_{0b}\mathrm{e}\mathrm{x}\mathrm{p}(-K_b{F}_i)$$

(10)

Where K_b is the extinction coefficient of direct light and is expressed by:

$$K_b=\mathrm{O}/\sin {\beta }_s$$

(11)

Whereas the sun elevation angle b_s can be calculated using the latitude and longitude of the study site. In this model, the attenuation of diffuse light within the canopy follows Lambert-Beer's law (Eq.1) [15] is assumed to decrease.

Photosynthetic rate in each layer: The photosynthetic rate per unit area in layer i (P_N,I, mmol CO₂/m²/s) is expressed by.

$$P_{N,i}=f_{sl,i}{P}_{sl,i}+{(1-f}_{sl,i})P_{sh,i}$$

(12)

The fraction of sun leaves, f_sl,i is expressed as follows.

$$f_{sl,i}=\exp \left(-K_b{F}_i\right)$$

(13)

Where P_sl,i and P_sh.i show each photosynthetic rate (mmol CO₂/m²/s) by sun leaves and sunshade leaves in layer I, respectively. They are calculated using the light response curve of photosynthesis approximated by non-rectangular hyperbola [23].

$$P_{sl,i}=\frac{\phi I_{sl,i}+P_m-\sqrt{{\left(\phi I_{sl,i}+P_m\right)}^2-4\phi I{P}_m\theta }}{2\theta }-r_D$$

(14)

Where f is the initial slope of the light response curve of photosynthesis, P_m is the photosynthetic rate per unit area under saturated light, q is the curvature, and r_D is dark respiration. P_sh,i is also calculated by inputting the light intensity received by shadow leaves: replacing I_sl,i by I_sh,i (Eq.6). By using Eq. 12, P_N,i can be calculated from all derived values$(f_{sl,i},P_{sl,i},P_{sh,i}$). By accumulating P_N in all layers i, the carbon dioxide fixed by photosynthesis for a given LAI can be estimated.

2.2. Experimental site

Two fields with alluvial cray loamy soil were used in the experimental paddy field of Tokyo University of Agriculture and Technology (N35.666, E139.471, 49 m height from sea level) (Figure 1). The area’s climatic conditions are mild and generally warm, with a mean annual temperature and precipitation of 15℃ and 1530 mm, respectively.

The experiment used Koshihikari (Japonica), one of the rice model cultivars and mostly cultivated in various parts of Japan. The experiments were established on three cropping seasons, Field A and B on June 3, 2020, and in Field B on May 20, 2021 and May 30, 2022. Three seedlings per hill were transplanted in every 30 cm (between rows) by 15 cm (between plants) in both Field A and B. Field A was divided into non-fertilized (0N) and fertilized (+N) areas for three replicates (total 6 plots) with plot sizes of 4.2 m by 3.6 m. Meanwhile, Field B has plot size of 7.8 m by 6.0 m with non-fertilized (0N) (Figure 1a).

Field A was used for UAV-based observation in 2020. Field B was used for ground-based observation using the tower with 3 m height (Figure 1b) and by hand from 2020 to 2022. With this set-up, we can make comparisons between fields under the same meteorological conditions within a year (Field A & B in 2020) and yearly difference from 3-year observations (Field B in 2020 to 2022).

2.3. Observation data

The observation data gathered in Field A and B are described as follows (Table 1).

2.3.1. Relative light intensity under rice canopy

The quantum sensors with memory (DEFI2-L, JFE Advantech) were set at five points in Field B (Figure 1b) and four points in plots of Field A (Figure 1a) to measure the transmitted light at the canopy bottom (I_t), and the incident light at the top of the canopy (I_0t) at the time of t. Each sensor was put between 15 cm of each plant.

For calibration, measurements were performed using all sensors under the same sunlight condition for about two weeks before and after the field measurement. The measurement interval was every 1 minute. The data measured from 7:00 to 18:00 was used to determine the daily relative light intensity (I/I₀). In general, the relative light intensity measured under cloudy weather has been used to determine the leaf extinction coefficient [16]. However, when using observation data under limited conditions, it might be difficult to automatically process the selection of usable data, or there would be no usable data depending on the weather conditions. In this study, we calculated the relative light intensity representing as daily value using the 660 values measured at 1-minute intervals under various sky conditions during the daytime using the following equation [24,25].

$$\frac{I}{I_0}=\frac{1}{n}\sum _{n=1}^{660}\frac{I_t}{I_{0t}}$$

(15)

2.3.2. Incident photosynthetic photon flux density

Two photosynthetic photon quantum sensors (LI-190SB, Licor Inc.) were installed at the top of tower in Field B (Figure 1b) to measure both of downward and upward photosynthetic photon flux density (PPFD, mmol/m²/s). The LI-190SB can scan every one second, and the data was saved as 10-minute mean at 10-minute intervals. These PPFDs were used as the input parameters of the incident light intensity at every 10-minute after dividing the direct and diffuse components (as described in 2.5).

2.3.3. Ground-based remote sensing

The single lens reflex camera (D5300; Nikon) with the wide-angle lens (EX/DC; SIGMA) in the waterproof case was installed at nadir direction on the tower of 3 m height in the Field B set-up (Figure 1b). The RGB image was captured by fixed aperture of f:5.6, auto sheeter-speed and auto white-balance. Images were saved as JPEG format at every 10 minutes during daytime in the periods from July 9 to September 9 in 2020 and from May 28 to September 6 in 2021. A total of 36 images per day captured from 9:00 to 15:00 were used to calculate Green Ratio (GR) and vegetation coverage in the region of interest (ROI) about 1.6 m by 1.6 m for each image (Figure 1b). The Green Ration was computed using the equation.

$$GR=\frac{G}{R+G+B}$$

(16)

Where R, G, and B represents the red, green, and blue bands. VC was calculated as the ratio of green leaves area. As for the detection of green leaves, we separated all pixels into the green area and the other area by using CIE L*a*b* which were computed from tristimulus values X, Y, and Z. The green area was segmented using threshold of a* value for each image. The calculated values of GR and VC from 36 images were averaged as daily representative values in Field B.

In addition, three bands of near infrared (NIR), red and green were captured using multi-spectral camera (HAWC; TETRACAM Inc.) at five points from about 1 m height by hand around 9 a.m. once a week rice growth survey. This multi-spectral camera has the incident light sensor. The image taken under the different sunlight condition can be calibrated as reflectance factor using the optional software, PexelWrench (TETRACAM Inc.), then the Normalized Difference Vegetation Index (NDVI) is calculated using the equation.

$$NDVI=\frac{NIR-R}{NIR+R}$$

(17)

We used the average measurements of NDVI in the five points as the weekly representative of the relative light-intensity measurement for Field B.

2.3.4. UAV-based remote sensing

To obtain the spatial and temporal information related to rice growth, the canopy surface model (CSM) and the ortho-mosaic image with spectral bands were generated using the overlapping images taken by UAV. Ground control points (GCPs) were set in Field A using the total station and auto level by conducting traverse surveying and levelling. The coordinates of the GCPs were referred to the Japan Geodetic Datum 2011/Plane Rectangular Coordinate System zone 9 as a map projection. These were fixed at the six points around the edge of the paddy field (Figure 1). A set of overlapped images covering the whole field was taken using an UAV (Inspire 2; DJI) equipped with an RGB camera (Zenmuse X4S; DJI) and a multispectral camera (RedEdge-MX; MicaSense). The flight altitude was fixed at 30 m above the rice canopy with a forward and lateral overlap rate of 85%.

The flight of the UAV was conducted between from 9:00 to 10:30 am, when the sun elevation angle was relatively high (around 50-70 degrees), and the sky was clear. Since the multispectral camera is equipped with an incident spectral-light sensor, the correction of reflectance factors can be performed in post-processing by taking a grayscale board for correction before the fight. Two kinds of aerial images with RGB and multi-spectral bands were taken every week from the beginning of June to the end of August in 2020 (10 - 13 flights used).

Generation of canopy surface model: The Metashape Professional (ver. 1.5.1, Agisoft), a Structure from Motion (SfM) software, was used to generate time-series digital surface model (DSM) by processing the UAV-acquired RGB images on days of 10 flights. First, tie points were automatically identified from the overlapped aerial images; then the tie points were used to calibrate the camera parameters such as the focal length of the lens, principal point positioning, and radial and tangential distortions. The parameters of external orientation (camera position and tilting angle) were estimated using the detected tie points and four installed GCPs, then, a 3D model was generated. This processing was performed to achieve a GCP accuracy within 1 pixel. A DSM with a spatial resolution of approximately 9 mm/pixel was developed. Using the DSM, CSM was calculated from the distance between the DSM of each observation day (DSM_n) and the first DSM after transplanting (DSM_1st, defined as the reference plane). The value of the CSM was defined as the canopy height (CH). This relation is mathematically expressed as

$${CSM}_n={DSM}_n-{DSM}_{1st},$$

(18)

where n represents the observation dates [26].

Generation of ortho-mosaic images: Time-series ortho-mosaic images with multispectral bands were generated from each set of multispectral images taken by total 10 flights at Field A using the same SfM software, Metashape. Each ortho-mosaic image was generated after creating DSM within 1 pixel (14mm) error using the same GCPs.

Calculation of canopy height (CH), NDVI and green ratio (GR): After generating CSMs and ortho-mosaic images, we created the polygons with the minimum size of 30 cm × 15 cm rectangle, which correspond to one hill of rice plant. Then, the mean value of the CH, NDVI and GR were calculated in each polygon.

2.3.5. Rice growth survey

In Field A, plant length (PL, m), LAI and above ground biomass (AGB, g/m²) were measured weekly for four hills in each plot as shown in gray of Figure 1b. LAI was measured using automatic area measurer (AAM-9; Hayashi Denko,). In Field B, PL was measured at five points for two hills, while AGB was measured at each week by two hills sampling from outside the ROI. The number of tillers for the target hills was also counted. The measuring and sampling areas in Fiend A and B are shown in Figure 1a-b.

The actual measured LAI was used to determine the extinction coefficient K_d by using relative light intensity (I/I₀), while the AGB was used to compare the daily growth production and accumulating biomass estimated by the canopy photosynthesis model (2.1.2).

2.4. LAI estimations

LAI can be calculated using relative light intensity (I/I₀) and K_d from Eq.19.

$$LAI=-\frac{1}{K_d}\mathrm{l}\mathrm{n}\left(\frac{I}{I_0}\right)$$

(19)

In this study, the estimation formulas of relative light intensity (I/I₀) were derived by comparing the variables of daily VC, daily and weekly GR, weekly NDVI, and weekly CH from ground- and UAV-based remote sensing.

2.5. Direct and diffuse components divided from incident global PPFD

According to the canopy photosynthetic model used in this study, the input light intensity needs to divide the direct and diffuse components. In the previous studies, it proposed to estimate the diffuse ratio (DR) proposed by using the clearness index (CI) which is the ratio of the solar radiation at the top of atmosphere to the global solar radiation measured at the ground level [27]. The direct and diffuse radiations basically have not observed as the public meteorological data, therefore, in the case of using the direct or diffuse components of the solar radiation or PPFD, it is necessary to estimate them using the global solar radiation e.g., [28,29,30], which commonly observed worldwide [31].

First, the formula was derived to estimate the diffuse ratio from CI by using the global and diffuse solar radiation observed at the campus of Tokyo University of Agriculture and Technology which is approximately 2 kilometers away from the experimental site.

$$\begin{gathered} DR=1(CI\le 0.250) \\ DR=7.511{CI}^3-10.894{CI}^3+3.112CI+0.782\left(0.250<CI\le 0.775\right) \\ DR=0.145(0.775<CI) \end{gathered}$$

(20)

Then, we integrated the amount of the extraterrestrial spectral irradiance (W/m²/mm) in the visible wavelength (0.4-0.7 mm) [32] by converting to the unit of photon flux density (mmol/m²/s) and obtained the extraterrestrial PPFD as the constant 2405 mmol/m ²/s (PPFD₀) [33]. The clearness index of PPFD (CI_ppfd) was calculated using the downward PPFD (PPFD_g) measured at Field B by

$${CI}_{ppfd}=\frac{{PPFD}_g}{{PPFD}_{toa}}$$

(21)

$${PPFD}_{toa}={\mathrm{P}\mathrm{P}\mathrm{F}\mathrm{D}}_0{\left(\frac{r_0}{r}\right)}^2\sin {\beta }_s,$$

(22)

where (r₀/r)² is the correction of the inverse square between r₀ (the mean Earth-Sun distance) and r (the Earth-Sun distance) on the observation day. Here, CI_ppfd was substituted for CI in Eq.20 to obtain DR.

After calculating DR, the downward PPFD was divided to diffuse PPFD and direct PPFD. In addition, the upward PPFD was multiplied by the ratio of diffuse and direct components. Each value subtracted from diffuse PPFD and direct PPFD was used as the diffuse and direct light intensity incident on the canopy (I_0d and I_0b) with 10 min-interval.

2.6. Estimation of daily growth biomass

Figure 2 shows the procedures to estimate daily growth biomass based on the canopy photosynthesis model using daily or weekly LAI estimated by ground- and UAV-based remote sensing (2.4) and incident light intensity with 10 min-interval (2.5).

As for the others input parameters, sun elevation angle b_s was calculated at 10 min-interval using the day of year and hour: minute at Japan standard time and the latitude and longitude at experimental field (N35.666, E139.471). Furthermore, the extinction coefficient K_d was determined by using the daily relative light intensity (I/I₀) (2.3.1) and actual measured LAI at Field A (2.3.5). As for the parameters of light response curve of photosynthesis, P_m: 21, f: 0.08, q: 0.8, and r_D: 0.5 were set by references [34,35] that describe about the results by measuring photosynthetic rate for Koshihikari.

As output, the photosynthetic rate (mmol CO₂/m²/s) at every 10min-interval was accumulated for one day (mol CO₂/m²/d), then converted to daily biomass (CH₂O g/m²/d) and obtained accumulated biomass (CH₂O g/m²). The photosynthesis rate, daily growth biomass and the accumulated biomass were compared with the incident light intensity and the measured AGB (g/m²) and discussed about these relationships.

3. Results and discussions

3.1. Rice growth survey

The results of the weekly rice growth survey, PL, AGB and LAI in Field A in 2020, and PL and AGB in Field B for 3 years are shown in Figure 3. There are four kinds of sampling areas, two non-fertilized plots (0N, #82, and #83) and two fertilized plots (+N, #84, and #85) in Field A (Figure 1b). Results show growth differences between 0N and +N starting from 35-40 days after transplanting (DAT) (Figure 3a). In Figure 3c, LAI declined after its peak around DAT60 in all the plots (fertilized and non-fertilized plots). This shows the same trends as the study conducted in 2018 and 2019 by Peprah, et al. [26]. In Field B, a total of 10 plants were sampled in the area where the 5 photon sensors were placed. The mean PL and AGB of the sample plants showed an increasing trend from 2020 to 2022 (Figure 3d-3f).

3.2. Weekly change of canopy height calculated from Canopy Surface Models

Figure 4a shows the results of canopy height calculated from aerial drone images taken in Field A. Differences in CH began to appear between the non-fertilized and fertilized plots from around 40 DAT which coincides with the actual PL measurement (Figure 3a). Although it was not obvious from the seasonal changes in PL (Figure 4a), CH begins to decrease slightly from around 75 DAT (heading stage). This can be explained by the method used to obtain the PL and CH measurements. PL is measured from the soil surface to the highest tip of the rice plant, whereas CH is the height in its natural state. Because the grain begins to fill and mature, the rice plants tend to bend toward harvesting.

The relationship between the measured PL and CH (Figure 4b), is expressed in a linear equation (CH = 0.993*PL - 0.286) with a slope nearly equal to 1.0. This indicates that a highly accurate canopy surface model (CSM) can be obtained using aerial images taken by UAV.

3.3. Daily change of relative light intensity (I/I₀) and the extinction coefficient (K_d)

The relative light intensity (I/I₀) is the light transmittance of leaves under the canopy and a parameter to estimate LAI (Eq.1) by non-destructive method. The daily average of relative light intensity (Eq.15) in the Field A in 2020 and Field B in 2020 and 2021 are shown in Figure 5. There are the clear differences between non-fertilized and fertilized plots in Figure 5a. The patterns of non-fertilized plots (#82 and #83) in Field B show different trends. Even on the same transplanting date (3rd June), differences in relative light intensity can be seen in different fields due to differences in growth during the tillering stage. In particular, the transmittance decreases from 3 weeks to 9 weeks after transplanting depending on whether the field is fertilized or not fertilized.

During the beginning of I/I₀ observation in 2021, DAT from 12 (30^th May) to 22 (10^th June) was linearly interpolated because the sensor was submerged in water and the data could not be used (Figure 5c). The error bar shows the standard deviation of 5 points of the quantum sensors measurement in Field B (Figure 5b-c). Comparing Field B in 2020 and 2021, variations in the 5 measurement points are seen especially in 2020. From around 70-75 DAT (heading stage), I/I₀ have been stable with less than 0.1 relative light intensity (Figure 5b, and 5c).

The extinction coefficient (K_d) was derived from the relationship between relative light intensity and measured LAI using data observed from early to late growth (Figure 5d). In this study, K_d value is 0.562, which was found to be very close to the value of existing studies for the extinction coefficient of rice e.g., [36,37]. Strictly speaking, it is known that the extinction coefficient differs depending on the growth stage [36]. However, in the canopy photosynthesis model used in this study, described in 2.1.2, the extinction coefficient of the leaves is assumed to be constant.

3.4. Relations of relative light intensity with the parameters by ground- and UAV-based observations

The relationship between relative light intensity (I/I0) and daily or weekly observation data (VC, GR, NDVI and CH) obtained using ground- and UAV-based remote sensing methods was determined to obtain the optimal parameters for estimating LAI (Figure 2).

3.4.1. Daily vegetation cover (VC) and green ratio (GR)

Figure 6 shows the negative relationships between daily VC, GR and I/I₀ at Field B in 2020 and 2021. At the early stage of growth (I/I₀: 0.95 to 0.8), the VC and GR increased quickly in 2020. Moreover, when the I/I₀ becomes 0.8 or less, VC (in 2020 and 2021) and GR (in 2020 only) increased with an almost constant slope. Comparing the two years, there are no major differences in the distribution for VC, however, the distributions of GR are significantly different between 2020 and 2021. It seems that other factors may be influencing GR [38,39,40]. When using VC as an estimation parameter, the image resolution must be high to calculate the percentage of pixels identified as leaves. Therefore, it is considered difficult to calculate VC using UAV images compared to ground-based images.

3.4.2. Weekly NDVI and CH at Field B in 2020 to 2022

Figure 7 shows the relationship of I/I₀ with NDVI and CH obtained at every week for three years from 2020 to 2022 during growth surveys in Field B. The CH was calculated from the measured PL in the growth survey using an estimation formula (CH = 0.993*PL - 0.286).

NDVI is a linear equation while CH is an exponential function, and a well-fitting regression equation has been derived. Comparing among years, although there is a slight difference in the intercept for NDVI, the slope is almost the same. This is because the multispectral camera (HARK) used is equipped with an incident light sensor. As a result, internal processing can correct the differences in solar radiation condition on the day and time the photo was taken. This result is significantly different from the GR result (Figure 6b) using the uncalibrated camera images [38,41].

Regarding CH, I/I₀ =1 at the top of the canopy can be theoretically regarded as the intercept, as defined by Lambert-Beer's law. There is almost no difference in the coefficient related to CH (2.86, 2.88 and 3.02), corresponding to the optical thickness within the rice canopy, over the three years of rice cultivation using the same variety and transplanting density. Since the physical quantity of height is used as a variable, there is a high possibility that plant height in conventional growth surveys can be used.

3.4.3. Weekly NDVI and CH by UAV-based observation at field A in 2020

The relationship between NDVI, GR, CH and relative light intensity based on UAV observation in Field A is shown in Figure 8. Slight differences in the coefficients can be seen on plots with or without fertilizer application. The difference in growth amount due to different fertilization is also shown in Figure 3. CH is an exponential function, and NDVI and GR show the linear regressions (Figure 7). The coefficients of three types of regressions are less in non-fertilized (0N) plots, indicating that even with the same NDVI, GR, and CH values, non-fertilization has high relative light intensity. In instance, the transmittance is different in the same leaf area and different leaf characteristics (ex. thickness). Although it also depends on the variety, the coefficient may be smaller in non-fertilized plots because the rice plants grow taller but with fewer leaves. Regarding such differences in the coefficients of regression equations, further investigation is required based on differences in varieties and fertilizer amounts.

3.4.4. Comparison with field A and field B

Non-fertilized plots in Field A and B were compared in 2020 observations (Figure 9). In terms of NDVI, there is a slight difference in the slope and a large difference in the intercept which largely due to the differences in the cameras used and spatial resolution (distance to the rice plants). In contrast, there is no difference in CH between Fields A and B, the coefficients are the same even in different fields when the fertilization conditions are the same.

Generation of DSM using overlapped images taken by UAV is relatively unaffected by the weather or time of day, therefore, the possibility of using time-series CSM and obtain spatially continuous information in the entire field is high. Furthermore, it is an effective method to observe the canopy height (corresponding to rice growth) from RGB images acquired by UAV, without using a calibrated multispectral camera [26]. However, satellite images (e.g., vegetation indexes etc.) can be used to expand to a wider regional scale [6,42,43], further knowledge of relationship between canopy height and satellite data is required [44,45].

3.5. Daily biomass estimation at the field scale

The daily LAI was estimated by using the observation results of daily relative light intensity (I/I ₀) with K_d as 0.562 (Eq.19). Then, the daily biomass (CH₂O g/m²/d) at Field B in 2020 and 2021 was assessed by using the incident diffuse and direct light intensity (I_0d and I_0b, mmol/m²/s) with 10 min-interval divided from the observation data of upward and downward PPFD (Eq.1-14). Figure 10 shows the sampled AGB, estimated daily biomass and accumulated biomass (Figure 10a-b), and the daily amounts of direct and diffuse PPFD (mol/m²/d) during the growing period (Figure 10c-d). Based on the results, the daily biomass amount varies depending on the light intensity. In addition, as LAI increases with growth, the daily biomass also increases before peaking at the early stage of heading (DAT60-70) in 2020 and 2021. After heading, the light intensity fluctuates depending on the daily weather. Due to unfavorable weather for one week around DAT 85-90 in 2021, the influence of daytime light intensity is found to decrease in daily biomass and accumulated biomass. Furthermore, the time-series changes in accumulated biomass and AGB measured by weekly growth surveys show closer values in 2020 and 2021.

To estimate daily biomass in the entire Field A, LAI with the mesh data of 15 cm by 30 cm size were estimated weekly by using the CSM generated by UAV-images from Eq.23.

$${LAI}_{ch}=-\frac{1}{K_d}\mathrm{l}\mathrm{n}\left(e^{-\tau *CH}\right)$$

(23)

Where LAI_ch is the LAI estimated from CH, and t is the optical thickness within the rice canopy. t was determined as 2.95 from the relationships between I/I₀ and CH (Figure 7c, Figure 8c, Figure 9b). The weekly estimated LAI_ch was interpolated to the daily LAI_ch by spline, then, the accumulated biomasses day by day from 14 (17^th June) to 97 (8^th Sep.) DAT were estimated spatially in Field A (Figure 11).

Daily changes of the biomass spatially, and the differences between non-fertilized and fertilized plots can be observed. As shown in the estimation results of time-series accumulated biomass at the quantum sensor points in plots #84 and #83 (Figure 12), daily growth and its accumulated amount fit well with AGB measured at the sampling area.

With the spread of UAVs, applied research on precision farming using UAV images is increasing. While there are many research examples using machine learning, there is a possibility that it cannot be applied to observed data due to different observation conditions (ex. differences in field or year) depending on the camera used. Therefore, engineers and researchers are required to understand the essence of observational data obtained through remote sensing, and then use their judgment to determine which images can be used as training data.

In this study, a method to estimate the relative light intensity under the rice canopy using a simple regression model based on limited observation data is proposed. From these results, it is possible to understand the daily growth amount by using LAI time series data estimated from the relative light intensity under the rice canopy based on the remotely sensed observation and the direct and diffuse PPFD data used for photosynthesis.

In recent years, the use of spatiotemporal data from ground, UAVs, and satellites remote sensing is expected to be utilized in smart agriculture as it has been used in an integrated manner with numerical models based on observational data [6,46]. However, since PPFD is not a frequent meteorological observation item, research on calculating direct and diffuse PPFD, utilizing global solar radiation, ground observation, and meteorological satellite data is required.

4. Conclusions

This study examined an efficient and nondestructive method for estimating LAI for precision farming by remote sensing techniques. The method was applied to a canopy photosynthesis model that calculates the photosynthetic rate considering temporal changes in sun incident light. The spatial estimation of daily biomass and accumulated biomass in the paddy field was also conducted.

The results show that daily rice plant growth can be quantitatively estimated, including spatial variation, by utilizing estimated LAI. It was estimated using weekly ground-based or UAV-based data and a canopy photosynthesis model with high temporal resolution PPFD 10-minute observations as input data. Conventional simulations based on numerical models have been limited to one-dimensional predictions representative of a farm field because spatio-temporal observation data were unavailable. With the recent global promotion of precision farming and the widespread use of UAVs, the possibility of near-real-time estimation and prediction has been demonstrated through the integrated use of observation data based on remote sensing. It is concluded that the near-real-time estimation of rice biomass by data assimilation with numerical models based on the field observation data is an effective method for the cultivation management of major crops.