Estimation of area- and mass-based leaf nitrogen contents of wheat and rice crops from water-removed spectra using continuous wavelet analysis

Abstract

Background

The visible and near infrared region has been widely used to estimate the leaf nitrogen (N) content based on the correlation of N with chlorophyll and deep absorption valleys of chlorophyll in this region. However, most absorption features related to N are located in the shortwave infrared (SWIR) region and the physical mechanism of leaf N estimation from fresh leaf reflectance spectra remains unclear. The use of SWIR region may help us reveal the underlying mechanism of casual relationships and better understand the spectral responses to N variation from fresh leaf reflectance spectra. This study combined continuous wavelet analysis (CWA) and water removal technique to improve the estimation of N content and leaf mass per area (LMA) by reducing the effect of water absorption and enhancing absorption signals in the SWIR region. The performance of the wavelet-based method was evaluated for estimating leaf N content and LMA of rice and wheat crops from fresh leaf reflectance spectra collected over a 2-year field experiment and compared with normalization difference (ND)-based spectral indices.

Results

The LMA and area-based N content (N_area) exhibited better correlations with the determined wavelet features derived from the water-removed (WR) spectra (LMA: R² = 0.71, N_area: R² = 0.77) than those from the measured reflectance (MR) spectra (LMA: R² = 0.62, N_area: R² = 0.64). The wavelet features performed remarkably better than the optimized ND indices for the estimations of LMA and N_area with MR spectra or WR spectra. Based on the best estimations of LMA and N_area with wavelet features from WR spectra, the mass-based N content (N_mass) could be retrieved with a high accuracy (R² = 0.82, RMSE = 0.32%) in the indirect way. This accuracy was higher than that for N_mass obtained in the direct use of a single wavelet feature (R² = 0.68, RMSE = 0.42%).

Conclusions

The enhancement of absorption features in the SWIR region through the CWA applied to water-removed (WR) spectra was able to improve the spectroscopic estimation of leaf N content and LMA as compared to that obtained with the reflectance spectra of fresh leaves. The success in estimating LMA and N with this method would advance the spectroscopic estimations of grain quality parameters for staple crops and individual dry matter constituents for various vegetation types.

Background

Rice and wheat are two major staple crops in the world and provide primary dietary calories and protein for the global population [1, 2]. Leaf nitrogen (N) content is an important indicator of crop photosynthetic capacity [3] and is needed by agronomists for making fertilization recommendations [4]. Quantification of leaf N content could provide valuable information for monitoring crop physiology [5], and practicing precise farming [6] so as to improve the use efficiency of nitrogen fertilizers. Remote sensing has been widely used as a non-destructive approach for estimating leaf N content in the past few decades [7,8,9]. The common practice is to establish linear or nonlinear relationships between leaf N content and spectral features derived from leaf reflectance spectra.

Curran [10] listed 42 absorption features caused by bending and stretching of foliar chemical bonds in the 400–2400 nm range and identified many of them related to N in the shortwave infrared region (SWIR) region (1000–2400 nm). Based on a number of N-sensitive wavelengths in the SWIR region, leaf N can be estimated accurately from the reflectance spectra of dried and ground leaves [11, 12]. However, these absorption features are masked by water absorption and hence not clearly visible in the SWIR reflectance spectra of fresh leaves, thereby leading to weaker signals of N in the entire spectra [13,14,15]. This has been proven by the findings that the accuracy of N estimation from the reflectance spectra of fresh leaves is lower than that of dry leaves or dried and ground leaves [11]. Therefore, the SWIR region in the reflectance spectra of fresh leaves has seldom been used for N estimation. Instead, the visible and near infrared (VNIR) region has been widely used for this purpose [8, 16, 17] because nitrogen and chlorophyll are well related [3, 18] and the latter has deep absorption valleys in this region [19, 20]. To date, most studies built N estimation models with the spectral information of fresh leaves in the VNIR region alone [21,22,23] or VNIR and SWIR regions [9, 24, 25]. A few studies made use of SWIR reflectance alone [12] but their focus was on the spectra of dried leaves or leaf powder, rather than fresh leaves that are easier to handle for reflectance measurements. How accurately the N content could be estimated from reflectance spectra of fresh leaves in the SWIR region alone is poorly understood.

As a group of N-containing compounds in fresh leaves, chlorophylls account for only a portion of the total nitrogen [3, 26]. The physical mechanism underlying the spectroscopic estimation of leaf N content remains unclear due to the use of indirect N-sensitive wavelengths in the VNIR region. The use of SWIR region may help us reveal the mechanism and better understand the spectral responses to N variation [10, 11]. A major problem constraining the satisfactory estimation of N content is the effect of leaf water absorption on reflectance spectra, which was claimed to be removed to within 10% by Kokaly and Clark [11]. Recently, a few studies adopted a water removal technique originally proposed by Gao and Goetz [27] to remove the effect of water absorption so as to improve the estimation of N content [15, 28] and N to phosphorous ratio [29]. These studies compared water removal to traditional spectral transformation techniques, but did not decompose the water removal process for understanding the underlying mechanism of this technique. To make use of the spectral information in the SWIR region, they applied traditional methods such as stepwise multiple linear regression and partial least-squares regression to the spectra after the water removal process. Based on these methods, some of the selected wavelengths were not related to the absorption features of the chemical being examined [12] and the regressions were often faced with model overfitting and indirect relationships [10]. In addition, how to enhance the causal absorption features in WR spectra for the improved estimation of leaf N content is poorly understood.

Recently, continuous wavelet analysis (CWA) has been widely used to estimate chlorophyll content [30,31,32], water content [33, 34], dry matter content [35], and leaf area index [36] from leaf and canopy reflectance spectra. After the application of continuous wavelet transform, a reflectance spectrum is decomposed into a number of scale components, which have the same length as the reflectance spectrum and are composed of wavelet features as a function of wavelength and scale [33]. Wavelet features have been proven to be superior to vegetation indices (VIs) in the characterization of absorption by foliar chemicals in reflectance spectra [33, 35, 37]. However, none of previous studies have used CWA and water removal techniques collectively and have investigated the application of CWA to WR spectra for examining dry matter and nitrogen related absorption characteristics. In contrast to the continuum removal operation used in Schlerf et al. [15] and Ramoelo et al. [28], the continuous wavelet transform (CWT) is a linear operation and enables us to decompose the CWT of WR spectra into that of spectral addition or subtraction. Combining CWA and water removal has the potential to improve the estimation of N content by reducing the effect of water absorption and enhancing N absorption signals.

Common measures for expressing leaf N content are either area-based (N_area, g/m²) or mass-based (N_mass, %). Because of its tight correlation with photosynthetic capacity [3] and the widespread use in fertilization management [38, 39], N_mass has been extensively studied and estimated from remotely sensed data [7, 23,24,25, 40]. Both N_area and N_mass can be directly obtained in a destructive way, but remote sensing usually works better for estimating N_area. The leaf biochemistry for dried ground samples is usually expressed as concentration (mass-based, mass per unit dry weight). For intact fresh leaves, the biochemistry is often expressed as content (area-based, mass per unit leaf area) [41]. The use of chemical content may be more suitable for remote sensing applications because it is a better representation of the interaction between matter and light per unit surface area. For example, the units of biochemical parameters used in PROSPECT model are all area-based [19, 42]. It may be proven by the higher correlations of many VIs with N_area than N_mass documented in Hansen and Schjoerring [22] and Jay et al. [43]. Although the interaction of leaves with light per unit leaf area is directly related to such area-based traits as N_area [41, 44], few studies paid attention to the estimation of N_area [13, 43, 45]. As pointed out by Wright et al. [46], N_mass and N_area are interconverted via leaf mass per area (LMA, g/m²). To the best of our knowledge, only one study made use of this connection and derived N_area indirectly from estimated N_mass and LMA [25]. Based on our understanding of N_mass and N_area, this study will investigate first estimating the two area-based factors (N_area and LMA) and then indirectly estimating N_mass instead. The SWIR region encompasses all of the major absorption features of dry matter and N and may provide sufficient spectral information to estimate LMA, N_area and N_mass.

Therefore, the objectives of this study were (1) to determine the wavelet features most sensitive to LMA and N_area from both leaf reflectance spectra and WR spectra in the SWIR region, (2) to evaluate the feasibility of improving LMA, N_area and N_mass estimations with the integration of CWA and the water removal technique, and (3) to compare the CWA and normalization difference (ND) index approaches in the performance of LMA, N_area and N_mass estimations.

Methods

Experimental design

Four experiments were conducted at the experimental station in Rugao, Jiangsu of eastern China (120°45′E and 32°16′N) with two for rice and two for wheat. Both crops were chosen as they were grown in rotation in this experimental area. We intended to develop robust models for both of them despite their differences in biochemical parameters and surface properties [47]. The pooled data represented a wider range of samples and were beneficial to test the stability of the method proposed in this study. The treatment in each experiment represented variations in cultivar type, nitrogen fertilization rate and planting density with three replications. These treatments were applied to create wide ranges for N_mass, N_area and LMA, which includes the extremely high and low values. There were a total of 36 plots with the same size of 5 m × 6 m for each experiment. Details about the experiments and sampling dates are shown in Table 1. Rice was transplanted and grown in a grid pattern with a plant spacing of 15 cm and two row spacings of 30 and 50 cm. Wheat was sowed in drill with different row spacings and the rows were oriented in a south-north direction. The 2-year experimental data for rice and wheat were combined to form RICE and WHEAT datasets, respectively.

Table 1 Summary of designs and ground samplings periods of the field plot experiments

Measurements of reflectance spectra and chemical constituents

Three leaves per plot were collected and each leaf made a sample. Their reflectance spectra were measured using an ASD FieldSpec Pro spectrometer (Analytical Spectral Devices, Boulder, CO, USA) assisted with a leaf clip accessory. Reflectance spectra were taken per leaf for three leaf positions and these spectra were averaged to represent the leaf sample spectrum. The spectrometer collects data at 1.4 nm and 2 nm sampling interval in the 350–1000 nm and 1000–2500 nm spectral regions, respectively. The spectral data were obtained at a 1 nm spectral interval.

After collection of reflectance spectra, the surface area (A, m²) and fresh weight (FW, g) were measured immediately for every leaf sample. The leaf area (A, m²) was determined with a HP G4050 scanner (HP Development Company, L.P., USA) and calculated as the product of pixel count per leaf and the actual area represented by each pixel. The leaf dry weight (DW, g) was measured after the leaves were dried in the oven at 80 °C for 48 h. Finally, equivalent water thickness (EWT, g/m²), leaf mass per area (LMA, g/m²) and area-based leaf nitrogen content (N_area, g/m²) were calculated with the equations below:

$${\text{EWT}} = \left( {{\text{FW}} - {\text{DW}}} \right)/{\text{A}}$$

(1)

$${\text{LMA}} = {\text{DW}}/{\text{A}}$$

(2)

$${\text{N}}_{\text{area}} = {\text{N}}_{\text{mass}} \times {\text{LMA}}$$

(3)

where mass-based leaf nitrogen content (N_mass, %) was determined using Kjeldahl method with SEAL AutoAnalyzer 3 HR (SEAL Analytical, Ltd., German). Note that the unit g/m² of EWT can be interconverted to cm (1 g/m² = 10⁻⁴ g/cm² = 10⁻⁴ cm) via the density of water (1 g/cm³: 1 g water = 1 cm³ water) [48].

The adjusted water-removal technique

The water removal technique was originally proposed by Gao and Goetz [27] with the assumption that there is a linear spectral background level and a nonlinear combination of a reflectance spectrum of leaf water and that of leaf dry matter for fresh leaves. However, the assumption of linear spectral background is only valid for small wavelength regions [27] and this technique is often performed on spectral segments, such as 1500–1780 nm and 2100–2200 nm [15, 27]. The use of narrow and discrete spectral segments may limit the utility of full-range spectral information. Therefore, this technique was adjusted in this study to obtain full-range water removed (WR) spectra by adopting PROSPECT-5B [42] for modeling leaf reflectance spectra. The input parameters of PROSPECT-5B include chlorophyll content (C_ab), carotenoid content (C_xc), EWT, LMA, and structural parameter (N_struc) [42]. To estimate the WR spectra, the first step was to find the best combination of input parameters for PROSPECT-5B by minimizing the following merit function:

$$min = \mathop \sum \limits_{\lambda = 1}^{n} \left( {R_{\lambda } - \widetilde{{R_{\lambda } }}} \right)^{2}$$

(4)

where R_λ and $\widetilde{{R_{\lambda } }}$ are the measured reflectance (MR) and modeled reflectance values of a spectrum at wavelength λ, respectively. The mathematical optimization was performed with the ‘CONSTRAINED_MIN’ function in IDL 8.3 (Exelis Visual Information Solutions, Boulder, CO, USA), which gives the best combination of parameters that minimizes the difference between modeled and measured reflectance for the specified spectral range. Since most nitrogen and dry matter related absorption features are located in the shortwave infrared (SWIR) region [10], the spectral range 1000–2400 nm was used in this study for implementing PROSPECT model inversion as well as the determination of WR spectra. In addition, the C_ab and C_xc were fixed at 30 and 3 μg/cm², respectively, because these pigments have no absorption in this spectral region [42, 47].

After the PROSPECT model inversion process, the second step was to model the leaf water reflectance spectra by inputting the inverted parameters into PROSPECT-5B with zero LMA values. Since pigments have no influences on the reflectance in this spectral region, the parameters C_ab and C_xc were fixed at 30 and 3 μg/cm², respectively [42, 47]. As water was the only absorbing chemical left in fresh leaves for this purpose, the modeled leaf water reflectance spectra were named as water-only (WO) spectra. The residual spectra, namely the WR spectra, were estimated as below:

$$WR_{\lambda } = \left( {R_{\lambda } - R_{\lambda }^{wo} } \right)/\bar{R}$$

(5)

where R ^wo _λ was the reflectance of a WO spectrum at wavelength λ. $\bar{R}$ was the mean reflectance of a MR spectrum over all wavelengths.

After the absorption of dry matter was removed, the amplitude of modeled leaf water reflectance increased for all wavelengths except for the strong water absorption region around 1900 nm (Fig. 1a). Many nitrogen-related absorption features that were not obvious in the MR spectrum became visible in the WR spectrum, such as the feature at 2180 nm (Fig. 1b). Note that the reflectance values have physical meanings and cannot be negative, but the amplitudes of WR spectra were relative values estimated from MR spectra and the negative MR values only meant the negative differences between a reflectance spectra and the corresponding water-only spectrum. Similar terminology is also seen in “derivative spectra”, which represents the derivative data of reflectance spectra and encompasses negative and positive values [49].

Fig. 1

a Represents a measured reflectance spectrum and the modeled water-only spectrum for a fresh leaf and b represents the corresponding water-removed spectrum. The vertical dashed lines correspond to the absorption features of nitrogen centered at 1510 nm, 1690 nm, 2060 nm, 2180 nm, and 2300 nm. In addition, three scaled and shifted mother wavelets (labeled as S3, S4, and S5) at scales 3–5 with the shifting factor of 2180 nm are presented in (b)

Continuous wavelet analysis (CWA)

Mathematically, continuous wavelet transform (CWT) is a linear operation that performs the convolution of reflectance spectrum with a scaled and shifted mother wavelet. The realization process is shown as below:

$$\psi_{a,b} \left( {{\uplambda }} \right) = \frac{1}{\sqrt a }\psi \left( {\frac{\lambda - b}{a}} \right)$$

(6)

$$W_{r} ( {a,b}) = \int\limits_{- \infty }^{+\infty } r(\lambda )\psi_{a,b} ( \uplambda)d\lambda$$

(7)

where $\psi \left( {{\uplambda }} \right)$ is the mother wavelet function and $\psi_{a,b} \left( {{\uplambda }} \right)$ is the scaled and shifted version of $\psi \left( {{\uplambda }} \right)$. W_r(a, b) is the wavelet coefficient (or wavelet feature, denoted as WF_b,a) for the scaling factor a and the shifting factor b. The scaling factor indicates the width of the scaled mother wavelet, which can be comparable with the width of an absorption feature. Narrow absorption features can be captured by a low scaling factor while broad features need high scaling factors. The scaling factor used in this study is at dyadic scales 2³ (scale 3), 2⁴ (scale 4), 2⁵ (scale 5), and 2⁶ (scale 6). The shifting factor determines the central wavelength of shifted mother wavelet, which is used to capture the peak or valley of an absorption feature. Taking the absorption feature of nitrogen at 2180 nm as an example, the three scaled and shifted mother wavelets (Scale 3, Scale 4, and Scale 5) are displayed in Fig. 1b. The spectral width of the corresponding wavelet feature increases with the scaling factor. Different with a spectral index calculated with two discrete wavelengths, a wavelet feature represents the information for a continuous spectral sub-region that determined by scaling and shifting factors. Although CWA is known to be resistant to signal noise, noise reduction is not the focus of this study. Rather, we took the advantage of CWA in capturing and enhancing the absorption features by nitrogen and dry matter. The multiscale property of CWA is well suited for characterizing the absorption changes that may occur in different widths over different wavelength ranges.

The one-dimensional reflectance spectra are converted into two-dimensional wavelet power (magnitude of wavelet coefficient) datasets after the application of CWT. Hence, a scalogram of coefficients of determination (R²) was obtained at all wavelengths and scales [33]. Moreover, the regions that compromised the top 1% R² values were extracted from the comprised scalograms and finally the wavelet feature with the highest R² was selected. In this study, the mother wavelet function was determined as the second derivative of the Gaussian function because of its similarity to the shapes of leaf absorption features [35]. CWT was conducted using the ‘WV_CWT’ function in IDL 8.3 (Exelis Visual Information Solutions, Boulder, CO, USA). Since the wavelet is linearly additive [37], the wavelet power derived from a WR spectrum equals the difference in wavelet power between the MR spectrum and the WO spectrum divided by the mean of reflectance over all wavelengths as follows:

$$WP(WR_{\lambda } ) = WP\left( {R_{\lambda } )/\bar{R} - WP(R_{\lambda }^{wo} } \right)/\bar{R}$$

(8)

where WP represents the wavelet power and other symbols are the same as in Eq. (5).

Evaluation of estimation accuracy

The widely used normalization difference (ND) indices were calculated in this study for comparison with the wavelet features in the evaluation of estimation accuracy. The ND index wavelengths were determined as the best band combination from R² scalograms except these combinations with difference between two wavelengths within 10 nm [50, 51]. Since the aim of this study was to build a robust model across rice and wheat crops, the most sensitive wavelet features and ND indices to LMA, N_area and N_mass were determined separately from pooled data. These spectral features were used to represent their regression relationships and generate direct estimates of those traits. Since area-based foliar traits are directly related to the interaction of foliar constituents with light and can be more easily estimated with spectral features than mass-based traits [41, 44], N_mass was also estimated in an indirect way as the ratio of estimated LMA and N_area values as below:

$${\text{Estimated }}\,{\text{N}}_{\text{mass}} = \frac{{{\text{Estimated }}\,{\text{N}}_{\text{area}} }}{{{\text{Estimate}}\, {\text{LMA}}}} \times 100{{\% }}$$

(9)

The indirect estimation of N_mass, based on two wavelet features separately sensitive to N_area and LMA, might outperform the direct estimation with a single wavelet feature.

The estimation accuracies were evaluated with the metrics R² and root mean square root (RMSE):

$$R^{2} = 1 - \frac{{\mathop \sum \nolimits_{i} \left( {y_{i} - y_{i}^{\prime } } \right)^{2} }}{{\mathop \sum \nolimits_{i} \left( {y_{i} - \bar{y}} \right)^{2} }}$$

(10)

$$RMSE = \sqrt {\frac{{\mathop \sum \nolimits_{i} \left( {y_{i} - y_{i}^{\prime } } \right)^{2} }}{n}}$$

(11)

where y_i and $y_{i}^{\prime }$ are the measured and estimated trait values for sample i. $\bar{y}$ is the arithmetic mean of trait and n is the number of samples. Due to the intrinsic discrepancy between rice and wheat leaves, the best feature for rice may not be exactly the same as that for wheat and vice versa. A generic model that would work best across rice and wheat samples was preferred over that specific to one crop alone. Therefore, the selected features were determined from the pooled data of rice and wheat. These evaluation metrics were calculated from two commonly used validation processes. One is a 10-fold cross validation procedure applied to the pooled data, and the other is the leave-one-out (LOO) method. In the LOO process, three sub-datasets were used for model calibration and the left one was used for validation. This validation process was repeated for the every one of the four sub-datasets (Rice 2015, Rice 2016, Wheat 2016, and Wheat 2017) in Table 1. Since many studies also used the VNIR region to calculate the Red-Edge Chlorophyll Index (CI_{red edge}) for leaf N content estimation [16, 52], a brief comparison between SWIR and VNIR regions was presented.

Results

Descriptive statistics of foliar traits

Differences in foliar traits between the two datasets could be observed in terms of range and magnitude (Fig. 2). Compared with RICE, the WHEAT dataset exhibited a higher mean and a wider range for EWT, N_area and N_mass (p < 0.01) but a lower mean and a narrower range for LMA (p < 0.01). In addition, N_area was positively correlated with N_mass for both RICE (r = 0.79, p < 0.01) and WHEAT (r = 0.92, p < 0.01) (Table 2). N_area was positively correlated with LMA for WHEAT (r = 0.65, p < 0.01) but not correlated with EWT for both datasets. N_mass was correlated with LMA for RICE (r = − 0.56, p < 0.01) and WHEAT (r = 0.31, p < 0.01) but correlated with EWT only for RICE (r = − 0.30, p < 0.01).

Fig. 2

Stacked frequency distributions of a EWT, b LMA, c N_area, and d N_mass for the two datasets RICE and WHEAT. The statistics before and after slash (/) in each distribution plot are shown for RICE and WHEAT, respectively

Table 2 Pearson’s correlation matrix for EWT, LMA, N_area and N_mass (top-right triangle for RICE and bottom-left triangle for WHEAT)

As one of the most important factors in this experiment, the effect of N treatment on leaf reflectance spectra is displayed in Fig. 3 by taking the WHEAT 2017 samples as an example. The N_mass, N_area, and LMA values for WHEAT 2017 increased with the N fertilization rate, and EWT was almost invariant for N1 and N2 (Fig. 3a–d). The leaf reflectance spectra under three N treatments exhibited marginal differences in the SWIR region (Fig. 3e). This phenomenon was expected for fresh leaves, which demonstrated the need of water removal technique to remove the water absorption and CWA to enhance the absorption features of nitrogen and dry matter.

Fig. 3

The boxplots of a N_mass, b N_area, c EWT and d LMA as well as e averaged leaf reflectance spectra under three nitrogen treatments (N0: 0 kg hm⁻², N1: 150 kg hm⁻², N2: 300 kg hm⁻²) for WHEAT 2017 samples

Relationships of N_area with optimal wavelet features and ND indices

The highest R² values between N_area and wavelet features for pooled data were 0.64 and 0.77 for MR and WR spectra, respectively (Fig. 4a, b). The two feature regions determined for MR spectra were centered at 2060 nm and 2180 nm, and both of them matched up with the absorptions of protein and nitrogen [10]. As for the WR spectra, a similar feature region with longer wavelengths was found, and a minor one exhibited much shorter wavelengths than the counterpart for the MR spectra.

Fig. 4

Coefficients of determination (R²) for the relationships of N_area with a MR spectra and b WR spectra derived wavelet features, and c MR spectra and d WR spectra derived ND indices. Regions highlighted in red on the scalograms represent the top 1% R² for wavelet features or ND indices. The downward arrows on the top of the figure indicate the wavelength locations of absorption features with both 2060 nm and 2180 nm for protein and nitrogen

The R² contour maps for the relationship between N_area and ND indices are shown in Fig. 4c, d for MR (maximum R² = 0.40) and WR (maximum R² = 0.37) spectra, respectively. Seven wavelength regions most sensitive to N_area were observed for the MR spectra and five of them were close to the diagonal line, which indicates the presence of similar wavelengths in an ND index. As for the WR spectra, the majority of wavelength combinations were located in the 1700–1800 nm region in which both λ₁ and λ₂ fall.

Figure 5 shows the relationships of N_area with four optimal spectral features (WF_2060,3, WF_2181,4, ND_2192,2202 and ND_1749,1773) determined from the top 1% correlations in Fig. 4. The two optimal wavelet features (WF_2060,3 and WF_2181,4) used the spectral information over two spectral sub-regions that matched well with the absorption features centered at 2060 nm and 2180 nm (Fig. 4a, b). Nevertheless, the two optimal ND indices (ND_2192,2202 and ND_1749,1773) used the spectral information at four individual wavelengths 1749 nm, 1773 nm, 2192 nm and 2202 nm (Fig. 4c, d). All of the relationships for WHEAT, RICE and pooled data were linear. After the application of the WR process, the N_area ~ WF correlation exhibited a pronounced improvement (R² increased from 0.71 to 0.83) over that from the MR spectra for the WHEAT dataset but a slight increase for RICE dataset (Fig. 5a, b). Overall, this caused an increase of R² from 0.64 to 0.77 for the pooled data. For N_area ~ ND correlations, the combination of an increase for RICE and a decrease for WHEAT resulted in an overall decrease in R² (Fig. 5c, d).

Fig. 5

Relationships of N_area with wavelet features (top row) and ND indices (bottom row) derived from MR spectra (left column) and WR spectra (right column). Red, blue and black lines represent linear fits on RICE, WHEAT and pooled data, respectively. All regressions are statistically significant (p < 0.001)

Because of the linearly additive property of CWT, the final regressions from the WR spectra could be decomposed into intermediate regressions to understand how the effect of water absorption was removed. Taking the N_area estimation as an example, Fig. 6 shows the intermediate and final regressions based on Eq. (8) to illustrate the advantages of the water removal technique. N_area was well correlated to the WF_2181,4 derived from the MR spectra (Fig. 6a) but not to the WF_2181,4 derived from the WO spectra for both RICE and WHEAT datasets (Fig. 6b). The WF_2181,4 ~ N_area relationship in the MR spectra was weaker than that in the WR spectra (Fig. 6c, d). After the removal of water absorption information in the WO spectra (Fig. 6e), the WR-derived WF_2181,4 that equaled the difference between MR-derived and WO-derived WF_2181,4 values contained no more water information (p < 0.001) (Fig. 6f). Clearly, the WR-derived WF_2181,4 showed improved relationships with N_area with increases of R² values by 0.23 and 0.18 for RICE and WHEAT datasets, respectively.

Fig. 6

Relationships between WF_2181,4 and (a, b, c) N_area and (d, e, f) EWT derived from MR spectra (left column), WO spectra (middle column) and WR spectra (right column). Particularly, the WF_2181,4 at right column was equal to the difference between WF_2181,4 at left column and middle column. Red and blue lines represent linear fits on RICE and WHEAT, respectively. All regressions are statistically significant (p < 0.001) except for that in (b, f)

Relationships of N_mass with optimal wavelet features and ND indices

The highest R² between N_mass and wavelet features for pooled data were 0.71 and 0.68 for MR and WR spectra, respectively (Fig. 7a, b). The three feature regions determined for MR spectra were centered at 2050 nm, 2120 nm and 2170 nm, all of which corresponded to the absorption wavelengths of protein and nitrogen [10] but showed an offset about 10 nm. As for the WR spectra, only one feature region sensitive to N_mass was found and its center wavelength (2110 nm) was 20 nm off the known absorption feature.

Fig. 7

Coefficients of determination (R²) for the relationships of N_mass with a MR spectra and b WR spectra derived wavelet features and c MR spectra and d WR spectra derived ND indices. Regions highlighted in red on the scalograms represent the top 1% R² for wavelet features or ND indices. The downward arrows on the top of the figure indicate the wavelength locations of absorption features with both 2060 nm and 2180 nm for protein and nitrogen and 2130 nm for protein

The R² contour maps for the relationships between N_mass and ND indices are shown in Fig. 7c, d for MR (maximum R² = 0.32) and WR (maximum R² = 0.37) spectra, respectively. Regardless of MR and WR spectra, wavelength regions sensitive to N_area were dispersed on the contour map and a majority of them were close to the diagonal line.

Figure 8 shows the linear relationships of N_mass with four optimal spectral features (WF_2053,3, WF_2113,4, ND_2281,2291 and ND_2147,2161) determined from the top 1% regions in Fig. 7. With the WR technique, the N_mass ~ WF correlation exhibited a decrease over that from the MR spectra for the individual datasets and for the pooled data (Fig. 8a, b). N_area ~ ND correlations increased from MR spectra to WR spectra but the R² values were still substantially lower than those for the N_mass ~ WF correlations (Fig. 8c, d).

Fig. 8

Relationships of N_mass with (a, b) wavelet features and (c, d) ND indices derived from MR spectra (left column) and WR spectra (right column). Red, blue and black lines represent linear fits on RICE, WHEAT and pooled data, respectively. All regressions are statistically significant (p < 0.001)

Relationships of LMA with optimal wavelet features and ND indices

The highest R² between LMA and wavelet features for pooled data were 0.62 and 0.71 for MR and WR spectra, respectively (Fig. 9a, b). The two feature regions determined for MR spectra were centered at 2210 nm and 2270 nm. The former did not coincide with any known absorption feature but the latter matched up with the absorption features of cellulose, sugar and starch [10]. After the water removal process, a similar feature region was found with longer wavelengths and a narrower range of scales. In addition, a new feature region with shorter wavelengths was found and exhibited a wavelength centered at 1580 nm that was associated with the absorptions of starch and sugar [10].

Fig. 9

Coefficients of determination (R²) for the relationships of LMA with a MR spectra and b WR spectra derived wavelet features and c MR spectra and d WR spectra derived ND indices. Regions highlighted in red on the scalograms represent the top 1% R² for wavelet features or ND indices. The downward arrows on the top of the figure indicate the wavelength locations of absorption features reported in [10] with 1580 nm for starch and sugar and 2270 nm for cellulose, sugar and starch

The R² contour maps for the correlations between LMA and ND indices for pooled data are shown in Fig. 9c, d for MR (maximum R² = 0.60) and WR (maximum R² = 0.60) spectra, respectively. Three wavelength regions most sensitive to LMA for the MR spectra were represented by center wavelengths λ₁ at 1450 nm, 1870 nm and 2230 nm, and λ₂ at 2240 nm. As for the WR spectra, the most significant wavelength combinations were located in the regions with both λ₁ and λ₂ shorter than 1800 nm. The ND index regions representing the top 1% R² with LMA were generally more scattered than the wavelet feature regions, and therefore more difficult to relate with major absorption features documented in the literature.

Figure 10 shows the relationships of LMA with four optimal spectral features (WF_2273,4, WF_1578,4, ND_2236,2247 and ND_1716,1727) determined as the best candidates from the top 1% correlations in Fig. 9. These relationships were linear for WHEAT, RICE and pooled data. With the WR technique, the LMA ~ WF correlation exhibited a pronounced improvement (R² increased from 0.45 to 0.66) over that from the MR spectra for the RICE dataset but did not change for WHEAT (Fig. 10a, b). This caused an increase of R² from 0.62 to 0.71 for the pooled data. For LMA ~ ND correlations, the combination of a decrease for RICE and an increase for WHEAT led to no change in R² for the pooled data (Fig. 10c, d).

Fig. 10

Relationships of LMA with (a, b) wavelet features and (c, d) ND indices derived from MR spectra (left column) and WR spectra (right column). Red, blue and black lines represent linear fits on RICE, WHEAT and pooled data, respectively. All regressions are statistically significant (p < 0.001)

Assessment for the estimations of LMA, N_area and N_mass

Table 3 shows the estimation accuracy for LMA, N_area and N_mass over the pooled data with the 10-fold cross-validation. In particular, N_mass was estimated with the determined wavelet features as shown in Fig. 10 (the direct way) and with the ratio of estimated N_mass and LMA values (the indirect way). In all cases, wavelet features outperformed ND indices substantially in the estimations of N_area and N_mass for both MR and WR spectra. For the estimation of N_area, the accuracy with the best correlated wavelet features was improved from MR spectra (R² = 0.64 and RMSE = 0.22 g/m²) to WR spectra (R² = 0.77 and RMSE = 0.18 g/m²), but that with ND indices even decreased slightly. For the direct estimation of N_mass with the best correlated wavelet features and ND indices, there was no significant difference in estimation accuracy before and after the application of water removal. As the ratio of wavelet-derived N_area and LMA, N_mass could be better estimated from WR spectra (R² = 0.82 and RMSE = 0.32%) than from MR spectra (R² = 0.71 and RMSE = 0.40%). As can be seen in Fig. 11, the scatter plot of estimated versus measured N_mass is closer to the 1:1 line, especially for samples in the WHEAT dataset. This indirect estimation of N_mass from WR spectra was also more accurate than the direct estimation. The improvement in the indirect estimation of N_mass from MR spectra to WR spectra was not seen for ND indices. The performance of these wavelet features derived from the SWIR region (R² = 0.77 and 0.82 for N_area and N_mass, respectively) was better than CI_{red edge} from the VNIR region (R² = 0.76 and 0.71 for N_area and N_mass, respectively), especially for N_mass.

Table 3 Assessment of LMA, N_area and N_mass estimations generated with the direct and indirect ways

Fig. 11

The scatter plots of measured and estimated N_mass derived as the ratio of LMA and N_area estimated from (a, b) wavelet features and (c, d) ND indices from MR spectra (leaf column) and WR spectra (right column). The LMA and N_area estimations were generated from the 10-fold cross validation process

In addition to the 10-fold cross validation, the LOO validation accuracies are listed in Table 4. Generally, the estimation accuracy from the LOO process was lower. However, both wavelet features and ND indices derived from the WR spectra still produced higher R² for LMA, N_area and the ratio of N_area to LMA (N_mass in an indirect way) than those from the MR spectra. As for N_mass estimation, the indirect way with wavelet features (R² = 0.79 and RMSE = 0.34%) still worked better than the direct way (R² = 0.61 and RMSE = 0.46%).

Table 4 Assessment of LMA, N_area and N_mass estimations generated with the direct ways and the indirect way

Discussion

Difference in sensitive spectral features derived from MR and WR spectra

When comparing the spectral features specifically sensitive to LMA, N_area and N_mass, we found that the difference in wavelet feature was less significant between MR and WR spectra than that in ND index. The highlighted feature regions for ND indices were generally more dispersed across the scalograms, which complicated the selection of the best ND indices and was also seen in the studies by le Maire et al. [51] and Wang et al. [53]. Two dispersed index regions with even similar R² values would lead to substantially different wavelength combinations, not to mention many of them as highlighted in Figs. 4, 7 and 9. The sensitive wavelet features were mostly concentrated in no more than two feature regions, from which it would be easier to select the best wavelet features in less than three scales and wavelength ranges. The mismatch of the wavelengths selected by ND indices with the absorption centers suggested the weakness of ND indices and the strength of wavelet features in capturing the changes in nitrogen absorption features. This might be attributed to the different calculation principles for ND indices [54] and wavelet features [33]. The former was determined by the reflectance values at two wavelengths and was sensitive to the amplitude of reflectance spectra, but the latter mainly captured the shape of a spectral segment [35]. Before and after water removal, the amplitude of spectra changed but the shapes of absorption features retained and became even more obvious. Therefore, common sensitive wavelet features can be observed between MR and WR spectra, but not for ND indices.

Necessity of removing the effect of water absorption with CWA

The absorption of leaf water is strong in the SWIR region [19, 42, 48], hence the subtle absorption features of such dry matter constituents as protein and nitrogen are masked in the reflectance spectra of fresh leaves (Fig. 1). Many studies have shown that this masking effect reduces the estimation accuracy of these chemical constituents [14, 15, 55]. This is supported by our findings that the estimations of LMA and N_area were significantly improved through the wavelet-based method after the effect of water absorption had been removed. Particularly, Kokaly and Clark [11] modeled the estimation accuracy of leaf chemistry under eight levels of water content in percentage of fresh weight and found that the accuracy decreased rapidly when leaf water content exceeded 30%. To date, the effect of water absorption on ND indices and wavelet features is still poorly understood. Here we present an exploration of their performance for N_area estimation under different water content levels in EWT. To avoid the confounding factors of leaf surface property and internal structure, the WHEAT dataset with a wider range of EWT was selected to generate three subsets with EWT levels of 130–140 g/m² (W1), 150–160 g/m² (W2), and 170–180 g/m² (W3).

The optimal index ND_2192,2202 was well correlated to N_area only at the lowest EWT level (Fig. 12a), but WF_2060,3 bore stronger relationships with N_area at every EWT level (Fig. 12b). In contrast to the close relationships of N_area and ND_2192,2202 under three EWT levels, the N_area ~ WF_2060,3 relationship was sensitive to EWT level. This is probably because WF_2060,3 made use of the spectral information in the 2020–2100 nm range [33], which is the right wing of the water absorption valley centered at 1950 nm [19] and whose shape changes as a function of water content [11]. In spite of the tight N_area ~ WF_2060,3 relationships, the sensitivity of WF_2060,3 to EWT reduced the model transferability, e.g., from subset W1 to subset W3. Therefore, it was essential to remove water absorption for building a robust model that can be directly applied to datasets with different EWT levels.

Fig. 12

Relationships of N_area with a ND_2192,2202 and b WF_2060,3 derived from the MR spectra at three levels of EWT (W1: 130–140 g/m², W2: 150–160 g/m² and W3: 170–180 g/m²) in the WHEAT dataset. Red, blue and green lines represent linear fits for the data points from W1, W2 and W3, respectively. All regressions are statistically significant (p < 0.001)

How did the water removal technique improve the estimations of N_area and N_mass?

The accuracies of LMA and N_area estimations using wavelet features were improved when the effect of water absorption was removed with the water removal technique. Some recent studies used WR spectra through spectral indices or partial least squares regression to improve the estimations of LMA and N_mass [15, 28, 29, 56], but they did not provide explanations as to why those improvements were obtained. This study reported the application of CWA to the WR spectra for estimating foliar traits by removing the effect of water absorption and enhancing the absorption features of foliar traits. The better WF_2181,4 ~ N_area relationships in the WR spectra than in the MR spectra (Fig. 6a, c) could be explained by the information on EWT in WF_2181,4, which was derived from the MR spectra (Fig. 6d) but not from the WR spectra (Fig. 6f).

With regard to the water removal technique for N_mass, the wavelet-based direct way did not improve the estimation accuracy but the indirect way exhibited a significant improvement (Table 3). This could be attributed to the improved estimations of both LMA and N_area. Kokaly and Clark [11] suggested any remote sensing algorithm for retrieving foliar chemicals from fresh leaf spectra should remove the influence of water absorption. Wang et al. [56] used the normalized dry matter index (NDMI) from WR spectra to estimate the LMA but only obtained slight improvement over MR spectra as for the use of ND indices in this study (Fig. 10c, d). The wavelet features determined for LMA estimation by Cheng et al. [35] from a wide range of plant species performed better than NDMI, but they might still suffer from the effect of water absorption and be hard to be applied to other datasets with different water content levels. Therefore, the difficulty in LMA estimation hampered the indirect estimation of N_mass due to the dominant role of LMA as the denominator in Eq. (9). For the first time, this study obtained improved estimations for LMA as well as N_area by applying CWA to the WR spectra and eventually increased the estimation accuracy for N_mass in the indirect way.

The close fits of data points to the linear regression and the 1:1 lines in Figs. 5 and 11 suggested the best models worked well for the whole ranges of N_area and N_mass. Although the rice and wheat plants in this experiment were deliberately treated in the field, the data were collected in natural light conditions and the models should be applicable to plants growing in the same condition and without any N treatment. To take advantage of the proposed strategy of combining CWA and water removal, it is imperative to have a dataset exhibiting weak covariance of EWT with LMA and N_area. Generally, this applies to leaf level but may not apply to canopy level where the traits are determined as the products of leaf traits and leaf area index and high correlations often exist between EWT_canopy and other traits. In that case, it is difficult to find a wavelet feature insensitive to EWT_canopy but sensitive to another canopy trait. Therefore, the water removal procedure would have limited contribution to the improvement of N_mass or N_area estimation. To test the suitability of our strategy at canopy level, one would need to avoid crop datasets with high correlations of canopy EWT and nitrogen related traits. Alternatively, the strategy could be applied to natural plant types such as trees [15] and grasses [28, 29]. Although this research was conducted at leaf level, it represented a significant advance in the mechanistic understanding of water removal technique for removing water signals and enhancing nitrogen and dry matter signals in fresh leaves. The leaf-level analysis concentrated on chemical signals and laid the foundation of water removal application for scaling up to the canopy level. To better understand the mechanism at canopy level, the strategy should be examined in detail by considering other external factors such as canopy structure [57, 58], soil background [59, 60], solar zenith angle [61, 62], view zenith angle [7, 43], and diffuse/direct light condition [63]. Finally, the signal to noise ratio should also be accounted for in case the nitrogen-related second overtones are suppressed by the spectral noise from instruments and atmospheric interference [10].

Physical interpretations of determined wavelet features

Low-scale wavelet features are able to capture the local shape of absorption features associated with biochemical parameters while high-scale ones characterize the broad continuum in leaf reflectance spectra [30, 33, 35]. The matching wavelet features could be found for N_area and LMA but not for N_mass, which might suggest that the area-based expressions (N_area and LMA) were better representations of the interaction between matter and light per unit surface area [41].

The absorption feature at 2270 nm caused by bending and stretching of chemical bonds in cellulose, sugar and starch does not obviously appear in the reflectance spectra (Figs. 1, 13a), but it is well represented in the wavelet power spectra of MR-derived WF_2273,4 (Fig. 13c). Since LMA is equivalent to the dry matter content by definition, we could relate spectral variation to dry matter absorption from the perspective of spectroscopy of foliar chemistry as done in previous studies [35]. The amplitude of reflectance in the 2200–2350 nm range increases with LMA (Fig. 13a), which contradicts the common sense that the reflectance should decrease with dry matter absorption. This might be attributed to the leaf internal scattering as well as the confounding effect of water absorption on leaf reflectance, which makes it difficult to derive the LMA with traditional techniques, such as ND index [55, 64]. In contrast, the depth of absorption center at 2270 nm in the wavelet power spectra follows the trend of LMA variation (Fig. 13c). Moreover, the absorption feature at 1580 nm by starch and cellulose are invisible in the MR spectra (Fig. 1), but it appears clearly in the WR spectra after water removal (Fig. 13b). The change in the amplitude at WR spectra-derived WF_2273,4 in the wavelet power spectra matches with the variation in LMA (Fig. 13d).

Fig. 13

Effect of LMA on spectral variation between 2200 and 2350 nm (leaf panel) and 1500–1650 nm (right panel). a, c Represent the MR spectra and the corresponding wavelet power spectra at scale 4. b, d Represent the WR spectra and the corresponding wavelet power at scale 4. The location and width of the Gaussian-shaped wavelets for the best LMA-sensitive wavelet features WF_2273,3 and WF_1578,4 from MR and WR spectra are presented in (a) and (b), respectively. The black arrows in (c) and (d) point to the wavelength locations of WF_2273,3 and WF_1578,4, respectively

Similar to dry matter, the absorption of nitrogen at 2060 nm and 2180 nm [10] could be well characterized by the optimal N_area sensitive wavelet features WF_2060,3 (MR spectra) and WF_2180,4 (WR spectra) (Fig. 14). The close matches in wavelength between the determined wavelet features and the documented absorption centers documented [10,11,12] convince us of the physical interpretations underlying the empirical WF-N_area relationships.

Fig. 14

Effect of N_area on spectral variation between 2000 and 2120 nm (leaf panel) and 2100–2250 nm (right panel). a, c Represent the MR spectra and the corresponding wavelet power spectra at scale 3. b, d Represent the WR spectra and the corresponding wavelet power at scale 4. The location and width of the Gaussian-shaped wavelets for the best N_area-sensitive wavelet features WF_2060,3 and WF_2181,4 from MR and WR spectra are presented in (a) and (b), respectively. The black arrows in (c) and (d) point to the wavelength location of the WF_2060,3 and WF_2181,4, respectively

Conclusions

The enhancement of absorption features of nitrogen in the SWIR region would help improve the spectroscopic estimation of leaf nitrogen content from the SWIR reflectance spectra of fresh leaves. This study evaluated the applicability of the water-removal technique for increasing the estimation accuracies of LMA, N_area and N_mass with CWA and spectral index methods. After removing water absorption, the estimation accuracies for LMA and N_area were significantly improved with the use of individual wavelet features but not with ND indices. The two wavelet features WF_1578,4 and WF_2181,4 derived from the WR spectra produced the best estimations for LMA (R² = 0.70, RMSE = 4.38 g/m²) and N_area (R² = 0.77, RMSE = 0.18 g/m²), respectively. Compared to that from the MR spectra, the estimation of N_mass from the WR spectra only changed marginally with individual wavelet features (MR spectra: R² = 0.71, RMSE = 0.40%; WR spectra: R² = 0.68, RMSE = 0.42%). However, a notable improvement was obtained when indirectly deriving N_mass as a ratio of N_area to LMA (MR spectra: R² = 0.71, RMSE = 0.40%; WR spectra: R² = 0.82, RMSE = 0.32%). Generally, the determined wavelet features performed remarkably better than the optimized ND indices for the estimations of all the three traits.

The results demonstrated the feasibility of improving the estimation of N_area and N_mass from the SWIR reflectance spectra by applying CWA to WR spectra. This procedure could suppress the effect of water absorption and enhance the absorption features of foliar chemicals. The characterization of nitrogen absorption variation in the SWIR region with wavelet features offers physical interpretations for the direct detection of leaf nitrogen content and complements the indirect detection of nitrogen content via the use of chlorophyll absorption bands in the red edge region. The integration of CWA and water removal provides a new insight for better understanding spectral responses to the variation in LMA and leaf nitrogen content and has the potential for applications to other crops and different vegetation types. The technique developed in this study could be useful for high throughput phenotyping of leaf nitrogen related traits and the non-destructive detection of nitrogen stress in crops.