Prediction of yield and quality in medicinal plant Ligusticum chuanxiong Hort. using uncrewed aerial vehicle multispectral measurement

Preparation of standard substance solutions

Ferulic acid (FA), Z-ligustilide (Z-L) and senkyunolide A (SA) standard substances were dissolved in methanol, respectively, as the standard solution mother solutions. Then, they were diluted into standard substance solutions with concentrations of 31.92 µg mL⁻¹, 275.00 µg mL⁻¹ and 160.00 µg mL⁻¹ for quantitative analysis, respectively.

Preparation of sample solutions

The dried L. chuanxiong rhizome samples were crushed into powder and collected through a sieve with a pore. About one g of the powder was accurately weighed and put into a conical bottle with a stopper, and then 50 mL of 75% methanol solution was accurately added. Samples were ultrasonically extracted (power: 360 W, frequency: 40 kHz) for 30 min. The supernatant was taken and filtered by 0.22 µm microporous membrane, and the filtrate was taken to obtain the sample solution.

Chromatographic conditions

Quantitative analysis was performed using an UltiMate 3000 UPLC system equipped with a DAD detector. The separation of the three active components was performed at 30 °C using an Agilent ZORBAX Eclipse Plus C18 column (2.1  × 50 mm, 1.8 µm). Use 0.03% phosphoric acid water (A)-methanol (B) as the gradient solvent system. The gradient elution conditions were as follows: 0∼3 min, 15%∼25% B; 3∼5 min, 25%∼28% B; 5∼9 min, 28%∼30% B; 10∼16 min, 48%∼50% B; 16∼18 min, 50%∼80% B. The flow rate and injection volume were 0.40 mL min⁻¹ and 1 µL, respectively. The detection wavelengths were 282 nm, 322 nm and 328 nm. The chromatographic peak separation results (Fig. 2) show that the chromatographic peak separation was good.

Multispectral image processing

The spectral images of the five bands obtained from each flight mission were imported into Pix4D Mapper 4.5.6 to complete the geometric correction and radiation correction with image information from the corresponding radiation calibration panels. Then, the spectral images were spliced into orthographic mosaic reflectance images. ENVI 5.3 software was used to overlay the orthographic mosaics of each band into a multispectral image. After classifying the canopy of L. chuanxiong and the background of soil and weeds by the Support Vector Machine class function of ENVI, a mask was constructed for segmentation, and the multispectral image of L. chuanxiong was obtained (Fig. 3). By selecting the location of each field sampling area and editing the corresponding region of interest (ROI), the VIs of each ROI were extracted and calculated to obtain 23 characteristic variables (Table 1).

VIs selection and models construction

The best model input variables were screened by Pearson correlation analysis between variables and the target values of FW, DW, FA, Z-L and SA in the corresponding samples. However, there may be a high degree of correlation and collinearity among variables. Screening only the characteristic variables with high Pearson correlation coefficient may lead to data redundancy and miss some effective information (Xu et al., 2023). Therefore, according to the order of the absolute value of the correlation coefficient between the variables and the target value from large to small, the characteristic variables with the correlation coefficient greater than 0.9 were eliminated sequentially. The remaining characteristic variables with significant correlation (p < 0.05) were used as model input variables to construct a prediction model corresponding to each target value.

Model construction and related data processing were performed in Python 3.11. Four machine learning ensemble algorithms, Random Forest (RF), adaptive boosting (AdaBoost), Gradient Boosting Decision Tree (GBDT) and extreme gradient boosting (XGBoost), were used to construct regression prediction models of FW, DW, FA, Z-L and SA contents, respectively. The quality discriminant models were based on the standard set by the Pharmacopoeia of the People’s Republic of China (2020) (ChP 2020) (Chinese National Pharmacopoeia Commission, 2020) that the FA content of chuanxiong rhizoma should not be less than 0.10% (1.0 mg g⁻¹). Considering that there may be a large difference between the sample sizes that meet the criteria or not, the Border line SMOTE was used to oversample the categories with a small sample size to obtain balanced data. Each data set was randomly divided into a training set and a testing set with a ratio of 7:3 for model building. In the process of model construction, each algorithm determined the optimal model parameters through grid search and five-fold cross validation.

Three indicators were used to evaluate the effectiveness and stability of the model. Coefficient of determination (R²) is used to measure the fitting degree of the regression model, and the value was 0∼1. The larger the value, the better the fitting degree of the model. Normalized root mean squared error (NRMSE) and mean absolute percentage error (MAPE) are indicators of the gap between the predicted value and the true value. The former is more sensitive to large errors because of the root mean square calculation, while the latter weights all errors equally. Both of them can compare the effects of models composed of different data sets. The performance of the model in the independent validation set was evaluated using NRMSE and MAPE. The calculation formula is as follows (Costa et al., 2022; Zhu et al., 2023): (1)${\displaystyle {\mathrm{R}}^2=\mathrm{1}-\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}}}_{\mathrm{i}}\right)}^2/\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}^2}$ (2)${\displaystyle \text{NRMSE}=\left(1/\overline{\mathrm{y}}\right)\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}}}_{\mathrm{i}}\right)}^2/\mathrm{n}}\cdot \mathrm{100}\text{\%}}$ (3)${\displaystyle \text{MAPE}=\left(1/\mathrm{n}\right)\sum _{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}}}_{\mathrm{i}}\right|/{\mathrm{y}}_{\mathrm{i}}\cdot \mathrm{100}\text{\%}}$

where n is the sample size of the input data, y_i is the true target value, $\overline{y}$ is the average value of the corresponding true target value, and is ${\hat{y}}_i$ the predicted target value calculated by the corresponding regression model.

The performance of the quality discriminant model was evaluated using the following indicators: Accuracy is used to evaluate the overall discriminant effect of the discriminant model. Precision, recall and F1-score derived from the confusion matrix can explain the discriminant effect in more detail. The area under curve (AUC) is the area under the receiver operating characteristic (ROC) curve, and the value is 0∼1. The larger the AUC, the better the discriminant effect of the model. The performance of the model in the independent validation set was evaluated using accuracy and AUC. The formula is as follows (Zhang et al., 2022a; Zhang et al., 2022b): (4)${\displaystyle \text{Accuracy}=\left(\mathrm{TP}+\mathrm{TN}\right)/\left(\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}\right)}$ (5)${\displaystyle \text{Precision}=\mathrm{TP}/\left(\mathrm{TP}+\mathrm{FP}\right)}$ (6)${\displaystyle \text{Recall}=\mathrm{TP}/\left(\mathrm{TP}+\mathrm{FN}\right)}$ (7)${\displaystyle \mathrm{F}\mathrm{1}-\text{score}=2\times \text{Precision}\times \text{Recall}/\left(\text{Precision}+\text{Recall}\right)}$

where TP is the number of samples whose instances are positive and predicted to be positive, TN is the number of samples whose instances are negative and predicted to be negative, FN is the number of samples whose instances are positive and predicted to be negative, FP is the number of samples whose instances are negative and predicted to be positive.

Yield prediction models

The R² (Table 2) of the four algorithm models in the FW training set and the test set data in both periods could reach more than 0.68. Overall, the model with the best prediction accuracy in both periods was XGBoost, the R² of the training set was close to 0.90, and the NRMSE and MAPE were less than 20.00%. The XGBoost over-fitting degree in D1 was lower than that in D2, whose R² on the testing set was 19.54% lower than that on the training set.

The difference between the DW prediction models (Table 3) was similar to that of FW, and the best prediction model was the XGBoost model in D1. The R² on the training and test sets were both greater than 0.80, and the NRMSE and MAPE were within the range from 19.00% to 21.00%. In the model in D2, the XGBoost model performed best in the training set with an R² of 0.82, but the over-fitting situation was still higher than other algorithm models.

Active component contents prediction models

The effects of FA, Z-L and SA content prediction models in D1 are shown in Table 4. Among the FA content prediction models in D1, the training set R² of GBDT was the largest, but the R² of the test set decreases greatly, up to 45.68%, while the AdaBoost model has a small degree of over-fitting while ensuring a low error. The degree of over-fitting of the Z-L content prediction models was high, and the AdaBoost had the lowest degree of over-fitting that the difference between the R² values of the training set and the test set was 0.24. The R² of all SA content prediction models in D1 performed well and the degree of over-fitting was low. Among them, the RF training set had the highest R² of 0.89, while the R² of the AdaBoost test set was higher and the degree of over-fitting was lower.

The effect of the three active component content prediction models in D2 is shown in Table 5. The FA and Z-L content prediction models of D2 had higher degree of over-fitting, except for AdaBoost. However, each training set R² of AdaBoost was lower, 0.35 and 0.51, respectively. The R² of the SA content prediction model in D2 was higher, and GBDT was the best, whose R² of the training set and the test set were 0.90 and 0.76, respectively. While the AdaBoost, with training set R² of 0.74 for both training and test sets, had the lowest degree of over-fitting.

Chuanxiong rhizoma quality discriminant models

ChP 2020 stipulates that the FA content in chuanxiong rhizoma should not be less than 0.10% (1.00 mg g⁻¹) (Chinese National Pharmacopoeia Commission, 2020). Therefore, we constructed quality discriminant models (Table 6) to evaluate whether the FA content reached the standard. The accuracy of XGBoost and GBDT models in D1 reached 0.9020, and the AUC was 0.9031 and 0.9095, respectively. The Accuracy of XGBoost, AdaBoost and RF constructed by D2 data was better and all of them were 0.8431. Among them, the AUC of RF was larger so the discriminative efficiency was the best.

Independent validation of yield prediction models

Since D2 and D3 are close to the harvesting periods of Shiyang Town and Aoping Town respectively, the VIs of D3 were used as the characteristic variable to be input according to the requirements of the yield prediction models in D2 to evaluate the performance of models on the independent validation set. As shown in Fig. 5, the relationship between the predicted value and the actual value is obtained. When applied to the independent validation set, each model underestimates FW and DW. Among them, RF was the best among the FW prediction models with NRMSE of 23.76%, and MAPE of 14.75%, while AdaBoost was the best among the DW prediction models with NRMSE of 34.60% and MAPE of 21.73%.

Independent validation of active component content prediction models

The prediction models of FA, Z-L and SA content in D2 were applied to predict the content of active components in the independent validation set, and the results were shown in Fig. 6. AdaBoost performs best in predicting FA content. Compared to the other models, the RF has a more concentrated distribution of the predicted values and smaller relative errors. The predicted Z-L content of each model was generally lower than the true values, among which the RF’s predicted values were the most concentrated with NRMSE of 34.35% and MAPE of 23.40%. The predicted SA content of each model was generally lower than the true values. Among them, XGBoost was the best at predicting with NRMSE of 45.26% and MAPE of 25.48%.

Independent validation of chuanxiong rhizoma quality discriminant models

To verify the effect of the constructed quality discriminant model based on the FA content standard of chuanxiong rhizoma stipulated in ChP 2020, the corresponding characteristic variables were entered into the quality discriminant model of D2, and the discriminant results were shown in Fig. 7. The best performing models were XGBoost and RF, with accuracy and AUC of 0.7083 and 0.7214, respectively, which had effective discriminative ability.