Exploring a best vegetation index to explain the seasonal variation of a forest photosynthesis using a hyper-spectral camera equipped with liquid crystal tunable filter

Liquid crystal tunable filter ( LCTF ) can change the transmissible wavelength by changing the applied voltage to the filter, which enables the drastic increase in the observable wavelength resolution in a small size system and is considered to be a powerful tool for the spectral earth observation from flying units or microsatellites. However, there is limited knowledge about its season-long application for the vegetation monitoring and the prediction of the ecosystem photosynthetic capacity. We compared the seasonal variation of spectral reflectance obtained by a LCTF camera with that obtained by a popular spectral radiometer in a cool-temperate young larch plantation in northern Hokkaido, Japan. Then we tried to find the best normalized difference spectral index ( NDSI ) to explain the seasonal variation of the ecosystem photosynthetic capacity using all pairs of two reflectances observed in the range of wavelength between 500 and 770 nm with 10-nm intervals ( 28 wavelength bands ) by the LCTF. The best NDSI among all combinations ( 28 × 27 ) of two reflectances was NDSI[770, 720] for the maximum gross primary production at light saturation and NDSI[530, 600] for the initial slope of the light-response curve, which reflect the red edge shift owing to the change in the chlorophyll content and relative strength of the light absorbance in the visible red wavelength region compared with that in the green wavelength region, respectively. Predicted daily gross primary production of the plantation using these NDSI agreed well with the observed values. NDSI[530, 600] was better to distinguish each vegetation type of the studied plantation.


Introduction
Terrestrial gross primary production (GPP) is the largest global CO 2 flux (Beer et al., 2010), and several vegetation indices (VIs) have been proposed to predict its temporal and spatial variations globally using remotely sensed surface reflectance (Turner et al., 2003(Turner et al., , 2006Running et al., 2004;Xiao et al., 2004;Sims et al., 2006b;Inoue et al., 2008). Several major VIs are expressed as the normalized difference of two spectral reflectances (Normalized Difference Spectral Index: NDSI), such as Normalized Difference Vegetation Index (NDVI), Green-Red Vegetation Index (GRVI), Photochemical Reflectance Index (PRI), and Chlorophyll Index (CI) ( Table 1), and have been used for estimation of ecosystem GPP or light use efficiency and their universal representativeness have been compared and discussed (Nakaji et al., 2008;Ide et al., 2010).
Recently microsatellites are expected as a powerful tool for earth observation because of their low cost for the development and operation ( Kurihara et al., 2018). A microsatellite RISING-2 had been developed for this purpose by a cooperative project of Tohoku University and Hokkaido University, and launched in 2014. This microsatellite first contains high precision telescope with liquid crystal tunable filter (LCTF), which enable multispectral earth surface monitoring with 5 m spatial resolution at ground surface ( Kurihara et al., 2018). The LCTF is a type of optical band pass filter whose center wavelength is electrically tunable and has been used for many applications such as, agriculture, healthcare, archaeology, and art (Gat, 2000). Shibayama et al. (2009) succeed in predicting leaf greenness and 1000-grain weight of ripening rice plots using the reflectance obtained by a LCTF imager. Ishida et al. (2018) and Kurihara et al. (2020) apply the LCTF camera to UAV (unmanned aerial vehicle)-based hyperspectral imaging for the surface vegetation classification. However, there is limited knowledge about its season-long application for the vegetation monitoring and the prediction of the ecosystem photosynthetic capacity using this imager.
The aim of this study is to confirm the performance of the LCTF camera by comparing the observed seasonal variation of spectral reflectance with that measured by a popular spectral radiometer in a young larch plantation in northern Hokkaido, and to find the best NDSI to explain the seasonal variation of the ecosystem photosynthetic capacity. We calculated seasonal variation of all NDSIs using every combination of two spectral reflectance having the waveband width of 10 nm and compared them with that of the canopy photosynthetic characteristics obtained from the CO 2 flux observation. We also confirmed that such NDSIs can distinguish each vegetation types in the plantation.

Study site
The study site is a young larch (Larix gmelinii (Rupr.) Kuzen. var. japonica (Maxim. ex Regel) Pilg. L. kaempferi (Lamb.) Carrière) plantation created after clearcutting of a mixed forest in January 2003 in the Teshio Experimental Forest of Hokkaido University (45 03'N, 142 06'E, 66 m a.s.l.) (Fig. 1). The larch plantation is mixed with naturally regenerated deciduous birch (Betula ermanii Cham. and B. platyphylla var. japonica (Miq.) Hara.) and oak (Quercus crispula Blume) trees, and the forest floor is covered with dense evergreen dwarf bamboos (Sasa senanensis Rehd. and Sasa kurilensis (Rupr.) Makino et Shibata). Maximum tree height was ca. 8 m but having large spatial heterogeneity. The plant area index (PAI; which includes the shade by stems, branches, and culms in addition by leaves) values for the canopy trees and the Sasa bamboos, measured using an LAI-2000 plant-canopy analyzer (Li-Cor, Lincoln, NE, USA), were 1.7 0.2 and 6.8 0.59 m 2 m -2 , respectively, at this parameter's seasonal maximum in 2016, showing major contribution of Sasa to the total plant area in this ecosystem.

Spectral radiation measurements
Custom-made liquid crystal tunable filter (LCTF) camera (Genesia Co., Tokyo, Japan) was set at a 16 m in height to observe southward vegetation canopy surrounding CO 2 flux monitoring tower with an angle of depletion at 30 (Figs. 1 and 2). This camera can capture 2 D spectral reflected radiation from 460 to 780 nm at 10-nm intervals (i.e. 33 wavelength bands) within 1 s for each wavelength band. The full width at half maximum is 6 nm at 460 nm and 21 nm at 780 nm and linearly increases with the wavelength. The angle of view and pixel number are 72 54 and 494 656 (Fig. 2). The 2D reflected radiation was obtained for each wavelength band every 15 min during 12 : 00 to 14 : 00 every day from 21 April to 10 November (snow-free period) in 2016. We regulated the gain of LCTF to get the digital number to prevent saturation at high radiation. Digital numbers (256 gradation, 10 bits) obtained by the LCTF camera were converted into the spectral radiation ( Kurihara et al., 2018), here the spatial heterogeneity in the sensitivity within a picture was calibrated for each pixel using coefficients obtained in a uniform light source (HELIOS USLR-D12LNMNN; Labsphere, Inc., North Sutton, NH, USA). Downward (solar) and upward (reflection from the canopy) spectral radiation was measured with two hemispherical spectroradiometers (MS700; EKO Instruments, Tokyo, Japan) at 30 m in height of the same tower. The sensor for measuring the downward (solar) spectral radiation was mounted on the top of the tower, and that for upward (reflection from the canopy) spectral radiation was mounted on the lower side at the end of the 1.6-m long horizontal boom which jutted out from the tower. The lower side sensor was equipped with a partial cover to exclude high reflection from the mounted aluminum tower ( Ide et al., 2016). The underestimation owing to the shade by the cover was corrected based on the ratio of the shade area of the sensor. The spectral solar and reflected radiation was measured every minute from 12 : 00 to 14 : 00 every day for 256 wavelength bands ranging from 350 nm to 1050 nm at 3.3-nm intervals, and the spectral reflected radiation was divided by the solar radiation for each wavelength band to obtain the whole ecosystem spectral reflectance.
The spectral reflected radiation in the region of interest (ROI) of 2D data from LCTF camera was averaged to represent each vegetation type (three replicates of deciduous larch, deciduous birch, and evergreen Sasa bamboo) or the whole ecosystem to compare with the spectral reflectance of MS700 or the ecosystem photosynthetic capacities obtained from the CO 2 flux observation. The spectral reflectance for each ROI was derived from the reflected radiation averaged for each ROI divided by the spectral solar radiation obtained by MS700 every 15 min during 12 : 00 to 14 : 00 and every 10 nm from 460 nm to 780 nm in the wavelength, then the daily value was obtained by averaging all the 15-min data for each wavelength band. Three spectral solar radiation data obtained by MS700 at 3.3-nm intervals in the wavelength bands were averaged to adjust to the wavelength band resolution of LCTF to calculate the reflectance. Data spikes in the spectral solar radiation caused by the cloud movements were removed before this averaging, in order to keep the coefficient of variation (standard variation divided by the average) of consecutive 5 data below 0.1. Spectral reflectance calculated for each of the three ROIs for each vegetation type was averaged to show the representative seasonal variation. We removed the daily spectral reflectance obtained during rainfall event (rainfall amount during 12 : 00 -14 : 00 > 0.5 mm).  Gitelson and Merzlyak (1994) R λ is the reflectance at wavelengths around λ nm, and subscripts nir, red, and grn indicate the reflectance at wavelength around near-infrared (740-900 nm), red (630-690 nm), and green (531-570 nm) bands, respectively.

CO 2 flux and micrometeorology measurements
During the measurement of spectral reflectance, CO 2 flux over the study site was evaluated using a closed-path eddy covariance system mounted at 7.4 m in height of another tower (Fig. 1), placed inside the LCTF camera imaging range. A sonic anemometer (DA600-3TV; Kaijo, Tokyo, Japan) and an infrared gas analyzer (IRGA, LI-7000; Li-Cor) were used to evaluate the fluxes. In this system, the separation distance between the  anemometer and the air intake was 20 cm, and the air was drawn into the IRGA at a flow rate of 10 L min -1 through a 1.0-µm filter and a 15-m Teflon tube. The IRGA reference cell was fed with N 2 gas (0 µmol mol -1 for the CO 2 concentration) at a flow rate of 10 to 20 mL min -1 . The CO 2 concentration was calibrated daily using two standard gases (320 and 420 ppm) between 23 : 30 and 24 : 00 h. The data were sampled at a frequency of 10 Hz with a digitizing data recorder (CR5000; Campbell Scientific, Logan, UT, USA) after low-pass filtering (with a cut-off frequency of 5 Hz). Takagi  Photosynthetically active radiation (PAR) was observed by a quantum sensor (SQ110; Apogee Instruments, Logan, UT, USA) at a height of 32 m, and air temperature (HMP45A; Vaisala, Helsinki, Finland) was observed at 3 m in height. These data were sampled every 5 s, and data were stored as 0.5-h means using two dataloggers (CR23Xs; Campbell Scientific).

Evaluation of CO 2 flux, gross primary production and photosynthetic parameters
Following the same procedure used in the previous studies at the site ( Takagi et al., 2009;Aguilos et al., 2014), we determined the daily sonic rotation angle for use in the planar fit rotation ( Wilczak et al., 2001) using the 30-min mean wind speeds in a 15-day moving window. We determined the fixed value of the sonic-tube lag time for CO 2 monthly by averaging the lag times obtained at 30-min intervals under turbulent conditions. Using these angles and lag times, half-hourly CO 2 fluxes (F c , µmol m -2 s -1 ) were calculated. Crosswind speed and water vapor concentration effects on the sensible heat flux were corrected using the methods of Kaimal and Gaynor (1991) and Hignett (1992), respectively, then we corrected the effect of air density fluctuations on the F c values ( Leuning and King, 1992). High-frequency losses for the sonic sensor's span and sensor separation were corrected using transfer functions related to the sources of signal damping (Moore, 1986), and losses for tube attenuation were corrected following the method of Kowalski et al. (2003), where correction factor was expressed as a linear function of the horizontal wind speed, and the coefficients of the function were fixed throughout the study period.
In checking the raw F c data, we applied the same procedure as Takagi et al. (2009), which includes quality control program developed by Vickers and Mahrt (1997) and Mano et al. (2007), instationarity ratio and integral turbulence characteristic tests (Foken and Wichura, 1996), and footprint evaluation of the observed F c using the model developed by Kormann and Meixner (2001) to remove the effect of fluxes from outside the plantation. Then u * (friction velocity) filtering was applied to the remaining F c data, setting the threshold value at 0.1 m s -1 . F c data passing these quality control procedures were used as the net ecosystem CO 2 exchange (NEE, µmol m -2 s -1 ).
Two photosynthetic parameters, (i.e. the maximum gross primary production (GPP) at light saturation (A max : μmol m -2 s -1 ) and initial slope of the light-response curve ( φ: mol mol -1 )), and daytime respiration (R d : μmol m -2 s -1 ) as the y-intercept of the light-response curve in equation (1), were determined every day by the least-squares method using daytime (PAR ≥ 10 µmol m -2 s -1 ) 30-min NEE and PAR data following Ueyama et al. (2012) and Takagi et al. (2015), but within a 7-day moving window: where θ is the convexity of the light-response curve and we used a fixed value (0.9) to apply the least-squares method. Seasonal variation of the two photosynthetic parameters, daily A max and φ are compared with that of NDSIs, however because of the limitation in the least-squares method (regression curve is drawn through the intermediate point of the high-low range of the NEE values for each PAR class), estimated A max from the obtained regression curve usually underestimate the observed maximum GPP. Thus, we used A max and φ of the light response curve of GPP passing through the upper limit of the 95 confidence intervals of the regression to compare with the NDSI values.
GPP (μmol m -2 s -1 ) was determined from observed NEE and micrometeorology data to validate our GPP simulation using the NDSI. Before the determination, we filled the gaps in NEE data mainly by using lookup tables (Falge et al., 2001) created every 30 days for each air temperature (2ºC interval) PAR (100 μmol m -2 s -1 ) class. A few data gaps not filled by the lookup tables were filled by the mean diurnal variation (MDV) approach ( Falge et al., 2001), in which missing NEE was replaced by the mean for that time based on the adjacent 9 days. Night-time NEE (PAR < 1 μmol m -2 s -1 ) was assumed to be equivalent to the ecosystem respiration (RE; μmol m -2 s -1 ), and the 30-min values were compiled for air temperature classes (2ºC interval), and the average NEE for each temperature class was related to the observed air temperature (T: K) by using an Arrhenius-type equation (Lloyd and Taylor, 1994): where R is the universal gas constant (8.134 JK -1 mol -1 ), and the constants R ref and E a are the respiration rate (μmol m -2 s -1 ) at the reference temperature (T ref = 283.16 K) and the apparent temperature sensitivity (or activation energy) (J mol -1 ), respectively, and fixed throughout the study period. Daytime RE was estimated every 30 min by using equation (2) and air temperature, and GPP was estimated as RE-NEE. Daily GPP was calculated by summing all 30-min GPP values for each day. j )) was calculated using all (33 32) pairs of two wavelength bands obtained by LCTF camera (only for ROI for whole ecosystem) every day from May 5 to October 30 (179 days), where R i and R j are the spectral reflectance of 10-nm width wavelength bands from 460 nm to 780 nm, but the NDSI using the spectral reflectance within the range from 500 nm to 770 nm (28 27 pairs) was used to obtain the coefficient of determination (r 2 ) between the seasonal variation of daily NDSI and that of two photosynthetic parameters (daily A max and φ), excluding uncertain spectral reflectance values discussed in the results and discussion section. For this comparison, we used daily NDSIs averaged from 3 days before to 3 days after that day to adjust to the averaging period for the photosynthetic parameters. We obtained two linear equations to predict daily A max and φ, using the best NDSI [ i, j ] as, daily A max = a daily NDSI [ i, j ] +b, and daily φ= c daily NDSI [ i, j ] +d, where a, b, c, d are constants for the linear regression. GPP (μmol m -2 s -1 ) was simulated every 30 min using the following equation, here simulated daily A max and φ was fixed through each day, and 30-min averages of observed PAR was used for the simulation.

Determining the best-fit vegetation index and GPP simulation
For the comparison with the observed GPP, daily GPP on a certain day was simulated by using the PAR on that day, and the two linear relationships to estimate daily A max and φ from NDSIs determined using all daily data during the study period except the 7days used to determine the daily A max , φ and NDSIs of that day.

Results and discussion
Spectral reflectance obtained by the LCTF camera tended to have less smooth spectral variation than that of MS700, additionally it showed decrease at 780 nm, and had slight protrusion around 480 nm (Fig. 3). This would mainly be caused by the smaller ROI and coarser spectral resolution of the LCTF. Additionally, we regulated the gain of LCTF to get the digital number to prevent saturation at high radiation. This might cause the black defects (thus, small overestimation) of radiation especially at < 500 nm in the wavelength, and careful usage would be required at this range. Despite these disadvantages, seasonal variation of the spectral reflectance obtained by the LCFT camera showed similar characteristics with that obtained by MS700. Coefficients of determination (r 2 ) between the seasonal variation of NDSI [ i, j ] of the LCTF camera and that Fig. 4. Matrix of the coefficient of determination (r 2 ) between the seasonal variation of daily A max or daily φ and that of daily where R i and R j is the spectral reflectance from 500 nm to 770 nm) obtained by LCTF camera with 10-nm for the waveband width.
The best NDSI to explain the seasonal variation of A max was the pair of reflectances at 770 nm and 720 nm, and the linear equation, daily A max = 85.76 daily NDSI [ 770, 720 ] + 2.00 (r 2 = 0.46) was obtained (Fig. 5). The best NDSI for φ was NDSI [ 530,600 ] , and the linear equation, daily φ= 0.2043 daily NDSI [ 530, 600 ] + 0.0295 (r 2 = 0.60) was obtained (Fig. 5). 2D images of the NDSI [ 770,720 ] and NDSI [ 530,600 ] captured the spatial distribution of the seasonal variation, and NDSI [ 530,600 ] seemed to distinguish each vegetation type more clearly (Fig. 6). Area-averaged NDSIs shows that the NDSI [ 530, 600 ] is better to distinguish each vegetation type (Fig. 7), although the seasonal variation of both NDSIs shows that birch and larch trees had higher NDSI from the leaf expansion in late May, then kept the highest values among the vegetations until late August. Evergreen Sasa expands current year's leaf from Mid-July, then the highest value was obtained from late August to Mid-September, hereafter Sasa NDSI became the highest, owing to the decreasing NDSI of deciduous birch and larch trees by the leaf yellowing. Whole ecosystem NDSI had intermediate value among the three species during the period.
Simulated GPP agreed well with the observed GPP (Fig. 8); accordingly NDSI [ 770,720 ] and NDSI [ 530,600 ] can be a good index to explain the seasonal variation of ecosystem photosynthesis. The slight gentle slope (0.87) of the regression equation was owing to the overestimation in the smaller GPP range, and this became 0.96 if we fixed the y-intercept at 0. If we simply use the regression parameters (A max and φ) without considering the upper boundary of the light response curve of GPP, this equation becomes y = 0.79 x + 0.49 (r 2 = 0.91, RMSE = 2.09), where y and x are simulated and observed daily GPP (gC m -2 d -1 ) respectively. Accordingly, our approach was better than standard regression to simulate observed GPP even without applying any stress functions for GPP estimation.
The range between 720 nm and 770 nm includes the red edge shift (drastic increase in the spectral reflectance at near infrared bands) around the wavelength of 720 nm, where the change in chlorophyll content is reflected (Sims et al., 2006a). This NDSI has similar features with Chlorophyll Index, CI ( Table 1) and Canopy Chlorophyll Index (CCI, D 720 /D 710 ) (Sims et al., 2006a), where D λ is the first derivative of reflectance at wavelengths around λ nm. Accordingly the seasonal variation in the maximum photosynthetic capacity, A max , reflected that in NDSI [ 770,720 ] . On the other hand, the best VI for φ was NDSI [ 530, 600 ] , which rather reflects the relative difference in the spectral reflectance between green (around 531 -570 nm) and red (around 630 -690 nm) indicating the relative strength of the light absorbance in the visible red wavelength region compared with that in the green wavelength region, and is similar with GRVI (Table 1). Inoue et al. (2008) applied a similar approach to explore the best NDSI to explain the seasonal variation of photosynthetic capacities or light use efficiencies for an irrigated rice field. Their proposed NDSI [ 750, 761 ] (r 2 = 0.67) for A max is similar with our result, but using narrower range between each of the two reflectance than our results. In contrast, they found the best NDSI for φ, using a pair of shorter and longer wavelength bands as NDSI [ 450,1330 ] (r 2 = 0.77), which are both outside of our measurement limitation.
Our data-mining approach found an alternative VI to explain the seasonal variation of photosynthetic capacities in a young larch plantation using an LCTF measurement, although the applicability is limited in the range of the wavelength between 500 and 770 nm. To confirm the observed difference between A max and φ in the relationship to NDSI, and the universal applicability of these NDSIs for predicting regional GPP by microsatellite or UAV observation, the relationships must be validated at other ecosystems and with longer-term observation.

Conclusions
The LCTF camera could capture the seasonal variations of ecosystem photosynthetic capacity and distinguish the vegetation type in a young larch plantation in northern Japan, although there is a limitation in the range and resolution of the wavelength. To explain the seasonal variation of ecosystem photosynthetic capacity, the best NDSI among all pairs of two reflectances was NDSI [ 770,720 ] for A max and NDSI [ 530,600 ] for φ, which reflect the red edge shift owing to the change in the chlorophyll content and relative strength of the light absorbance in the visible red wavelength region compared with that in the green wavelength region, respectively. Predicted GPP using these NDSIs agreed well with the observed values. Universal applicability of these NDSIs and the relationship with the photosynthetic capacity must be validated at other ecosystems and with longer-term observation for predicting regional GPP by microsatellite or UAV observation.