• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Evaluation of UAV-derived multimodal remote sensing data for biomass prediction and drought tolerance assessment in bioenergy sorghum

    2022-10-12 09:30:56JitingLiDnielSchchtmnCoyCreechLinWngYufengGeYeyinShi
    The Crop Journal 2022年5期

    Jiting Li ,Dniel P.Schchtmn ,Coy F.Creech ,Lin Wng ,Yufeng Ge,c ,Yeyin Shi,*

    a Department of Biological Systems Engineering,University of Nebraska-Lincoln,Lincoln,NE 68583,USA

    b Department of Agronomy and Horticulture,Center for Plant Science Innovation,University of Nebraska-Lincoln,Lincoln,NE 68588,USA

    c Center for Plant Science Innovation,University of Nebraska-Lincoln,Lincoln,NE 68588,USA

    d Department of Agronomy and Horticulture,Panhandle Research and Extension Center,University of Nebraska-Lincoln,Scottsbluff,NE 69361,USA

    Keywords:Crop phenotyping Unmanned Aerial System (UAS) Thermal Machine learning Drought stress

    ABSTRACT Screening for drought tolerance is critical to ensure high biomass production of bioenergy sorghum in arid or semi-arid environments.The bottleneck in drought tolerance selection is the challenge of accurately predicting biomass for a large number of genotypes.Although biomass prediction by lowaltitude remote sensing has been widely investigated on various crops,the performance of the predictions are not consistent,especially when applied in a breeding context with hundreds of genotypes.In some cases,biomass prediction of a large group of genotypes benefited from multimodal remote sensing data;while in other cases,the benefits were not obvious.In this study,we evaluated the performance of single and multimodal data (thermal,RGB,and multispectral) derived from an unmanned aerial vehicle(UAV) for biomass prediction for drought tolerance assessments within a context of bioenergy sorghum breeding.The biomass of 360 sorghum genotypes grown under well-watered and water-stressed regimes was predicted with a series of UAV-derived canopy features,including canopy structure,spectral reflectance,and thermal radiation features.Biomass predictions using canopy features derived from the multimodal data showed comparable performance with the best results obtained with the single modal data with coefficients of determination(R2)ranging from 0.40 to 0.53 under water-stressed environment and 0.11 to 0.35 under well-watered environment.The significance in biomass prediction was highest with multispectral followed by RGB and lowest with the thermal sensor.Finally,two well-recognized yieldbased drought tolerance indices were calculated from ground truth biomass data and UAV predicted biomass,respectively.Results showed that the geometric mean productivity index outperformed the yield stability index in terms of the potential for reliable predictions by the remotely sensed data.Collectively,this study demonstrated a promising strategy for the use of different UAV-based imaging sensors to quantify yield-based drought tolerance.

    1.Introduction

    Bioenergy sorghum is a photoperiod-sensitive sorghum that either flowers late or not at all in the northern hemisphere.Owing to its enhanced photoperiod sensitivity and prolonged vegetative growth duration,bioenergy sorghum has been found to produce more than twice the biomass of grain sorghum [1,2] and is therefore being developed as an important feedstock for bioenergy production.It grows in marginal areas,such as arid regions;hence,drought tolerance is one of the most important traits for the production of bioenergy sorghum.Although the prolonged vegetative growth stage is one factor that enhances its drought-tolerance ability [3],selecting the most drought tolerant sorghum genotypes is an important breeding goal to maximize biomass yields in marginal environments.

    Over many decades,researchers have used direct or indirect selection criteria for drought tolerance.As a direct selection criterion,primary traits such as grain yield or biomass are of great interest to researchers.Based on primary traits,numerous drought tolerance indices have also been developed to evaluate the drought-adaptive performance of different plant genotypes.Some commonly used indices include geometric mean productivity(GMP),mean productivity(MP),harmonic mean(HM),stress tolerance index (STI),stress susceptibility index (SSI),tolerance index(TOL),and yield stability index(YSI)[4,5].These indices can be largely grouped into two categories.The first group includes GMP,MP,HM,and STI,which have been found to be more reliable for identifying genotypes with both high yield potential and stress tolerance potential[6,7].YSI and SSI fall into the second group that is useful for discriminating genotypes with higher stability across multiple environments[8,9].The indices within each group generally have strong correlations with each other[5,9].In this study,we were interested in the genotypes with both high biomass and yield stability.Therefore,we identified two drought tolerance indices,i.e.,GMP and YSI,as the drought tolerance indicators for selecting the best drought tolerant bioenergy sorghum genotypes.

    In addition to the direct selection criterion,numerous secondary traits have also been evaluated and used as indirect selection criteria for drought tolerance[10].Secondary traits are related to primary traits but provide additional information on crop growth.Compared to primary traits that are usually measured at the end of the season,the secondary traits are easier and faster to measure during the season,are highly heritable,and may have good correlations with the primary traits [11,12].For sorghum,a suite of secondary traits has been evaluated in drought tolerance breeding,such as leaf rolling,leaf water potential,stay-green,stomatal conductance,plant height,and canopy temperature[13,14].Although these traits have shown to be useful in drought tolerance selection,phenotyping them in a high-throughput manner in the field remains very labor-intensive and slow.A promising solution is to replace the traditionally labor-intensive manual measurements with new advanced remote sensing techniques.

    Remote sensing is capable of collecting data in a faster,nondestructive,and more cost-effective manner [15].Remote sensing data can be further used to estimate the above primary and secondary traits.For example,one of the most well-characterized drought adaptative traits in sorghum,stay-green [16],has traditionally been evaluated using visual scores based on the fraction of senesced leaf area[17]or the fraction of green leaves[18].Using remote sensing tools,it is possible to estimate stay-green nondestructively from spectral indices,such as the normalized difference red edge (NDRE) [19] and the normalized difference vegetation index (NDVI) [20].Additionally,another important drought adaptative trait,canopy temperature,has been well-documented in sorghum [21],wheat [22],potato [23],and other crops [24-26].It is usually measured by remote sensors,such as thermal infrared imagery and infrared thermometry.Early in 1989,researchers used handheld infrared thermometers to measure canopy temperature and reported that the drought-resistant sorghum genotypes had lower canopy temperatures than the drought susceptible ones[27].These successful applications of remote sensors in measuring drought responsive secondary traits indicate the possibility of applying such techniques to a drought tolerance breeding program.

    Among various remote sensing platforms,the unmanned aerial vehicle (UAV) has gained increasing attention since it provides higher spatial resolution than satellite platforms and allows faster data collection than ground platforms.Moreover,with the miniaturization of optical cameras,a large number of successful applications of UAV systems have been reported in agricultural fields[28-31].Many canopy features can be derived from the sensors attached to UAVs that may be directly or indirectly related to drought responsive traits,such as canopy temperature,vegetation indices,and plant height.This has led to the use of the UAV system for evaluating the drought tolerance in a bioenergy sorghum association panel in this study.

    Only a limited number of studies have applied UAV systems to drought tolerance screening [5,32-34].UAV thermal imaging was used to evaluate the drought responses of 503 poplar tree genotypes [5].One-fourth of the population was identified as candidates for drought-tolerant genotypes using a stress susceptibility index calculated with canopy temperature.In addition to UAVacquired thermal imaging,RGB images were also added [32,33]to assess the drought tolerance of 390 genotypes of forage grasses.In that study,the drought tolerance was rated by visual scoring of the amount of green biomass using the human eye.Multiple RGBbased indices and thermal-based indices correlated with the breeder scores.The results showed that the thermal indices had weaker correlations than the RGB indices with the visual scores.Given the intrinsic complexity of plant drought adaptation mechanisms[10],different genotypes may have a unique set of physiological and biochemical responses.Therefore,many studies have emphasized the importance of combining multiple selection criteria for selecting drought-tolerant genotypes[11,22].Using UAV technology,the utilization of data sets from multiple sensors such as thermal,RGB,and multispectral cameras makes this approach more feasible.Image features extracted from multiple UAV-based cameras were used together with a support vector machine classifier to distinguish between two groups of genotypes with differential wilting times with an average classification accuracy of 80% [34].While earlier results suggested that UAVs provide a rapid,objective,and labor-saving method for distinguishing drought tolerance within large populations of different plant genotypes,further studies are needed to quantitatively evaluate the potential of UAVbased multimodal images (thermal,RGB,and multispectral) for yield-based drought tolerance selection.

    In our study,the quantitative evaluation was focused on the estimation of biomass since this was the primary trait for the selection of drought tolerant bioenergy sorghum genotypes.Although UAV-based remote sensing data have been used to estimate biomass for various crops,such as maize [35],wheat [36],pea [37],and sorghum [38,39],the estimation performance varies amongst studies on the same crop species.For example,Masjedi et al.[38]conducted experiments for sorghum biomass prediction using UAV-based RGB,hyperspectral and LiDAR sensors.The number of genotypes evaluated ranged from 4 to 840.Results showed that the R2between estimated and measured end-of-season biomass ranged from 0.64 to 0.89.In another study that targeted biomass estimation of ten sorghum varieties using vegetation indices derived from UAV-based remote sensing data,the maximum R2achieved was 0.91[39].Differences in the estimation performance among these studies may be attributed to variation in experimental design,types of sensors used,modeling approach,number of genotypes,sample size,and other factors.In contrast to other studies highlighted above,our study sought to assess the use of UAV-based single and multimodal data (thermal,RGB,and multispectral) to quantify biomass production and to evaluate drought tolerance in a large bioenergy sorghum association panel [40].A linear support vector regression model was first established to predict biomass using the canopy features extracted from UAV multimodal images.Thereafter,we assessed the reliability of using the UAV predicted biomass to derive the yield-based drought tolerance indices (i.e.,GMP and YSI).

    2.Materials and methods

    2.1.Field experiment and biomass sampling

    The experiment was conducted at the Panhandle Research and Extension Center of the University of Nebraska-Lincoln located in Scottsbluff,NE,USA (Fig.1).The soil was relatively uniform throughout the field.Four randomized blocks were arranged in the field,including two well-watered (WW) replicates and two water-stressed (WS) replicates in 30 plots × 12 ranges grid.The size of each block was 75 m long and 52 m wide.A total of 360 sorghum genotypes were planted in each block,consisting of 357 accessions of the Bioenergy Association Panel (BAP) population[40],one Grassl [40],one BTx623 [41],and one PI_92264.The entries in each block were fully randomized.Individual plots were 4.5 m long and 1.5 m wide,consisting of two rows spaced 76 cm apart.Each plot was planted with 100 seeds,with 50 seeds per row.

    On April 16,2020,112 kg ha-1of nitrogen urea was applied to the field and incorporated using tillage.Planting occurred on June 1,2020 for all 1440 plots using a precision vacuum planter.For the WW replicates,29.35 cm of supplemental irrigation was applied on a weekly basis.For the WS replicates,3.05 cm of supplemental irrigation was applied initially to establish the plants in the field.Fresh and dry biomass was harvested on September 2-4,using a tractor pulled a forage harvester with a weigh hopper.The entire plot was harvested and chopped.Whole plot weights were recorded by load cells below the weigh hopper to the nearest 0.23 kg.A subsample of the wet chopped-up plant material was taken from each plot,and its fresh biomass was obtained.This subsample was then put into a small mesh drawstring bag and placed in drying ovens until sample weights stopped decreasing.Dried plant material was taken out of the bag and weighed.Both fresh and dry biomass units were finally converted to kg m-2for each plot.

    2.2.Statistical analysis for fresh and dry biomass

    Two-way ANOVA(analysis of variance)was carried out to evaluate the effects of water treatment,genotype,and the interaction between genotype and treatment (G × T) on the observed fresh and dry biomass.The difference was considered to be statistically significant at α=0.001.From the ANOVA output,the expected mean of squares (EMS) can be estimated by ‘mean of squares’,which were further used to calculate genotypic varianceG×T variance,and residual variancefollowing Equation (1-3) [11].Thereafter,the broad-sense heritabilitywas estimated over combined treatments according to Equation(4)[5].

    Where EMSG,EMSG×T,and EMSεare the mean of squares from ANOVA output for genotype,G×T,and residual error components,respectively.The number of replicates for a genotype within each treatment is r,and t is the total number of treatments.The total phenotypic variance is

    Additionally,the impact of water deficit on biomass was evaluated by computing mean,percent reduction (PR) in WS as compared to WW (Equation (5)) [42],and coefficient of variation (CV,Equation(6)).Pearson correlation coefficient(r)was used to examine the correlations between biomass and UAV-derived individual canopy feature.

    Where μwwand μwsare the mean value under well-water and water-stressed conditions.SD and μ are the standard deviation and mean of the target variable under well-water or waterstressed conditions.

    2.3.UAV systems and flight missions

    Two unmanned aerial systems were used for data collection,as depicted in Fig.1.The six-rotor UAV is the DJI Matrice 600 Pro(DJI,Shenzhen,Guangdong,China),equipped with the Zenmuse X5R RGB camera (DJI,Shenzhen,Guangdong,China) and the five-band RedEdge multispectral camera (MicaSense,Inc.Seattle,WA,USA).RGB and multispectral images were obtained using this system.Detailed specifications of the RGB and the multispectral cameras were provided in Table S1.The four-rotor aerial vehicle,DJI Matrice 210 RTK (DJI,Shenzhen,Guangdong,China),was used to collect thermal images by the carried Zenmuse XT2 camera (DJI,Shenzhen,Guangdong,China).This camera is a dual-sensor system: a radiometric thermal sensor made up of uncooled Vox microbolometers and a visible sensor.In this study,only the images from the thermal sensor were used.Each thermal pixel records temperature value in Celsius degrees.More specifications of this thermal sensor are listed in Table S1.

    Two UAV data collections were conducted in 2020.The corresponding drone flight settings and weather conditions were listed in Table S2.The weather during the flights was mostly sunny with mild wind.All the images were collected within the period from 11:30 to 13:30 local time,except on August 5 when the thermal images were collected at 15:50 local time due to unexpected system malfunction.As suggested by Kelly et al.[43],the thermal camera was turned on and stabilized for about 30 mins before the flight.The purpose was to alleviate any temperature dependency effect,which is the effect of the thermal sensor’s temperature on the thermal responsive signal [44].

    2.4.Image processing and canopy feature extraction

    Pix4Dmapper (Pix4D,Lausanne,Switzerland) software was used to generate orthomosaic from raw aerial images.During the map generation,geometric correction was performed in Pix4Dmapper to correct geometric distortion and to ensure that maps from different sensors or different times were well-aligned.The correction utilized five ground control points,i.e.,checkered wooden boards,that were distributed evenly in the field.The geo-locations of the center of those boards were surveyed in the field with a Real-Time Kinematic GPS receiver(Topcon Positioning Systems,Inc.,Tokyo,Japan) with centimeter-level accuracy.The radiometric calibration for the multispectral images involved two sources of calibration references: (1) the MicaSense’s Calibrated Reflectance Panel provided by the manufacturer;(2) three standard reflectance tarps laid out in the field.Details on the multispectral radiometric calibration are described in Li et al.[45].The thermal calibration description can be found in Fig.S1.Specifically,three reference targets were laid out in the middle of the field,including a cold reference (a cooler containing water and ice mixture),a plant sample representing medium temperature,and a hot reference (the darkest radiometric correction tarp).The surface temperature of these targets was measured concurrently with the flights,using a handheld infrared radiometer (Apogee Instruments,Inc.,Logan,UT,USA).For each data collection,a linear calibration function was established between the measured reference temperatures of the targets and their corresponding pixel values averaged from the thermal map (Fig.S1).The linear calibration function was then used to calibrate the original thermal map.From the corrected orthomosaic(i.e.,thermal,RGB,and multispectral maps),six canopy features were extracted: normalized relative canopy temperature (NRCT),canopy cover,plant height,and three vegetation indices.

    2.4.1.Normalized relative canopy temperature (NRCT) and canopy cover

    The NRCT and canopy cover were extracted together,using the RGB co-registration method to remove soil pixels from the RGB and thermal map [46].The co-registration method makes use of the higher spatial resolution of the RGB map and could partially alleviate the effect of the mixed pixels in the thermal map with lower spatial resolution.The mixed pixel is the pixel that covers both vegetation and soil [47].The RGB co-registration workflow is shown in Fig.S2.First,the RGB map was converted to the HSV(i.e.,Hue,Saturation,Value) color space.Based on the literature[48] and manual adjustment,a pixel was considered to be plant pixel when its value in the Hue channel ranged from 20 to 120,in the Saturation channel ranged from 50 to 255,and in the Value channel ranged from 80 to 255.Using this approach,a binary plant mask was created,where one represents plant pixel and zero represents non-plant pixel.Second,the original thermal map was resampled to the same size as the plant mask using the nearest neighbor interpolation method.The resampling allows us to overlay the plant mask onto the thermal map to generate the segmented thermal map,where the non-plant pixels were set to Null value.

    From the segmented thermal map,mean canopy temperature(°C)was derived from each plot.To normalize across different data collection events,the canopy temperature was converted to NRCT,according to Equation(7)[49,50].In other studies canopy-air temperature difference was used,but because of the speed with which the UAV can cover the entire field(approximately 10 mins),the air temperature did not change significantly so that it was not necessary to factor the air temperature into our data analysis.The canopy cover was also derived from the segmented thermal map,as the ratio between the number of plant pixels and the number of all pixels falling into each plot.

    Where Tcanopyis the mean canopy temperature (°C) from the soil-removed thermal map;Tminand Tmaxrepresent the minimum and maximum temperature (°C) within all the pixels in the soilremoved thermal map.

    2.4.2.Plant height

    To obtain the plant height map,the digital surfaced model(DSM)was generated from the RGB images using the Pix4Dmapper software.The DSM represents the elevation above sea level of the surface,including the surface of natural or artificial features.Moreover,a digital elevation model (DEM) indicating the elevation above sea level of the bare soil was also required.In this study,the DEM was created by interpolating soil points sampled from the DSM map,using the Kriging interpolation tool in ArcMap 10.5.1(Esri Inc.,Redlands,CA,USA).The plant height map was then derived as the difference between the DSM and the DEM maps.Then,for each plot,the trimmed mean value was extracted to represent plot-level average plant height: the mean value of all the pixels falling into the range of 25th to 75th percentiles of plant height pixel values within one plot.

    2.4.3.Vegetation indices

    Finally,three vegetation indices were computed from the multispectral bands: NDVI,NDRE,and RDVI (i.e.,renormalized difference vegetation index).These indices have already been widely used and were reported to be efficient in estimating crop water stress [51].The formulae for these indices were given in Equation(8-10).

    Where RNIR,RRED,RRedEdgeare reflectance value of the nearinfrared,red,and red edge spectral band,respectively.

    2.5.Linear support vector regressor for biomass prediction

    To estimate the biomass using the extracted canopy features,the support vector regressor with the linear kernel (Linear-SVR)was utilized.The support vector regressor is the regression equivalent of the support vector machine(SVM),a classical and popular supervised machine learning algorithm [52].Different nonlinear kernel functions have been integrated into SVM to solve nonlinear problems,such as Gaussian radial basis function,polynomial function,and sigmoid function[53].While in this study,the linear kernel was found to be the most efficient and was thus selected.The objective of the SVR model is to find an optimal hyperplane that fits the data.Besides,a flexible boundary is formed symmetrically around the hyperplane.Within this boundary,the estimation errors are tolerant during model training [54].The width of this boundary can be predefined or optimized by the analyst,referring to the hyperparameter ε.Conversely,the estimations outside this boundary are penalized by a regularization term C.In this study,both ε and C were optimized for each model.

    Regression models were established separately for WS and WW conditions.Under each water regime,70% of the 360 genotypes were randomly selected as the training set (i.e.,252 genotypes),and 30% of the genotypes were used as the testing set (i.e.,108 genotypes).Given that each treatment was replicated twice,the resulting training sample size was 504,and the testing sample size was 216.The genotypes assigned to the testing set were consistent between two water regions.Model hyperparameters were optimized through a ten-fold cross-validation nested in a grid search process,using the training samples only.Then,the model with its optimized hyperparameters was further evaluated on the stand-alone testing set.Three evaluation criteria,i.e.,the coefficient of determination (R2) (Equation (11)) and the root mean squared error (RMSE) (Equation (12)),and relative RMSE (RRMSE,Equation (13)) were computed between the model predicted and the measured biomass.

    Where yiandare the actual and predicted values for ith sample;is the mean value of n total number of samples.

    2.6.Yield-based drought tolerance indices

    Selection criteria for drought-tolerant genotypes vary from case to case.In this study,we selected GMP and YSI as drought tolerance indices,computed by equations(14)[6]and(15)[55],respectively.These formulae were originally established based on grain yield production.In this study,however,as biomass production is the primary trait of bioenergy sorghum,grain yield terms were replaced with biomass.Accordingly,we renamed these two indices as ‘biomass-GMP’ and ‘biomass-YSI’,respectively.We understand that the final decision on whether a genotype is drought tolerant or not depends on individual breeders and breeding purposes in each breeding program.In the present study,we considered the top-ten genotypes as drought tolerant ones by ranking the genotypes based on biomass-GMP and biomass-YSI indices.

    In addition,to evaluate the biomass predicted from UAVderived canopy features,biomass-GMP and biomass-YSI computed using the predicted biomass(i.e.,all the canopy features were used as predictors) were compared with those computed using ground truth biomass.

    Where BiomassWWis the biomass under well-watered conditions;BiomassWSis the biomass under water-stressed condition.

    3.Results

    3.1.Statistical analysis for fresh and dry biomass

    The biomass was statistically different (α=0.001) due to treatment,genotype,but not the interaction between treatment and genotype (Table 1).The effects of water deficit on fresh and dry biomass were further compared using mean values and CV within each water regime and the percent reduction (PR) under WS condition relative to WW condition(Table 1).In general,biomass production was reduced under the water-stressed environment.The reduction in the fresh biomass (32.99% PR) was larger than that of the dry biomass (10.81% PR),which probably indicates that the water deficit had a significant impact on plant water content.Using the output of ANOVA,we computed the variance components and the broad-sense heritability of fresh and dry biomass(Table 1).Results showed that the heritability of fresh biomass was higher than that of dry biomass.

    Table 1 Mean value,coefficient of variation (CV),percent reduction (PR),effect of source of variable in ANOVA,and variance component analysis for fresh and dry biomass.

    3.2.Correlation between individual UAV derived canopy feature and biomass

    The correlations (i.e.,correlation coefficient r) between UAVderived canopy features and biomass under each water regime were explored (Fig.2).In general,the correlations under WS conditions were stronger than those under WW conditions,regardless of fresh or dry biomass.For the three vegetation indices,the correlations on August 5 were all stronger than those on August 28,regardless of the water regimes or biomass components.Among all canopy features derived from two dates of UAV data collections,RDVI under WS conditions on August 5 had the strongest correlation with fresh biomass(r=0.70)and with dry biomass(r=0.69).On the other hand,the plant height under WW conditions on August 5 and NRCT under WW conditions on August 28 showed nonsignificant correlations with fresh biomass (α=0.001).For dry biomass,NRCT under WW conditions on both August 5 and 28,as well as the canopy cover from the WS block on August 5,showed nonsignificant correlations.

    3.3.Biomass predicted by combined canopy features

    To predict fresh and dry biomass from UAV-based canopy features,the linear-SVR model was trained with 70%of the 360 genotypes and was tested on the remaining 30% genotypes under two water regimes.Prediction results were given in Fig.3 for fresh biomass and in Fig.4 for dry biomass.According to the dates of data collection,three sets of input variables were compared: (1) combined six canopy features(NRCT,plant height,canopy cover,NDVI,RDVI,and NDRE)from August 5(Fig.3A,D,and 4A,D);(2)combined six canopy features from August 28 (Fig.3B,E,and 4B,E);(3) integrating (1) and (2) (Fig.3C,F,and 4C,F).Overall,the model performed better under WS (R2=0.47-0.53 for fresh biomass in Fig.3D-F;R2=0.40-0.47 for dry biomass in Fig.4D-F) than WW environment (R2=0.28-0.35 for fresh biomass Fig.3A-C;R2=0.11-0.29 for dry biomass in Fig.4A-C).Among the three sets of inputs,the one that combined the canopy features from two UAV data collections improved the estimation of the fresh biomass but had no significant improvement for dry biomass prediction.

    To further compare the contribution of different sources of aerial images,we grouped these canopy features based on their source image.The source images were: (1) thermal;(2) RGB;(3) multispectral images.The corresponding canopy features were: (1)NRCT from 5 and 28 August;(2) plant height and canopy cover from two dates;(3) NDVI,RDVI,and NDRE from two dates.The model was trained on each group of features and compared with the model that used the canopy features from all three image sources.As shown in Fig.5,except for the dry biomass prediction under WW condition,the multispectral images had the lowest RRMSE and highest R2values among the three images sources.On the other hand,the thermal images mostly performed the worst among the three image sources.Nevertheless,the model trained with all canopy features performed the best for estimating the fresh and dry biomass,which suggests that utilizing multimodal data for yield-based drought tolerance evaluation is optimal.

    Fig.2.Correlations between UAV-derived canopy features and biomass.(A)Fresh biomass.(B)Dry biomass.Canopy features were derived from August 5,2020 and August 28,2020.Blue circles represent well-watered conditions;orange circles represent water-stressed conditions.*indicates correlation is statistically significant given α=0.001.NRCT,normalized relative canopy temperature;PH,plant height;CC,canopy cover;NDVI,NDRE,and RDVI are three vegetation indices.

    Fig.3.Predicted versus measured fresh biomass from the testing dataset(n=216)under(A-C)well-watered condition and(D-F)water-stressed condition.Canopy features input into each model were derived from (A,D) August 5,2020,(B,E) August 28,2020,and (C,F) August 5 and 28,2020.

    Fig.4.Predicted versus measured dry biomass from the testing dataset (n=216) under(A-C) well-watered condition and (D-F) water-stressed condition.Canopy features input into each model were derived from (A,D) August 5,2020,(B,E) August 28,2020,and (C,F) August 5 and 28,2020.

    Fig.5.Performance of different image sources(i.e.,thermal,RGB,and multispectral)for(A,C)fresh biomass and(B,D)dry biomass prediction.R2 and RRMSE values are based on the testing dataset (n=216).Multimodal indicates the combination of all three types of images.

    3.4.Yield-based drought tolerance indices

    In this section,we extended the biomass prediction to the drought tolerance assessment using two well-recognized yieldbased drought tolerance indices,i.e.,the biomass-GMP and the biomass-YSI indices.As shown in Fig.6,biomass-GMP and biomass-YSI were computed using predicted biomass and compared with those computed using manually measured biomass.A significantly better correlation was found for the biomass-GMP index than that for the biomass-YSI index.

    In Fig.7,we highlighted the top ten genotypes by ranking the biomass-GMP and biomass-YSI calculated from both measured and predicted biomass.An assumption we made here was that the higher biomass-GMP or biomass-YSI value corresponds to the better drought tolerance.Results showed that using the biomass-GMP index,four out of the top ten genotypes were selected by the index derived by predicted biomass and the one derived by manually measured biomass were in common (Fig.7A,B);while two or three out of ten genotypes were in common between biomass-YSI indices calculated from predicted and manually measured biomass(Fig.7C,D).These findings(Figs.6,7)indicated that the biomass-GMP drought tolerance index might be the one that can be estimated more reliably from remotely sense data.

    4.Discussion

    4.1.Capability and limitation of UAV-based remote sensing data for biomass prediction and drought tolerance assessment of bioenergy sorghum

    In this study,we investigated the potential of UAV-based multimodal remote sensing data (including thermal,RGB,and multispectral data) for fresh and dry biomass predictions and drought tolerance evaluation over a large group of bioenergy sorghum genotypes and systematically compared the results with those obtained from single data sources.Multimodal data has been used previously for yield predictions for other crops such as soybean[56] and maize [57],but to the best of our knowledge,this study is the first for a large number of sorghum genotypes.In our study,multimodal data did not result in noticeably better predictions than the single-source data for both fresh and dry biomass.Results showed that the multimodal data had comparable or slightly better performance than the multispectral data in this study (Fig.5).Yet,in the study on soybean yield prediction [56],the multimodal data had consistently superior prediction performance than single modal data.The better performance obtained using multimodal data might be due to the fewer soybean varieties (i.e.,three varieties),or different algorithms employed,or the addition of texture features to spectral,structure,and thermal features.We noticed that the correlation between plant height and sorghum fresh biomass was 0.41 in the current study,which is not as high as in a similar study with only 24 bioenergy sorghum genotypes(r=0.79) [58].A major challenge in our study was the large amount of genotypic diversity(360 genotypes).On the other hand,we published a similar finding regarding the performance ranking of three UAV-based sensors in the soybean study[56].Namely,the multispectral sensor contributed the most to the estimation of sorghum biomass.RGB was also important,but its performance in modeling was inferior to multispectral images in most cases.The thermal sensor had the worst performance.Nevertheless,with relatively better performance for fresh biomass prediction identified using multiple image sources (R2=0.53 and RRMSE=23.57%),and given the inconsistency among different studies,collecting multimodal data to start with the investigation may still be a more reliable way to ensure better biomass prediction accuracy for yield-based drought tolerance assessment in breeding studies.

    Fig.6.Drought tolerance indices calculated from predicted versus measured FB (fresh biomass) and DB (dry biomass).(A,B) Biomass-GMP.(C,D) Biomass-YSI.

    Fig.7.Top ten genotypes selected by (A,B) biomass-GMP and (C,D) biomass-YSI calculated from measured versus predicted biomass.The x-axis label is the pseudo-ID number for each genotype.The length of the vertical color bar represents the magnitude of biomass-GMP or biomass-YSI value.FB is fresh biomass;DB is dry biomass.

    In the applications that a drought tolerance index is preferred,this study demonstrated the possibility of predicting the biomass-GMP index from the remotely sensed data (Fig.6).The UAV predicted fresh and dry biomass was used to compute two drought tolerance indices,namely,biomass-GMP and biomass-YSI.The results (Figs.6,7) suggested that the UAV predicted biomass was more reliable when using biomass-GMP as an index of drought tolerance.The capability to obtain a rapid and reliable estimation of biomass-GMP from the UAV multimodal data is encouraging because biomass-GMP is one of the most widely used indexes to identify elite genotypes in drought tolerance screening programs[7-9,59-61].On the other hand,the biomass predicted in the current study appeared to be less reliable for derivation of biomass-YSI.Future UAV-based remote sensing studies will be needed to improve the biomass prediction for better estimation of biomass-YSI,but the multimodal approach used in this study shows great promise for the estimation of biomass-GMP.

    Admittedly,the UAV-based remote sensing data applied in this study also had certain limitations.As the results shown in Figs.3,4,the overall performance under WS condition was better than that under WW condition.This coincides with the generally stronger correlations between canopy features and biomass under the WS conditions than those under the WW conditions(Fig.2).One factor that possibly impacts the UAV data’s performance under WW environment is the saturation of some canopy features,such as the NDVI.The canopy under WW becomes dense quicker than that under WS.However,the broadband reflectance that was sensed in this study mainly comes from the top canopy layer.As the biomass kept accumulating under water sufficient environment,some of the remotely sensed canopy features would reach a plateau.By looking at the distribution of each canopy feature (Fig.S3),values of NDVI under WW condition had less variation and almost reached one.The genotypic variation in NDVI was not enough to explain the genotypic variation in the biomass.Previous studies have suggested that using vegetation indices based on the red edge band,such as NDRE,is better than NDVI in terms of saturation issue [62].In this study,however,the correlations between NDRE and fresh/dry biomass were not as good as those between NDVI and fresh/dry biomass,even though NDRE showed less saturation(Fig.2).This may be due to the complication of unknown factors related to genotypic differences.

    4.2.Perspectives on the use of UAV-equipped thermal camera in yieldbased drought tolerance screening

    The utility of UAV-equipped miniaturized thermal cameras in agriculture has been demonstrated in a number of studies that have estimated plant water status to detect water deficit [28-31,46,50,63,64].Most of these studies used thermal data to quantify instantaneous plant water status,usually with a single or a few genotypes.Based on these studies,the use of thermal indices appears to be useful for measuring short-term responses,such as stomatal conductance and transpiration.For example,the crown temperature within a citrus orchard was found to be well correlated with manually measured stomatal conductance (r=0.88)and plant water potential (r=0.58) [30].In a study of cotton[63],the crop water stress index (CWSI) was derived from a UAV thermal map and correlated well with stomatal conductance(r=0.81).However,the good correlations between canopy temperature and instantaneous traits may not be replicable in a drought tolerance screening program that targets final yield or biomass.In a study assessing genetic variation of 120 recombinant inbred lines of Setaria grass under drought stress,correlations ranging from-0.32 to-0.49 were observed between canopy temperature and total biomass [65].In our study,negative but relatively weak correlations were found between canopy temperature and biomass,with absolute r value ranged from 0.09 to 0.44(Fig.2).Moreover,the comparative results among the three image sources from this study show that the thermal images performed the worst for modeling biomass.

    Under drought,the most common early response of plants is stomatal closure [66].Stomatal closure also may vary according to a diurnal cycle with a midday reduction in transpiration occurring due to elevated temperatures[67].This may partly explain the lower correlation between canopy temperature and biomass for the drought stressed genotypes as compared to the stronger correlations with NDVI,NDRE,and RDVI.However,canopy temperature is not solely regulated by stomatal conductance and may not be perfectly linked to transpiration rates in sorghum under drought[67].Canopy temperature is also affected by air temperature,solar radiation,wind speed,etc.[68].In addition,drought responses are not only crop species specific but also vary amongst genotypes of a single species [69].Although midday reductions in stomatal conductance occur and variation in genotype responses were expected,even with a large number of sorghum genotypes there was a significantly negative correlation between canopy temperature and biomass under drought conditions(Fig.2).An earlier study evaluated 300 sorghum genotypes under irrigated conditions[21] and showed no strong correlation between the leaf temperature and grain yield,which is similar to our results under well-watered conditions.The earlier study [21] also validated the simulation models and the hypothesis [67] that higher canopy temperatures in sorghum may not be tightly linked to lower yields due to stomatal closure that leads to conservation of water resulting in higher biomass.Therefore,the significant negative correlations observed in our dataset between canopy temperature and biomass may have been less than that observed for the spectral indices because some high biomass genotypes conserved water by stomatal closure as hypothesized [67] and shown previously[21].In this regard,we suggest that the use of thermal sensors alone may not be sufficient to identify the most drought tolerant sorghum lines in a screening program.Additional sources of phenotypic information are needed to fully elucidate the complexities in genotypic responses for drought tolerance in sorghum.

    From the sensing perspective,another possible factor that limits the performance of the thermal index(i.e.,NRCT)in estimating sorghum biomass is the sensitivity and accuracy of the thermal camera.In our study,the thermal camera we used was an uncooled camera,which is suitable for UAV applications due to its compact size and lightweight.However,compared to the cooled thermal cameras that use cooling systems to reduce measurement noise,the uncooled type is more likely to be affected by the temperature of the camera body,as well as the environmental conditions such as air temperature,relative air humidity,wind speed,cloud cover,etc.[70,71].Furthermore,according to the manufacturer,the estimated temperature accuracy of the thermal camera used in this study was around ± 5 °C under ideal environmental conditions.This accuracy may not be sufficient in a screening or breeding program that uses small plot designs,where the required accuracy may be greater than 1 °C [72].The measurement sensitivity to weather conditions,along with the relatively low sensing accuracy,may have reduced the performance of the thermal sensor used in this study.Future studies will be required to investigate whether the major limitation is because of different genotypic response rates to canopy temperature and biomass production or because of the measurement accuracy of this thermal camera.For the latter question,a thermal camera with higher measurement accuracy,such as ICI 8640P thermal camera that has an accuracy of ± 1 °C[72],maybe an alternative for the UAV system.

    4.3.Future improvements

    In our study,images including six spectral bands were used for sorghum biomass prediction and drought tolerance assessment:blue,green,red,red edge,near-infrared,and thermal infrared.The combination of this information provided acceptable performance for evaluating drought tolerance which was benchmarked with biomass accumulation.However,given the complexity in genotypic responses under drought stress,future work is required to facilitate the UAV-based remote sensing in drought tolerance selection.First,additional improvement will likely be gained if more sources of data could be included,such as the information from the hyperspectral imagery that has hundreds of continuous spectral bands.With the development of sensing technologies,UAV systems can now carry miniaturized spectrometers [73] that have been used for yield or biomass estimation [74,75],crop disease detection [76],and water stress detection [30].In the future,UAV-based point spectrometers or hyperspectral cameras will likely improve our ability to more precisely evaluate the drought tolerance of a large population of plant genotypes.Another future improvement is the collection of time-series data.In the present study,the canopy temperature index was measured at only two time points and thus provided only a snapshot in time.Therefore,it does not provide an integrated canopy assessment that reveals a long-term response to water deficit (i.e.,biomass) [33].Instead of single or a few snapshots of crop traits,multitemporal crop traits may be more beneficial for predicting biomass which is an accumulated product over the growth season.In addition,intensive data collection over the course of crop growth allows the researcher to pinpoint the critical growth stages for final biomass estimation,which will be helpful when frequent data collection is not practical.

    5.Conclusions

    This study examined the use of UAV-based multimodal images(i.e.,thermal,RGB,and multispectral) for estimating biomass and assessing drought tolerance of bioenergy sorghum.Biomass predictions using canopy features derived from the multimodal data showed comparable performance with the best results obtained from the single modal data with R2up to 0.57 under the WS condition where the canopy was less dense.Vegetation index RDVI derived from the multispectral data had the highest correlations with biomass and had a major contribution in the predictions when multimodal data were used.RGB contributed less than the multispectral images but was also important.Thermal images alone had the lowest prediction power,indicating the limitation of using this thermal camera alone for this particular drought tolerance screening project.These findings were similar to some previous studies but different than some others;hence,it indicated the advantage of starting with multimodal data under unknown conditions.Finally,for drought tolerance assessments using yield-based indices,it turned out that the biomass-GMP index outperformed the biomass-YSI index in terms of the potential to be directly and more reliably predicted by the remotely sensed data.Follow-up improvement in biomass prediction within a breeding context based on low-altitude remote sensing is promising,by enriching the available data sources through the use of timeseries data,hyperspectral imagery,or by incorporating a higher accuracy thermal camera.Ultimately,UAV systems will be an essential tool for drought tolerance screening programs.In this case,instead of destructively sampling biomass or yield production from all experimental plots,subsampling can be done together with UAV mapping to greatly reduce cost and decrease time commitment.

    Data availability

    The UAV image data has been published on Dryad and is available at https://datadryad.org/stash/share/fxAffx2rXRHVdIyGQR 9RUV2GMq0epPMuW6wqHtxsPqc.

    Declaration of competing interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    CRediT authorship contribution statement

    Jiating Li:Conceptualization,Data curation,Formal analysis,Investigation,Methodology,Software,Visualization,Writing -original draft,Writing -review &editing.Daniel P.Schachtman:Conceptualization,Funding acquisition,Data curation,Investigation,Methodology,Project administration,Resources,Writing -review &editing.Cody F.Creech:Resources,Data curation,Writing -review &editing.Lin Wang:Investigation,Writing -review&editing.Yufeng Ge:Methodology,Writing -review &editing.Yeyin Shi:Conceptualization,Funding acquisition,Methodology,Project administration,Resources,Supervision,Writing -review&editing.

    Acknowledgments

    This research was funded by US Department of Energy,BER(DE-SC0014395 to DPS),a USDA-NIFA Grant (2021-67021-34417),and the Nebraska Agricultural Experiment Station through the Hatch Act Capacity Funding Program(1011130)from the USDA National Institute of Food and Agriculture.We would like to thank Allyn Pella for her help in coordinating the data,thank Stephanie Futrell for help in ground truth data collection and lab analysis,and thank Dr.Xin Qiao and his team for their help during UAV data collections.

    Appendix A.Supplementary data

    Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.04.005.

    精品久久久久久久久亚洲| 成人二区视频| 国产伦理片在线播放av一区| 在线观看免费日韩欧美大片 | 只有这里有精品99| 色哟哟·www| 国产在视频线精品| 丝瓜视频免费看黄片| 日韩国内少妇激情av| 精品人妻偷拍中文字幕| 777米奇影视久久| 18禁裸乳无遮挡免费网站照片| 晚上一个人看的免费电影| 国产伦理片在线播放av一区| 国产精品无大码| 黑人猛操日本美女一级片| 一二三四中文在线观看免费高清| 九九爱精品视频在线观看| 免费看光身美女| 2022亚洲国产成人精品| 在线观看国产h片| 国产精品99久久99久久久不卡 | 建设人人有责人人尽责人人享有的 | 内射极品少妇av片p| 日产精品乱码卡一卡2卡三| 国产亚洲欧美精品永久| 国产片特级美女逼逼视频| 老女人水多毛片| 国产精品一区二区在线观看99| 免费观看a级毛片全部| 久久ye,这里只有精品| 下体分泌物呈黄色| 国产又色又爽无遮挡免| 只有这里有精品99| av在线播放精品| 热re99久久精品国产66热6| 亚洲精华国产精华液的使用体验| 自拍偷自拍亚洲精品老妇| 国产亚洲最大av| 日韩 亚洲 欧美在线| 日韩视频在线欧美| 男女啪啪激烈高潮av片| 国产精品精品国产色婷婷| 国产黄色视频一区二区在线观看| 国产精品不卡视频一区二区| 亚洲熟女精品中文字幕| av国产精品久久久久影院| 色视频www国产| 尾随美女入室| 美女主播在线视频| 久久 成人 亚洲| 午夜精品国产一区二区电影| 欧美日韩综合久久久久久| 最新中文字幕久久久久| 亚洲国产最新在线播放| 亚洲国产毛片av蜜桃av| 两个人的视频大全免费| 国产高清国产精品国产三级 | 国产精品秋霞免费鲁丝片| 久久人人爽av亚洲精品天堂 | 亚洲av.av天堂| 嫩草影院入口| 亚洲av成人精品一区久久| 久久人妻熟女aⅴ| 18禁裸乳无遮挡动漫免费视频| 高清欧美精品videossex| 又黄又爽又刺激的免费视频.| 你懂的网址亚洲精品在线观看| 日韩大片免费观看网站| 久久久色成人| 国产精品人妻久久久久久| 亚洲精品亚洲一区二区| 91在线精品国自产拍蜜月| 国产v大片淫在线免费观看| 岛国毛片在线播放| 亚洲av男天堂| 国产精品99久久99久久久不卡 | 成人一区二区视频在线观看| 午夜福利在线在线| 久久人人爽av亚洲精品天堂 | 成人国产av品久久久| 成人漫画全彩无遮挡| 少妇熟女欧美另类| 黄色配什么色好看| 久久久久视频综合| 麻豆乱淫一区二区| 大香蕉久久网| 一本色道久久久久久精品综合| 春色校园在线视频观看| a级一级毛片免费在线观看| 久久热精品热| 日本av手机在线免费观看| 干丝袜人妻中文字幕| 亚洲国产日韩一区二区| 日韩三级伦理在线观看| 麻豆乱淫一区二区| 亚洲欧美日韩东京热| 精品国产乱码久久久久久小说| 免费观看无遮挡的男女| 国产日韩欧美亚洲二区| 午夜福利网站1000一区二区三区| 少妇裸体淫交视频免费看高清| 一级av片app| 少妇高潮的动态图| 亚洲国产欧美在线一区| 性色av一级| 免费看光身美女| 成人国产麻豆网| 插阴视频在线观看视频| 国产色婷婷99| 久久国产乱子免费精品| 婷婷色综合www| 国产日韩欧美在线精品| 亚洲欧美精品专区久久| 高清不卡的av网站| 国产成人a区在线观看| 日韩在线高清观看一区二区三区| 日韩中字成人| 男人和女人高潮做爰伦理| av专区在线播放| av在线观看视频网站免费| 在线观看美女被高潮喷水网站| 免费看日本二区| 久久女婷五月综合色啪小说| 免费av中文字幕在线| 久久精品久久久久久久性| 九草在线视频观看| 中文精品一卡2卡3卡4更新| 国产人妻一区二区三区在| kizo精华| 国内揄拍国产精品人妻在线| 插阴视频在线观看视频| 国产伦理片在线播放av一区| 中国三级夫妇交换| 成年免费大片在线观看| 国产有黄有色有爽视频| 麻豆国产97在线/欧美| 午夜激情福利司机影院| 国产黄片美女视频| 精品久久国产蜜桃| 午夜精品国产一区二区电影| 欧美精品人与动牲交sv欧美| 最近中文字幕高清免费大全6| 国产深夜福利视频在线观看| 男女边摸边吃奶| 综合色丁香网| 免费黄频网站在线观看国产| 日产精品乱码卡一卡2卡三| 亚洲欧美日韩卡通动漫| 精品一区二区三卡| 黄片wwwwww| 夜夜骑夜夜射夜夜干| 国产在线男女| 黄片wwwwww| 久久国内精品自在自线图片| 久久av网站| 日韩精品有码人妻一区| 亚洲国产精品成人久久小说| 不卡视频在线观看欧美| 中文精品一卡2卡3卡4更新| 建设人人有责人人尽责人人享有的 | 观看av在线不卡| 久久午夜福利片| 久久婷婷青草| 精品一区二区免费观看| 麻豆乱淫一区二区| 日本-黄色视频高清免费观看| 国产精品爽爽va在线观看网站| 纯流量卡能插随身wifi吗| 国产精品一区二区在线观看99| 日韩av免费高清视频| 国产亚洲一区二区精品| 肉色欧美久久久久久久蜜桃| 少妇裸体淫交视频免费看高清| 国产日韩欧美亚洲二区| 哪个播放器可以免费观看大片| 欧美成人一区二区免费高清观看| 人妻一区二区av| 亚洲国产av新网站| 国产精品一区www在线观看| 99国产精品免费福利视频| 国产成人91sexporn| 狂野欧美白嫩少妇大欣赏| 精品酒店卫生间| 国产男人的电影天堂91| 国产欧美日韩一区二区三区在线 | h日本视频在线播放| 777米奇影视久久| 国产精品久久久久久精品古装| 老司机影院毛片| 97在线视频观看| 中文字幕制服av| 天堂中文最新版在线下载| 丰满迷人的少妇在线观看| 亚洲图色成人| 国产午夜精品久久久久久一区二区三区| 99久久综合免费| 久久这里有精品视频免费| 欧美日韩精品成人综合77777| 国产精品嫩草影院av在线观看| 亚洲性久久影院| 亚洲精品乱久久久久久| 内地一区二区视频在线| 午夜福利高清视频| 少妇猛男粗大的猛烈进出视频| 亚洲三级黄色毛片| 久久久成人免费电影| 亚洲人与动物交配视频| 久久久久久久久久久免费av| 国产综合精华液| 亚洲人成网站高清观看| 国内揄拍国产精品人妻在线| 色视频www国产| 亚洲电影在线观看av| 偷拍熟女少妇极品色| 久久久a久久爽久久v久久| 国产精品成人在线| 成年人午夜在线观看视频| 岛国毛片在线播放| 欧美一区二区亚洲| 精品一品国产午夜福利视频| 欧美bdsm另类| 人妻少妇偷人精品九色| 精品国产露脸久久av麻豆| 国产精品女同一区二区软件| 日本欧美视频一区| 少妇被粗大猛烈的视频| 人妻系列 视频| 最近中文字幕2019免费版| 亚洲精品国产色婷婷电影| 老熟女久久久| 黄片wwwwww| 亚洲人与动物交配视频| 在线看a的网站| 高清不卡的av网站| 高清视频免费观看一区二区| 精品久久久久久久久av| 免费黄频网站在线观看国产| 日韩制服骚丝袜av| 国产一级毛片在线| 亚洲精品视频女| 日韩欧美精品免费久久| 欧美激情国产日韩精品一区| 搡老乐熟女国产| 丝袜喷水一区| 少妇人妻久久综合中文| 成年av动漫网址| 亚洲激情五月婷婷啪啪| 我要看黄色一级片免费的| 伊人久久精品亚洲午夜| 欧美一区二区亚洲| 18+在线观看网站| 国产精品国产三级国产av玫瑰| 美女高潮的动态| 大话2 男鬼变身卡| 在线观看免费高清a一片| 精品一区二区三卡| 一本色道久久久久久精品综合| 香蕉精品网在线| 日韩av不卡免费在线播放| 丰满少妇做爰视频| 成人亚洲精品一区在线观看 | 联通29元200g的流量卡| 我的女老师完整版在线观看| 亚洲av在线观看美女高潮| 免费高清在线观看视频在线观看| 国产精品久久久久久精品电影小说 | 久久久久人妻精品一区果冻| 秋霞伦理黄片| 又粗又硬又长又爽又黄的视频| 午夜精品国产一区二区电影| 日本午夜av视频| 乱系列少妇在线播放| 男人和女人高潮做爰伦理| 国产成人精品福利久久| 在线观看美女被高潮喷水网站| 成年女人在线观看亚洲视频| 国产男人的电影天堂91| 久久久久久久亚洲中文字幕| 日韩精品有码人妻一区| 中文字幕制服av| 亚洲真实伦在线观看| 国内精品宾馆在线| av.在线天堂| 国产成人a∨麻豆精品| 亚洲精华国产精华液的使用体验| 97在线视频观看| 99热这里只有是精品50| 十分钟在线观看高清视频www | 国产日韩欧美在线精品| 一本色道久久久久久精品综合| 精品99又大又爽又粗少妇毛片| 久久99精品国语久久久| 久久久久国产精品人妻一区二区| 交换朋友夫妻互换小说| 最近手机中文字幕大全| 黄色日韩在线| 国产av国产精品国产| 国产伦理片在线播放av一区| 精品久久久久久久久av| 精品视频人人做人人爽| 大码成人一级视频| 99热这里只有是精品50| 国产精品人妻久久久影院| 亚洲人成网站在线播| 久久99蜜桃精品久久| 欧美xxxx性猛交bbbb| 日本黄大片高清| 午夜老司机福利剧场| 岛国毛片在线播放| 日日摸夜夜添夜夜爱| 人人妻人人爽人人添夜夜欢视频 | 国产免费视频播放在线视频| 国产一区二区三区综合在线观看 | 亚洲成人av在线免费| 校园人妻丝袜中文字幕| 一个人看视频在线观看www免费| 国产亚洲午夜精品一区二区久久| 成人一区二区视频在线观看| 国产男女内射视频| 成人一区二区视频在线观看| 大码成人一级视频| 男人添女人高潮全过程视频| 国产男人的电影天堂91| 大又大粗又爽又黄少妇毛片口| 国产成人91sexporn| 久久久久久人妻| 亚洲熟女精品中文字幕| 欧美日韩亚洲高清精品| 少妇丰满av| 欧美3d第一页| 亚洲天堂av无毛| 国产69精品久久久久777片| 亚洲精品自拍成人| 亚洲精品成人av观看孕妇| 精品人妻偷拍中文字幕| 亚洲高清免费不卡视频| 国产午夜精品一二区理论片| 伦理电影免费视频| 亚洲av成人精品一区久久| 日产精品乱码卡一卡2卡三| 国产精品熟女久久久久浪| 国产熟女欧美一区二区| 六月丁香七月| 国产av精品麻豆| 一区二区三区免费毛片| 少妇高潮的动态图| 亚洲国产欧美人成| 成年av动漫网址| 亚洲欧美日韩东京热| 欧美人与善性xxx| 中国三级夫妇交换| 日韩一区二区三区影片| 精品国产一区二区三区久久久樱花 | 欧美老熟妇乱子伦牲交| 交换朋友夫妻互换小说| 日韩中文字幕视频在线看片 | 欧美+日韩+精品| 蜜臀久久99精品久久宅男| 国产国拍精品亚洲av在线观看| 最近手机中文字幕大全| 久久精品国产a三级三级三级| 国产欧美亚洲国产| 亚洲美女黄色视频免费看| 在线观看国产h片| 在线免费观看不下载黄p国产| 夜夜爽夜夜爽视频| 天堂中文最新版在线下载| 国产成人免费观看mmmm| 午夜福利视频精品| 狠狠精品人妻久久久久久综合| 不卡视频在线观看欧美| 国产真实伦视频高清在线观看| 亚洲精品国产成人久久av| 色婷婷久久久亚洲欧美| tube8黄色片| 国产av码专区亚洲av| 色吧在线观看| 亚洲精品国产av成人精品| 涩涩av久久男人的天堂| 大陆偷拍与自拍| 久久99热这里只有精品18| 毛片一级片免费看久久久久| 久久久久久人妻| 亚洲av二区三区四区| 在线观看一区二区三区| 久久精品久久精品一区二区三区| 美女高潮的动态| 国产伦在线观看视频一区| 汤姆久久久久久久影院中文字幕| 日韩强制内射视频| 少妇人妻精品综合一区二区| 国产欧美日韩精品一区二区| 91久久精品电影网| 久久久国产一区二区| 国产 精品1| 免费看不卡的av| 免费人妻精品一区二区三区视频| 国产精品爽爽va在线观看网站| 99热6这里只有精品| 国产av一区二区精品久久 | 国产精品免费大片| 国产乱人偷精品视频| 国产极品天堂在线| 亚洲自偷自拍三级| 丰满少妇做爰视频| 精品亚洲成国产av| 成年美女黄网站色视频大全免费 | 伊人久久国产一区二区| 人妻制服诱惑在线中文字幕| 中文欧美无线码| 久久久久人妻精品一区果冻| 一级黄片播放器| 国产精品人妻久久久久久| 男女啪啪激烈高潮av片| 国产精品欧美亚洲77777| 国产成人aa在线观看| a级一级毛片免费在线观看| 久久精品夜色国产| 亚洲国产高清在线一区二区三| 91在线精品国自产拍蜜月| 成人午夜精彩视频在线观看| 我的女老师完整版在线观看| 国产一级毛片在线| 免费大片黄手机在线观看| 免费播放大片免费观看视频在线观看| av卡一久久| 亚洲国产毛片av蜜桃av| 男的添女的下面高潮视频| 亚洲精品国产色婷婷电影| 国产精品久久久久久精品古装| 色综合色国产| 一级毛片 在线播放| 丝袜脚勾引网站| 色5月婷婷丁香| 男女下面进入的视频免费午夜| 亚洲高清免费不卡视频| 亚洲第一av免费看| 天天躁夜夜躁狠狠久久av| 超碰av人人做人人爽久久| 熟妇人妻不卡中文字幕| 欧美日韩一区二区视频在线观看视频在线| 亚洲国产欧美在线一区| 少妇 在线观看| 99re6热这里在线精品视频| 亚洲美女视频黄频| 在现免费观看毛片| 小蜜桃在线观看免费完整版高清| 国产一区二区三区av在线| 99久久综合免费| 美女中出高潮动态图| 一区二区三区精品91| 欧美成人a在线观看| 在线观看一区二区三区激情| 国产永久视频网站| 啦啦啦中文免费视频观看日本| 亚洲不卡免费看| 我的老师免费观看完整版| 少妇的逼水好多| 亚洲精品视频女| 免费观看无遮挡的男女| 观看免费一级毛片| 亚洲av中文av极速乱| 国产av国产精品国产| 99热这里只有是精品50| 亚洲真实伦在线观看| 亚洲高清免费不卡视频| 免费人成在线观看视频色| 久久99精品国语久久久| 欧美成人一区二区免费高清观看| 中文字幕久久专区| 亚洲欧美精品专区久久| 欧美日韩国产mv在线观看视频 | 中国三级夫妇交换| 亚洲精品成人av观看孕妇| 三级经典国产精品| 国产国拍精品亚洲av在线观看| 热re99久久精品国产66热6| 99视频精品全部免费 在线| 亚洲精品国产成人久久av| 高清毛片免费看| 日韩一本色道免费dvd| 少妇猛男粗大的猛烈进出视频| 成人一区二区视频在线观看| 久久国产乱子免费精品| av不卡在线播放| 在线精品无人区一区二区三 | 亚洲国产精品国产精品| 亚洲激情五月婷婷啪啪| 王馨瑶露胸无遮挡在线观看| 久热这里只有精品99| 国内揄拍国产精品人妻在线| 国产人妻一区二区三区在| 伦理电影大哥的女人| 国产精品欧美亚洲77777| 国产久久久一区二区三区| 亚洲一区二区三区欧美精品| 欧美xxxx黑人xx丫x性爽| 亚洲欧美成人精品一区二区| 亚洲av男天堂| 国产无遮挡羞羞视频在线观看| 国产精品久久久久久久电影| 欧美一区二区亚洲| 久久久久久久亚洲中文字幕| 一区在线观看完整版| a 毛片基地| 波野结衣二区三区在线| 激情 狠狠 欧美| tube8黄色片| 五月伊人婷婷丁香| 人妻少妇偷人精品九色| 老司机影院毛片| 国产精品女同一区二区软件| 国产免费又黄又爽又色| 欧美bdsm另类| 亚洲欧美日韩另类电影网站 | 最新中文字幕久久久久| 亚洲精品成人av观看孕妇| 一级毛片黄色毛片免费观看视频| 天堂8中文在线网| 日韩电影二区| 少妇人妻久久综合中文| 色视频www国产| 亚洲av欧美aⅴ国产| 性色avwww在线观看| 在线精品无人区一区二区三 | av在线老鸭窝| 久久精品久久久久久噜噜老黄| 韩国av在线不卡| 国产免费一级a男人的天堂| 六月丁香七月| 久久精品国产亚洲av涩爱| 夫妻午夜视频| 亚洲av成人精品一区久久| 久久精品久久久久久噜噜老黄| 久久精品国产亚洲网站| 国产精品一区二区性色av| 亚洲国产色片| 欧美最新免费一区二区三区| 成人漫画全彩无遮挡| 亚洲三级黄色毛片| 美女国产视频在线观看| 极品教师在线视频| 男人舔奶头视频| 久久久久网色| 大陆偷拍与自拍| 啦啦啦视频在线资源免费观看| 97超碰精品成人国产| 国产精品人妻久久久影院| 国产在线一区二区三区精| 国产精品偷伦视频观看了| 中文乱码字字幕精品一区二区三区| av天堂中文字幕网| av国产精品久久久久影院| 欧美日本视频| 一级a做视频免费观看| 国产成人a区在线观看| 亚洲伊人久久精品综合| 亚洲国产成人一精品久久久| 欧美精品一区二区免费开放| 亚洲精品成人av观看孕妇| 欧美精品一区二区大全| 干丝袜人妻中文字幕| 日日啪夜夜撸| 最近2019中文字幕mv第一页| 国产免费福利视频在线观看| 大片免费播放器 马上看| 黄色配什么色好看| 18禁在线播放成人免费| 亚洲人成网站在线观看播放| 亚洲内射少妇av| 国产色婷婷99| 国产欧美另类精品又又久久亚洲欧美| 99久久人妻综合| 在线观看一区二区三区| 久热久热在线精品观看| 街头女战士在线观看网站| 日韩一区二区三区影片| 亚洲av中文av极速乱| 狂野欧美白嫩少妇大欣赏| 欧美精品一区二区大全| 在现免费观看毛片| 午夜福利在线观看免费完整高清在| 午夜老司机福利剧场| 亚洲成人av在线免费| 久久午夜福利片| av天堂中文字幕网| 美女内射精品一级片tv| 久久女婷五月综合色啪小说| 卡戴珊不雅视频在线播放| 九草在线视频观看| 啦啦啦视频在线资源免费观看| 日韩,欧美,国产一区二区三区| 青青草视频在线视频观看| 久久av网站| 综合色丁香网| 国产成人freesex在线| 美女xxoo啪啪120秒动态图| 在线播放无遮挡| 99九九线精品视频在线观看视频| 小蜜桃在线观看免费完整版高清| 亚洲国产欧美人成| 小蜜桃在线观看免费完整版高清| 欧美少妇被猛烈插入视频| 成人午夜精彩视频在线观看| av一本久久久久| 欧美日本视频| 国产综合精华液| 成人高潮视频无遮挡免费网站| 国产av国产精品国产| 女性生殖器流出的白浆| 美女内射精品一级片tv| 亚洲成人av在线免费| 99久久精品国产国产毛片| av在线蜜桃|