Huangqi state recognition method and system

Abstract

The application discloses a method and system for identifying the state of Astragalus membranaceus. The method comprises the following steps: acquiring a surface image of healthy Astragalus membranaceus and a surface image of Astragalus membranaceus to be identified, extracting a local chroma contrast vector from a healthy area and fitting a bivariate Gaussian mixture model to establish a chroma contrast distribution model of the current production batch; performing multi-scale local chroma contrast coding and directional gradient coding on the surface image of Astragalus membranaceus to be identified, counting the joint occurrence of the chroma contrast mode and the directional gradient mode and performing multi-scale fusion to obtain a joint feature vector, and determining the state of Astragalus membranaceus according to the joint feature vector. The joint feature vector can also be used to form a visual fingerprint and be associated with production batch data to perform causal diagram learning, intervention inference and counterfactual prediction.

Description

Technical Field

[0001] This invention relates to the field of machine vision inspection technology for Chinese medicinal materials, specifically to a method and system for recognizing the state of Astragalus membranaceus. Background Technology

[0002] In the processing, sorting, and quality control of Astragalus membranaceus, it is usually necessary to identify surface defects such as black spots and cracks, and further complete grading and batch recording. Manual visual inspection is easily affected by the experience of the inspectors, the work cycle, and on-site conditions. Therefore, using machine vision to acquire surface images of Astragalus membranaceus and perform automatic condition determination is a technical approach to improve sorting consistency and processing efficiency.

[0003] Existing methods for inspecting the appearance of agricultural products often employ artificial features such as color histograms, local binary patterns, oriented gradient histograms, and gray-level co-occurrence matrices, or classification models built based on training samples. Local binary patterns mainly describe the relative relationships between gray levels, while color histograms mainly describe the global color distribution. When using only one of these features, it is difficult to simultaneously express local color variations and texture orientation structures.

[0004] For Astragalus membranaceus, healthy samples from different origins, seasons, and production batches may exhibit natural color variations. When the color quantification boundary or statistical benchmark remains fixed after system deployment, normal color variations between different production batches may be mistakenly identified as defective color variations, or the distinction between actual defects and healthy colors may be reduced. Therefore, it is necessary to adapt the statistical basis of color coding to the actual color distribution of the healthy region in the current production batch.

[0005] On the other hand, anomalies such as black spots, fungal patches, and cracks on the surface of Astragalus exhibit significant differences in size and directional structure. Smaller anomalies may occupy only a few pixels, while longer cracks show obvious directionality. If features are extracted only at a single spatial scale, it is difficult to take into account anomalies of different sizes; if color and gradient features are processed separately without local joint statistics, it is also difficult to describe the co-occurrence relationship between color changes and directional structure in the same local region.

[0006] During the cultivation and processing of Astragalus membranaceus, production batch data such as environmental humidity, temperature, light duration, soil pH, watering frequency, and quality defects can be generated. Existing traceability systems are typically used to record and query batch information, with further data analysis often limited to correlation statistics. When it is necessary to assess the actual effect of a certain adjustable production factor on quality defects, or to analyze the potential quality outcomes of a batch that has already experienced defects under another intervention condition, relying solely on correlation statistics is insufficient to provide corresponding conclusions. Summary of the Invention

[0007] The main technical problem solved by this invention is: given the natural color differences among different production batches of Astragalus membranaceus, how to use the current production batch of healthy Astragalus membranaceus as a statistical benchmark to jointly characterize the local color changes and directional structures of the Astragalus membranaceus surface at multiple spatial scales to obtain stable joint features for state determination; on this basis, the visual detection results can also be correlated with the production batch data to provide a data foundation for subsequent causal graph learning, intervention inference, and counterfactual prediction.

[0008] In a first aspect, the present invention provides a method for identifying the state of Astragalus membranaceus. The method involves acquiring surface images of healthy Astragalus membranaceus from the current production batch and surface images of Astragalus membranaceus to be identified; extracting a set of local chromaticity contrast vectors from the healthy regions of the healthy Astragalus membranaceus surface images; fitting this set of vectors to a bivariate Gaussian mixture model to obtain a chromaticity contrast distribution model characterizing the chromaticity contrast distribution of the healthy regions in the current production batch; performing local chromaticity contrast encoding and local directional gradient encoding on the surface images of the Astragalus membranaceus to be identified at multiple spatial scales; obtaining a chromaticity contrast pattern based on the chromaticity contrast distribution model and a directional gradient pattern based on the local directional gradient encoding; statistically analyzing the joint occurrence of the two patterns at each spatial scale and performing multi-scale fusion to obtain a joint feature vector; and determining the state determination result of the Astragalus membranaceus to be identified based on the joint feature vector.

[0009] In one implementation, the surface image of Astragalus membranaceus to be identified is decomposed into Gaussian pyramids to generate multiple scale-level images, and chromaticity contrast encoding and directional gradient encoding are performed in the local neighborhood of each scale-level image.

[0010] In one implementation, the scale-level image is converted from the RGB color space to the HSV color space, and a local chromaticity contrast vector is constructed based on the hue and saturation differences between neighboring pixels and the center pixel. The chromaticity contrast mode corresponding to the local chromaticity contrast vector is determined using a bivariate Gaussian mixture model fitted to the healthy region of the current production batch.

[0011] In one implementation, the set of local chromaticity contrast vectors of the healthy region of the current production batch is used as training data, and the mixing coefficients, mean vector and covariance matrix of the bivariate Gaussian mixture model are fitted online using the expectation-maximization algorithm; after the production batch is switched, the bivariate Gaussian mixture model is refitted according to the healthy region of the switched production batch.

[0012] In one implementation, gradient magnitude and gradient direction are calculated on images at various scale levels, the gradient direction is quantized into multiple direction intervals, a gradient direction co-occurrence matrix is ​​constructed based on the adjacency relationship of different direction intervals in the local neighborhood, and the directional gradient pattern is determined based on the gradient direction co-occurrence matrix.

[0013] In one implementation, the chroma contrast mode and directional gradient mode are used as two-dimensional statistical dimensions, and their joint occurrence in the corresponding spatial scale is statistically analyzed to obtain a two-dimensional joint histogram. The two-dimensional joint histograms at different spatial scales are normalized and concatenated in scale order to obtain a joint feature vector.

[0014] In one implementation, a hash operation is performed on the joint feature vector, and the resulting hash value is used as the visual fingerprint of the corresponding detection record. The visual fingerprint is then associated with and stored in conjunction with the status determination result, the origin environmental factor data, the cultivation management data, and the quality defect data to form production batch traceability data.

[0015] In one implementation, environmental variables, cultivation management variables, and quality defect indicators are determined based on traceability data from multiple production batches; the PC algorithm is used to determine the adjacency relationship between variables based on conditional independence tests, and the directed relationship is determined by combining domain knowledge priors; for non-root nodes in the causal graph, the set of parent nodes of the non-root node is used as input, and Gaussian process regression is used to estimate the nonlinear structural equation.

[0016] In one implementation, intervention inference is performed on the modifiable variable and the target quality defect index. By deleting directed edges pointing to the modifiable variable and setting an intervention value, sampling is performed using a nonlinear structural equation to determine the average causal effect. For the target production batch, the posterior distribution of the exogenous variable is estimated using the Markov chain Monte Carlo method, a counterfactual intervention value is set, and the counterfactual prediction result of the target quality defect index is obtained using the structural equation after intervention.

[0017] Secondly, the present invention provides an Astragalus membranaceus state recognition system, including an image acquisition module, a distribution model construction module, a feature encoding module, and a state determination module. Each module is used to perform the above-mentioned image acquisition, current batch color contrast distribution model construction, multi-scale joint feature encoding, and state determination processing, respectively.

[0018] This invention does not directly process different production batches by using fixed color judgment boundaries. Instead, it uses the local color contrast vector of healthy Astragalus membranaceus in the current production batch to establish a color contrast distribution model online, so that the color pattern encoding of the Astragalus membranaceus to be identified corresponds to the actual color contrast distribution of the healthy area in the current batch.

[0019] This invention performs chromaticity contrast encoding and directional gradient encoding at multiple spatial scales, statistically analyzes the joint occurrence of the two types of patterns, and then fuses the joint features from different spatial scales so that the joint feature vector simultaneously contains local chromaticity changes, local directional structure, and scale information. This provides a unified state judgment input for surface anomalies of Astragalus membranaceus with different sizes and directional features.

[0020] Furthermore, by associating the hash value of the joint feature vector with production batch environmental data, cultivation management data, and quality defect data, a correspondence can be established between visual inspection records and production data, thereby providing a data foundation for performing causal structure learning, intervention inference, and counterfactual prediction based on data from multiple production batches. Attached Figure Description

[0021] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments of this application and their descriptions are used to explain this application and do not constitute an undue limitation of this application.

[0022] Figure 1 This is a schematic diagram of the data flow for the joint processing of multi-scale local chromaticity contrast mode and directional gradient mode in one embodiment of the present invention.

[0023] Figure 2 This is a schematic diagram of the intervention inference process for an interventionable variable in one embodiment of the present invention.

[0024] Figure 3 This is a schematic diagram of the counterfactual inference process for a target production batch according to one embodiment of the present invention. Detailed Implementation

[0025] The present invention will be further described below with reference to specific embodiments. The device parameters, sample size, model parameters, and verification data in the following embodiments are used to illustrate the specific implementation of the technical solution and do not constitute a limitation on the scope of protection. Unless otherwise specified, image acquisition, region cropping, color space conversion, gradient calculation, model training, and data storage can all be implemented using conventional methods or devices.

[0026] In one specific implementation, an industrial camera is used to acquire RGB images of the Astragalus membranaceus surface, and the region containing the Astragalus membranaceus is cropped to obtain a Region of Interest (ROI) image. The industrial camera can use a 1920×1080 resolution and a global shutter, along with a 16mm focal length lens and a ring LED light source; the color temperature of the light source can be set to 5000K and the illuminance can be set to 5000 lux.

[0027] The ROI image size is normalized to 256×256 pixels, and a three-level Gaussian pyramid decomposition is performed on the ROI image. Adjacent scale levels are obtained by smoothing the previous scale level using Gaussian smoothing and then downsampling by a factor of 2. In a specific parameter configuration, the standard deviation of the Gaussian smoothing kernel is 1.0, resulting in three scale levels: I0, I1, and I2, with sizes of 256×256 pixels, 128×128 pixels, and 64×64 pixels, respectively.

[0028] The purpose of employing multiple scale levels is to enable the same Astragalus membranaceus surface image to form local feature representations at different spatial resolutions. Higher resolution scales retain more subtle local variations, while lower resolution scales form a larger equivalent receptive range for longer structures, and subsequently, local chromaticity and orientation features are independently formed at each scale.

[0029] Healthy Astragalus membranaceus (Huang Qi) samples, manually confirmed to be free of any defects to be detected, are selected from the current production batch, and their healthy regions are determined. For each scale level image, the RGB color space is converted to the HSV color space. Let p be the center pixel, N(p) be the k×k pixel neighborhood centered at p, and q be the neighboring pixels in N(p). A local chromaticity contrast vector is constructed based on the chromaticity difference between the neighboring pixel q and the center pixel p. Specifically, the hue component difference and saturation component difference between the neighboring pixel q and the center pixel p are calculated respectively, and a two-dimensional local chromaticity contrast vector is formed in the order of hue component difference and saturation component difference. The hue components of the neighboring pixel q and the center pixel p are denoted as Hq and Hp, respectively, and the saturation components are denoted as Sq and Sp, respectively. The local chromaticity contrast vectors obtained within the healthy regions of the current production batch are summarized to form the healthy chromaticity contrast vector set DB corresponding to the current production batch B.

[0030] When calculating the hue component difference, the minimum ring difference between two hue values ​​can be determined based on the periodic characteristics of the HSV hue components, so as to avoid similar colors on both sides of the hue period boundary being characterized as large hue differences.

[0031] A bivariate Gaussian mixture model with M Gaussian components is used to model the healthy chromaticity contrast vector set DB. For any local chromaticity contrast vector, the probability density of the vector under the bivariate Gaussian distribution corresponding to each Gaussian component is calculated. Each probability density is multiplied by the mixing coefficient of the corresponding Gaussian component, and the product results corresponding to the M Gaussian components are summed. The result is used as the overall probability density of the local chromaticity contrast vector under the bivariate Gaussian mixture model.

[0032] The parameters of the bivariate Gaussian mixture model include the mixing coefficients of each Gaussian component, a two-dimensional mean vector, and a covariance matrix. Each mixing coefficient is non-negative, and the sum of all mixing coefficients is 1. The two-dimensional mean vector represents the central position of the corresponding Gaussian component in the hue difference and saturation difference dimensions. The covariance matrix represents the degree of dispersion in the two dimensions and the correlation between them.

[0033] Using the healthy chromaticity contrast vector set DB as input, the mixing coefficients, mean vector, and covariance matrix of each Gaussian component are estimated online through the expectation-maximization algorithm. The expectation step calculates the posterior probability of each local chromaticity contrast vector to each Gaussian component based on the current model parameters; the maximization step updates each model parameter based on the posterior probability; the expectation step and the maximization step are executed alternately until the log-likelihood reaches a preset convergence condition or the number of iterations reaches a preset upper limit.

[0034] In one optional implementation, 200 healthy Astragalus membranaceus roots are randomly selected from the current production batch. The value k is set to 3, meaning the chromaticity contrast vector of the 8 neighboring pixels surrounding the central pixel is calculated using a 3×3 pixel neighborhood. The number of Gaussian components M is set to 8, and the initial parameters of each Gaussian component are initialized using the k-means++ method. The convergence condition for the expectation-maximization algorithm is set to a log-likelihood increment of less than 10^-6 or an iteration count of 500. These parameters are for illustrative purposes only.

[0035] For each local chromaticity contrast vector of Astragalus membranaceus to be identified, its posterior probability of belonging to each Gaussian component is calculated. Specifically, the mixing coefficient of the corresponding Gaussian component is multiplied by the probability density of the local chromaticity contrast vector under that Gaussian component, and then normalized using the overall probability density of the local chromaticity contrast vector under all Gaussian components to obtain the posterior probability of the local chromaticity contrast vector belonging to the corresponding Gaussian component. The Gaussian component with the highest posterior probability is determined as the chromaticity contrast pattern corresponding to the local chromaticity contrast vector. For multiple chromaticity contrast vectors obtained from the same local neighborhood, chromaticity contrast pattern codes can be formed respectively, and the joint occurrence of the chromaticity pattern and the corresponding local directional pattern is statistically analyzed in the subsequent joint statistical process.

[0036] When a production batch is switched, the surface image and healthy area of ​​healthy Astragalus membranaceus are re-acquired from the switched production batch to form a new set of healthy color contrast vectors. The online fitting of the bivariate Gaussian mixture model is then re-executed so that the subsequent color encoding is based on the healthy color distribution of the switched production batch.

[0037] In an optional implementation, the Mahalanobis distance of the local chromaticity contrast vector to be detected relative to each Gaussian component can also be used to characterize the degree of deviation of the local chromaticity contrast vector from the healthy chromaticity distribution of the current production batch. For any Gaussian component, firstly, the difference vector between the local chromaticity contrast vector to be detected and the two-dimensional mean vector of the Gaussian component is calculated. Then, the difference vector is normalized by covariance scaling using the inverse matrix of the covariance matrix of the Gaussian component. The inner product between the difference vector and the normalization result is calculated and the square root is taken to obtain the corresponding Mahalanobis distance.

[0038] Calculate the Mahalanobis distance from the local chromaticity contrast vector to each Gaussian component, and take the minimum value as the deviation of the local chromaticity contrast vector from the healthy chromaticity distribution of the current production batch.

[0039] Furthermore, for each scale level of the image, the Sobel operator is used to calculate the horizontal and vertical gradients respectively. The gradient magnitude of the corresponding pixel is determined by the square root of the sum of the squares of the horizontal and vertical gradients, and the arctangent operation, which can distinguish the four quadrants, is used to determine the gradient direction of the corresponding pixel based on the vertical and horizontal gradients.

[0040] The gradient direction θ is quantized into Nθ directional intervals. In one specific implementation, Nθ = 8, and each directional interval covers 45°. Within the local pixel neighborhood, an Nθ × Nθ gradient direction co-occurrence matrix C is constructed based on the adjacency relationship of different directional intervals, where C(i,j) represents the frequency of adjacent occurrence of pixels in directional interval i and pixels in directional interval j within the local neighborhood.

[0041] The distribution patterns of gradient orientation co-occurrence matrices corresponding to multiple local neighborhoods are clustered to form directional gradient patterns. In one optional implementation, the k-means algorithm is used to cluster the distribution patterns of the co-occurrence matrices, dividing them into 8 categories of directional gradient patterns. Each local neighborhood to be identified obtains its directional gradient pattern encoding based on the clustering result corresponding to its gradient orientation co-occurrence matrix.

[0042] To reduce the impact of posture rotation on directional features during transport or imaging of Astragalus membranaceus, rotation normalization based on the principal direction can be performed. In one specific implementation, neighboring pixels in a 3×3 neighborhood are weighted according to gradient magnitude to form a directional histogram with 8 directional intervals, and the peak direction is taken as the principal direction of the local neighborhood. The gradient direction within the neighborhood is then angularly shifted relative to the principal direction and requantized. Based on the normalized directional quantization result, a directional gradient pattern for subsequent joint statistics is formed.

[0043] At each spatial scale, a two-dimensional joint histogram is constructed using chromaticity contrast pattern and directional gradient pattern as two statistical dimensions. For the combination of chromaticity contrast pattern and directional gradient pattern appearing at a corresponding local location, the corresponding statistical units in the two-dimensional joint histogram are accumulated, thereby obtaining the joint occurrence statistics of local chromaticity change and directional structure at the current scale.

[0044] For an implementation with M chroma contrast modes and Nθ directional gradient modes, the two-dimensional joint histogram at each spatial scale has M×Nθ statistical units. In one specific parameter configuration, M=8 and Nθ=8, thus each scale corresponds to a 64-dimensional joint histogram. L2 normalization is performed on the two-dimensional joint histograms corresponding to each spatial scale, and the normalized two-dimensional joint histograms are concatenated end-to-end according to a preset spatial scale order to form the final joint feature vector. In one specific parameter configuration, the number of chroma contrast modes is 8, and the number of directional gradient modes is 8, thus each spatial scale corresponds to a 64-dimensional joint histogram; when the number of spatial scales is 3, the final joint feature vector is 192-dimensional. Normalizing at different spatial scales can reduce the impact of differences in the number of pixels at different scales on the statistical magnitude.

[0045] The joint processing flow of the above-mentioned multi-scale local chromaticity contrast mode and directional gradient mode is as follows: Figure 1 As shown.

[0046] Furthermore, the obtained joint feature vector is input into the state determination model, which outputs the state determination result of the Astragalus membranaceus to be identified. The state determination model can be a tree model classifier trained based on labeled Astragalus membranaceus samples. In an optional specific implementation, the XGBoost classification architecture can be used, with the 192-dimensional joint feature vector as the main input. It first distinguishes between good state and defective state, and can further output the black spot state, crack state, or mixed defect state according to the detection task.

[0047] During the training of the state determination model, the training samples undergo the same image preprocessing, Gaussian pyramid decomposition, local color contrast encoding, directional gradient encoding, and multi-scale joint statistical procedures as the samples to be identified. Optionally, a 5-fold cross-validation method can be used to train the XGBoost classifier, and the learning rate, maximum depth, subsampling ratio, column sampling ratio, and regularization parameters can be determined through Bayesian optimization over 50 iterations. The above training method pertains to a specific implementation.

[0048] In an extended implementation, the status determination result can also be combined with the defect degree, size parameter, and fullness parameter to form a grading result. The grading module receives the status score, defect score, size score, and fullness score, and determines the comprehensive score in the following manner: assign corresponding weights to the above four types of scores, multiply each type of score by its corresponding weight, and then add the resulting products to obtain the comprehensive score, wherein the sum of the four weights is 1.

[0049] In one alternative implementation, the weights corresponding to the status score, defect score, size score, and fullness score can be set to 0.30, 0.25, 0.25, and 0.20, respectively.

[0050] The defect score can be determined by interval classification based on the dark spot area ratio (Ddark) and the maximum crack length (Lmax). In one optional example, the defect score is 95 points when Ddark is not greater than 0.02 and Lmax is not greater than 1.0 mm; the defect score is 50 points when Ddark is greater than 0.10 or Lmax is greater than 3.0 mm, and the intermediate interval can be determined according to a preset grading rule. The size score can be determined based on the relative deviation of the measured size parameter relative to the corresponding standard size. Specifically, the absolute value of the difference between the measured size parameter and the standard size is calculated, and the absolute value is divided by the standard size to obtain the size relative deviation; the size relative deviation is subtracted from 1, and if the result is less than 0, it is treated as 0. The result is then converted to a scoring scale of 0 to 100 to obtain the size score.

[0051] Fullness can be characterized by the roundness of the Astragalus membranaceus cross-section. Cross-sectional roundness is determined by the area of ​​the cross-sectional region and the perimeter of the cross-sectional outline. Specifically, the area of ​​the cross-sectional region is multiplied by four times pi, and then divided by the square of the perimeter of the cross-sectional outline to obtain the corresponding cross-sectional roundness. The closer the cross-sectional roundness is to 1, the closer the corresponding cross-section is to a circle. The obtained cross-sectional roundness can be further converted to a scoring scale from 0 to 100 to obtain a fullness score.

[0052] After obtaining the joint feature vector and state determination result, a hash operation can be performed on the joint feature vector, and the hash value can be used as the visual fingerprint of the corresponding detection record. The visual fingerprint is used to establish the association between the visual detection record and the production batch data, and it is not required that the same hash value be obtained when re-photographing the same Astragalus membranaceus at different times.

[0053] The system links and stores visual fingerprints, status determination results, optional grading results, production environment factor data, cultivation management data, and quality defect data. This allows for the location of corresponding inspection records via visual fingerprints and the aggregation of environmental, management, and quality defect information by production batch, providing a dataset for subsequent cross-batch analysis.

[0054] In one optional implementation, a set of visual recognition verification data for Astragalus membranaceus includes 5200 Astragalus membranaceus roots, of which 3120 are good samples, 1040 are black spot samples, and 1040 are cracked samples. The training set and the test set are divided in a 7:3 ratio, and the average value is obtained from 5 independent repeated experiments.

[0055] In the verification, the F1 score of LBP plus color histogram was 89.2%±1.2%, HOG was 85.7%±1.5%, GLCM was 81.3%±1.8%, Opponent-LBP was 91.3%±1.0%, CLBP was 90.8%±1.1%, and the ResNet-50 fine-tuning scheme was 93.5%±0.8%. The multi-scale local chromaticity-direction joint feature scheme with a fixed BGMM baseline was 94.1%±0.6%, and the scheme using the online fitting of the BGMM baseline of the current production batch was 97.6%±0.4%.

[0056] In the detailed results, the recall rate for black spots was 98.5% ± 0.5%, the recall rate for minor cracks was 96.2% ± 0.7%, and the recall rate for good condition was 99.2% ± 0.3%. Furthermore, when using spring samples from Longxi as the training set and conducting cross-batch testing with autumn samples from Hubei, the F1-score for the fixed BGMM scheme was 89.2% ± 1.8%, while the online BGMM scheme achieved 96.8% ± 0.9%. These data illustrate the effectiveness of the current batch-healthy chromaticity contrast distribution model under cross-batch color variation conditions.

[0057] After accumulating traceability data from multiple production batches, environmental variables, cultivation management variables, and quality defect indicators are selected from the traceability data to form a set of variables for causal analysis. One set of production data includes records from 3 consecutive years, 12 planting bases, and 364 complete growth cycles.

[0058] Each growth cycle record can include 12-month environmental humidity, temperature, and light duration sequences, quarterly soil pH records, and watering frequency sequences; quality defect indicators such as the area ratio of dark spots (Ddark) and the number of cracks (Dcrack) at harvest. It can also include cultivation management information such as fertilization type. For implementations using conditional independence tests of continuous variables based on partial correlation, continuous environmental variables, continuous management variables, and continuous / counting defect indicators can be used as direct analysis variables in the PC algorithm. Other management information can be used for batch records, domain knowledge priors, or preprocessed using methods suitable for its data type before being included in the analysis.

[0059] In one data processing method, 12 months of time-series data are aggregated by quarter, and corresponding statistical values ​​are calculated for each quarter, so that environmental variables and cultivation management variables form time-series observations according to growth stages; the data are divided into training data and test data in an 8:2 ratio for causal graph learning and structural equation verification.

[0060] Furthermore, the PC algorithm is used to determine the adjacency relationships between variables from the observed data. For testing conditional independence between continuous variables, the Fisher's Z-test based on the partial correlation coefficient can be used. Specifically, given a set of conditional variables, the partial correlation coefficient between the two variables to be tested is calculated. Based on this partial correlation coefficient, the Fisher transform value is calculated, and the Fisher transform value is scaled according to the sample size used for this test and the number of variables in the set of conditional variables to obtain the corresponding test statistic. Based on the test statistic and a preset significance level, it is determined whether the two variables to be tested satisfy a conditional independence relationship given the set of conditional variables.

[0061] When the PC algorithm determines candidate adjacency relationships between variables, domain knowledge is introduced to restrict connections that do not conform to production timelines and causal directions. For example, quality defect indicators are not used as parent nodes of early-stage environmental variables; later-stage variables are not used as causes of earlier-stage variables; for variables with clear physical coupling relationships that need to be preserved, the adjacency relationships can be maintained during the skeleton learning stage, and the allowed directions can be determined by combining the results of conditional independence and time order constraints during the direction determination stage.

[0062] For each non-root node in the causal graph structure, the set of all parent nodes of that non-root node is used as input, and Gaussian process regression is employed to learn the nonlinear mapping relationship between the set of parent nodes and the non-root node. The value of the non-root node is determined by the nonlinear prediction result obtained from the set of parent nodes and the corresponding random noise.

[0063] Gaussian process regression can use a radial basis function kernel to measure the similarity between input samples from different parent nodes. Specifically, the squared Euclidean distance between two parent node input samples is calculated, and the squared distance is normalized according to a preset length scale parameter so that the similarity between samples decreases exponentially as the input distance increases. At the same time, the overall amplitude of the kernel function is adjusted using a signal variance parameter.

[0064] The signal variance, length scale, and noise variance parameters of Gaussian process regression can be determined by maximizing the log-marginal likelihood corresponding to the training data. During optimization, a kernel matrix is ​​formed based on the input samples, and the noise variance is incorporated into this kernel matrix. Considering both the degree of matching between the observed values ​​and the model's predicted distribution, as well as the model's complexity, the log-marginal likelihood under the current parameter conditions is calculated. Then, the L-BFGS numerical optimization algorithm is used to iteratively adjust each hyperparameter until the log-marginal likelihood meets the preset convergence condition.

[0065] Furthermore, for the modifiable variable X and the target quality defect index Y in the causal graph structure, an intervention value x is set. Graph operations are performed on the causal graph to delete all directed edges pointing to X, and X is fixed as the intervention value x. Monte Carlo sampling is then performed on the causal graph after intervention by calling the corresponding nonlinear structural equations in topological order to obtain the predicted distribution of the target quality defect index Y under the intervention state.

[0066] Let the reference value corresponding to the interventionable variable X under the baseline state be the baseline intervention value. Calculate the expected value of the target quality defect indicator Y when the interventionable variable X is fixed to the target intervention value, and the expected value of the target quality defect indicator Y when the interventionable variable X is fixed to the baseline intervention value. The difference between the former expected value and the latter baseline expected value is taken as the average causal effect of the interventionable variable X on the target quality defect indicator Y. Performing the above process for multiple interventionable variables separately yields the average causal effect of different interventionable variables on the same target quality defect indicator. The intervention inference process for the above interventionable variables is as follows: Figure 2 As shown.

[0067] In one optional specific example, taking the impact of the average monthly humidity (HQ2) in the second quarter on the dark spot area ratio (Ddark) as an example, HQ2 was intervened from 78% to 62%. A Gaussian process regression structural equation was used to perform 10,000 Monte Carlo samplings on the intervention plot, resulting in an ACE of -0.102. A 95% confidence interval [-0.128, -0.076] was obtained through 1,000 Bootstrap resampling operations. These values ​​are used to illustrate the implementation process of the intervention inference.

[0068] Furthermore, for the target production batch that has already occurred, let the set of endogenous variables in the causal model be V and the set of exogenous variables be U. Extract the actual observed variable values ​​(vobs) of the production batch, and estimate the posterior distribution of the exogenous variables under the actual observation conditions based on the actual observed variable values ​​of the target production batch.

[0069] In one specific implementation, the Markov chain Monte Carlo method is used to estimate the posterior distribution of exogenous variables, for example, using Metropolis-Hastings sampling. Then, an interventionist variable X is selected, and a counterfactual intervention value xcf is set. The interventionist variable X is fixed to the set counterfactual intervention value, forming the corresponding counterfactual intervention condition. An exogenous variable sample is drawn from the posterior distribution of the exogenous variable, and this sample is input into the post-intervention structural equation and forward propagated according to the causal graph to obtain the predicted distribution of the target quality defect index under the counterfactual condition. The counterfactual inference process for the above target production batch is as follows: Figure 3 As shown.

[0070] In one optional example, batch B3 from the Longxi base is used as the subject. The sample size of this batch is 1200 roots, the actual dark spot area ratio (Ddark) is 0.15, and the average monthly humidity (HQ2) in the second quarter is 78%; the counterfactual assumption is to set HQ2 to 62%. The Metropolis-Hastings algorithm is used to perform 10,000 samplings, with the first 2000 samplings used as burn-in to estimate the posterior distribution of exogenous variables under the condition that the actual observed variable values ​​of the target production batch are known. Then, 1000 sets of exogenous variables are extracted from the posterior distribution, and forward calculation is performed under the counterfactual intervention condition that the average monthly humidity in the second quarter is fixed at 62%, resulting in the counterfactual Ddark median of 0.048 and the 95% confidence interval [0.031, 0.072].

[0071] In a corresponding field validation: during the next planting cycle, the average monthly humidity in the second quarter for the 1200 samples in the experimental group was controlled at 62%±3%, while the average monthly humidity in the second quarter for the 1200 samples in the control group was 76%±5%. After harvest, the mean black spot area ratio in the experimental group was 0.052±0.008, and in the control group it was 0.143±0.015. This validation data is used to illustrate that the counterfactual prediction results can be compared with subsequent actual production results.

[0072] In further validation at three planting sites, the sample size at each site was no less than 800 and the planting cycle was repeated for two consecutive years. After adopting the humidity control recommendations, the black spot rate decreased by an average of 52% in the following season, with a 95% confidence interval of 38% to 66%. The above validation data are specific implementation results and should not be used as numerical limits for average causal effects or counterfactual prediction results.

[0073] This embodiment also provides an Astragalus membranaceus state recognition system for implementing the above method, including an image acquisition module, a distribution model construction module, a feature encoding module, and a state determination module. The image acquisition module acquires surface images of healthy Astragalus membranaceus from the current production batch and surface images of Astragalus membranaceus to be identified; the distribution model construction module extracts a set of local chromaticity contrast vectors from the healthy region and uses the expectation-maximization algorithm to fit a bivariate Gaussian mixture model of the current production batch.

[0074] It should be noted that after switching production batches and refitting the bivariate Gaussian mixture model, a unified component arrangement order is determined based on the two-dimensional mean vector of each Gaussian component; or a corresponding match is made based on the mean vector and covariance characteristics of the Gaussian components before and after the switch, so that Gaussian components with corresponding chromaticity distribution characteristics in different production batches use the same chromaticity contrast mode number.

[0075] The feature encoding module performs Gaussian pyramid decomposition, local chroma contrast encoding, local directional gradient encoding, chroma-direction joint statistics, and multi-scale joint feature concatenation; the state determination module inputs the joint feature vector into the state determination model and outputs state determination results such as good, defective, or further black spots, cracks, etc.

[0076] The system can also be configured with a grading module, a traceability data association module, and a causal analysis module. The grading module can form grades based on status, defect severity, size, and fullness in the manner described above; the traceability data association module associates visual fingerprints with production batch data; and the causal analysis module performs the aforementioned causal graph learning, intervention inference, and counterfactual prediction based on cross-batch data.

[0077] The modules described above can be implemented by software functional units in industrial computers, servers, or other data processing equipment, or by a processor executing computer programs stored in memory. Images, bivariate Gaussian mixture model parameters, joint feature vectors, state determination results, visual fingerprints, and production batch data can be stored in a local database or a server database.

[0078] In one optional example, partial results from six months of continuous operation of the system under trial production conditions are presented: approximately 3.6 million Astragalus roots were processed cumulatively, with a monthly overall F1-score of 97.6% ± 0.3%; during production, four online BGMM updates were triggered due to seasonal or geographical location changes, with the F1-score remaining within the range of 97.2% to 97.8% after each update; the end-to-end latency from image acquisition to traceability label printing was 220ms ± 15ms, with 10,000 measurements performed. This data is used to illustrate the feasibility of the system's production operation.

[0079] The above embodiments are used to illustrate the technical solution of the present invention. For those skilled in the art, without changing the main technical line of modeling the healthy chromaticity distribution of the current production batch, constructing multi-scale chromaticity-direction joint features, and determining the state based on the joint features, corresponding adjustments can be made according to the detection equipment, image resolution, number of spatial scales, number of Gaussian components, number of gradient direction intervals, and parameters of the state determination model.

Claims

1. A method for identifying the state of Astragalus membranaceus, characterized in that, include: Acquire surface images of healthy Astragalus membranaceus from the current production batch and surface images of Astragalus membranaceus to be identified; A set of local chromaticity contrast vectors is extracted from the healthy region of the surface image of the healthy Astragalus membranaceus. A bivariate Gaussian mixture model is fitted with the set of local chromaticity contrast vectors to obtain a chromaticity contrast distribution model that characterizes the chromaticity contrast distribution of the healthy region in the current production batch. Local chroma contrast encoding and local directional gradient encoding are performed on the surface image of Astragalus membranaceus to be identified at multiple spatial scales. The chroma contrast pattern is obtained according to the chroma contrast distribution model, and the directional gradient pattern is obtained according to the local directional gradient encoding. The joint occurrence of the chromaticity contrast mode and the directional gradient mode at each spatial scale is statistically analyzed to obtain the joint features corresponding to each spatial scale. The joint features of multiple spatial scales are then fused to obtain a joint feature vector. The state determination result of the Astragalus membranaceus to be identified is determined based on the joint feature vector.

2. The method for identifying the state of Astragalus membranaceus according to claim 1, characterized in that, Performing the local chromaticity contrast encoding and the local directional gradient encoding at multiple spatial scales includes: performing Gaussian pyramid decomposition on the surface image of the Astragalus membranaceus to be identified to obtain images at multiple scale levels; and performing the local chromaticity contrast encoding and the local directional gradient encoding in the pixel neighborhood of each scale level image.

3. The method for identifying the state of Astragalus membranaceus according to claim 2, characterized in that, The local chromaticity contrast encoding includes: converting the scale-level image to the HSV color space; determining a local chromaticity contrast vector for a pixel neighborhood in the scale-level image based on the hue and saturation differences between the neighboring pixels and the center pixel; and determining the chromaticity contrast mode corresponding to the local chromaticity contrast vector based on the chromaticity contrast distribution model.

4. The method for identifying the state of Astragalus membranaceus according to claim 3, characterized in that, Using the set of local chromaticity contrast vectors corresponding to the healthy region of the current production batch as training data, the mixing coefficients, mean vector, and covariance matrix of the bivariate Gaussian mixture model are fitted online using the expectation-maximization algorithm. After the production batch is switched, the bivariate Gaussian mixture model is refitted based on the set of local chromaticity contrast vectors corresponding to the healthy region of the switched production batch.

5. The method for identifying the state of Astragalus membranaceus according to any one of claims 2 to 4, characterized in that, The local directional gradient encoding includes: calculating the gradient magnitude and gradient direction of the scale-level image, performing directional interval quantization on the gradient direction; constructing a gradient direction co-occurrence matrix based on the adjacent occurrence relationship of different gradient directions in the pixel neighborhood, and determining the corresponding directional gradient mode based on the gradient direction co-occurrence matrix.

6. The method for identifying the state of Astragalus membranaceus according to claim 5, characterized in that, The joint occurrence of the chroma contrast mode and the directional gradient mode at each spatial scale is statistically analyzed, including: for each spatial scale, constructing a two-dimensional joint histogram based on the combination of different chroma contrast modes and different directional gradient modes, and counting the number of occurrences of each combination in the corresponding scale level image; normalizing the two-dimensional joint histograms corresponding to each spatial scale, and splicing them according to the spatial scale order to obtain the joint feature vector.

7. The method for identifying the state of Astragalus membranaceus according to any one of claims 1 to 6, characterized in that, After obtaining the state determination result, the method further includes: performing a hash operation on the joint feature vector, using the resulting hash value as the visual fingerprint of the corresponding detection record; and associating and storing the visual fingerprint with the state determination result and the corresponding production batch's place of origin environmental factor data, cultivation management data, and quality defect data to form traceability data.

8. The method for identifying the state of Astragalus membranaceus according to claim 7, characterized in that, Also includes: Environmental variables, cultivation management variables, and quality defect indicators are determined based on the traceability data corresponding to multiple production batches. The PC algorithm is used to determine the adjacency relationships among the environmental variables, cultivation management variables, and quality defect indicators based on conditional independence tests, and the directed relationships are determined by combining domain knowledge priors to obtain a causal graph structure. The domain knowledge priors include at least the time order constraints of the variables and the directional constraints that the quality defect indicators are not used as parent nodes of the prior environmental variables. For the non-root nodes in the causal graph structure, the corresponding nonlinear structural equations are estimated using the set of parent nodes of the corresponding non-root nodes as input and Gaussian process regression.

9. The method for identifying the state of Astragalus membranaceus according to claim 8, characterized in that, Also includes: For the modifiable variable and the target quality defect index in the causal graph structure, the directed edges pointing to the modifiable variable are deleted and an intervention value is set. The nonlinear structural equation is used for sampling to obtain the expected value of the target quality defect index under the intervention state. The average causal effect of the modifiable variable on the target quality defect index is determined based on the difference between the expected value and the expected value of the target quality defect index under the baseline state. For a target production batch, based on the observed data of the target production batch, the posterior distribution of the exogenous variables is estimated using the Markov chain Monte Carlo method, a counterfactual intervention value is set, and the counterfactual prediction result of the target quality defect index is obtained using the posterior distribution of the exogenous variables and the nonlinear structural equation after intervention.

10. A system for recognizing the state of Astragalus membranaceus, characterized in that, It includes an image acquisition module, a distribution model construction module, a feature encoding module, and a state determination module; The image acquisition module is used to acquire surface images of healthy Astragalus membranaceus from the current production batch and surface images of Astragalus membranaceus to be identified. The distribution model construction module is used to extract a set of local chromaticity contrast vectors from the healthy region of the surface image of the healthy Astragalus membranaceus, and fit a bivariate Gaussian mixture model with the set of local chromaticity contrast vectors to obtain the chromaticity contrast distribution model of the current production batch. The feature encoding module is used to perform local chroma contrast encoding and local directional gradient encoding on the surface image of Astragalus membranaceus to be identified at multiple spatial scales. It obtains the chroma contrast pattern according to the chroma contrast distribution model, obtains the directional gradient pattern according to the local directional gradient encoding, and statistically analyzes the joint occurrence of the chroma contrast pattern and the directional gradient pattern and performs multi-scale fusion to obtain a joint feature vector. The state determination module is used to determine the state determination result of the Astragalus membranaceus to be identified based on the joint feature vector.

Images