A method for breeding and screening of spirulina algae with high antioxidant activity

Abstract

This invention discloses a method for breeding and screening Spirulina strains with high antioxidant activity, relating to the field of microalgae biotechnology. The method involves culturing different candidate Spirulina strains in multiple culture units, installing multispectral sensor arrays in each culture unit to collect real-time spectral data and environmental parameter data during the Spirulina growth process, and obtaining spectral feature vectors related to antioxidant activity. This method enables independent identification and comparative analysis of the optimal activity state of each strain, effectively avoiding missing the activity peak window, significantly improving the screening efficiency and stability of superior high antioxidant activity strains, and providing reliable technical support for the targeted breeding and large-scale application of Spirulina strains.

Description

Technical Field

[0001] This invention relates to the field of microalgae biotechnology, specifically to a method for breeding and screening Spirulina strains with high antioxidant activity. Background Technology

[0002] In the field of microalgae resource development and functional foods, spirulina has attracted widespread attention due to its rich content of various antioxidant active ingredients. Its antioxidant capacity is often used as an important indicator for evaluating the quality and application value of algal strains. In existing technologies, the selection and screening of spirulina strains with high antioxidant activity typically involves batch culture of different algal strains followed by sampling and testing. The antioxidant activity of each algal strain at specific time points is determined through chemical analysis or in vitro antioxidant experiments, and the test results are used as the basis for screening.

[0003] However, the growth process of Spirulina exhibits distinct stages and dynamic characteristics, with the synthesis and accumulation of its antioxidants influenced by multiple factors, including growth stage, nutritional status, and environmental conditions. The trends in antioxidant content within the algae differ at different growth stages, and the test results for the same algal strain at different time points within the culture cycle may vary significantly. Current technologies largely rely on manual experience or fixed culture cycles for sampling and testing, making it difficult to accurately reflect the dynamic changes in antioxidant activity of each algal strain throughout the complete culture cycle. Furthermore, under parallel culture conditions of multiple algal strains, the growth rate, metabolic intensity, and the timing of peak antioxidant accumulation often differ among the strains. If uniform testing is conducted only at predetermined time points, some algal strains may be evaluated before reaching their peak antioxidant activity or have missed their optimal activity stage, thus affecting the objectivity and accuracy of the screening results. This mismatch in testing time points introduces significant uncertainty into the screening results for algal strains with high antioxidant activity. Meanwhile, existing screening methods mostly rely on a single indicator or a small amount of detection data as the evaluation basis, lacking a comprehensive analysis of the peak level of antioxidant activity, the time of peak occurrence, and the trend of changes throughout the entire cycle. This makes it difficult to establish a systematic and quantitative evaluation standard, resulting in a highly subjective screening process for superior algal species with insufficient repeatability and stability. Summary of the Invention

[0004] The purpose of this invention is to provide a method for breeding and screening Spirulina strains with high antioxidant activity, thereby solving the problems existing in the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for breeding and screening Spirulina strains with high antioxidant activity, comprising: S1, culturing different candidate Spirulina strains in multiple culture units, installing a multispectral sensor array in each culture unit to collect real-time spectral data and environmental parameter data during the growth process of Spirulina, and obtaining spectral feature vectors related to antioxidant activity; S2, inputting the obtained spectral feature vectors related to antioxidant activity into a pre-trained support vector machine model for classification, determining the current growth stage of the strain, and calculating the dynamic change trend of the growth stage by fusing environmental parameter data after obtaining the classification results; S3, if the determined growth stage dynamic change... If the trend shows that biomass growth is slowing down, then the random forest algorithm is used to perform regression prediction on the spectral feature vector to obtain the predicted antioxidant concentration values, and obtain continuous time series data of the concentration values; S4, the sliding window method is used to analyze the obtained continuous time series data of the concentration values, calculate the concentration gradient change within each window, determine whether the gradient change exceeds the preset threshold, and obtain a set of candidate points for potential peak moments; S5, the obtained set of candidate points for potential peak moments is clustered, and the K-means algorithm is used to group similar points to obtain clustered peak moment groups. The point with the highest concentration value is selected from the group to determine the peak antioxidant content of the corresponding algal strain during its culture cycle.

[0006] Preferably, step S1 includes acquiring real-time spectral data and environmental data during the growth process of different candidate Spirulina strains using a multispectral sensor array; processing the real-time spectral data and environmental data using a wavelet denoising algorithm to obtain a smoothed spectral sequence and a standard environmental sequence; processing the smoothed spectral sequence using Fourier transform infrared spectroscopy analysis to obtain frequency domain spectral features; constructing a multidimensional feature tensor based on the frequency domain spectral features and the standard environmental sequence; and extracting spectral feature vectors related to antioxidant activity from the multidimensional feature tensor if the multidimensional feature tensor meets a preset activity threshold condition.

[0007] Preferably, step S2 includes normalizing the spectral feature vectors related to antioxidant activity to obtain a standard feature matrix; inputting the standard feature matrix into a pre-trained support vector machine model for classification and recognition to obtain the current algal growth stage; retrieving historical environmental parameter data within the corresponding time window based on the current algal growth stage to obtain an environmental time series set; performing multi-dimensional spatial mapping between the environmental time series set and the current algal growth stage to obtain stage environmental correlation coordinates; if the stage environmental correlation coordinates deviate from the preset steady-state center point, calculating the projection displacement of the stage environmental correlation coordinates on the time axis to obtain the dynamic change trend of the growth stage.

[0008] Preferably, step S3 includes acquiring the algal density collected by the sensor and calculating the biomass growth rate; if the biomass growth rate is lower than a preset dynamic threshold range, retrieving the original spectral data at the corresponding time point; performing principal component analysis on the original spectral data to extract dimensionality-reduced feature vectors; inputting the dimensionality-reduced feature vectors into a pre-trained random forest regression model to obtain predicted concentration values; constructing a concentration time series based on the distribution pattern of the predicted concentration values ​​on the time axis; and using the concentration time series to fit the material accumulation trajectory and determine the concentration values ​​of antioxidant substances.

[0009] Preferably, step S4 includes acquiring continuous time series data of antioxidant concentration values; performing sliding segmentation on the continuous time series data according to a preset window step size to obtain multiple overlapping local concentration intervals; performing linear fitting on sampling points within the local concentration intervals to calculate concentration gradient changes; if the positive or negative sign of the concentration gradient change is reversed and the absolute value exceeds a preset threshold, then extracting the slope extreme points within the local concentration intervals as candidate points for potential peak moments; performing clustering processing based on the distribution density of candidate points on the time axis, removing outliers caused by time series fluctuations to determine the final candidate point set; and using the sampling frequency corresponding to the candidate point set to perform local data smoothing to obtain a candidate point set that can reflect the peak of antioxidant accumulation.

[0010] Preferably, step S5 includes obtaining the time coordinates in the candidate point set and constructing a one-dimensional feature vector space, calculating the Euclidean distance between different time coordinates in the one-dimensional feature vector space; iteratively dividing the candidate point set according to the Euclidean distance using the K-means algorithm to obtain peak time groups where the cluster center no longer shifts and has cluster labels; retrieving the concentration values ​​within the peak time groups and extracting the sampling time with the largest concentration value as a local extremum; mapping the local extremum to the time axis of the algal cycle to identify the growth stage, and determining the peak value of antioxidant content of the corresponding algal strain in its culture cycle by comparing the local extremums of different growth stages.

[0011] Preferably, the method further includes S6: constructing an antioxidant activity evaluation index system based on the peak value of antioxidant substances, the time of peak occurrence, and the dynamic change trajectory of each candidate Spirulina strain throughout its entire culture cycle. Specifically, this includes obtaining the time series sequence of antioxidant substance content for each candidate Spirulina strain throughout its entire culture cycle; obtaining a set of content change rates by performing a first-order difference operation on the time series sequence of antioxidant substance content; determining the peak value of antioxidant substance content and its corresponding peak occurrence time using the set of content change rates; performing a normalization mapping based on the peak value of antioxidant substance content and the peak occurrence time; importing the mapped values ​​into a preset growth trajectory fitting model; and obtaining a morphological feature vector reflecting the dynamic change trajectory throughout the entire culture cycle by extracting features from the growth trajectory fitting model.

[0012] Preferably, step S6 further includes calculating the trajectory deviation of each candidate Spirulina strain relative to the benchmark growth model using morphological feature vectors. If the trajectory deviation exceeds a preset stability threshold, the instantaneous slope of the candidate Spirulina strain in the deviation interval is extracted as an activity fluctuation variable. By performing a weighted matrix operation on the activity fluctuation variable, the peak value of antioxidant content, and the peak occurrence time, a comprehensive activity score for each candidate Spirulina strain is obtained, and an antioxidant activity evaluation index system is constructed using the comprehensive activity score.

[0013] Preferably, the method further includes S7: ranking and screening multiple candidate Spirulina strains according to the antioxidant activity evaluation index system, and determining the superior algae strain with the highest antioxidant activity as the target algae strain for breeding. Specifically, this includes arranging multiple candidate Spirulina strains in descending order according to the comprehensive activity score to obtain an initial sequence set; constructing a fluctuation discrete matrix using the activity fluctuation variables of each candidate Spirulina strain in the initial sequence set; and extracting the eigenvalue vector reflecting growth stability by performing singular value decomposition on the fluctuation discrete matrix.

[0014] Preferably, step S7 further includes weighting and correcting the initial sequence set according to the feature root vector to obtain a corrected activity ranking sequence; if the trajectory deviation of the first candidate Spirulina strain in the activity ranking sequence is lower than a preset stability threshold, the candidate Spirulina strain is determined as an excellent algae species; by verifying the extreme value of the peak value of antioxidant content in the dynamic change trajectory of the excellent algae species throughout the entire cycle, the excellent algae species with the highest antioxidant activity is determined as the target algae species for breeding.

[0015] As can be seen from the above technical solution, the present invention has the following beneficial effects: This method for breeding and screening Spirulina strains with high antioxidant activity achieves continuous tracking and precise peak identification of the accumulation process of antioxidant substances in each strain throughout the complete culture cycle by simultaneously culturing different candidate Spirulina strains in multiple culture units and combining multispectral real-time monitoring with dynamic modeling and analysis of antioxidant activity. This avoids the evaluation bias problems caused by traditional fixed-time-point sampling and testing. By comprehensively analyzing the peak content of antioxidant substances, the time of peak occurrence, and the trajectory of changes throughout the entire cycle, a systematic antioxidant activity evaluation index system is constructed, transforming the screening process from single-index judgment to multi-dimensional quantitative evaluation, thus improving the objectivity and accuracy of the screening results. At the same time, considering the differences in growth stages and the asynchronous accumulation of antioxidant substances among different strains, the method enables independent identification and comparative analysis of the optimal activity state of each strain, thereby effectively avoiding missing the activity peak window and significantly improving the screening efficiency and stability of superior strains with high antioxidant activity. This provides reliable technical support for the targeted breeding and large-scale application of Spirulina strains. Attached Figure Description

[0016] Figure 1This is a flowchart of the method of the present invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] like Figure 1 As shown, the present invention provides a technical solution: a method for breeding and screening Spirulina strains with high antioxidant activity, comprising: S1. By culturing different candidate Spirulina strains in multiple culture units, and installing a multispectral sensor array in each culture unit to collect real-time spectral data and environmental parameter data during the growth process of Spirulina, the spectral feature vectors related to antioxidant activity are obtained. S2. Based on the obtained spectral feature vectors related to antioxidant activity, input them into a pre-trained support vector machine model for classification, determine the current growth stage of the algal strain, and then integrate environmental parameter data to calculate the dynamic change trend of the growth stage after obtaining the classification results. S3. If the dynamic change trend of the determined growth stage shows that the biomass growth is slowing down, then the random forest algorithm is used to perform regression prediction on the spectral feature vector to obtain the predicted antioxidant concentration value and obtain continuous time series data of the concentration value. S4. The continuous time series data of the obtained concentration values ​​are analyzed using the sliding window method. The concentration gradient change within each window is calculated, and it is determined whether the gradient change exceeds the preset threshold to obtain a set of candidate points for potential peak times. S5. Cluster the obtained potential peak moment candidate point set, use the K-means algorithm to group similar points to obtain the clustered peak moment group, select the point with the highest concentration value from the group, and determine the peak content of antioxidant substances of the corresponding algal strain in its culture cycle. S6. Based on the peak content of antioxidant substances, the time of peak occurrence, and the dynamic change trajectory of the entire cycle of each candidate Spirulina strain, an antioxidant activity evaluation index system was constructed. S7. Based on the antioxidant activity evaluation index system, multiple candidate Spirulina strains were sorted and screened to determine the superior algae strain with the highest antioxidant activity as the target algae strain for breeding.

[0019] In the above scheme, each culture unit is an independent and controllable photobioreactor, equipped with a light source control module, a temperature control module, a pH adjustment module, and a stirring and aeration module to ensure that different candidate algal strains grow under relatively consistent or adjustable gradient culture conditions. A multispectral sensor array, including visible light and near-infrared sensors, is installed on the sidewall or top of the culture unit to collect reflectance or transmission spectral signals in the 400nm–900nm range. Environmental parameters include temperature, pH, dissolved oxygen concentration, light intensity, and carbon dioxide concentration.

[0020] After denoising, standardization, and principal component analysis (PCA) dimensionality reduction, the collected raw spectral data were used to extract feature bands related to antioxidants such as phycocyanin, chlorophyll a, and carotenoids, constructing spectral feature vectors associated with antioxidant activity. A pre-trained support vector machine (SVM) model was used to establish a classification hyperplane based on labeled growth stage samples to identify the lag phase, exponential growth phase, stationary phase, and decline phase. The classification results were then fused with environmental parameter data, and a weighted time series model or a Kalman filter model was used to calculate the dynamic trends of growth stages to determine changes in biomass growth rate.

[0021] When the biomass growth rate shows a significant downward trend, it indicates that the algal strain may be entering a metabolic redistribution phase, which is often accompanied by the accumulation of antioxidants. At this time, a random forest regression model is used to perform nonlinear mapping prediction on the spectral feature vector, outputting the concentration values ​​of antioxidants. The time sampling interval can be set to 10 min to 60 min to form continuous time series data.

[0022] The sliding window method uses a fixed time length (e.g., 6 hours or 12 hours) as the window width and progressively slides to calculate the rate of change of concentration gradient within the window. When the rate of change of gradient exceeds a preset threshold, that time point is marked as a potential peak candidate point. Subsequently, the K-means algorithm is used to perform cluster analysis on the candidate points. The number of clusters is adaptively determined based on the silhouette coefficient, grouping candidate points that are close in time and have similar concentration change trends into the same group. The time point with the largest predicted concentration value in each group is selected as the peak time of antioxidant content.

[0023] Based on peak concentration, peak occurrence time, and the dynamic trajectory of concentration changes throughout the entire cycle, a multidimensional antioxidant activity evaluation index system was constructed, including peak intensity, peak lead time, stability, and cumulative antioxidant yield. Multiple candidate algal strains were comprehensively ranked using a weighted scoring model or analytic hierarchy process (AHP) to screen out superior algal strains with the highest and most stable antioxidant activity.

[0024] S1 includes acquiring real-time spectral data and environmental data during the growth process of different candidate Spirulina strains through a multispectral sensor array; processing the real-time spectral data and environmental data using a wavelet denoising algorithm to obtain a smoothed spectral sequence and a standard environmental sequence; processing the smoothed spectral sequence using Fourier transform infrared spectroscopy analysis to obtain frequency domain spectral features; constructing a multidimensional feature tensor based on the frequency domain spectral features and the standard environmental sequence; and extracting spectral feature vectors related to antioxidant activity from the multidimensional feature tensor if the multidimensional feature tensor meets the preset activity threshold condition.

[0025] In this embodiment, before the culture is started, multiple culture units are assigned to different candidate Spirulina strains. The initial volume of the algal solution in each culture unit is set to 3 liters, and the initial inoculation density is set to 0.2 g per liter. The initial volume of 3 liters is determined by matching the sensor's field of view coverage area with the optical path length. It is required that the sensor's observation area falls completely within the effective liquid surface range and the optical path fluctuation does not exceed 1 mm. The initial inoculation density of 0.2 g per liter is determined by converting the dry weight measurement results of the pre-cultured algal solution. During the conversion, three weighing results are taken for the same batch of pre-cultured algal solution, and the average value is taken as the density conversion benchmark to avoid the error of a single weighing being amplified to the online model input.

[0026] A multispectral sensor array is arranged at the side wall light-transmitting window and top observation port of each culture unit. The sensor coverage wavelength range is set to 400 nm to 900 nm. The wavelength range of 400 nm to 900 nm is determined by the coverage requirements of the main pigment absorption and scattering bands of Spirulina, requiring the inclusion of the chlorophyll absorption peak neighborhood, the phycocyanin characteristic band neighborhood, and the near-infrared scattering sensitive region, so that the spectral characteristics are simultaneously sensitive to changes in biomass and antioxidant-related components. The sensor sampling period is set to 10 minutes. The 10-minute sampling period is determined by balancing the rate of spectral change during growth with the computational load of data processing. First, during a continuous 24-hour pre-run phase, the mean and variance of the spectral differences between adjacent sampling points are statistically analyzed. If increasing the sampling period causes the mean difference to change by more than 0.2 times the original mean, the time resolution is considered insufficient, and the 10-minute setting is maintained. If decreasing the sampling period causes the data volume per unit time to cause the single-machine processing delay to exceed the sampling period, the computing power is considered insufficient, and the 10-minute setting is maintained. Five environmental data points—temperature, pH, dissolved oxygen, light intensity, and carbon dioxide concentration—are collected synchronously. The sampling timestamps are consistent with the spectral sampling timestamps and are aligned using the same clock source to avoid false correlations introduced by timing misalignments.

[0027] Real-time spectral data and environmental data are processed by a denoising and standardization module. Wavelet denoising is applied to the spectral sequences, with a decomposition layer of 5. The number of decomposition layers (5) is determined by the noise band separation effect. First, spectra collected from the same culture unit in an algae-free culture medium are used as noise samples. The energy proportion of high-frequency components is calculated at different decomposition layer numbers. The decomposition layer number is gradually increased from 3. When adding one more layer results in a decrease in the high-frequency energy proportion of less than 0.05 times and the change in the width of the main peak of the reconstructed spectrum exceeds 0.02 times the original width, over-decomposition is considered to have occurred, and the previous layer number is taken as the decomposition layer number, resulting in 5 layers. A symmetric orthogonal wavelet basis is selected, based on the highest correlation between the reconstructed spectrum and the original smooth trend, and the lowest residual energy of high-frequency noise. Specifically, decomposition and reconstruction are performed on three candidate wavelet bases, the correlation between the reconstructed spectrum and the original sequence is calculated, and the energy of the high-frequency residual after reconstruction is calculated. The wavelet base with the highest correlation and the lowest residual energy is selected as the one used. The denoising threshold adopts a soft thresholding method, and the threshold value is determined by noise sample estimation. The threshold determination process involves extracting the absolute value sequence of the highest frequency component, taking the median as the noise scale indicator, and then selecting an amplification coefficient related to the number of sampling points. The amplification coefficient ranges from 0.8 to 1.6, traversing it with a step size of 0.1, and performing denoising reconstruction one by one. The signal-to-noise ratio improvement after denoising is compared with the peak distortion. The peak distortion is measured by the change in the full width at half maximum (FWHM) of the main peak and the shift in the peak position. The coefficient corresponding to the largest signal-to-noise ratio improvement and the smallest peak distortion is selected to obtain the final threshold. Environmental data denoising adopts the same wavelet denoising process, and then a standard environmental sequence is formed. The standardization baseline period is set to 6 hours. The 6-hour baseline period is determined by the natural fluctuation period of environmental parameters. First, the autocorrelation of the environmental data in the pre-running stage is calculated, and the time interval corresponding to the first significant peak of the autocorrelation is taken as the main fluctuation period. If the main fluctuation period falls between 5 hours and 7 hours, 6 hours is taken. The standardization process involves calculating the mean and standard deviation of each environmental parameter within each 6-hour sliding segment. Then, the mean is subtracted from each sampled value within that segment, and the result is divided by the standard deviation to obtain a dimensionless standard environmental sequence, thereby eliminating the weight bias caused by different dimensions.

[0028] The smoothed spectral sequence is then fed into the frequency domain feature extraction module, where it is mapped to the frequency domain using Fourier transform infrared spectroscopy. To ensure consistent sampling intervals during the transform, the wavelength axis sampling points are resampled with a step size of 2 nanometers. This 2-nanometer step size is determined by the sensor channel bandwidth and feature resolution requirements, ensuring it is no greater than 0.5 times the half-bandwidth of the narrowest channel to avoid spectral aliasing caused by inter-channel interpolation. Linear interpolation is used for the resampling method, with the interpolation interval strictly limited to the original wavelength coverage to avoid extrapolation. A window function is then applied to the resampled spectrum at each time point to suppress boundary effects. The Hanning window is chosen based on its low sidelobe leakage and moderate main lobe width, which can reduce spectral pseudo-peaks while preserving peak position information. The transform length is set to 256 points. This length is determined by taking the nearest power of 2 from the number of wavelength points after resampling, with the rounding rule being that it should not be less than the actual number of points and the difference should be minimized, thereby reducing computational load and maintaining stable frequency resolution. After the transformation is completed, frequency domain spectral features are extracted. The features include the position of the main peak of the amplitude spectrum, the amplitude of the main peak, the half-peak width of the main peak, the proportion of low-frequency energy, the proportion of mid-frequency energy, the proportion of high-frequency energy, and the spectral area of ​​the specified frequency band. The boundaries of each frequency band are determined through historical calibration samples. The calibration process involves calculating the spectral energy distribution of the samples that have completed offline antioxidant substance determination, finding the three energy accumulation regions with the highest correlation with the measurement results, and defining them as the boundaries of the low-frequency, mid-frequency, and high-frequency bands to maximize the correlation.

[0029] The frequency domain spectral features were aligned with the standard environmental sequence on the time axis, with an alignment tolerance of 30 seconds. This 30-second tolerance was determined statistically by analyzing the maximum delay between sensor acquisition and data transmission. The statistical method involved recording the delay from each acquisition trigger to data entry during 12 consecutive hours of monitoring, taking the maximum delay, and adding a 10-second redundancy to obtain the 30-second tolerance. After alignment, a multidimensional feature tensor was constructed, containing four dimensions: algal strain dimension, time dimension, frequency domain feature dimension, and environmental dimension. The length of the algal strain dimension was taken from the number of candidate algal strains, the length of the time dimension from the number of consecutive sampling points, the length of the frequency domain feature dimension from the number of the aforementioned frequency domain feature items, and the length of the environmental dimension from five environmental parameters. To ensure consistency in the magnitude of features among different algal strains, frequency domain features were normalized to the same scale before tensor construction. The normalization benchmark was set as the mean and standard deviation of the features of the same algal strain in the middle of the exponential growth phase. The middle of the exponential growth phase was determined based on growth stage label data or offline biomass curves. The middle interval was formed by extending 12 hours before and after the period with the highest biomass growth rate. The mean and standard deviation were calculated within this interval. Then, the full-cycle features were processed by subtracting the mean and dividing by the standard deviation.

[0030] After constructing the multidimensional feature tensor, an activity threshold condition is determined. The activity threshold condition is based on an activity score, which is obtained by weighting frequency domain features and environmental features. The weights are determined through fitting the calibration samples. The fitting process uses the offline measured concentration of antioxidants as the target value, takes the frequency domain features of the calibration samples and the standard environmental sequence as input, and uses the minimum error criterion to calculate the contribution of each feature to the target value. The greater the contribution, the higher the weight. Then, all weights are scaled until the sum of weights is 1, forming a stable summarization rule. The activity threshold is set at an activity score of 0.7. The activity threshold of 0.7 is determined through screening performance verification. First, the calibration samples are divided into training and validation segments according to time sequence. In the validation segment, the activity threshold is traversed from 0.5 to 0.9 with a step size of 0.02. The number of correctly detected high-activity samples and the number of falsely detected low-activity samples are counted respectively. The threshold corresponding to the largest number of correct detections and the smallest number of false detections is selected, resulting in 0.7. The determination process involves calculating the activity score for each algal strain at each sampling time. If the activity score is not lower than 0.7 for three consecutive sampling times, the activity threshold condition is considered to be met for that period. The setting of three consecutive sampling times is determined by the requirement of suppressing occasional noise. The duration of single-point anomalies is usually no more than 20 minutes. Taking three sampling points to cover 30 minutes can reduce the interference of single-point spikes.

[0031] After meeting the activity threshold condition, spectral feature vectors related to antioxidant activity are extracted from the multidimensional feature tensor. Feature vector extraction employs a two-stage screening process: Stage 1 involves correlation screening within the frequency domain feature dimension. Correlation calculations are based on calibrated samples, calculating the correlation strength between each frequency domain feature and the offline concentration, retaining features with a correlation strength not lower than 0.6. The correlation strength threshold of 0.6 is determined through erroneous deletion risk control. During the validation phase, the threshold is iterated from 0.4 to 0.8 with a step size of 0.05, comparing the prediction error with the number of features, selecting the threshold with the smallest prediction error and no more than 12 features, resulting in 0.6. Stage 2 involves redundancy removal. Redundancy determination is based on the pairwise correlation strength between features. If the correlation strength between two features is higher than 0.85, the feature with the higher correlation strength to the offline concentration is retained, and the other is removed. The redundancy threshold of 0.85 is determined through multicollinearity control. During the validation phase, the correlation matrix of the model input features is observed; when it exceeds 0.85, the model coefficient fluctuations significantly increase and the generalization error increases, setting 0.85 as the removal threshold. After two-stage screening, continuous time segments meeting the activity threshold are extracted along the time dimension. The frequency domain features retained within these segments are then concatenated with a standard environmental sequence from the same time point in a predetermined order to form the final spectral feature vector. The concatenation order is frequency domain first, followed by environment, based on the fact that subsequent classification models are more sensitive to frequency domain features, and arranging frequency domain features first helps maintain the consistency of the input structure. Before outputting the feature vector, an integrity check is performed. The check rule is that if any feature is missing for more than one sampling period, the data at that time point is considered invalid. Missing data is filled using linear interpolation of two adjacent valid sampling points. The interpolation trigger condition is that the missing length does not exceed two sampling periods. The upper limit of two sampling periods for missing length is determined by interpolation distortion control. If the missing length exceeds two sampling periods, the interpolation causes a significant trend shift, and the input for that time period is directly discarded.

[0032] S2 includes normalizing the spectral feature vectors related to antioxidant activity to obtain a standard feature matrix; inputting the standard feature matrix into a pre-trained support vector machine model for classification and recognition to obtain the current algal growth stage; retrieving historical environmental parameter data within the corresponding time window based on the current algal growth stage to obtain an environmental time series set; performing multi-dimensional spatial mapping between the environmental time series set and the current algal growth stage to obtain stage environmental correlation coordinates; if the stage environmental correlation coordinates deviate from the preset steady-state center point, calculating the projection displacement of the stage environmental correlation coordinates on the time axis to obtain the dynamic change trend of the growth stage.

[0033] In this embodiment, for the spectral feature vectors related to antioxidant activity, normalization processing is first performed to form a standard feature matrix. The normalization processing starts with feature dimension consistency: for each spectral feature vector formed by the same algal strain at the current time, the presence of sudden increases or decreases is checked dimension by dimension. The determination of sudden increases or decreases is completed using the stable fluctuation boundary of the training samples. The stable fluctuation boundary is obtained by summarizing all steady-state sample points on the same feature dimension during the training phase. First, the mean and standard deviation of all steady-state sample points in this dimension are calculated. Then, the upper boundary is set as the mean plus 3 times the standard deviation, and the lower boundary is set as the mean minus 3 times the standard deviation. The process of determining 3 is as follows: in the steady-state training samples, the proportion of points exceeding the boundary is counted. 1, 2, 3, and 4 are tested step by step. The value that makes the number of points exceeding the boundary not exceed 1 / 100 of the total number of points and has the smallest boundary width is selected to obtain 3. When the real-time feature value is higher than the upper boundary, the feature value is replaced with the upper boundary. When the real-time feature value is lower than the lower boundary, the feature value is replaced with the lower boundary, thereby eliminating the pull of single-point spikes on subsequent classification and discrimination values. After completing the boundary constraints, dimensional normalization is performed: For each feature dimension, the minimum and maximum values ​​of that dimension during the training phase are retrieved. The process for determining the minimum and maximum values ​​is to first process the training samples according to the above boundary constraints, and then take the minimum and maximum values ​​of each dimension in the processed training samples. The normalization calculation process is to first subtract the minimum value of that dimension from the current feature value to obtain the difference relative to the starting point, and then divide the difference by the span between the maximum and minimum values ​​of that dimension, so that the output of that dimension falls within the range of 0 to 1. 0 and 1 are determined by the normalized output range definition. The construction of the standard feature matrix is ​​completed by stacking in chronological order: Using the normalized feature vectors of the six most recent consecutive sampling times as input segments, the normalized vectors of each sampling time are arranged in chronological order as a row to obtain the standard feature matrix. The process for determining the six consecutive sampling times is to test the validation data at sampling times 3, 4, 5, 6, 7, and 8 respectively, compare the number of stage label jumps with the stage delay, and select the length with the fewest jumps and a stage delay of no more than one sampling period, resulting in 6.

[0034] The standard feature matrix is ​​fed into the support vector machine (SVM) classifier to complete the identification of growth stages. During the training phase, the SVM classifier has already learned the boundaries of the stage samples. The training samples consist of standard feature matrix fragments labeled with growth stages. Stage labeling is generated based on a joint rule derived from changes in biomass growth rate and changes in the key spectral bandwidth. The joint rule generation process is as follows: for each culture batch, the growth amount of adjacent sampling intervals of the biomass curve is calculated; then, intervals of continuous decrease and continuous increase are calculated for the growth amount sequence; these intervals are matched with the amplitude change direction of the main absorption band in the spectrum; when a match is found, the intervals are labeled as the same stage. The classifier parameters include penalty strength parameters and kernel scale parameters. The candidate sets for penalty strength parameters are set to 0.1, 0.3, 1, 3, and 10, and the candidate sets for kernel scale parameters are set to 0.1, 0.3, 1, 3, and 10. The process of determining the candidate sets is to expand outwards from 1 at intervals of approximately 3 times, so that the parameters cover three ranges: low penalty, moderate penalty, and high penalty. The parameter selection process is as follows: for each set of candidate parameters, after boundary learning is completed on the training data, the number of misclassifications and sensitivity to noise perturbation are calculated segment by segment on the validation data. The sensitivity calculation process is to superimpose the perturbation amplitude obtained by the sensor noise statistics on each dimension of the validation data. The perturbation amplitude is taken as the standard deviation of the same dimension of the feature in the blank culture medium collection sequence. After superimposing in both the forward and reverse directions, the data is reclassified. The number of times the label changes during the stage is counted. The parameter set with the fewest changes and the fewest misclassifications is selected as the final parameters. The classification and recognition process is as follows: feature vectors are extracted row by row from the standard feature matrix and fed into the classifier in the same dimension order of the training phase. The classifier outputs a phase discrimination result for each row and performs a consistency decision on the discrimination results of 6 consecutive rows. The consistency decision is calculated by counting the phase label that appears most frequently in the 6 rows of discrimination results and taking the phase label that appears most frequently as the current phase label. If the phase label that appears most frequently appears less than 4 times, the phase label of the previous time step is kept unchanged. The process of determining 4 is to count the phase false jump events on the validation data, test the thresholds 3, 4 and 5 step by step, and select the threshold with the fewest false jump events and a phase update delay of no more than 2 sampling periods to obtain 4.

[0035] After obtaining the current algal growth stage, the system retrieves historical environmental parameter data within the corresponding time window to form an environmental time series set. The time window length is set to 360 minutes. The process of determining 360 minutes involves identifying environmental mutation events in the training data and measuring the stage response lag: the identification of environmental mutation events is completed using the change rate threshold of environmental parameters. The change rate threshold is obtained by adding 3 times the standard deviation to the mean change rate of steady-state samples. The method for determining 3 is consistent with the aforementioned boundary determination. When the change rate of any environmental parameter exceeds the threshold, it is marked as the starting point of the mutation event. At the same time, the time when the stage label first changes is recorded. The time difference between the two is used as a lag sample. After summing all lag samples, the maximum lag is taken and a 60-minute safety margin is added to obtain 360 minutes. The process of determining the 60 minutes involves statistically analyzing the fluctuation range of the maximum value of the lag sample between adjacent batches and taking the upper limit of this fluctuation range as the margin. The process of constructing the environmental time series set involves looking back 360 minutes from the current moment, extracting five data points—temperature, pH, dissolved oxygen, light intensity, and carbon dioxide concentration—in sequence according to the sampling timestamp to form a multivariate time series, and then performing standardization processing on each environmental parameter. The standardization process involves calculating the training mean and training standard deviation for each environmental parameter in the training data. The real-time value is first subtracted from the training mean, and then the difference is divided by the training standard deviation to transform different dimensions into the same scale. The handling of missing environmental measurement points employs length-constrained interpolation: the upper limit of the missing measurement length is set to 2 sampling periods. The process of determining 2 involves manually creating missing measurements of different lengths in the training data and interpolating them. Then, the interpolated environmental sequence is input into the stage environment mapping module, and the mapping coordinate deviation is statistically analyzed. The minimum missing measurement length whose deviation exceeds the steady-state distance threshold is taken as the upper limit, resulting in 2. When the missing measurement length does not exceed 2 sampling periods, the interpolation adopts linear interpolation of the nearest valid points on both sides. The linear interpolation calculation process is to allocate the weights of the valid values ​​on both sides according to the time ratio based on the time position of the missing measurement point between the two valid points and sum them to ensure the continuity of the interpolation sequence. When the missing measurement length exceeds 2 sampling periods, the window is judged as an invalid window and the mapping calculation is stopped to avoid the interpolation error being transmitted to the trend output.

[0036] Subsequently, a multi-dimensional spatial mapping is performed to obtain the stage-related coordinates of the environment. The mapping input consists of the temporal representation features of the environmental time series set and the current stage indication features. The calculation process of the temporal representation features involves calculating four types of quantities for each environmental parameter within a 360-minute window: the first type is the window mean, obtained by summing all sampled values ​​within the window and dividing by the number of sample points; the second type is the window change rate, obtained by subtracting the window start value from the window end value and then dividing by the window duration; the third type is the window fluctuation amplitude, obtained by taking the difference between the maximum and minimum values ​​within the window; and the fourth type is short-term consistency, obtained by dividing the window into six consecutive 60-minute segments, calculating the mean of each segment, and then calculating the maximum difference between the six segment means. The smaller the maximum difference, the stronger the consistency. The six segments and 60 minutes are obtained by equally dividing the 360-minute window. The stage indication features are expressed using a stage number vector. The construction process of the stage number vector involves establishing a set of position identifiers for a preset number of stages. The position corresponding to the current stage is set to 1, and the other positions are set to 0. The values ​​of 1 and 0 are determined by the definition of the indication meaning. After the splicing is completed, it is input into the mapping model. The mapping model outputs 2D coordinates as the stage environment association coordinates. The process of determining the 2D coordinates is to test 2D, 3D and 4D on the training data respectively, calculate the clustering degree of steady-state samples in the same stage in the coordinate space and the separation degree of samples in different stages. The clustering degree is represented by the average distance from the sample to the center in the same stage, and the separation degree is represented by the minimum distance between the centers in different stages. The dimension with the minimum clustering degree and the maximum separation degree is selected to obtain the 2D coordinates. The training objective of the mapping model is achieved using distance constraints: During the training phase, a steady-state sample set and a transfer sample set are formed for each stage. The selection criteria for the steady-state sample set are that the stage label remains consistent across 12 consecutive sampling points and the window fluctuation amplitudes of the five environmental parameters are all below their corresponding steady-state fluctuation thresholds. The process of determining the 12 consecutive sampling points is to count the shortest duration of the steady-state segment on the training data, and take the number of sampling points corresponding to the shortest duration as the lower limit, resulting in 12. The steady-state fluctuation threshold is obtained by adding twice the standard deviation to the mean of the fluctuation amplitude of the environmental parameter in the steady-state sample window. The process of determining 2 is to test 1, 2, and 3 in the training steady-state samples, and select the value with the fewest number of steady-state windows misclassified as transfer and the smallest threshold, resulting in 2. During the mapping training process, a set of mapping coefficients is calculated to minimize the average distance between steady-state samples in the same stage in the output coordinates, maximize the distance between the steady-state centers of different stages, and make the transfer samples unidirectionally discrete relative to the steady-state centers of their corresponding stages.

[0037] After generating the phase-related coordinates, steady-state deviation judgment is performed and dynamic change trends are calculated. The steady-state center point is determined separately for each phase: the average coordinates of the steady-state sample coordinates for that phase are first calculated as the initial center, and then the distance from each steady-state sample to the initial center is calculated. The samples with the largest distance are sorted by distance and the top 1 / 20 are removed. The process of determining 1 / 20 is to test 1 / 50, 1 / 20, and 1 / 10 on the training steady-state samples, compare the sensitivity of the center point to outliers with the center point offset, and select the proportion with the smallest center point offset and the smallest removal ratio to obtain 1 / 20. After removal, the center is recalculated and repeated once. The number of repetitions is set to 2. The process of determining 2 times is to test 1, 2, and 3 times, compare the change in center point, and terminate when the change in change is less than 1 / 5 of the previous round's decrease after 2 times, thus obtaining 2 times. The deviation threshold is represented by a distance threshold. The process of determining the distance threshold is as follows: calculate the distance from each sample to the steady-state center point in the training steady-state samples to obtain a distance sequence. Take the mean of the distance sequence plus 3 times the standard deviation as the distance threshold. The process of determining 3 is as follows: count the number of false triggers on the training steady-state samples, test 2, 3, and 4 step by step, and select the value with the smallest number of false triggers that does not exceed 1 / 200 of the total number of windows and the smallest threshold value to obtain 3. When the distance from the real-time environment associated coordinates to the steady-state center point exceeds the distance threshold, it is determined that a deviation has occurred and trend calculation is initiated.

[0038] Trend calculation uses time-axis projected displacement output. First, a coordinate trajectory is formed: the stage environment-related coordinates of the most recent 18 sampling points are collected over continuous time. The 18 is determined by statistically analyzing the average number of continuous sampling points from the occurrence of deviation to the actual change of the stage label on the training data, and then rounding this average up to obtain 18. Second, instantaneous displacement is calculated: the lateral and longitudinal differences are calculated for each adjacent coordinate point, and then the two differences are combined into a displacement length. The displacement length is combined by squaring each point separately, adding them, and then taking the square root to represent the planar distance. The squaring and square root operations are performed by the numerical calculation module according to standard operational rules. Third, cumulative displacement is calculated: the displacement length at each moment is accumulated in chronological order to obtain the cumulative displacement sequence from the trajectory starting point to the current moment. The projection direction is then determined: the projection direction is determined by the principal migration direction of this stage. The process of determining the principal migration direction is as follows: collect sample trajectories of migration from this stage to the next stage in the training data, calculate the direction vector from the steady-state center point of the starting point to the migration center point of the ending point of each trajectory, take the average value of the migration sample coordinates of the migration center point, and then perform normalization on the direction vector. The normalization process is as follows: first calculate the length of the direction vector, then divide each component of the direction vector by the length to make the length of the direction vector become 1, where 1 is determined by the definition of unit length; average all normalized direction vectors to obtain the principal direction, and then perform normalization on the principal direction to obtain the final projection direction. Finally, the projected displacement is calculated: for the cumulative displacement vector at each moment, the projected length is calculated along the projection direction. The calculation process of the projected length is to sum the components of the cumulative displacement in the same direction along the projection direction to obtain a scalar result. A positive scalar indicates that the deviation along the migration direction is enhanced, and a negative scalar indicates that the deviation is regressed to the steady state direction. The scalar is output in time order to form a dynamic change trend sequence of the growth stage, and short-term smoothing is performed on the trend sequence to eliminate single-point fluctuations. The smoothing window length is set to 3 sampling points. The process of determining 3 is to test 2, 3, 4, and 5 on the validation data, compare the number of peaks of the trend curve with the response delay, and select the length with the fewest peaks and a response delay of no more than 1 sampling period to obtain 3.

[0039] S3 includes acquiring algal density collected by sensors and calculating biomass growth rate; if the biomass growth rate is lower than a preset dynamic threshold range, retrieving the original spectral data at the corresponding time point; performing principal component analysis on the original spectral data to extract dimensionality-reduced feature vectors; inputting the dimensionality-reduced feature vectors into a pre-trained random forest regression model to obtain predicted concentration values; constructing a concentration time series based on the distribution pattern of the predicted concentration values ​​on the time axis; and using the concentration time series to fit the material accumulation trajectory and determine the concentration values ​​of antioxidant substances.

[0040] In this embodiment, algal density is acquired by an online density sensor within the culture unit. The density sampling period is set to 10 minutes, determined based on the combined results of the density signal's rate of change during the slowdown phase and data processing delay constraints. During model training, the same batch of historical data is replayed at sampling periods of 5 minutes, 10 minutes, and 15 minutes, respectively. The fluctuation amplitude of biomass growth rate and the advance amount of triggering the slowdown determination are statistically analyzed. The fluctuation amplitude is represented by the average absolute value of the difference between two adjacent rate values, and the advance amount is represented by the time difference between the triggering time and the offline biomass inflection point. The 10-minute interval corresponds to the lowest fluctuation amplitude and an advance amount no later than the 15-minute interval. The density sensor output first completes baseline correction. Baseline correction involves collecting culture medium reference values ​​at the same time each day to form the daily baseline. The baseline update cycle is set to 1 day, determined based on the statistical results of the cumulative rate of sensor zero-point drift. The intraday drift amplitude and interday drift amplitude are calculated from 7 consecutive days of offline operation data. When the interday drift is significantly higher than the intraday drift, a 1-day update is adopted. The net density value is obtained by subtracting the daily baseline from the real-time density value, and the net density value is then used for biomass conversion.

[0041] Biomass conversion was performed using the density-to-dry weight mapping relationship established during the training phase. This mapping relationship was fitted using a second-order polynomial. The determination of the second-order polynomial involved establishing first-order, second-order, and third-order fitting relationships for the same training data, calculating the conversion error on the validation data one by one. The conversion error was characterized by the mean absolute difference between the online converted dry weight and the offline weighed dry weight. Simultaneously, the monotonicity of the fitting relationship was checked across the entire range; the second-order relationship had the smallest error and met the monotonicity requirements. The training samples used for fitting came from paired data at multiple time points under the same culture conditions. Each paired point consisted of the net density value and the offline dry weight value at the same time point. The offline dry weight value was obtained by filtering, drying, and weighing a fixed volume of algal solution. The fixed volume was set at 50 ml, determined based on the repeatability of offline weighing and the constraint on disturbance to the culture system: sampling volumes of 30 ml, 50 ml, and 70 ml were tested respectively. The 50 ml sample volume showed the lowest relative fluctuation in weighing, and the proportion of a single sampling volume to the total culture volume was less than 1 / 100. To suppress the amplification of instantaneous density noise on the conversion, the net density value is subjected to 3-point weighted smoothing before conversion. The 3 points correspond to a 30-minute time span. The process of determining the 3 points is to test 3-point, 5-point, and 7-point smoothing on the training data respectively, and compare the number of peaks and trend lag of the dry weight sequence after conversion. The number of peaks is counted as the number of times the absolute value of the difference between adjacent points exceeds the upper limit of the training steady-state difference. The trend lag is counted as the delay of the inflection point of the converted dry weight relative to the inflection point of the original net density. The number of peaks corresponding to the 3 points is significantly reduced and the lag does not exceed 10 minutes.

[0042] The biomass growth rate was calculated based on a smoothed converted biomass sequence, with a calculation window set to 60 minutes. The 60-minute window was determined to provide robustness for slowdown identification. The rate was calculated using 30-minute, 60-minute, and 90-minute windows on the training data. The rate was calculated by subtracting the biomass at the beginning of the window from the biomass at the end of the window to obtain the net growth. The net growth was then converted into growth per minute according to the window duration. Subsequently, the number of false triggers and the slowdown identification delay were statistically analyzed. The number of false triggers was counted as the number of times the slowdown determination was triggered within the offline non-slowdown interval. The identification delay was counted as the lag time between the trigger time and the offline biomass inflection point. The 60-minute window had the lowest number of false triggers and the identification delay did not exceed 30 minutes. To avoid single-point anomalies affecting the window endpoint difference, the instantaneous value of the window endpoint biomass is replaced by the median within the window. The median is selected by sorting all sampling points within the window and taking the middle value. The median replacement trigger condition is that the deviation of the endpoint value relative to the window median exceeds the training steady-state deviation upper limit. The steady-state deviation upper limit is obtained by adding 3 times the standard deviation to the mean of the endpoint deviations of the training steady-state interval. The process of determining 3 is to test 2, 3, and 4 and count the proportion of the steady-state interval that is replaced. The replacement proportion does not exceed 1 / 200 and the minimum upper limit corresponds to 3.

[0043] When biomass growth rate is used to trigger slowdown determination, a dynamic threshold range is employed. The dynamic threshold range consists of an upper threshold and a lower threshold. The upper threshold indicates entry into the slowdown warning zone, and the lower threshold indicates successful slowdown. Both the upper and lower thresholds are determined separately for each growth stage. The determination process involves summarizing all growth rate samples for each growth stage during the training phase, calculating the stage mean and stage standard deviation. The upper threshold is set as the stage mean minus one stage standard deviation, and the lower threshold is set as the stage mean minus two stage standard deviations. The determination of thresholds 1 and 2 involves testing the upper threshold multiples of 0.5, 1, and 1.5 and the lower threshold multiples of 1, 2, and 3 on the validation data, respectively, and counting the number of false triggers and missed triggers. The number of false triggers is defined as the number of times slowdown is successfully triggered in the non-slowdown zone, and the number of missed triggers is defined as the number of times slowdown is not successfully triggered in the slowdown zone. An upper threshold of 1 and a lower threshold of 2 correspond to the simultaneous minimum of both the number of false triggers and the number of missed triggers. To reflect dynamic updates, the stage standard deviation is corrected during online operation based on the fluctuation range of the rate sequence over the most recent 24 hours. The 24-hour period is determined based on the coverage requirement of daily environmental disturbances on growth rate fluctuations: rate autocorrelation is calculated in the training data to locate the main period, and 24 hours is selected when the main period falls within the 20-26 hour range. The correction method is to calculate the ratio of the rate standard deviation of the most recent 24 hours to the standard deviation of that stage during the training phase, and multiply this ratio by the current threshold interval width to obtain the online threshold interval width. The slowdown trigger adopts a continuity constraint, with the number of consecutive points set to 3. The process of determining 3 is to test consecutive points of 2, 3, and 4 on the validation data and compare the number of false triggers with the recognition delay. 3 corresponds to a significant decrease in the number of false triggers and a delay of no more than 20 minutes. The triggering rule is that the growth rate is below the upper limit threshold for 3 consecutive sampling points, and at least one of the sampling points is below the lower limit threshold. After meeting these conditions, the spectral regression prediction link is entered.

[0044] After the trigger is established, the system retrieves the original spectral data at the trigger time, with a retrieval window length set to 120 minutes. The 120-minute window is determined based on the need for short-term spectral evolution information in principal component extraction: spectral sample matrices are constructed on the training data with windows of 60 minutes, 120 minutes, and 180 minutes respectively, and a homogeneous regression model is trained. The validation error and computational delay are compared. The validation error is represented by the mean absolute difference between the predicted concentration and the offline concentration, and the computational delay is represented by the processing time of a single window. 120 minutes corresponds to the lowest error and a delay not exceeding the sampling period. After the original spectral data is retrieved, a consistency processing is performed, including dark current subtraction and reference spectrum correction: the dark current is obtained from the sensor output under shading conditions, with an update cycle set to 12 hours. The 12-hour period is determined based on the semi-diurnal scale statistical results of dark current drift; the reference spectrum is obtained from the same optical path under algae-free culture conditions and is used to correct the gain differences in each band. The correction method is to divide the real-time intensity by the reference intensity for each band to obtain the relative intensity, and then scale the relative intensity to a uniform range according to the same dimensions during the training phase. Subsequently, a quality threshold screening was performed. The saturation threshold was set to 0.98 of the sensor's full scale. The process of determining 0.98 was to statistically analyze the spectral distortion probability when the sensor was close to full scale in the training data. When the inflection point where the distortion probability increased significantly fell in the range of 0.97 to 0.99, 0.98 was selected. If any key band exceeded the saturation threshold, the spectrum at that moment was judged as invalid and removed.

[0045] Principal component extraction (PCE) constructs a sample matrix using the effective spectra within a window. The sample matrix is ​​arranged chronologically, with time sampling points as rows and band channels as columns. The sample matrix is ​​first centered by channel. This centering process involves averaging all sampled values ​​within the window for each band channel to obtain the channel mean. Then, each sampled value is subtracted from the corresponding channel mean to make the channel mean zero, thus eliminating the bias of the channel's absolute magnitude on the principal direction. Next, the co-variation structure between channels is calculated. This is done by taking the centered sequences of any two channels, multiplying them point-by-point, and averaging them to form a characterization of the correlation strength between channel pairs. All channel pairs are then summed to form the co-variation structure. This co-variation structure is decomposed to obtain a set of principal directions sorted by contribution. The contribution is calculated by calculating the explanatory power of each principal direction and sorting them from largest to smallest. The number of principal directions is set to 10. The process of determining 10 involves training a random forest regression on the training data with principal directions of 6, 8, 10, 12, and 14 respectively, and then statistically analyzing the prediction error and output fluctuation on the validation data. Output fluctuation is represented by the standard deviation of the difference between two adjacent prediction points. A value of 10 corresponds to the lowest prediction error and output fluctuation no higher than schemes with 12 and 14. The calculation process for the dimensionality reduction feature vector involves projecting the centered spectrum at each time step along the 10 principal directions to obtain 10 projection coefficients. These projection coefficients are sorted by principal direction to form a vector, and dimensionality standardization is performed on the 10 coefficients. The standardization parameters are the mean and standard deviation of the projection coefficients in that dimension during the training phase, which are calculated by summing the training projection coefficients.

[0046] The dimensionality-reduced feature vectors are fed into the random forest regression model to output predicted concentration values. The training samples for the random forest regression model consist of paired dimensionality-reduced feature vectors and offline concentrations. The offline concentrations are derived from chemical analysis results at corresponding times within the same batch. The range of the offline concentration is determined by the minimum and maximum values ​​of the training samples, and this range is used for subsequent quality discrimination threshold normalization. Random forest parameters are determined during the training phase through discrete search. The number of trees is set to 300. The determination of 300 involves training models with 100, 200, 300, 400, and 500 trees respectively and calculating the validation error. When the number of trees increases from 300 to 400, the decrease in validation error is less than 1 / 5 of the decrease when increasing from 200 to 300, and the prediction time increases by more than 30%, so 300 is selected. The maximum depth of a single tree is set to 12. The determination of 12 involves testing depths of 8, 10, 12, 14, and 16 and comparing... Comparing validation error and output fluctuation, scheme 12 corresponds to the lowest error with fluctuation no higher than schemes 14 and 16. The minimum sample size for leaf nodes is set to 5. The process of determining 5 is to test 3, 5, and 7 and compare the validation error with the number of jitters near the peak. The number of jitters is counted based on the number of times the sign of adjacent points near the peak changes repeatedly. Scheme 5 corresponds to the lowest number of jitters. The number of features participating in the comparison for each split is set to 4. The process of determining 4 is to test 2, 3, 4, and 5 with a principal direction of 10 and compare the validation error. Scheme 4 corresponds to the lowest validation error. During training, sampling for each tree is performed using sampling with replacement. The sampling size is set to 1 times the total number of training samples. The determination of 1 times is based on the balance between ensemble variance and bias: 0.7 times, 1 times, and 1.3 times are tested on the training data respectively. Scheme 1 times corresponds to the lowest validation error. During online prediction, 300 trees output concentration prediction values, and the final predicted concentration is the arithmetic mean of the 300 outputs. Simultaneously, the standard deviation of the tree outputs is calculated as the inter-tree divergence value. The inter-tree divergence threshold is set to 0.03 times the training concentration range. The determination of 0.03 involves testing 0.02, 0.03, and 0.04 on validation data and counting the number of false positives and false negatives for unstable markers. False positives are defined as the number of times the inter-tree divergence exceeds the threshold but the prediction error remains within an acceptable range, while false negatives are defined as the number of times the inter-tree divergence does not exceed the threshold but the prediction error increases significantly. 0.03 corresponds to the lowest number of both false positives and false negatives. When the inter-tree divergence exceeds the threshold, the prediction value at that moment is considered invalid and is not included in the time series fitting.

[0047] Predicted concentration values ​​are used to construct a concentration time series in chronological order, with the time interval consistent with the sampling period of 10 minutes. The time series is first processed for missing values, with an upper limit of 2 sampling periods allowed. This 2 is determined by artificially creating missing segments on the training data and performing subsequent trajectory fitting. The fitting error is statistically analyzed as the missing length changes; a significant jump in fitting error occurs when the missing length reaches 3 sampling periods, thus setting the upper limit to 2. For missing lengths not exceeding 2 sampling periods, linear interpolation is used to restore continuity. The linear interpolation process involves locating the two nearest valid concentration values ​​before and after the missing point and generating interpolated values ​​proportional to the time interval between the missing point and the preceding and following valid points. This generation process ensures the interpolated sequence is continuous at both ends of the missing segment. For missing lengths exceeding 2 sampling periods, this segment of the time series is not included in trajectory fitting and is marked as invalid. The time series then undergoes consistency screening. The jump threshold for consistency screening is set as the upper limit of single-step concentration change during the training phase. The upper limit of single-step concentration change is calculated from the absolute value sequence of the concentration difference between adjacent valid points in the training sample. The mean of the sequence is taken and three times the standard deviation is added to form the threshold. The process of determining the threshold is to test 2, 3, and 4 and count the retention rate of the true peak and the rejection rate of abnormal jumps. The retention rate and rejection rate are both at their highest corresponding to the threshold. If the absolute value of the concentration difference between a certain time and the concentration at the previous time exceeds the threshold, the concentration at that time is judged as abnormal and processed according to the missing point rule.

[0048] The concentration time series is then fitted to the material accumulation trajectory. The fitting objective is to reduce short-term fluctuations and retain true turning points without altering the time sequence. Piecewise constrained smoothing is employed: first, a 5-point weighted smoothing is performed on the entire sequence to obtain the initial trajectory. The 5 points are determined by testing points 3, 5, and 7 on the training data and comparing peak position shifts with noise residue. Peak position shifts are represented by the time difference between the smoothed peak moment and the original predicted peak moment, while noise residue is represented by the standard deviation of the difference between adjacent points after smoothing. The peak position shift corresponding to the 5 points should not exceed 10 minutes, and noise residue should be minimized. Turning point localization is then performed based on the stability of continuous difference sign: the sign difference between adjacent two points is calculated point by point on the initial trajectory. When the sign difference changes from positive to negative and is maintained for 3 consecutive sampling points, this moment is marked as the peak turning point. The number of consecutively maintained points is set to 3. The determination of 3 involves statistically analyzing the duration of brief sign reversals caused by noise on the training data. When the maximum duration falls within 2 sampling points, 3 is chosen to exclude noise reversals. The sequence is divided into growth and decline segments based on peak inflection points. The growth segment is subject to monotonic constraints: if the current point is less than the previous valid point within the growth segment, the current point is replaced by the previous valid point. The decline segment is subject to reverse monotonic constraints: if the current point is greater than the previous valid point within the decline segment, the current point is replaced by the previous valid point. A change magnitude check is added to the replacement trigger condition, with the check threshold using the aforementioned upper limit of single-step concentration change to avoid over-constraining real rapid changes. After segmentation constraints are completed, local smoothing is performed at the connection between the end of the growth segment and the beginning of the decline segment. The local smoothing window length is set to 3 sampling points. The process of determining 3 involves testing 3 and 5 on the training data and comparing the number of peaks at the connection point with the peak time offset. 3 corresponds to the lowest number of peaks and an offset not exceeding 10 minutes.

[0049] The fitted trajectory is output as a sequence of antioxidant concentration values. The output includes the fitted concentration value and key moment indicators for each sampling time. The key moment indicators include peak concentration and peak time. The peak concentration is determined by the maximum value of the fitted trajectory, and the peak time is determined by the timestamp corresponding to the maximum value. At the same time, the peak time of the accumulation rate in the interval from 0 minutes to 120 minutes after the slowdown trigger is also output. The accumulation rate is calculated by calculating the concentration difference between adjacent sampling points of the fitted trajectory and converting it to an increment per minute in 10-minute increments. Then, the time corresponding to the maximum increment in this interval is taken as the peak time of the accumulation rate.

[0050] S4 includes acquiring continuous time series data of antioxidant concentration values, performing sliding segmentation on the continuous time series data according to a preset window step size to obtain multiple overlapping local concentration intervals; performing linear fitting on sampling points within the local concentration intervals to calculate concentration gradient changes; if the positive or negative sign of the concentration gradient change is reversed and the absolute value exceeds a preset threshold, then extracting the slope extreme points within the local concentration interval as candidate points for potential peak moments; performing clustering processing based on the distribution density of candidate points on the time axis, removing outliers caused by time series fluctuations to determine the final candidate point set; and using the sampling frequency corresponding to the candidate point set to perform local data smoothing to obtain a candidate point set that can reflect the peak of antioxidant accumulation.

[0051] In this embodiment, continuous time series data of antioxidant concentration values ​​are first acquired, and time scale consistency processing is performed on the time series to ensure consistency between the time intervals used for subsequent segmentation and fitting. This time series is generated from the upstream concentration prediction output, with a sampling interval of 10 minutes. The 10-minute interval is determined by the combined constraints of the output beat of the concentration prediction link and the peak inflection resolution requirement: The steepest change segment of the concentration curve within a 2-hour range before and after the actual peak is statistically analyzed in the training samples. In most samples, the duration is not less than 40 minutes. Setting the sampling interval to 10 minutes ensures that at least four sampling points are formed within this change segment, thereby ensuring sufficient resolution for slope sign reversal. The time scale processing is implemented as follows: the time difference between two adjacent sampling times is calculated point by point, and the time difference with the highest frequency is taken as the master time difference. When a time difference deviates from the master time difference by more than 2 minutes, the sampling point is marked as a time series outlier and removed from subsequent fitting samples. The 2-minute interval is determined by adding a 1-minute margin to the statistical upper bound of the acquisition link clock jitter. This margin is used to cover short-term fluctuations in communication delay. The removed gaps are handled according to the missing point rules: when the missing length is no more than 2 sampling periods, linear interpolation is performed; when the missing length exceeds 2 sampling periods, the segment is marked as an unusable segment and the output of candidate points in the segment is stopped. The determination of 2 sampling periods is based on constructing missing segments in the training samples and observing the effect of linear interpolation on the peak time offset. When the missing length reaches 3 sampling periods, the peak time offset exceeds 20 minutes, and the upper limit of the missing length is set to 2 sampling periods to control the offset within 20 minutes.

[0052] After time series processing, the continuous time series is divided into multiple overlapping local concentration intervals based on a preset window length and window step. The window length is set to 12 sampling points, corresponding to 120 minutes. The 120-minute window is determined based on the peak width statistics of the training samples: in the training samples, with the offline peak time as the center, a continuous duration with a concentration not less than 0.9 times the peak concentration is searched forward and backward, and this continuous duration is used as the peak width sample; the median width of all peak width samples is calculated to obtain approximately 240 minutes, and the window length is set to half of the median width to obtain 120 minutes. This ensures that the window covers the main area of ​​peak transition while avoiding the slope from being canceled out by different trends before and after the peak due to an excessively long window. The window step size is set to 2 sampling points, corresponding to 20 minutes. The 20-minute step size is determined based on the candidate point localization error control requirements: in the training samples, the segmentation process is replayed with different step sizes, and the absolute value of the offset between the candidate point time and the offline peak time is counted. When the step size is 20 minutes, the proportion of samples with an absolute offset value not exceeding 20 minutes is not less than 19 / 20, and the number of false detections is minimized. Therefore, the step size is set to 20 minutes to form highly overlapping intervals and improve the localization density near the peak. The sliding segmentation process is as follows: the starting point of the first local concentration interval is placed at the beginning of the time series, and 12 sampling points are taken consecutively to form the interval; the starting point of the next interval is moved back 20 minutes along the time axis, and 12 sampling points are taken consecutively again to form the interval; this process is repeated until the end of the interval covers the end of the time series. If there are fewer than 12 sampling points at the end, the insufficient segment at the end is marked as an unusable segment and does not participate in the fitting.

[0053] For each local concentration interval, linear fitting is performed to calculate the concentration gradient change within that interval, and this gradient is used as the interval slope for subsequent sign inversion determination. To avoid single-point anomalous jumps affecting the fitting results, the concentration difference between adjacent sampling points within the interval is first screened: the concentration difference between each pair of adjacent sampling points within the interval is calculated, and its magnitude is recorded; if the magnitude exceeds the single-step change upper limit, the latter sampling point is marked as an outlier and removed from the fitted sample. The single-step change upper limit is determined based on the true rate of change boundary of the training sample: adjacent concentration differences are extracted from the training sample from the intervals from 2 hours before the peak to 2 hours after the peak to obtain a single-step change sample set. The mean and standard deviation of this set are calculated, and the mean plus 3 times the standard deviation is set as the single-step change upper limit; the determination of 3 is based on the error deletion control requirements. In the training sample, 2, 3, and 4 are used as multiples respectively, and the number of points removed in the neighborhood of the true peak is counted. When the multiple is 3, the proportion of removed points is no higher than 1 / 200, and the single-step change upper limit is the smallest, so 3 is used. After outlier removal, each sampling point within the interval is converted into a time sequence number, incrementing from 0 to the end of the interval. Then, the linear fitting slope is calculated: first, the average of the time sequence numbers and the average of the concentration values ​​within the interval are calculated; then, for each sampling point within the interval, the difference between the time sequence number and its average, and the difference between the concentration value and its average are calculated; the two differences for each sampling point are multiplied and summed to obtain the cumulative multiplicative measure; the time sequence number difference for each sampling point is multiplied by itself and summed to obtain the cumulative squared measure; the cumulative multiplicative measure is divided by the cumulative squared measure to obtain the interval slope. The sign of the interval slope indicates whether the interval is rising or falling overall, and the magnitude of the interval slope indicates the strength of the rise or fall. To ensure comparability of slopes across different algal strains, the slope is scaled uniformly according to the sampling interval: the slope is converted to a change per minute by dividing the interval slope by 10 minutes, ensuring a uniform unit of slope across all samples.

[0054] After obtaining the slopes of all local concentration intervals, a slope sequence is formed in chronological order, and a joint determination of sign flip and amplitude threshold is performed on the slopes of adjacent intervals. The execution process of sign flip determination is as follows: the slopes of two adjacent intervals are taken in sequence. When the slope of the previous interval is positive and the slope of the next interval is negative, a flip from rising to falling is determined; when the slope of the previous interval is negative and the slope of the next interval is positive, a flip from falling to rising is determined. The amplitude threshold is used to exclude false flips caused by low-amplitude fluctuations. The amplitude threshold is set to be no less than twice the upper limit of the training noise slope per minute. The upper limit of the training noise slope is obtained by statistical analysis of the non-peak segments of the training samples: the training samples are removed from the 3-hour segments before and after the offline peak time, and the remaining segments are taken as non-peak segments; the interval slope of the non-peak segments is calculated and converted into the change per minute using the aforementioned method, and its mean and standard deviation are calculated. The mean plus three times the standard deviation is set as the upper limit of the training noise slope; the multiplier of 3 is determined based on the false trigger control requirements. In the non-peak segments of the training samples, the upper limit of the noise slope is calculated using 2, 3, and 4 respectively, and the number of false flip triggers is counted. When the multiplier is 3, the number of false flip triggers is no higher than 1 / 500 and the upper limit of the noise slope is the smallest. Therefore, 3 is used. The amplitude threshold is set to twice the upper limit of the noise slope. The determination of 2 is based on the requirement of preserving true flips: in the training samples, the flip slope values ​​in the neighborhood of the true peak are statistically analyzed. The proportion of samples with values ​​lower than twice the upper limit of the noise slope is no more than 1 / 50. Setting the threshold to twice can suppress false flips while maintaining a true flip retention rate of no less than 49 / 50. The joint determination trigger rule is: when the slopes of two adjacent intervals flip, and the values ​​of the slopes of the two intervals before and after the flip are not lower than the amplitude threshold, the potential peak neighborhood corresponding to the flip is determined and entered into the candidate point extraction.

[0055] After extracting the candidate trigger points, the slope extreme points are located within the overlapping time range of the two intervals that trigger the flip and are used as potential peak moment candidate points. The location of the slope extreme points adopts the criterion of maximizing the local slope change. The specific implementation process is as follows: within the overlapping range, the local differential slope is calculated point by point according to the sampling points. The local differential slope is taken as the concentration of the current sampling point minus the concentration of the previous sampling point and converted into the change per minute. Then, the local differential slope is smoothed by a short window. The short window length is set to 3 sampling points, corresponding to 30 minutes. The 30 minutes is determined based on the statistical results of the main fluctuation period of the concentration prediction noise in the training samples: the autocorrelation of the local differential slope in the non-peak segment is calculated and the main peak is located. When the period corresponding to the main peak falls in the range of 20 minutes to 40 minutes, 30 minutes is taken as the smoothing span to suppress the peak reversal within the noise period. The smoothing process involves calculating a weighted average of three consecutive local difference slopes. The weights are set to 3, 2, and 1 respectively, from closest to furthest in time. The weights are determined based on the principle that the center point contributes the most and the contribution of the edge points decreases, ensuring that the slope retains its inflection point after smoothing. After smoothing, the change between two adjacent smoothed slopes is calculated point by point. The change is calculated as the difference between the previous and subsequent smoothed slopes, and the magnitude of the change is recorded. Within the overlapping range, the sampling point with the largest change is found, and this sampling point is selected as a candidate for the slope extremum. If multiple points have the same maximum change, the concentration value is used as the criterion, and the moment with the largest concentration value is selected as the candidate for the slope extremum. The concentration value criterion ensures that the candidate point is closer to the accumulated peak rather than local slope abrupt noise. To prevent boundary effects from causing extreme points to fall at both ends of the overlapping range, additional boundary constraints are added: the time interval between the candidate point and the boundary of the overlapping range is not less than 20 minutes. The 20-minute interval is determined by a window step size of 20 minutes, so that the candidate point contains at least one complete step of information before and after. If the candidate point violates the boundary constraints, the sampling point with the second largest change is selected and the decision is repeated until the boundary constraints are met.

[0056] After forming the initial candidate point set, clustering is performed based on the distribution density of candidate points on the time axis, and outliers are removed to obtain the final candidate point set. The distribution density is calculated using neighborhood counting: the time neighborhood radius is set to 30 minutes. The 30-minute timeframe is determined based on the statistical results of the time spread width of candidate points in the real peak neighborhood in the training samples. In the training samples, the proportion of candidate points falling within the 30-minute range before and after the offline peak time is counted. This proportion is not less than 19 / 20. The 30-minute timeframe is set as the neighborhood radius to cover the main candidate points in the real peak neighborhood. For each candidate point, the number of candidate points contained within its 30-minute range before and after is counted as the neighborhood count of that candidate point. The minimum neighborhood count is set to 3. The 3 is determined based on noise outlier suppression: the same candidate point extraction process is run in the non-peak segment of the training samples. The number of events in which any point in the non-peak segment forms 3 candidate points within its 30-minute neighborhood is counted. This number of events is not higher than 1 / 500. The minimum neighborhood count is set to 3 to make it difficult for outliers to pass the density threshold. The clustering and removal process is as follows: First, candidate points with a neighborhood count less than 3 are removed. The remaining candidate points are sorted by time. If the time interval between two adjacent candidate points does not exceed 30 minutes, they are grouped into the same cluster; if the time interval exceeds 30 minutes, a new cluster is created. For each cluster, the cluster center time is calculated, taking the average of the times of all candidate points within that cluster. Then, the time distance from each candidate point within the cluster to the cluster center time is calculated, and the allowable radius within the cluster is set to 20 minutes. The 20-minute radius is determined by adding twice the standard deviation to the mean of the time distances within clusters in the training samples, forming an upper bound. The multiple of 2 is used to reduce false deletions while maintaining cluster compactness. Candidate points exceeding 20 minutes are removed from the cluster. If the number of remaining points in a cluster is less than 3 after removal, the entire cluster is removed.

[0057] Local data smoothing and consistency verification matching the sampling frequency are performed on the final candidate point set to output a candidate point set that reflects the peak of antioxidant accumulation. Local smoothing extracts a neighborhood window centered on the candidate point, with a window span of 5 sampling points, corresponding to 50 minutes. The 50-minute timeframe is determined based on the statistical results of the curve bending radius near the peak: the intensity of the secondary change in the concentration curve within a 25-minute range before and after the actual peak time in the training samples is calculated. This range can cover the segment with the most obvious bending. Setting the window span to 50 minutes allows the smoothing to focus on the bending segment without introducing distant trends. The smoothing calculation process is as follows: Two sampling points are taken before and after the candidate point to form a 5-point sequence with the candidate point itself. These 5 concentration values ​​are assigned weights decreasing from the center outwards, with weights set to 5, 4, 3, 2, and 1 respectively. The weights are set based on the principle that the center point contributes the most and the contributions of the edge points gradually decrease. The concentration of each point is multiplied by its weight and summed. The sum is then divided by the total weights to obtain the smoothed concentration value of the candidate point. This smoothing calculation is repeated for each sampling point in the neighborhood to obtain a smoothed sequence. After smoothing, a consistency check is performed: Using the candidate point as the center, the local difference slope is recalculated for the two sampling points to the left and two to the right of the candidate point and converted into a change per minute. The local difference slope on the left must be positive and its value must not be lower than the aforementioned amplitude threshold; the local difference slope on the right must be negative and its value must not be lower than the aforementioned amplitude threshold. Simultaneously, the smoothed concentration value of the candidate point must be the maximum value in the 5-point smoothed sequence. Candidate points that do not meet any of these conditions are directly eliminated.

[0058] S5 includes obtaining the time coordinates in the candidate point set and constructing a one-dimensional feature vector space, calculating the Euclidean distance between different time coordinates in the one-dimensional feature vector space; using the K-means algorithm based on the Euclidean distance to iteratively divide the candidate point set, obtaining peak time groups where the cluster center no longer shifts and has cluster labels; retrieving the concentration values ​​within the peak time groups and extracting the sampling time with the largest concentration value as the local extremum; mapping the local extremum to the time axis of the algal cycle to identify the growth stage, and determining the peak content of antioxidant substances of the corresponding algal strain in its culture cycle by comparing the local extremums of different growth stages.

[0059] In this embodiment, based on the obtained candidate point set, the time coordinates of each candidate point are extracted and a one-dimensional feature vector space is constructed. The time coordinates are uniformly expressed using relative time values, calculated in minutes. The calculation process is as follows: read the timestamp of the candidate point and the timestamp of the algal strain's incubation start time; subtract the incubation start timestamp from the candidate point's timestamp to obtain the relative time difference; then convert the relative time difference into minutes to form a relative time value. The determination of using minutes as the unit is based on the combined result of the sampling interval and cluster resolution constraints: when the sampling interval is 10 minutes, using minutes can obtain integer or half-integer level resolution when calculating the center displacement and inter-cluster interval, avoiding the amplification of discretization errors caused by using hours. After completing the relative time value conversion, the relative time values ​​of each candidate point are arranged in the order of appearance to form a one-dimensional feature vector set, thereby forming a one-dimensional feature vector space.

[0060] In a one-dimensional feature vector space, the Euclidean distance between different time coordinates is calculated using the absolute value of the time difference. The specific calculation process is as follows: Two candidate points are randomly selected, their relative time values ​​are read, and the difference between the two is obtained. The absolute value of this time difference is then taken to obtain the Euclidean distance between the two points. This calculation is repeated for all candidate point pairs within the candidate point set, forming the basis of the distance metric. Using absolute values ​​as the basis for distance is consistent with the one-dimensional Euclidean distance definition and the monotonicity of the time axis. The distance results correspond one-to-one with the degree of separation of candidate points on the time axis, and can be directly used for subsequent iterative partitioning.

[0061] Subsequently, the candidate point set is iteratively divided using the mean clustering algorithm based on Euclidean distance to obtain peak time groups. Before iterative division, the number of clusters is determined and the cluster centers are initialized. The number of clusters is determined using a time interval segmentation rule. The calculation process is as follows: the relative time values ​​of candidate points are sorted in ascending order, and the time interval between two adjacent candidate points is calculated by subtracting the previous time value from the next time value and taking the positive value, resulting in a time interval sequence. Large interval positions are identified in the time interval sequence, and the determination of large interval positions is based on an interval threshold. The interval threshold is set to 30 minutes. The determination of 30 minutes is based on the statistical results of the time diffusion width of candidate points in the neighborhood of a single real peak in the training data: the occurrence range of candidate points is counted in the training data with the offline peak time as the center, the difference between the maximum and minimum time values ​​of candidate points is calculated to obtain the diffusion width, and an upper bound is calculated for the diffusion width of all training batches. When the upper bound falls on 60 minutes, the interval threshold is set to half of the diffusion width, resulting in 30 minutes, so that the interval between adjacent candidate points in the same peak neighborhood does not exceed 30 minutes, while the interval across peaks or across fluctuation segments exceeds 30 minutes. The number of clusters is calculated as "the number of large-interval locations plus 1", which makes each segment of consecutive candidate points divided by large intervals form a cluster.

[0062] The initialization of cluster centers adopts a uniform selection rule covering the entire time range. The calculation process is as follows: The minimum and maximum relative time values ​​of candidate points are read, and the full range span is calculated as the maximum value minus the minimum value. The span is then divided equally by the number of clusters to obtain an equal step size, calculated as the span divided by the number of clusters. Starting from the minimum value, initial center time values ​​with the same number of clusters are generated sequentially according to the equal step size, and these initial centers are used as the initial cluster centers. The basis for using uniform selection is that the initial centers cover the entire range, which can reduce the risk of empty clusters caused by centers being concentrated in local areas in the early stages of iteration.

[0063] The iterative partitioning of the mean clustering algorithm is performed in the order of "allocation, update, displacement evaluation, and convergence determination". The allocation step is as follows: for each candidate point, calculate the Euclidean distance from that candidate point to each cluster center, using the absolute value of the time difference; assign the candidate point to the cluster corresponding to the center with the smallest distance, and record a cluster label for the candidate point, using integer numbers starting from 1. The update step is as follows: for each cluster, summarize the relative time values ​​of all candidate points within the cluster and calculate the average value. The average value is calculated by summing the values ​​and then dividing by the number of candidate points. This average value is used as the new cluster center for that cluster. When a cluster does not contain any candidate points after one allocation, the cluster center is reset based on the "time value of the candidate point farthest from all candidate points". The reset calculation process involves calculating the distance from each candidate point to the nearest center, selecting the candidate point with the largest distance as the reset center, so that the empty cluster can attract objects in the next iteration. The calculation process for the displacement assessment step is as follows: calculate the absolute value of the time difference between the center before and after the update for each cluster center as the center displacement, and record the maximum value of all center displacements.

[0064] The convergence determination employs a dual-condition constraint: the first condition is that the center displacement does not exceed the convergence threshold, and the second condition is that the cluster label assignment remains consistent for two consecutive rounds. The convergence threshold is set to 5 minutes, determined based on half the sampling interval and consistent with the center positioning resolution. When the sampling interval is 10 minutes, a center displacement less than or equal to 5 minutes indicates that the center change is less than a single sampling interval, and the center position change will not cause a physical time-span shift in the peak group to which the candidate point belongs. In the training data, the clustering process is replayed with different thresholds, and the offset of the final center position relative to the offline peak time is statistically analyzed. The absolute value of the offset corresponding to 5 minutes remains stable, and the number of iterations does not increase. The calculation process for cluster label consistency determination is as follows: the cluster label of each candidate point generated in the current iteration is compared point-by-point with the corresponding candidate point cluster label from the previous round. If all are consistent, it is recorded as consistent; otherwise, it is recorded as inconsistent. The second condition is satisfied when both rounds record consistency. To ensure algorithm termination, a maximum number of iterations is set to 25. This 25 iterations are determined based on the upper bound of convergence iterations in the training data: the same initialization rule is run on the training data, and the number of iterations required for each batch to achieve biconditional convergence is recorded. The maximum number of iterations is 18. Seven redundancies are added to this maximum to obtain 25 iterations, ensuring that extreme distributions can still converge or terminate at the upper limit. The iteration ends when biconditional convergence is achieved or the number of iterations reaches 25. The algorithm outputs the peak time groups with cluster labels and the corresponding cluster center time values.

[0065] After obtaining the peak time groups, the concentration values ​​corresponding to candidate points are retrieved within each group, and the sampling time with the largest concentration value is extracted as the local extremum. The retrieval process uses the candidate point timestamp as an index, extracting the concentration values ​​corresponding to the same timestamp from the concentration time series to form a "candidate point-concentration value" pairing set. The calculation process for local extremum extraction is as follows: the concentration values ​​of all candidate points within the same group are sorted from largest to smallest, and the candidate point time corresponding to the first position in the sort is determined as the local extremum time, and the concentration value corresponding to the first position in the sort is determined as the local extremum concentration value. When there is a tie for the largest concentration value, a tie-breaking decision is performed. The tie-breaking decision is completed in the following order: first, the time distance from the tied candidate point to the cluster center of the group is calculated. The time distance is calculated as the relative time value of the candidate point minus the center time value and the absolute value is taken, and the one with the smallest time distance is selected; if the time distance is still tied, the one with the smaller relative time value is selected so that the local extremum time is output earlier; if they are still tied, the peak sharpness is calculated and the candidate point with the greater sharpness is selected. The sharpness of the peak is calculated using the variation of a local window. The length of the local window is set to 5 sampling points, corresponding to 50 minutes. The 50-minute interval is determined based on the coverage requirement of the peak transition zone: when the sampling interval is 10 minutes, the 5 sampling points cover 20 minutes before and after the peak and include the peak itself, thus covering the main variation range of the rising and falling segments of the transition. The sharpness calculation process is as follows: taking the first two and last two sampling points as the center of the candidate point, the rise amplitude is calculated as the candidate point concentration minus the smaller concentration value of the first two sampling points, and the fall amplitude is calculated as the candidate point concentration minus the smaller concentration value of the last two sampling points. The rise amplitude and fall amplitude are then added together to obtain the sharpness, and the candidate point with the greater sharpness is selected as the local extreme moment.

[0066] After determining the local extrema, they are mapped to the algal culture cycle timeline to identify growth stages. The local extrema are then compared across different growth stages to determine the peak value for the entire cycle. The mapping process uses timestamp consistency matching: the timestamp of the local extremum is read, and the corresponding stage label is located in the growth stage label sequence. If the local extremum is located at the stage label switching boundary, causing instability in a single label, neighborhood consistency adjudication is used to determine the stage affiliation. The neighborhood length is set to three sampling points before and after, corresponding to 30 minutes before and after. The 30-minute period is determined based on the statistical results of the stage label consistency constraint window: the duration of erroneous stage label transitions is counted in the training data, with a maximum duration of 20 minutes. Setting the neighborhood coverage to 30 minutes ensures complete coverage of erroneous transitions and achieves majority adjudication. The calculation process for neighborhood consistency adjudication is as follows: seven stage labels are taken: the three sampling points before the local extremum, the current time, and the three sampling points after it. The frequency of each stage label is counted, and the one with the highest frequency is determined as the stage affiliation corresponding to the local extremum. If the frequency is tied, the stage label with the same frequency as the local extremum is selected.

[0067] When determining the peak value for the entire cycle, all local extreme values ​​are first grouped according to stage affiliation. Then, the local extreme value with the largest concentration is selected as the peak value of antioxidant content within the culture cycle of that algal strain. The comparison process uses a primary criterion and a tie-breaking method: the primary criterion is the magnitude of the local extreme value concentration; the local extreme value with the larger concentration is directly determined as the peak value for the entire cycle. When there are tied maximum concentration values, the peak times are first compared, and the one with the smaller relative time value is selected, indicating that the peak appeared earlier. If they are still tied, the sharpness of the peak shape is compared, and the one with greater sharpness is determined as the peak value for the entire cycle. The output of the peak value for the entire cycle includes the peak concentration value, the peak time timestamp, the relative time value, and the growth stage label corresponding to the peak value.

[0068] S6 includes obtaining the time-series sequence of antioxidant content of each candidate Spirulina strain during the complete culture cycle; obtaining the content change rate set by performing a first-order difference operation on the antioxidant content time-series sequence; determining the peak value of antioxidant content and its corresponding peak occurrence time using the content change rate set; performing normalization mapping based on the peak value of antioxidant content and the peak occurrence time; importing the mapped values ​​into a preset growth trajectory fitting model; obtaining morphological feature vectors reflecting the dynamic change trajectory throughout the entire cycle by extracting features from the growth trajectory fitting model; calculating the trajectory deviation of each candidate Spirulina strain relative to the benchmark growth model using the morphological feature vectors; if the trajectory deviation exceeds a preset stability threshold, extracting the instantaneous slope of the candidate Spirulina strain in the deviation interval as the activity fluctuation variable; obtaining the comprehensive activity score of each candidate Spirulina strain by performing a weighted matrix operation on the activity fluctuation variable, the peak value of antioxidant content, and the peak occurrence time; and constructing an antioxidant activity evaluation index system using the comprehensive activity score.

[0069] In this embodiment, for each candidate Spirulina strain, the time-series sequence of antioxidant content throughout its complete culture cycle is first obtained, and this sequence is then standardized to the same time scale before being used in the calculation. The time scale is standardized with the start time of culture as the zero point, and each sampling point in the sequence includes a timestamp and the corresponding content value. The sampling interval is set to 10 minutes, which is determined based on the combined constraints of the upstream concentration output beat and the peak transition recognition resolution: the duration of the steepest change segment of the content curve within 2 hours before and after the actual peak in the training data is statistically analyzed, and the minimum duration is not less than 40 minutes. The 10-minute interval can form no less than 4 sampling points within this segment, so that the differential rate and transition discrimination have sufficient resolution. The time consistency verification process is to calculate the time difference between two adjacent timestamps point by point, and to count the time difference that occurs most frequently as the master time difference; when a time difference deviates from the master time difference by more than 2 minutes, the corresponding sampling point is marked as a time series anomaly and removed from the differential calculation input. The 2-minute time difference is determined by adding a 1-minute margin to the statistical upper bound of the acquisition link clock jitter, and the margin is used to cover short-term fluctuations in communication delay. The gaps resulting from the removal process are processed using missing data handling. The upper limit for the missing data length is set to two sampling periods. The determination of two sampling periods is based on the impact assessment of the missing data simulation on the peak time shift: different missing data lengths are constructed in the training data and restored using linear interpolation before peak localization is performed. When the missing data length reaches three sampling periods, the peak time shift exceeds 20 minutes, so the upper limit is set to two sampling periods to control the shift within 20 minutes. When the missing data length does not exceed two sampling periods, linear interpolation is performed. The interpolation calculation process involves locating the two nearest valid content points before and after the missing segment, and generating interpolated values ​​according to the time interval ratio between the missing point and the valid points on both sides, so that the interpolated sequence is continuous at both ends of the missing segment. When the missing data length exceeds two sampling periods, the segment is marked as an unusable segment, and differencing and fitting calculations are stopped in that segment.

[0070] After time series processing, a first-order difference operation is performed on the time series of antioxidant content to obtain a set of content change rates. The difference operation is performed as follows: two adjacent sampling points are taken in chronological order, and the content value of the latter sampling point is subtracted from the content value of the former sampling point to obtain the content change. The content change is then converted into the change per unit time according to the main sampling interval. The conversion process is to divide the content change by 10 minutes and convert it into the change per minute to form the content change rate of the adjacent interval. The above operation is repeated for all adjacent intervals in the entire period to obtain a set of content change rates covering the entire period. To suppress rate spikes caused by single-point noise, the rate set enters a smoothing and anomaly removal process. The smoothing window length is set to 3 rate points, corresponding to 30 minutes. The 30-minute window is determined based on the statistical results of the main fluctuation period of rate noise in the non-peak segment of the training data: the autocorrelation of the rate series in the non-peak segment is calculated and the main peak is located. When the period corresponding to the main peak falls in the interval of 20 minutes to 40 minutes, 30 minutes is taken to cover one noise fluctuation period. The smoothing calculation process involves taking a weighted average of three consecutive rate values, with the center value having the highest weight. The weights are set to 3, 2, and 1 respectively. The weights are set based on the principle that the center point contributes the most and the contribution of the edge points decreases, so that the smoothed rate retains the turning points. The abnormal rate removal threshold is set as the upper limit of the training steady-state rate. The process of determining the upper limit of the training steady-state rate is as follows: remove the three-hour intervals before and after the offline peak time from the training data, and the remaining intervals are regarded as steady-state and non-peak fluctuation intervals; calculate the rate within this interval using the aforementioned difference and reduction method to obtain the steady-state rate sample set; calculate the mean and standard deviation of this set, and add 3 times the standard deviation to the mean as the upper limit of the steady-state rate. The determination of 3 is based on the error deletion control requirements: in the training steady-state interval, 2, 3, and 4 are used as multiples to count the proportion of the true trend segment that is error-deleted. When the multiple is 3, the error deletion ratio is no higher than 1 / 200 and the upper limit is the smallest. 3 is used to balance noise suppression and fidelity preservation. When the absolute value of the real-time rate exceeds the upper limit of the steady-state rate, the original adjacent point pair corresponding to that rate is marked as an abnormal pair and removed from the rate consistency statistics confirmed by the peak.

[0071] When determining the peak content and peak occurrence time of antioxidant substances using the set of content change rates, peak candidate localization is performed first, followed by rate consistency confirmation. The calculation process for peak candidate localization is as follows: traverse the content values ​​of each sampling point in the full-cycle content sequence, record the current maximum content value and its corresponding timestamp, and after traversal, obtain the maximum content value of the entire cycle as the peak candidate value, and obtain the corresponding timestamp as the peak occurrence time candidate. Rate consistency confirmation is used to eliminate misjudgments caused by plateau segments and local noise spikes. Its implementation process involves setting a verification window around the peak occurrence time candidate and statistically analyzing sign consistency. The verification window length is set to 12 sampling points, corresponding to 120 minutes. The 120-minute timeframe is determined based on the coverage requirements of the main change segments before and after the peak transition in the training data: statistically analyze the length of the continuous interval where the rate changes from positive to negative before and after the peak transition in the training data. The upper bound of the interval length falls between 90 and 110 minutes. The window is set to 120 minutes to cover the complete transition interval and leave a margin of 10 to 30 minutes. The validation window is formed by taking 6 sampling intervals forward and 6 sampling intervals backward from the candidate time point to create a first half-window and a second half-window. The calculation process for symbol consistency statistics is as follows: in the rate set corresponding to the first half-window, the number of positive rates is counted and divided by the number of rate points in the first half-window to obtain the rising ratio of the first half-window; in the rate set corresponding to the second half-window, the number of negative rates is counted and divided by the number of rate points in the second half-window to obtain the falling ratio of the second half-window. The thresholds for both the rising and falling ratios are set to 0.7. The determination of 0.7 is based on the consistency distribution of the true peak neighborhood in the training data and the misjudgment control: in the training data, the rising and falling ratios are calculated for the true peak neighborhood, and the same ratio is calculated for non-peak fluctuation segments. The thresholds of 0.6, 0.7, and 0.8 are traversed, and the true peak retention rate and non-peak false entry rate are counted. 0.7 corresponds to a retention rate of no less than 49 / 50 and a false entry rate of no more than 1 / 100. When the rising proportion of the first half-window is not less than 0.7 and the falling proportion of the second half-window is not less than 0.7, the candidate time is confirmed as the peak occurrence time, and the content value at that time is confirmed as the peak content of antioxidant substances. If the confirmation condition is not met, the selected candidate points are temporarily excluded from the content sequence, and the second largest content point is searched as a new peak candidate, and the confirmation process is repeated until the confirmation condition is met or the number of candidate points is exhausted. This process ensures that the peak output simultaneously satisfies the structural characteristics of "maximum content" and "rate changing from positive to negative".

[0072] After obtaining the peak antioxidant content and peak occurrence time for each algal strain, a normalization mapping is performed to eliminate dimensional differences across algal strains and provide a unified scale input for trajectory fitting and subsequent weighted scoring. The calculation process for peak normalization mapping is as follows: the peak set of all candidate algal strains is summarized, and the minimum and maximum peak values ​​in the set are calculated; for a peak value of a certain algal strain, the relative increment is obtained by subtracting all minimum peak values ​​from the peak value, and then the normalized peak value is obtained by dividing the relative increment by the difference between all maximum peak values ​​and all minimum peak values. The normalized peak value falls within the range of 0 to 1, where 0 and 1 are determined by the normalization output range definition. The calculation process for the peak occurrence time normalization mapping is as follows: The peak occurrence time of each algal strain is converted into a relative time value in minutes. The relative time value is obtained by subtracting the culture start time stamp from the peak timestamp and then converting it to minutes. The set of relative peak times for all algal strains is summarized, and the earliest and latest relative peak times are calculated. For a given algal strain's relative peak time, the earliest relative peak time is subtracted to obtain the relative time increment. This increment is then divided by the difference between the latest and earliest relative peak times to obtain the normalized peak time. To avoid abnormal algal strains stretching the mapping scale and causing a decrease in resolution, the mapping endpoints are determined only using algal strain samples that have passed quality control. The quality control condition is set to ensure that the proportion of effective sampling points for the entire algal strain's lifecycle is not less than 19 / 20. The determination of 19 / 20 is based on the fitting stability requirement: reducing the proportion of effective points in the training data and observing the fluctuation of the fitting residuals shows that when the proportion of effective points is lower than 19 / 20, the fitting fluctuation increases significantly. Boundary truncation is performed on peak or peak time inputs that exceed the endpoint range. The truncation rule is that inputs exceeding the upper endpoint are taken as the upper endpoint value, and inputs below the lower endpoint are taken as the lower endpoint value, so that the normalized output remains stable and avoids extreme values ​​from dominating the subsequent model.

[0073] Normalized peak value and normalized peak value time are entered into the growth trajectory fitting model. Simultaneously, the full-cycle content time series serves as the main input for fitting, used to obtain a fitted trajectory reflecting the dynamic changes throughout the entire cycle and extract morphological feature vectors. The fitting model employs a constrained smoothing fitting framework, with the fitting process divided into three stages: initial smoothing, peak constraint correction, and full-cycle consistency verification. Initial smoothing uses sliding weighted smoothing to suppress high-frequency fluctuations. The smoothing window length is set to 5 sampling points, corresponding to 50 minutes. The 50-minute window is determined based on the coverage requirement of the peak inflection curve's bending radius: when the sampling interval is 10 minutes, the 5-point window covers 20 minutes before and after the peak and includes the peak point, suppressing noise while preserving the peak shape. The weighted smoothing calculation process involves taking 5 points (2 points before and after each center point, plus the center point itself), assigning weights decreasing from the center outwards (5, 4, 3, 2, 1 respectively). The content of each point is multiplied by its weight, summed, and then divided by the total weight to obtain the smoothed content of that center point. This calculation is repeated along the time axis for all center points to obtain the initial smoothed trajectory. The peak constraint correction process involves establishing a constraint window near the time position corresponding to the normalized peak time. The constraint window length is set to 3 sampling points, corresponding to 30 minutes. The 30-minute period is determined based on peak positioning error control: the preceding peak confirmation process ensures that the peak time deviation is stable and does not exceed 20 minutes, and 30 minutes is used to cover this error range. The correction rule is to find the maximum point of the initial smooth trajectory within the constraint window. If this maximum point is inconsistent with the peak candidate point, the trajectory within the constraint window is locally raised or lowered so that the content of the maximum point within the window is consistent with the content level corresponding to the normalized peak. The raising or lowering amplitude is achieved by difference compensation, where the difference is the difference between the target peak content and the current maximum point content. This difference is distributed in decreasing increments according to the time distance from the peak point to maintain the continuity of the peak shape. Full-cycle consistency verification is used to prevent the introduction of non-physiological oscillations during correction. The process involves recalculating the first-order difference rate sequence of the fitted trajectory and checking that the number of rate sign flips matches the training baseline range. The upper limit for the number of sign flips is set to 4, determined by the upper bound statistical analysis of multi-peak interference in the true-content trajectory of the training data. Flips exceeding 4 are often introduced by noise or overfitting. If this is exceeded, the smoothing weight decay is increased and correction is repeated until the number of sign flips does not exceed 4 or the iteration limit is reached. The iteration limit is set to 10, determined by the maximum number of iterations required for convergence in the training data (7) plus 3 redundancies.

[0074] After obtaining the fitted trajectory, morphological feature vectors are extracted to characterize the dynamic change trajectory throughout the entire cycle. Each dimension of the morphological feature vector is generated according to explicit calculation rules, and all dimensions undergo scaling to ensure comparability across algal strains. The starting point of the rising segment is set at the position where the first three consecutive sampling points of the fitted trajectory have positive rates and the rate values ​​are not less than twice the upper limit of the noise rate. The determination of the three consecutive sampling points is based on suppressing single-point noise flips, and the determination of twice the upper limit of the noise rate is based on selecting the lower limit of retaining the true rising segment after statistically analyzing the rate distribution in the non-peak region. The ending point of the rising segment is set at the peak occurrence time. The calculation process for the rising segment duration is the relative time value of the rising segment ending point minus the relative time value of the rising segment starting point. The calculation process for the average rate of the rising segment is the net increment obtained by subtracting the peak content from the content at the starting point of the rising segment, and then dividing the net increment by the rising segment duration to obtain the average growth rate per unit time. The peak platform width is calculated using 0.95 times the peak content as the platform judgment line. The 0.95 is determined based on the fidelity requirements of the peak platform in the training data: platform widths are calculated in the training data using different judgment lines of 0.9, 0.95, and 0.97, and compared with the offline platform annotations; 0.95 corresponds to the highest consistency. The platform width calculation process involves searching backwards from the peak time, taking the moment when the fitted trajectory first falls below the platform judgment line as the platform start point, and then searching backwards from the peak time, taking the moment when it first falls below the platform judgment line as the platform end point. The platform width is the time difference between the platform end point and the platform start point. The descent segment start point is set to the peak occurrence time, and the descent segment end point is set to the position where the absolute value of the rate of the fitted trajectory for the first six consecutive sampling points after the peak does not exceed the upper limit of the noise rate. Six consecutive points correspond to 60 minutes, determined based on the requirement that post-peak stabilization identification needs to cover one noise cycle and suppress short-term rebounds. The descent segment duration is the time difference between the descent segment end point and the descent segment start point. The average rate of the descent segment is calculated by subtracting the peak value from the value at the end of the descent segment to obtain the net change, and then dividing the net change by the duration of the descent segment to obtain the average rate of decrease per unit time. The curvature intensity near the peak is characterized by the local differential change formed by two sampling points before and after the peak. This is calculated by calculating the difference between the rates of the two segments before and after the peak, summing the difference amplitudes to obtain the curvature intensity index, which characterizes the sharpness of the transition. The full-cycle fluctuation energy is characterized by the cumulative deviation between the fitted trajectory and the initial smoothed trajectory. This is calculated by calculating the difference between the fitted value and the initial smoothed value at each sampling point throughout the entire cycle, taking the absolute value, and then summing the absolute values ​​over the entire cycle to obtain the fluctuation energy. A larger fluctuation energy indicates that more corrections are needed to meet the constraints during the fitting process, reflecting the degree of trajectory instability. The stability of the post-peak decline is characterized by the fluctuation range of the content within the stable period after the peak. The calculation process is as follows: after the end of the decline, 12 consecutive sampling points corresponding to 120 minutes are taken as the stable period. The determination of 120 minutes is based on the coverage requirements of short-term stability in production evaluation. Within the stable period, the difference between the maximum content and the minimum content is calculated as the fluctuation range. The smaller the fluctuation range, the higher the stability.After completing the above feature calculations, scale adjustment is performed on each feature dimension. Scale adjustment uses the minimum and maximum values ​​across algal strains to construct a scale. The calculation process is to summarize the minimum and maximum values ​​of the same dimension feature of all algal strains, and map the feature of this dimension of a single algal strain to the range of 0 to 1 by "subtracting the minimum value and then dividing by the difference between the maximum and minimum values", so that each dimension of the morphological feature vector enters the same numerical scale.

[0075] The baseline growth model is used to define a stable reference trajectory, and it is established based on a baseline sample set. The selection criteria for the baseline sample set are: a high comprehensive activity score after offline verification, a valid sampling point ratio of no less than 19 / 20 throughout the entire cycle, and a peak confirmation consistency that meets the aforementioned 0.7 ratio requirement. For each algal strain in the baseline sample set, a morphological feature vector is calculated, and then the baseline morphological vector is obtained by averaging the results dimension by dimension. The average is calculated by summing the results dimension by dimension and then dividing by the number of baseline samples. Simultaneously, the dispersion is calculated dimension by dimension using the standard deviation. The standard deviation is calculated by averaging the squared differences between each sample and the baseline mean, and then taking the square root to obtain the standard deviation, which describes the natural fluctuation range of that dimension in the stable samples. The trajectory deviation is used to quantify the deviation strength of candidate algal strains relative to the baseline growth model. Its calculation process is as follows: the absolute values ​​of the differences between the morphological feature vectors of candidate algal strains and the baseline morphological vectors are calculated dimension by dimension to obtain the dimension-wise deviation; then, the dimension-wise deviations are weighted and summed according to their importance weights to obtain a single deviation value. The weights are determined based on the sensitivity of the screening target and the consistency of the ranking. The determination process involves setting a candidate set of weights in the training data and iterating through the combinations. For each combination, the candidate algae strain ranking is calculated and compared with the offline evaluation ranking. Ranking consistency is represented by the number of overlaps in the top 10 of the two rankings. The larger the overlap, the higher the consistency. When the consistency reaches its maximum, the ranking stability is further compared. Stability is represented by the average difference in ranking fluctuations under different batches of training subsets. The smaller the average difference in ranking, the more stable the system. The weight combination with the highest consistency and the smallest average difference in ranking is selected as the final weight, so that the deviation remains sensitive to the screening target and the output is stable.

[0076] The stability threshold is used to determine whether trajectory deviation reaches a level that affects repeatability. The threshold is determined by the deviation distribution of the baseline samples. The calculation process for determining the threshold is as follows: For each baseline sample, calculate its deviation relative to the baseline morphological vector to obtain a baseline deviation set; calculate the mean and standard deviation of this set, and add 3 times the standard deviation to the mean as the stability threshold. The determination of 3 is based on misclassification control: in the training data, 2, 3, and 4 are used as multiples respectively to count the number of baseline samples judged as unstable. When the multiple is 3, the number of misclassified baseline samples does not exceed 1 / 200 and the threshold is the smallest. Using 3 ensures that the threshold is strict and misclassification is limited. When the trajectory deviation of a candidate algal strain exceeds the stability threshold, the process of deviation interval identification and activity fluctuation variable extraction begins.

[0077] Deviation interval identification uses the local difference threshold after time phase alignment as the criterion. The local difference threshold is determined by the upper limit of the fluctuation of the benchmark sample in the same time phase. Time phase alignment uses the normalized peak time as the alignment anchor point. The alignment process involves shifting the time axis of each algal strain according to the "peak time alignment" method, so that the peaks of different algal strains fall at the same phase position, and then comparing the fitting content at the same phase point. After alignment, the point-by-point difference between the fitting trajectory of the candidate algal strain and the benchmark trajectory is calculated. The difference is calculated by subtracting the benchmark content from the candidate fitting content and taking the absolute value to form a difference sequence. The local difference threshold is set as the mean of the benchmark difference plus 2 standard deviations. The determination of 2 is based on the requirement that local differences need to be more sensitive than the global threshold: in the training data, 1, 2, and 3 are used as multiples respectively, and the detection rate and false detection rate of deviation intervals are statistically analyzed. When the multiple is 2, the detection rate is not less than 49 / 50 and the false detection rate is not more than 1 / 100. The process of determining the deviation interval involves finding consecutive sampling point segments in the difference sequence that exceed the local difference threshold. The minimum length of the consecutive segment is set to 3 sampling points, corresponding to 30 minutes. The 30-minute timeframe is determined based on suppressing single-point noise difference spikes: the maximum duration of noise difference spikes in the training data is 20 minutes, and the 30-minute timeframe covers this length and achieves filtering. The consecutive segments that meet the minimum length requirement are determined as the deviation interval.

[0078] The active fluctuation variable is extracted based on the instantaneous slope within the deviation interval. The instantaneous slope is calculated by taking the change in content between two adjacent sampling points within the deviation interval and converting it to a change per minute over 10 minutes, forming an instantaneous slope sequence. The active fluctuation variable is calculated based on this instantaneous slope sequence, and consists of three sub-quantities. The first sub-quantity is the absolute value of the maximum instantaneous slope, calculated by taking the absolute value of each point in the instantaneous slope sequence and finding the maximum value, representing the strongest fluctuation amplitude. The second sub-quantity is the mean of the absolute values ​​of the instantaneous slopes, calculated by taking the absolute values ​​of each point in the instantaneous slope sequence, summing them, and dividing by the number of points, representing the average fluctuation level. The third sub-quantity is the number of instantaneous slope sign flips, calculated by comparing the signs of two adjacent instantaneous slopes in chronological order, incrementing the count by 1 when the sign changes from positive to negative or vice versa, accumulating the flip count, representing the frequency of fluctuations. The above three sub-quantities are entered into scaling. The scaling scale is determined by the minimum and maximum values ​​of the same sub-quantities in the training stable samples. The calculation process is to summarize the set of stable sample sub-quantities, find the minimum and maximum values, and map the candidate algal strain sub-quantities to the range of 0 to 1 by "subtracting the minimum value and then dividing by the difference between the maximum and minimum values", so that the activity fluctuation variable can be compared between different algal strains.

[0079] The comprehensive activity score is obtained by performing a weighted matrix operation on the activity fluctuation variable, the peak value of antioxidant content, and the peak occurrence time. Numerically, this operation is represented as a process of "multiplying by weights and then summing." The scoring input vector consists of three categories of indicators in sequence: the first category is the normalized peak value, the second category is the normalized peak occurrence time, and the third category is the normalized values ​​of the three sub-quantities of the activity fluctuation variable, forming a total of five input components. The weight matrix consists of five corresponding weight coefficients. The weight coefficients are determined based on offline evaluation consistency and ranking stability: a candidate set of weights is set in the training data, and the weight combinations are traversed. For each combination, a comprehensive score is calculated and a ranking is generated. The ranking consistency is compared with the offline evaluation ranking. The ranking consistency is represented by the number of overlaps in the top 10 and the average difference between all rankings. The larger the number of overlaps and the smaller the average difference between rankings, the higher the consistency. In the combination with the highest consistency, the interpretation direction of the weights is further constrained: the peak weight is positive to increase the intensity contribution, the peak occurrence time weight is negative to reflect that earlier is better, and the activity fluctuation variable weight is negative to reflect the fluctuation penalty, thereby ensuring that the scoring logic is consistent with the selection target. The calculation process of the weighted matrix operation is as follows: each component of the input vector is multiplied by its corresponding weight to obtain five weighted components; the five weighted components are summed to obtain the comprehensive activity score; before outputting the comprehensive activity score, interval sorting is performed. Interval sorting uses the minimum and maximum values ​​of the comprehensive scores in the training data to construct a scale, mapping the score to the range of 0 to 1, so that the score results of different batches can be directly compared. When constructing the antioxidant activity evaluation index system from the comprehensive activity score, the comprehensive activity score is used as the core evaluation quantity, while retaining the peak intensity score, peak time series score, and stability penalty score as interpretable sub-items. The sub-item scores are derived from the cumulative result of the product of the corresponding input component and the weight.

[0080] S7 includes: sorting multiple candidate Spirulina strains in descending order based on their comprehensive activity scores to obtain an initial sequence set; constructing a fluctuation discrete matrix using the activity fluctuation variables of each candidate Spirulina strain in the initial sequence set; extracting eigenvectors reflecting growth stability by performing singular value decomposition on the fluctuation discrete matrix; weighting and correcting the initial sequence set based on the eigenvectors to obtain a corrected activity ranking sequence; if the trajectory deviation of the first candidate Spirulina strain in the activity ranking sequence is lower than a preset stability threshold, the candidate Spirulina strain is identified as a superior algae species; and verifying the extreme values ​​of the peak antioxidant content of the superior algae species in the dynamic change trajectory throughout the entire cycle to determine the superior algae species with the highest antioxidant activity as the target algae species for breeding.

[0081] In this embodiment, the comprehensive activity scores of all candidate Spirulina strains are first read, and an initial sequence set is generated accordingly. The generation process of the initial sequence set is as follows: the comprehensive activity score of each strain is paired with its unique identifier to form a score list; the score list is sorted from largest to smallest according to the comprehensive activity score, and the sorting is completed by a comparison exchange rule, that is, the scores of any two strains are compared one by one, and the strain with the larger score is ranked first, and the strain with the smaller score is ranked last, until all strains meet the condition that there is no reverse order between adjacent scores. In order to avoid uncertainty in the ranking due to tied scores, a tie adjudication order is set. The tie adjudication compares the trajectory deviation, peak occurrence time, and peak content of antioxidant substances in sequence: when the comprehensive activity scores are the same, the trajectory deviation is compared first, and the strain with the smaller trajectory deviation is ranked first; if the trajectory deviation is still the same, the peak occurrence time is compared, and the strain with the earlier peak occurrence time is ranked first; if the peak occurrence time is still the same, the peak content of antioxidant substances is compared, and the strain with the larger peak content is ranked first; after the tie adjudication is completed, the initial sequence set is output, which contains the candidate strain number sequence and the corresponding comprehensive activity score sequence arranged in order.

[0082] Based on the initial sequence set, a fluctuation discrete matrix is ​​constructed using the activity fluctuation variables of each candidate algal strain. The activity fluctuation variable consists of three components: the maximum absolute value of the instantaneous slope in the deviation interval, the mean absolute value of the instantaneous slope, and the number of sign flips of the instantaneous slope. All three components are derived from the results of previous deviation interval identification and instantaneous slope calculation. Scale adjustment is performed before constructing the fluctuation discrete matrix to prevent the matrix direction from being dominated by a single component due to differences in the numerical ranges of the three components. The scaling process is as follows: For all candidate algal strains, the same component is aggregated to form a component set. After removing samples that fail quality control, the minimum and maximum values ​​of that component are calculated. The quality control condition is set to ensure that the number of valid sampling points for the entire algal strain throughout the entire cycle is no less than the total number of sampling points minus two sampling cycles. The determination of two sampling cycles is based on the impact assessment of missing data simulation on the fitting residual. When the missing length reaches three sampling cycles, the fitting residual significantly increases. Therefore, the upper limit of allowed missing data is set to two sampling cycles to ensure that the endpoints of the scaling process are not contaminated by missing data. For a given algal strain, the minimum value of that component is subtracted to obtain the relative increment. Then, the relative increment is divided by the difference between the maximum and minimum values ​​of that component to obtain the scaling result in the range of 0 to 1. 0 and 1 are determined by the scaling output range definition. When the value of a component of an algal strain exceeds the endpoint range, values ​​exceeding the upper endpoint are processed according to the upper endpoint, and values ​​below the lower endpoint are processed according to the lower endpoint, thus keeping the scaling output stable. After scaling the three components, the three components of each algal strain are arranged into row vectors in a predetermined order and stacked row by row according to the initial sequence set to form a fluctuating discrete matrix. The number of rows in the matrix is ​​equal to the number of candidate algal strains, and the number of columns in the matrix is ​​3.

[0083] To eliminate overall bias, the fluctuation discrete matrix is ​​centered column-wise. The centering process is as follows: the average of all rows in the first column is calculated to obtain the column mean. Then, the column mean is subtracted from the values ​​of each row in that column, making the column mean zero. The same calculation is repeated for the second and third columns to obtain the centered fluctuation discrete matrix. The column mean is calculated by summing the values ​​and then dividing by the number of rows, which is the total number of rows in the centered matrix, ensuring that the centering process treats different candidate algal strains equally.

[0084] Singular value decomposition (SVD) is performed on the centered fluctuation discrete matrix to extract eigenvectors reflecting growth stability. The SVD calculation process adopts a chain of "correlation matrix construction, principal direction iterative solution, and vector normalization output" to avoid relying on opaque black-box solutions. The correlation matrix construction process is as follows: a 3x3 correlation matrix is ​​constructed using the centered matrix as input; the value in the first row and first column of the correlation matrix is ​​obtained by multiplying each value in the first column of the centered matrix by itself and summing the results, then dividing the sum by the number of candidate algal strains to obtain the average product; the value in the first row and second column is obtained by multiplying the values ​​in the first and second columns of the centered matrix in the same row, summing the results, and dividing by the number of candidate algal strains; the first row and third column are calculated in the same way as the first and third columns; the second and third rows are calculated column-wise to obtain the complete 3x3 correlation matrix. Subsequently, the principal direction vector is solved from this correlation matrix, and the principal direction vector corresponds to the direction of maximum stability difference. The main direction iterative solution is implemented using power iteration. The calculation process of power iteration is as follows: Initialize a 3-dimensional vector, setting all three components of the initial vector to 1. The initial value of 1 is to ensure that the initial value is not 0 and that the three components are unbiased. Multiply the correlation matrix with the current vector to obtain a new vector. The product operation is implemented by multiplying and adding the components of the current vector in rows and columns. That is, the first row of the correlation matrix is ​​multiplied by the three components of the current vector and summed to obtain the first component of the new vector. The second row is multiplied by the current vector and summed to obtain the second component of the new vector. The third row is multiplied by the current vector and summed to obtain the third component of the new vector. Normalize the new vector. The normalization calculation process is to first calculate the sum of squares of the three components of the new vector, then take the square root of the sum of squares to obtain the vector length. Then, divide each component of the new vector by this length to make the normalized vector length 1. 1 is determined by the definition of unit length. Replace the current vector with the normalized vector to enter the next iteration. Iterative convergence is determined using a component change threshold, set to 0.0001. This threshold is based on the stability verification of the principal direction on the training data: when the threshold is set to 0.001, fluctuations in the principal direction component across different batches of training subsets still cause ranking fluctuations; when the threshold is set to 0.0001, changes in the principal direction component no longer have a significant impact on the ranking results, and the number of iterations does not exceed 30. After each iteration, the absolute value of the difference between the current vector and the previous vector in three components is calculated. If the absolute value of the difference in all three components does not exceed 0.0001, convergence is determined, and the current vector is output as the feature root vector. To ensure termination, the maximum number of iterations is set to 30. This 30 is determined by the maximum number of iterations required for convergence in the training data (22) plus 8 redundancies. The redundancy is used to cover changes in the condition number of the correlation matrix caused by variations in the number of candidate algal strains.

[0085] After obtaining the eigenvectors, the initial sequence set is weighted and corrected based on the eigenvectors to obtain the corrected activity ranking sequence. The weighted correction first calculates the stability projection score for each algal strain. The calculation process is as follows: read the 3D activity fluctuation row vector and eigenvector of the algal strain, multiply each corresponding component and sum them. The sum is the projection score. The multiplication and summation are performed sequentially by the numerical computation module to ensure reproducibility. Subsequently, the projection scores are scaled from 0 to 1. The scale endpoints are determined by the minimum and maximum projection scores of all candidate algal strains. The scaling process uses the rule of "subtracting the minimum value and then dividing by the difference between the maximum and minimum values." Values ​​exceeding the endpoints are treated as endpoints to maintain comparability between different batches of candidate sets. After the stability projection scores are scaled, a corrected score is constructed. The corrected score is obtained by fusing the comprehensive activity score and the stability projection score. The fusion uses a structure of "main score minus penalty term" to reflect the fluctuation penalty. The penalty term is equal to the penalty intensity multiplied by the stability projection score. The penalty intensity is determined through a traversal during the training phase: The candidate penalty intensity set is set to 0.1, 0.2, 0.3, 0.4, and 0.5, with an interval of 0.1. The 0.1 setting is based on the fact that the consistency of the ranking changes with the intensity at this step size, exhibiting a single-peak structure with a clear peak position. A corrected score is calculated and ranked for each candidate intensity. The ranking results are compared with the offline verification ranking, which is a synthesis of offline antioxidant indicators and stability manual verification results. Ranking consistency is selected based on the criterion of "maximum number of overlaps in the top 10 and minimum average difference in overall rankings." The average difference in rankings is obtained by summing the absolute values ​​of the differences between the two rankings for each algal strain and dividing by the number of algal strains. The candidate intensity that meets the criteria is determined as the penalty intensity. Subsequently, the candidate algal strains are re-ranked from largest to smallest according to the corrected score, forming a corrected activity ranking sequence. The decision on tied rankings still follows the order of "smaller trajectory deviation first, earlier peak occurrence time first, larger peak value first," ensuring a single and definite output sequence.

[0086] Stability threshold determination is performed on the corrected activity ranking sequence to identify superior algal species. Threshold determination is achieved by comparing trajectory deviation with a stability threshold. Trajectory deviation is obtained by a weighted sum of the dimension-wise deviations of the morphological feature vectors relative to the baseline growth model; this deviation was generated in previous steps and organized using the same scale. The stability threshold is determined based on the upper bound of the baseline sample deviation distribution. The determination process is as follows: a baseline sample set is selected, with the selection criteria being that the comprehensive activity score ranks among the top 10 in the candidate set and the number of valid sampling points throughout the entire cycle meets the aforementioned upper limit constraint on missing data. For each algal species in the baseline sample set, its trajectory deviation is calculated to form a baseline deviation set. The mean and standard deviation of the baseline deviation set are calculated. The stability threshold is set as the mean plus three times the standard deviation, where the determination of 3 is based on misclassification control: in the training data, 2, 3, and 4 are used as multiples to statistically count the number of baseline samples judged as unstable. When the multiple is 3, the number of misclassified baseline samples does not exceed 1 out of 200 samples, and the threshold is minimized, thus maintaining a strict threshold while constraining misclassification. The threshold determination process is as follows: starting from the first position of the corrected activity ranking sequence, the trajectory deviation of the algal strain is read and compared with the stability threshold. If the trajectory deviation is less than the stability threshold, the algal strain is identified as a superior algal species and the threshold traversal is stopped. If the trajectory deviation is not less than the stability threshold, the comparison is repeated for the next algal strain until an algal strain that meets the conditions is found or the sequence traversal ends.

[0087] After identifying superior algal species, extreme value verification is performed on the peak content of antioxidants in their full-cycle dynamic change trajectory to confirm that the peak value belongs to an effective extreme value in the entire cycle and satisfies the consistency of the transition structure. The input for extreme value verification is the full-cycle content time series sequence of the superior algal species and the first-order difference rate set. The rate set is generated according to the aforementioned process of "converting adjacent content differences into changes per minute and smoothing". Extreme value verification first performs full-cycle maximum value verification. The verification process involves traversing the full-cycle content sequence and recording the maximum content value and its corresponding time. After the traversal, the recorded maximum content value is compared with the previously determined peak content. If they are the same, the maximum value verification is passed. If they are different, the peak time is replaced by the time corresponding to the full-cycle maximum content value and enters the subsequent transition verification. If there are multiple parallel maximum content points in the entire cycle, the transition structure adjudication is performed, and the adjudication is based on the consistency of rate signs and the minimum plateau width. The rate sign consistency calculation process is as follows: A 120-minute validation window is set centered on the candidate peak time. The 120-minute window is determined based on the coverage requirement of the main change segment before and after the peak transition in the training data. The upper limit of the continuous segment length where the peak transition rate changes from positive to negative falls within 110 minutes, so the window is set to 120 minutes to cover the complete segment. The number of positive rate points in the first half of the candidate peak time is counted and divided by the number of rate points in the first half of the window to obtain the increase ratio. The number of negative rate points in the second half of the candidate peak time is counted and divided by the number of rate points in the second half of the window to obtain the decrease ratio. Both the increase ratio threshold and the decrease ratio threshold are set to 0.7. The 0.7 threshold is determined based on the joint optimality of the true peak neighborhood retention rate and the non-peak false entry rate in the training data. When the threshold is 0.7, the ratio of the number of true peak retained to the number of true peak samples is not less than 49:50, and the ratio of the number of non-peak false entries to the number of non-peak windows is not more than 1:100. Parallel peak points that meet the requirements of an increase ratio of not less than 0.7 and a decrease ratio of not less than 0.7 are compared in terms of plateau width. The platform width calculation uses 0.95 times the peak content as the platform judgment line. The determination of 0.95 is based on the consistency verification of platform identification in the training data: the platform width is calculated using 0.9, 0.95, and 0.97 as judgment lines respectively, and compared with the offline platform annotation. 0.95 corresponds to the highest consistency and is not sensitive to noise spikes. The platform width calculation process is as follows: searching backward from the candidate peak time, the moment when the content first falls below the platform judgment line is taken as the platform start time; searching backward from the candidate peak time, the moment when the content first falls below the platform judgment line is taken as the platform end time. The platform width is the time difference between the platform end time and the platform start time. The candidate peak point with the smaller platform width is determined as the final peak time. If the platform widths are still the same, the earlier peak time is selected as the final peak time to ensure the uniqueness of the output. After completing the extreme value verification, the verified superior algal species are output as the target algal species for breeding. The output content includes the algal strain number, the corrected activity ranking, the comprehensive activity score, the corrected score, the trajectory deviation, the peak content, and the peak occurrence time.

[0088] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for breeding and screening Spirulina strains with high antioxidant activity, characterized in that, include: S1. By culturing different candidate Spirulina strains in multiple culture units, and installing a multispectral sensor array in each culture unit to collect real-time spectral data and environmental parameter data during the growth process of Spirulina, the spectral feature vectors related to antioxidant activity are obtained. S2. Based on the obtained spectral feature vectors related to antioxidant activity, input them into a pre-trained support vector machine model for classification, determine the current growth stage of the algal strain, and then integrate environmental parameter data to calculate the dynamic change trend of the growth stage after obtaining the classification results. S3. If the dynamic change trend of the determined growth stage shows that the biomass growth is slowing down, then the random forest algorithm is used to perform regression prediction on the spectral feature vector to obtain the predicted antioxidant concentration value and obtain continuous time series data of the concentration value. S4. The continuous time series data of the obtained concentration values ​​are analyzed using the sliding window method. The concentration gradient change within each window is calculated, and it is determined whether the gradient change exceeds the preset threshold to obtain a set of candidate points for potential peak times. S5. Cluster the obtained potential peak moment candidate point set, use the K-means algorithm to group similar points to obtain the clustered peak moment group, select the point with the highest concentration value from the group, and determine the peak content of antioxidant substances of the corresponding algal strain in its culture cycle.

2. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 1, characterized in that: S1 includes: Real-time spectral data and environmental data during the growth process of different candidate Spirulina strains were acquired using a multispectral sensor array. Wavelet denoising algorithm is used to process real-time spectral data and environmental data to obtain smooth spectral sequences and standard environmental sequences; Frequency domain spectral characteristics were obtained by processing smooth spectral sequences using Fourier transform infrared spectroscopy analysis. Construct a multidimensional feature tensor based on frequency domain spectral characteristics and standard environmental sequences; If the multidimensional feature tensor meets the preset activity threshold condition, then the spectral feature vector related to antioxidant activity is extracted from the multidimensional feature tensor.

3. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 1, characterized in that: S2 includes: The standard feature matrix is ​​obtained by normalizing the spectral feature vectors related to antioxidant activity; The standard feature matrix is ​​input into a pre-trained support vector machine model for classification and identification to obtain the current growth stage of the algae strain. Based on the current growth stage of the algae, historical environmental parameter data within the corresponding time window are retrieved to obtain an environmental time series set. By performing multi-dimensional spatial mapping between the environmental time series set and the current algal growth stage, the stage-environment-related coordinates are obtained; If the stage environment-related coordinates deviate from the preset steady-state center point, the dynamic change trend of the growth stage is obtained by calculating the projected displacement of the stage environment-related coordinates on the time axis.

4. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 1, characterized in that: S3 includes: Acquire algal density collected by sensors and calculate biomass growth rate; If the biomass growth rate is lower than the preset dynamic threshold range, the original spectral data of the corresponding time point will be retrieved. Principal component analysis is performed on the original spectral data to extract dimensionality-reduced feature vectors. The dimensionality-reduced feature vectors are input into a pre-trained random forest regression model to obtain the predicted concentration values; Concentration time series are constructed based on the distribution pattern of predicted concentration values ​​on the time axis; The concentration-time series was used to fit the accumulation trajectory of substances and determine the concentration values ​​of antioxidants.

5. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 1, characterized in that: S4 includes: Continuous time series data of antioxidant concentration values ​​are obtained, and the continuous time series data is slid segmented according to a preset window step size to obtain multiple overlapping local concentration intervals. Linear fitting is performed on sampling points within a local concentration range to calculate the concentration gradient change; If the sign of the concentration gradient change is reversed and the absolute value exceeds a preset threshold, the extreme point of the slope within the local concentration range is extracted as a candidate point for the potential peak moment. Clustering is performed based on the distribution density of candidate points on the time axis, and outliers caused by time fluctuations are removed to determine the final candidate point set. By using the sampling frequency corresponding to the candidate point set to perform local data smoothing, a candidate point set that can reflect the potential peak moment of antioxidant accumulation is obtained.

6. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 1, characterized in that: S5 includes: Obtain the time coordinates in the candidate point set and construct a one-dimensional feature vector space, then calculate the Euclidean distance between different time coordinates in the one-dimensional feature vector space. The candidate point set is iteratively divided using the K-means algorithm based on Euclidean distance to obtain the peak time group where the cluster center no longer shifts and has cluster label. Retrieve the concentration values ​​within the peak time group and extract the sampling time with the highest concentration value as the local extremum; Local extreme values ​​are mapped to the time axis of the algal cell cycle to identify growth stages. By comparing the local extreme values ​​of different growth stages, the peak value of antioxidant content of the corresponding algal cell during its culture cycle is determined.

7. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 1, characterized in that, It also includes S6, which constructs an antioxidant activity evaluation index system based on the peak content of antioxidant substances of each candidate Spirulina strain, the time of peak occurrence, and their dynamic changes throughout the entire cycle, specifically including: The antioxidant content time series of each candidate Spirulina strain was obtained during the complete culture cycle. The content change rate set was obtained by performing a first-order difference operation on the antioxidant content time series. The antioxidant content peak and its corresponding peak occurrence time were determined by using the content change rate set. Normalization mapping is performed based on the peak content of antioxidants and the time of peak occurrence. The mapped values ​​are then imported into a preset growth trajectory fitting model. By extracting features from the growth trajectory fitting model, a morphological feature vector reflecting the dynamic change trajectory throughout the entire cycle is obtained.

8. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 7, characterized in that: S6 further includes: The trajectory deviation of each candidate Spirulina strain relative to the baseline growth model is calculated using morphological feature vectors. If the trajectory deviation exceeds the preset stability threshold, the instantaneous slope of the candidate Spirulina strain in the deviation interval is extracted as the activity fluctuation variable. By performing a weighted matrix operation on the activity fluctuation variable, the peak content of antioxidant substances, and the time of peak occurrence, a comprehensive activity score for each candidate Spirulina strain was obtained, and an antioxidant activity evaluation index system was constructed using the comprehensive activity score.

9. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 7, characterized in that, This also includes S7, which involves ranking and screening multiple candidate Spirulina strains based on an antioxidant activity evaluation index system to determine the superior algal strain with the highest antioxidant activity as the target algal strain for breeding. Specifically, this includes: The initial sequence set was obtained by sorting multiple candidate Spirulina strains in descending order based on their comprehensive activity scores. A fluctuation discrete matrix is ​​constructed using the activity fluctuation variables of each candidate Spirulina strain in the initial sequence set; By performing singular value decomposition on the fluctuation discrete matrix, eigenvectors reflecting growth stability are extracted.

10. The method for breeding and screening Spirulina strains with high antioxidant activity according to claim 9, characterized in that: The S7 also includes: The initial sequence of positions is weighted and corrected based on the eigenvalue vectors to obtain the corrected activity ranking sequence. If the trajectory deviation of the first candidate Spirulina strain in the activity ranking sequence is lower than the preset stability threshold, then the candidate Spirulina strain is identified as a superior algal species. By verifying the extreme values ​​of the peak antioxidant content in the dynamic change trajectory of superior algal species throughout the entire cycle, the superior algal species with the highest antioxidant activity were identified as the target algal species for breeding.

Images