Chemical fiber raw material impurity spectrum intelligent analysis system

Abstract

The invention relates to the technical field of big data service, in particular to a chemical fiber raw material impurity spectrum intelligent analysis system which comprises a spectrum feature deconstruction module, an impurity feature modeling module, a periodic drift analysis module and a component quantification output module. According to the method, noise is suppressed through moving window scanning in combination with Savitzky-Golay filtering, peak shape features are reserved, the signal-to-noise ratio is increased, feature distortion caused by mean filtering is avoided, segments are divided through extreme value density, weak impurity recognition is enhanced in combination with a dynamic threshold value, and a multi-dimensional fusion model is established through peak width change rate and baseline offset convolution. Linear dimension reduction limitation is broken through to improve discrimination precision, dynamic time warping is used for aligning peak height difference and symmetry degree time sequence, drift and interference influences are eliminated, periodic error accumulation is solved, Kalman filtering is used for carrying out recursive optimization on compensated concentration, model lag is reduced, and a closed-loop optimization system is constructed through multi-stage feature decoupling and dynamic compensation. The sensitivity is improved; and the misjudgment rate is reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of big data services, and particularly to an intelligent spectral analysis system for impurities in chemical fiber raw materials. Background Art

[0002] The technical field of big data services includes a full-process technical system for data acquisition, storage, processing, analysis, and visualization. Its core lies in realizing the value mining of massive heterogeneous data through a distributed computing framework and machine learning algorithms. This field involves spectral signal feature extraction technology, high-dimensional data dimensionality reduction methods, real-time stream data processing architectures, and a quality evaluation system based on statistical models, focusing on solving the problems of standardized processing and knowledge discovery of multi-source heterogeneous data in industrial detection scenarios.

[0003] Among them, the intelligent spectral analysis system for impurities in chemical fiber raw materials refers to a data compression method based on near-infrared spectral characteristic wavelength selection technology combined with principal component analysis. By establishing a mapping relationship model between the spectral absorption peaks of raw materials and impurity content, using a sliding window mean filter to eliminate spectral noise, and using the partial least squares regression algorithm to construct a quantitative analysis model, relying on the historical data of the spectral database to achieve impurity type pattern matching, and finally forming a technical solution for raw material quality grading determination rules.

[0004] The sliding window mean filter uses a fixed window size and neighborhood mean calculation. The residual high-frequency noise causes weak signals to be masked in low signal-to-noise ratio scenarios, and the blurring of the edges of characteristic peaks affects the positioning accuracy. The linear dimensionality reduction characteristic of principal component analysis cannot capture local non-linear features in high-dimensional data, and some impurity components with low correlation are weakened during the compression process. The partial least squares regression model does not consider the baseline drift caused by temperature and humidity changes, and environmental factors during long-term operation cause the prediction deviation to be amplified. The pattern matching driven by the historical data of the spectral database lacks a dynamic update mechanism, and the threshold needs to be frequently calibrated manually when there are minor changes in raw material batches or processes. The real-time stream processing architecture does not design a time series drift compensation module, and the errors in continuous detection scenarios accumulate with time, affecting the long-term reliability of the quality grading results. For example, the spectral fluctuations of sampling points in the same batch respond sluggishly, and equipment vibration noise is easily misjudged as a real impurity signal. Summary of the Invention

[0005] The purpose of the present invention is to solve the disadvantages existing in the prior art and propose an intelligent spectral analysis system for impurities in chemical fiber raw materials.

[0006] To achieve the above purpose, the present invention adopts the following technical solutions: The intelligent spectral analysis system for impurities in chemical fiber raw materials includes: A spectral feature deconstruction module, which is used to scan the full-band spectrum through a moving window, extract peak width and peak height parameters by using the Savitzky-Golay filtering algorithm, divide characteristic sections based on extreme value density, generate an impurity feature marking matrix, and transfer the impurity feature marking matrix to the impurity feature modeling module; An impurity feature modeling module, which is used to receive the impurity feature marking matrix, perform a convolution operation on the characteristic section for the peak width change rate and the baseline offset, generate a feature-component weight vector, and transfer the feature-component weight vector to the periodic drift analysis module; A periodic drift analysis module, which is used to call the feature-component weight vector, perform dynamic time warping matching on the peak height difference and symmetry of continuous peak clusters, construct a warping matching matrix through path bending cost calculation and morphological difference accumulation, generate a component-related drift coefficient, and transfer the component-related drift coefficient to the component quantification output module; A component quantification output module, which is used to perform drift compensation on the initial content based on the component-related drift coefficient, update the concentration parameter by using the Kalman filtering algorithm, and output an impurity quantification analysis report.

[0007] As a further solution of the present invention, the impurity feature marking matrix specifically includes peak width parameters, peak height parameters, and extreme value density partitions. The feature-component weight vector includes peak width change rate weights, baseline offset weights, and convolution operation coefficients. The component-related drift coefficient specifically refers to peak spacing variation coefficients, fusion degree attenuation parameters, and dynamic warping matching degrees. The impurity quantification analysis report includes xylene derivative concentration, ester by-product concentration, and Kalman gain parameters.

[0008] As a further solution of the present invention, the combined feature extraction method of the Savitzky-Golay filtering algorithm and extreme value density division enhances the peak shape recognition accuracy through second derivative transformation, and improves the efficiency of characteristic section division by combining extreme point density statistics; The Kalman filtering algorithm uses a state equation and an observation equation , where A is the state transition matrix taking values in the range of 0.95 - 1.05, H is the observation matrix taking the identity matrix, and the process noise covariance Q and the observation noise covariance R are determined through a spectrometer calibration experiment.

[0009] As a further solution of the present invention, the spectral feature deconstruction module includes: A window scanning sub-module acquires full-band spectral data, uses a fixed step size to move a rectangular window to cover the spectral band, records the spectral intensity sequence within each window, and linearly interpolates and splices the overlapping regions of adjacent windows to generate a window spectral sequence covering the full band; Based on the window spectral sequence, the peak shape parameter extraction sub-module performs second derivative transformation on each window spectrum using the Savitzky-Golay filtering algorithm, calculates the spectral curvature change rate within the window, locates the coordinates of local extreme points, extracts the left half-peak width, right half-peak width, and peak height value on the left side of multiple extreme points, and generates a peak shape parameter set; The characteristic section division sub-module calls the peak shape parameter set, counts the number of extreme points within the unit wavelength interval as the density reference, divides the continuous band according to the density threshold optimized by the gradient descent method, merges the boundaries of adjacent high-density intervals, calculates the mean peak width and peak height coefficient of variation within multiple independent sections, and generates an impurity characteristic marking matrix.

[0010] As a further solution of the present invention, the half-peak width measurement uses a nanometer wavelength unit and is subjected to dimension normalization processing by a standard wavelength calibrator.

[0011] As a further solution of the present invention, the impurity characteristic modeling module includes: The peak width dynamic analysis sub-module extracts the characteristic sections marked as impurities and their corresponding sampling points based on the impurity characteristic marking matrix, establishes a measurement mechanism for the distance between adjacent sampling points, calculates the ratio of the peak width change amount in the current time window to the change amount in the previous time window, and performs median filtering on the ratios of five consecutive time windows using the sliding window method to generate a dynamic gradient matrix; The baseline offset quantization sub-module calls the baseline trajectory data of the characteristic section, constructs a moving average model of the reference baseline value, calculates the absolute deviation amount between the baseline values of multiple sampling points and the corresponding reference values, and performs arithmetic mean operation on the deviation amounts of three consecutive sampling points using a three-point sliding window to establish an offset intensity tensor; The convolution weight fusion sub-module splits the dynamic gradient matrix by row vectors, performs dot product operations with the column vectors of the offset intensity tensor, and uses the formula: ; Performs normalization processing on the dot product result, accumulates multi-vector components along the time dimension, and generates a feature-component weight vector; Among them, represents the weight vector of the i-th characteristic section and the j-th component, is the dimensionless peak width change rate, representing the ratio of the peak width change amount at the k-th sampling point in the i-th characteristic section to the change amount in the previous time window. i represents the characteristic section serial number, and k represents the sampling point serial number. is the baseline offset intensity, representing the absolute deviation amount between the baseline value at the k-th sampling point in the j-th characteristic section and the moving average value. j represents the characteristic section serial number, n represents the total number of sampling points, max(ΔP) represents the maximum peak width gradient of the current characteristic section, and max(ΔB) represents the maximum baseline offset intensity of the current characteristic section.

[0012] As a further solution of the present invention, the periodic drift analysis module includes: The peak cluster alignment sub-module calls the feature-component weight vector, locates the horizontal and vertical coordinates of the peak tops of continuous peak clusters, calculates the morphological difference between adjacent peak clusters using the peak height difference and the symmetry axis offset, screens effective peak cluster pairs based on the peak height difference threshold determined by spectral standard sample testing, constructs a peak spacing and symmetry difference sequence in time series, and generates a peak cluster alignment sequence; The dynamic programming sub-module, based on the peak cluster alignment sequence, uses the dynamic time warping algorithm to match the peak cluster spacing sequences in the differential time dimension, calculates the cumulative path cost of the peak height difference sequence and the symmetry sequence, adjusts the cumulative amount of morphological difference through the bending path weight coefficient, and generates a warping matching matrix; The drift coefficient generation sub-module calls the warping matching matrix, extracts the product factor of the path bending cost and the cumulative amount of morphological difference, performs range normalization processing in combination with the feature-component weight vector, calculates the periodic drift influence weights corresponding to multiple components, and generates component-related drift coefficients.

[0013] As a further solution of the present invention, the component quantification output module includes: The drift compensation sub-module calls the component-related drift coefficients, performs a Hadamard product operation on the initial content value and the drift coefficients, compensates for the concentration offset between adjacent time windows using the piecewise linear interpolation method, and adjusts the interpolation weight according to the compensation intensity coefficient of the floating-point number in the range of 0.5-1.2 to generate a drift compensation sequence; The parameter update sub-module, based on the drift compensation sequence, constructs a Kalman filter state equation and an observation equation, calculates the prior estimation covariance matrix through the prediction step, and calculates the Kalman gain coefficient by fusing the observation noise covariance matrix in the update step, and iteratively outputs an updated concentration parameter set; The quantification analysis sub-module calls the updated concentration parameter set, calculates the kurtosis coefficient and skewness coefficient of the concentration values of multiple components, constructs a concentration probability density function in combination with the impurity determination threshold determined by three parallel experiments, and outputs an impurity quantification analysis report.

[0014] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In the present invention, spectral noise suppression is optimized by combining moving window scanning with the Savitzky-Golay filtering algorithm, which improves the signal-to-noise ratio while preserving the peak shape characteristics and reduces the feature smoothing distortion of traditional mean filtering. The extreme density division feature section enhances the ability to identify weak impurity signals through a dynamic threshold adjustment mechanism, achieving precise positioning under complex spectral backgrounds. The convolution operation of the peak width change rate and the baseline offset establishes a multi-dimensional feature fusion model, breaking through the limitation of linear dimensionality reduction and improving the discrimination accuracy of impurity types. The dynamic time warping matching technology performs temporal alignment on the peak height difference and symmetry, eliminating the influence of instrument drift and environmental interference on continuous detection and solving the problem of periodic error accumulation. The Kalman filtering algorithm recursively optimizes the concentration parameters after drift compensation, reducing the model hysteresis and making the quantization result fit the actual production fluctuations. The multi-stage feature decoupling and dynamic compensation form a closed-loop optimization system, achieving an improvement in impurity detection sensitivity and a decrease in the false positive rate. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 is the system flow chart of the present invention; Figure 2 is the flow chart of the spectral feature deconstruction module of the present invention; Figure 3 is the flow chart of the impurity feature modeling module of the present invention; Figure 4 is the flow chart of the periodic drift analysis module of the present invention; Figure 5 is the flow chart of the component quantization output module of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0016] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0017] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the accompanying drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus should not be construed as limiting the present invention. In addition, in the description of the present invention, the meaning of "a plurality of" is two or more unless otherwise specifically defined.

[0018] Embodiment 1:

[0019] Please refer to Figure 1, the intelligent spectral analysis system for chemical fiber raw material impurities includes: A spectral feature deconstruction module, which is used to scan the full-band spectrum through a moving window, extract peak width and peak height parameters using the Savitzky-Golay filtering algorithm, divide the characteristic sections based on the extreme value density, generate an impurity characteristic marking matrix, and transfer the impurity characteristic marking matrix to the impurity characteristic modeling module; An impurity characteristic modeling module, which is used to receive the impurity characteristic marking matrix, perform a convolution operation on the characteristic section for the peak width change rate and the baseline offset amount, generate a feature-component weight vector, and transfer the feature-component weight vector to the periodic drift analysis module; A periodic drift analysis module, which is used to call the feature-component weight vector, perform dynamic time warping matching on the continuous peak clusters for the peak height difference and symmetry, construct a warping matching matrix through path bending cost calculation and morphological difference cumulative amount, generate a component-related drift coefficient, and transfer the component-related drift coefficient to the component quantification output module; A component quantification output module, which is used to perform drift compensation on the initial content based on the component-related drift coefficient, update the concentration parameter using the Kalman filtering algorithm, and output an impurity quantification analysis report.

[0020] The impurity characteristic marking matrix specifically includes peak width parameters, peak height parameters, and extreme value density partitioning. The feature-component weight vector includes peak width change rate weight, baseline offset amount weight, and convolution operation coefficient. The component-related drift coefficient specifically refers to peak spacing variation coefficient, fusion degree attenuation parameter, and dynamic warping matching degree. The impurity quantification analysis report includes xylene derivative concentration, ester by-product concentration, and Kalman gain parameter.

[0021] The combined feature extraction method of Savitzky-Golay filtering algorithm and extreme value density partitioning enhances the peak shape recognition accuracy through second derivative transformation and improves the feature section division efficiency by combining extreme point density statistics; The Kalman filtering algorithm uses the state equation and the observation equation , where A is the state transition matrix taking values in the range of 0.95 - 1.05, H is the observation matrix taking the identity matrix, and the process noise covariance Q and the observation noise covariance R are determined through the spectrometer calibration experiment.

[0022] Please refer to Figure 2 , the spectral feature deconstruction module includes: The window scanning sub-module obtains the full-band spectral data, uses a fixed step size to move a rectangular window to cover the spectral band, records the spectral intensity sequence within each window, and linearly interpolates and stitches the overlapping regions of adjacent windows to generate a window spectral sequence covering the full band; Full-band spectral data, specifically in an application for on-line monitoring of industrial product purity, are the spectral information continuously collected by a detection device in the wavelength range of 200 nm to 1100 nm. This information is stored in the form of a data matrix. Each row of the matrix corresponds to one sampling, the columns represent different wavelength points, and the element value is the light intensity at that wavelength point. For example, at 200.0 nm in a certain sampling, the light intensity is recorded as 0.051 units, at 200.1 nm it is 0.053 units, and until 1100.0 nm it is 0.012 units, forming a sequence containing data points.

[0023] Subsequently, a rectangular window is used to scan the spectral data. The width of the window is set to 50 nm. This width is determined based on the analysis of the spectral data of historical batch samples, by statistically calculating the average full width at half maximum (FWHM) of the known impurity characteristic peaks, and ensuring that the window can completely cover the vast majority of independent impurity peaks and have the ability to resolve adjacent peaks. Specifically, through spectral analysis of 150 samples containing the target impurities, the FWHM of the impurity peaks is mainly distributed between 15 nm and 25 nm, with an average value of 20 nm and a standard deviation of 4 nm. Considering the baseline fluctuation and noise effects, and to ensure the complete capture of features, the window width is set to twice the margin of the average FWHM plus three times the standard deviation, that is nm. After rounding and adding a safety margin, it is set to 50 nm to cope with wider characteristic peaks. The fixed step size of the spectral scan is set to 5 nm. The selection of this step size comprehensively considers the highest hardware resolution of the spectrometer (such as 0.1 nm), the minimum expected size of the detected features, and the data processing efficiency, and is selected as 50 times the hardware resolution, that is nm. This setting can achieve a balance between effectively capturing spectral details and controlling data redundancy.

[0024] The window starts moving from the starting wavelength of 200 nm of the spectrum. The first window covers the wavelength band of 200 nm to 250 nm. The light intensity values corresponding to all wavelength points within this window are completely recorded, constituting the spectral intensity sequence of the first window. For example, its data points are [0.051, 0.053, …, 0.155], with a total of intensity values. After completing the recording of the first window, the window moves 5 nm in the direction of increasing wavelength. The wavelength band covered by the second window is 205 nm to 255 nm, and the spectral intensity sequence in this interval is also recorded. This process is repeated until the window moves to the end of the spectral data. The last window will cover, for example, the wavelength band of 1050 nm to 1100 nm.

[0025] For any two adjacent windows, such as the first window (200 - 250 nm) and the second window (205 - 255 nm), there is an overlapping region of 45 nm (205 nm to 250 nm) between them. Within this overlapping region, the spectral intensities recorded by the two windows need to be smoothly transitioned, and specifically, linear interpolation is used for splicing. For any wavelength point within the overlapping region , its updated spectral intensity value is calculated by the following method: , where is the intensity at this wavelength point in the previous window, is the intensity at this wavelength point in the next window. The weighting factor linearly decreases from 1 at the starting point (205 nm) of the overlapping region to 0 at the ending point (250 nm) of the overlapping region, while the weighting factor linearly increases from 0 at the starting point to 1 at the ending point. Taking the wavelength of 210 nm as an example, this point is the 101st data point in the first window (counting from 200 nm at 0.1 nm intervals) and also the 51st data point in the second window (counting from 205 nm), and its position ratio in the overlapping region is . Therefore, , . If while , then . By performing such linear interpolation splicing on the overlapping regions of all adjacent windows, a window spectral sequence that covers the entire 200 nm to 1100 nm band and has a smooth transition at the window connections is finally integrated.

[0026] Based on the window spectral sequence, the peak shape parameter extraction sub-module uses the Savitzky-Golay filtering algorithm to perform a second derivative transformation on each window spectrum, calculates the spectral curvature change rate within the window, locates the coordinates of the local extreme points, extracts the left half-width, right half-width, and peak height values of multiple extreme points, and generates a set of peak shape parameters; The window spectral sequence generated by the previous sub-module is a series of spliced window spectral segments with a width of 50 nm that cover the entire band. For each window spectral sequence among them, such as the aforementioned window segment from 200 nm to 250 nm, the sequence of spectral intensity data it contains is [0.051, 0.053,..., 0.155]. Perform a second derivative transformation operation on the data points in this sequence. This transformation does not directly name a specific algorithm but is described as: for each data point ( is the index of the data point within this window), its second derivative value is obtained by examining its five neighboring data points (i.e., ), the spectral intensity values, use the least squares method to fit a quadratic polynomial to these five points, and then analytically calculate the second derivative of the polynomial at the point obtained, and the specific calculation is , where is the sampling interval of the spectral data, that is, 0.1 nanometer. Through this calculation, the spectral curvature change rate values at each wavelength point within the window are obtained. A large negative (such as -0.5 unit / nm²) indicates that the original spectrum is concave downward and has a significant curvature at this point, while a positive value (such as 0.2 unit / nm²) indicates an upward convexity.

[0027] In the calculated second derivative spectral sequence, local extreme points are located. The location process is as follows: Identify the zero points where the second derivative value transitions from positive to negative (corresponding to the left shoulder of the absorption peak or the bulge caused by noise in the original spectrum), and the zero points where it transitions from negative to positive (corresponding to the bottom or peak position of the absorption peak in the original spectrum). At the same time, pay attention to the minimum points of the second derivative itself (usually with large negative values), which precisely correspond to the peak positions in the original spectrum. To exclude noise interference, a threshold for the absolute value of the second derivative is set. The method for determining this threshold is: Collect the spectral data of 10 groups of pure solvents (or blank substrates), calculate their respective second derivative spectra, and then statistically calculate the standard deviation of these 10 groups of second derivative spectra at all wavelength points, and take the average standard deviation. For example, if 0.015 unit / nm² is obtained, the noise threshold is set to three times of it, that is unit / nm². Any point where the sign of the second derivative value changes and the absolute values before and after the change are both greater than 0.045 unit / nm², or the minimum point of the second derivative with an absolute value greater than 0.045 unit / nm² itself, the wavelength coordinate within the window is recorded. For example, at 225.3 nm, the second derivative value is -0.152 unit / nm², which is a significant minimum point, so 225.3 nm is identified as a peak candidate position.

[0028] For each extreme point located as a peak, such as the peak at 225.3 nm, its original spectral intensity value is 0.120 unit. First, determine the local baseline of this peak. Draw a straight line by connecting the two nearest valley points (or preset baseline anchor points) on the left and right sides of this peak. The intensity value of this straight line at 225.3 nm is the local baseline height, which is set to 0.030 unit. Then the peak height value of this peak is the original intensity minus the baseline height, that is unit. Then, search to the left (short wavelength direction) of this peak until a wavelength point is found whose original spectral intensity value is equal to half of the peak height plus the baseline height, that is unit, and record the wavelength of this point as , set to 224.8 nm. Then the left half-peak width is the difference between the peak wavelength and this wavelength, i.e., nm. Similarly, looking to the right (long-wavelength direction) of the peak for the point with an intensity of 0.075 units, record its wavelength as , set to 225.9 nm. Then the right half-peak width is nm. Collect the wavelength coordinates of all the identified extreme points within this window, the calculated left half-peak width, right half-peak width, and peak height values to form a parameter list, such as [{wavelength: 225.3 nm, left half-peak width: 0.5 nm, right half-peak width: 0.6 nm, peak height: 0.090 units}, {wavelength: 240.1 nm, left half-peak width: 0.4 nm, right half-peak width: 0.5 nm, peak height: 0.110 units}]. Repeat this process for all window spectral sequences and finally converge into a complete set of peak shape parameters.

[0029] The characteristic section division sub-module calls the set of peak shape parameters, counts the number of extreme points within the unit wavelength interval as the density reference, divides the continuous band according to the density threshold optimized by the gradient descent method, merges the boundaries of adjacent high-density intervals, calculates the average peak width and the coefficient of variation of peak height within multiple independent sections, and generates an impurity characteristic marking matrix; The set of peak shape parameters is the result obtained after the previous sub-module processes the full-band spectrum, containing the detailed parameters of all the identified extreme points. Its specific form is a structured list, for example: [{wavelength: 225.3 nm, left half-peak width: 0.5 nm, right half-peak width: 0.6 nm, peak height: 0.090 units}, {wavelength: 228.1 nm, left half-peak width: 0.4 nm, right half-peak width: 0.5 nm, peak height: 0.110 units}, {wavelength: 230.5 nm, left half-peak width: 0.3 nm, right half-peak width: 0.3 nm, peak height: 0.065 units},..., {wavelength: 1050.0 nm, left half-peak width: 0.8 nm, right half-peak width: 0.7 nm, peak height: 0.040 units}].

[0030] First, it is necessary to count the number of extreme points within a unit wavelength interval, which serves as the benchmark for spectral feature density. The width of the "unit wavelength interval" here is set to 10 nanometers. The basis for this width selection is as follows: Through spectral analysis of a mixed standard sample containing multiple known impurities, it is observed that the characteristic peak clusters of different impurities usually show an aggregation phenomenon within the range of 10 - 30 nanometers. Selecting 10 nanometers as the statistical window can achieve a balance between effectively distinguishing the dense peak region and the sparse peak region, while avoiding the loss of details caused by excessive smoothing. The specific operation is to divide the entire spectral measurement range (200 nanometers to 1100 nanometers) into units of 10 nanometers, and calculate the number of extreme points (from the above-mentioned peak shape parameter set) falling within each 10 - nanometer interval. For example, within the interval from 220 nanometers to 230 nanometers, if the peak shape parameter set contains two extreme points at 225.3 nanometers and 228.1 nanometers, then the number of extreme points in this interval is 2. Within the interval from 230 nanometers to 240 nanometers, if it only contains one extreme point at 230.5 nanometers, then the number is 1. Traversing the entire spectral range in this way, a density sequence is obtained, such as [..., 0, 1, 2, 1, 0, 3, 2, 1,...], where each value corresponds to the extreme point count within a 10 - nanometer - wide band.

[0031] Next, based on an optimized density threshold, continuous bands are segmented to identify high - density feature regions. The optimization process of this density threshold is as follows: First, calculate the average value of the number of extreme points in all 10 - nanometer intervals within the entire spectral range, denoted as , and calculate its standard deviation . For example, through the analysis of the spectra of 100 samples under different process conditions, the average number of extreme points in all 10 - nanometer intervals is 0.8, and the standard deviation is 0.5. The initial density threshold can be set to extreme points / 10 nanometers. Then, use this threshold to perform a preliminary segmentation of the spectrum, and compare the segmented high - density region with a set of known impurity feature regions manually labeled by experienced analysts, and calculate their overlap degree (such as using the Dice coefficient). If the overlap degree is lower than the preset target (such as 0.85), then adjust the threshold: If there are many known impurity regions not recognized (high false negative rate), then appropriately lower the threshold, for example, adjust it to ; if a large number of background noise regions are misidentified as high - density regions (high false positive rate), then appropriately increase the threshold, for example, adjust it to Repeat this iterative adjustment and verification process until the coincidence degree between the segmentation result and the expert annotation reaches an acceptable level, for example, the Dice coefficient is above 0.90. Through this process, the final density threshold is determined to be 1.1 extreme points / 10 nm. Any continuous band whose extreme point density in each 10-nm sub-interval inside is greater than or equal to 1.1 is initially marked as a high-density interval. For example, if the sequence is [..., 0.5 (low), 1.0 (low), 2.0 (high), 1.5 (high), 0.8 (low), 3.0 (high), 2.5 (high), 1.2 (high), 0.6 (low),...] (the density has been normalized to per 10 nm), then the bands corresponding to 2.0 and 1.5 form a high-density section, and the bands corresponding to 3.0, 2.5, and 1.2 form another high-density section.

[0032] Subsequently, boundary merging is performed on these initially marked high-density intervals. If two adjacent high-density intervals are separated by one or more continuous low-density intervals, but the total length of these low-density intervals is less than a preset merging length threshold, then these two high-density intervals will be merged into a longer characteristic section. The merging length threshold is set to 20 nm, based on the fact that the width of typical impurity peak clusters usually does not exceed this range. If the interval is too small, it may belong to different characteristic peaks of the same component. For example, if 220 - 240 nm is a high-density section and 250 - 280 nm is another high-density section, and they are separated by a 10-nm low-density section from 240 - 250 nm, since 10 nm is less than the 20-nm merging length threshold, these two high-density sections will be merged into a single characteristic section from 220 - 280 nm.

[0033] For each such formed independent characteristic section, such as section one (from 220 nm to 280 nm), it is necessary to calculate the mean peak width and the coefficient of variation of peak height for all the original extreme points inside this section. The peak width uses the full width at half maximum of each extreme point, that is, the sum of its left half-width and right half-width. Suppose the extreme points and their full widths included in section one are: (225.3 nm, 1.1 nm), (228.1 nm, 0.9 nm), (240.1 nm, 1.2 nm), (255.5 nm, 1.0 nm), (270.3 nm, 1.3 nm). Then the mean peak width of this section is nm. At the same time, calculate the coefficient of variation of the peak heights (which have been calculated in the peak shape parameter extraction sub-module as the net height after baseline removal) of these extreme points. The coefficient of variation is defined as the standard deviation of the peak height divided by the average peak height, and then multiplied by 100% to be expressed in percentage form. Suppose the peak heights of the above five extreme points are 0.090, 0.110, 0.085, 0.100, 0.120 (unit). Their average peak height is unit. Their standard deviation is unit. Then the coefficient of variation of peak height is .

[0034] Finally, organize the identification number, start wavelength, end wavelength, number of internal extreme points, average peak width, coefficient of variation of peak height, and a preliminary impurity correlation mark (e.g., if the average peak width is within a specific range of 0.5 - 2.0 nm and the coefficient of variation of peak height is less than 30%, then mark as "potential impurity characteristic region") of each independent characteristic section to form a structured impurity characteristic marking matrix. Each row of this matrix represents a characteristic section and its attributes.

[0035] Table 1 Example of impurity characteristic marking matrix:

[0036] As shown in Table 1, some of the divided characteristic sections and their calculated parameters are listed, and these parameters will be used for subsequent dynamic analysis and modeling.

[0037] The full width at half maximum is measured in nanometer wavelength units and is dimensionally normalized by a standard wavelength calibrator.

[0038] Please refer to Figure 3 , the impurity characteristic modeling module includes: Based on the impurity characteristic marking matrix, the peak width dynamic analysis sub-module extracts the characteristic sections marked as impurities and their corresponding sampling points, establishes a measurement mechanism for the spacing between adjacent sampling points, calculates the ratio of the peak width change in the current time window to the change in the previous time window, and uses the sliding window method to perform median filtering on the ratios of five consecutive time windows to generate a dynamic gradient matrix; The impurity characteristic marking matrix (as shown in Table 1) is the input of this sub-module, which contains multiple sections initially identified as "potential impurity characteristic regions" and their static parameters. This module focuses on analyzing the dynamic changes of spectral characteristics within these sections over time. Assume that the sample is continuously monitored, and full-band spectral data is collected seven times in the time series . For the characteristic sections marked as potential impurities, such as section 1 (220 - 280 nm) in Table 1, it is necessary to extract the changes in the peak width of representative sampling points within this section at different time points. Here, the "representative sampling point" can be the peak top wavelength of the peak with the strongest signal within the section, or some average or representative value of all extreme point parameters within the section. For simplicity of explanation, we select the peak at 225.3 nm within section 1 as the monitoring point, and its full peak width (sum of the left and right half peak widths) measured at different times is as follows: nm, nm, nm, nm, nm, nanometer, nanometer.

[0039] For each such monitoring point, within each time window (for ), calculate the ratio of the peak width change in the current time window to the peak width change in the previous time window. This ratio is defined as the dimensionless peak width change rate (where represents the serial number of the characteristic section, represents the serial number of the monitoring point within this section. In this example ). Its calculation formula is: , and the denominator is not zero.

[0040] Taking the monitoring point at 225.3 nanometers as an example: At : . . At : . . At : . . At : . . At : . . Thus, a peak width change rate sequence is obtained: [2.50, -0.60, -2.00, 0.33, -2.00], corresponding to the time points to .

[0041] Next, perform median filtering on this change rate sequence using the sliding window method. The window size is set to five consecutive time windows. This means that the length of the filtered sequence will be shortened. Since our change rate sequence has only 5 points, the filtering operation will act on this complete sequence. Sort these five values: [-2.00, -2.00, -0.60, 0.33, 2.50]. Their median is the third value in the sequence, which is -0.60. In the case of a longer sequence, for example, for a sequence containing change rates, take five consecutive change rate values, find the median of these five values, and use it as the output corresponding to the central time point of these five values after filtering. Move the sliding window forward by one time point and repeat this process. If the original sequence is , the first median filtering result is corresponding to the time (assuming corresponds to ), the second is corresponding to the time , and so on.

[0042] For all monitoring points within all characteristic segments, this peak width change rate calculation and median filtering operation are performed in all time windows. Organize these filtered dimensionless peak width change rate ratios to form a dynamic gradient matrix. The rows of this matrix represent different characteristic segments (or the monitoring points within them), and the columns represent different time windows. For example, for two characteristic segments (monitoring point 1 of segment 1, monitoring point 1 of segment 2), at three effective time points obtained after filtering (assuming as ), the dynamic gradient matrix has the form: Fill in specific values (here, use the unique median -0.60 calculated by us as an example. Assume this is the value of segment 1 at a certain effective time point, and supplement other assumed values to form the matrix): Each element in this matrix is a smoothed dimensionless peak width change rate.

[0043] The baseline offset quantization sub-module calls the baseline trajectory data of the characteristic segment, constructs a moving average model of the reference baseline value, calculates the absolute deviation amount between the baseline values of multiple sampling points and the corresponding reference values, and performs an arithmetic average operation on the deviation amounts of three consecutive sampling points using a three-point sliding window to establish an offset intensity tensor; This sub-module processes the same time series spectral data as the "peak width dynamic analysis sub-module". The focus is on the baseline behavior of the characteristic segments (such as segment 1, 220 - 280 nm) identified in the impurity feature marking matrix (as shown in Table 1). The baseline trajectory data of the characteristic segment refers to the baseline intensity values of the selected representative wavelength points (or the average baseline level of the entire segment) at each time point within this segment. The estimation of the baseline can be achieved by applying a baseline subtraction algorithm to the spectral data at each time point. For example, apply the asymmetric least squares fitting to the spectrum of each characteristic segment to obtain the baseline value of this segment.

[0044] Assume that for characteristic segment 1 (220 - 280 nm), its baseline intensity values at the central wavelength of 250 nm in the time series are respectively: , , , , , , (unit: spectral intensity unit).

[0045] Construct a moving average model of the reference baseline value to smooth short-term fluctuations and obtain a more stable baseline reference. The window length of the moving average Set to 5 time points. This length is chosen to balance effective noise smoothing and quick response to baseline trend changes, typically selecting a window that covers several typical measurement cycles but is not too long. For the time points (where ), its reference baseline value is calculated as the arithmetic mean of the observed baseline values from to for these time points: . For example, at ( , so the summation is from to , that is, using the data of ): units. At : units. At : units. Obtain the reference baseline sequence (starting from ): [0.0336, 0.0356, 0.0364].

[0046] Calculate the absolute deviation amount between the observed baseline value at each sampling point (here referring to Feature Section 1) and the corresponding reference baseline value at each defined reference value time point . At : units. At : units. At : units. Obtain the absolute deviation amount sequence (starting from ): [0.0044, 0.0044, 0.0004].

[0047] For this deviation amount sequence, perform an arithmetic mean operation using a three-point sliding window to further smooth the data. Here, "three-point" refers to three consecutive deviation amount values in the time series. For the smoothed deviation amount at time point is calculated as . Since our deviation amount sequence has only three points (corresponding to ), applying the three-point sliding average will only produce one valid output point, corresponding to the central point . units.

[0048] Organize the smoothed baseline absolute deviation amounts of all feature segments (taking segment 1 as an example here) in all time windows (after moving average processing) to form an offset intensity tensor. The dimensions of this tensor are (number of feature segments number of time windows). If there are multiple feature segments and multiple valid time points after filtering, this tensor will be more complete. Taking the above calculation as an example, if this is the value of segment 1 at a certain valid time point and other assumed values are supplemented, the offset intensity tensor is in the following form (assuming two feature segments and three valid time points ): Fill in specific values (using the calculated 0.00307 and supplementing other assumed values): where each element is the smoothed baseline offset intensity of the corresponding feature segment at the corresponding time point, in spectral intensity units.

[0049] The convolution weight fusion sub-module splits the dynamic gradient matrix into row vectors and performs a dot product operation with the column vectors of the offset intensity tensor, using the formula: ; Normalize the dot product result, accumulate the multi-vector components along the time dimension, and generate a feature-component weight vector; where, represents the weight vector of the i-th feature segment and the j-th component, is the dimensionless peak width change rate, representing the ratio of the peak width change amount at the k-th sampling point of the i-th feature segment to the change amount in the previous time window. i represents the feature segment serial number, k represents the sampling point serial number, is the baseline offset intensity, representing the absolute deviation amount between the baseline value at the k-th sampling point of the j-th feature segment and the moving average value. j represents the feature segment serial number, n represents the total number of sampling points, max(ΔP) represents the maximum peak width gradient of the current feature segment, and max(ΔB) represents the maximum baseline offset intensity of the current feature segment.

[0050] The dynamic gradient matrix is generated by the "peak width dynamic analysis sub-module", and its element represents the th feature segment (or its internal monitoring point , simplified here as representing the segment, representing the time point / sampling point serial number) of the dimensionless peak width change rate. The offset intensity tensor is generated by the "baseline offset quantization sub-module", and its element represents the th feature segment at the The baseline offset intensity at each time point / sampling point (in spectral intensity units). The formula is: Where is the th characteristic section and the th characteristic section (or the characteristic section representing the influence of the th component). is the peak width change rate of the th characteristic section at the th sampling point (or time point). is the baseline offset intensity of the th characteristic section at the th sampling point (or time point). is the total number of sampling points / time points for summation. is the th characteristic section's maximum absolute value in the sequence (although the absolute value is not specified in the original text here, for symmetry with the denominator term and to avoid the influence of negative values, usually the maximum absolute value or directly the maximum value is taken. We follow the original text to take the maximum value, but note that if the sequence contains negative numbers). is the th characteristic section's maximum value in the sequence.

[0051] To ensure the comparability of each term in the formula, especially the summation terms in the denominator, parameters with different physical meanings need to be uniformly represented. is already a dimensionless ratio. And has spectral intensity units. To make contribute reasonably in the formula, especially when added to the denominator and it needs to be converted to a dimensionless form. The conversion rule is: divide each by the average baseline intensity of its corresponding characteristic section during the observation period . That is . The basis for this conversion rule is that the relative magnitude of the baseline offset can better reflect its significance than its absolute value. For example, for section 1, its average baseline intensity units. Then the previously calculated units, and its dimensionless form is . All are processed in this way to obtain . The and in the formula should be understood as using these dimensionless values.

[0052] Now for an example calculation, calculate the characteristic section 1 ( ) and the characteristic section 1 ( ) Its own fusion weight , assuming there are valid time points . Starting from , take the first row as sequence: . Starting from , take the first row as the original sequence: (unit). Nondimensionalize the sequence (assuming the unit is approximately valid for all these points): So, the nondimensional sequence is: .

[0053] Calculate and : . (If the absolute value is the largest, it is 0.70. Here, follow the original formula symbol) .

[0054] Calculate the summation term : Summation ; Calculate : (If is taken as , then . This is more in line with the conventional positive value characteristic of the weight. We adopt this modified understanding.) Assume refers to . Then . .

[0055] This calculated is a quantitative fusion weight for the correlation between the peak width change of feature section 1 and its own (nondimensionalized) baseline shift. Perform this calculation for all feature section pairs to be examined ( ). If the original text "accumulate multi-vector components along the time dimension" means calculating its with all other components / features to form a weight vector , then what is calculated here is an element of this vector. Pool the such weight vectors of all feature sections to generate a set of feature-component weight vectors. This result ​Indicates that the synchronous correlation degree (after normalization and dimensionless processing) between the peak width change of feature section 1 and its own baseline shift is approximately 0.20678. If this value is higher than a preset correlation threshold (e.g., 0.15, which is set by analyzing the distribution, and a value higher than this indicates a strong correlation), it means that the peak shape change of this feature section has a strong coupling with baseline instability. This weight will be used for subsequent cycle drift analysis.

[0056] Please refer to Figure 4 , the cycle drift analysis module includes: The peak cluster alignment sub-module calls the feature-component weight vector to locate the horizontal and vertical coordinates of the peak tops of consecutive peak clusters, calculates the morphological differences between adjacent peak clusters using the peak height difference and the symmetry axis offset, screens valid peak cluster pairs based on the peak height difference threshold determined by spectral standard sample tests, constructs a sequence of peak spacing and symmetry degree differences in time series, and generates a peak cluster alignment sequence; The set of feature-component weight vectors is generated by the previous sub-module, where quantifies the correlation strength between different feature sections or between a feature section and a specific influencing factor. This module uses these weights and combines the original or processed spectral data to locate those (e.g., ) related to the main components (i.e., having a higher or value) with other main components in the time series. For these peak clusters, extract the wavelength (abscissa) of their peak tops and the net peak height (ordinate) after baseline subtraction.

[0057] Suppose at two consecutive time points and , a pair of peaks considered to correspond in time evolution is identified. Peak 1 is at wavelength at nm with a peak height of units; Peak 2 is at wavelength at nm with a peak height of units. Calculate the morphological differences between this pair of peaks: the peak height difference units; the symmetry axis (peak wavelength) offset nm.

[0058] In order to screen out valid peak cluster pairs that truly belong to the same source evolution, a peak height difference threshold needs to be set. The determination of this threshold is based on continuous spectral tests of known-concentration pure standard samples (ideally without significant drift) and standard samples known to have slight concentration or morphological changes.

[0059] Table 2 Experimental data table for determining the peak height difference threshold:

[0060] As shown in Table 2, for a stable standard sample, the difference in peak height between consecutive measurements is usually very small (e.g., not exceeding 0.008 units at most). When considering instrument noise and minor actual fluctuations, a reasonable threshold for peak height difference can be set slightly higher than the maximum difference observed for the stable standard sample, while being much lower than the difference for significant drift or peaks from different sources. Based on the experimental data, 0.010 units is selected as the peak height difference threshold, that is pairs of peaks with a height difference of units less than the threshold units are initially considered valid. For the aforementioned example peak pairs, their height difference

[0061] For all consecutive valid peak cluster pairs in the time series, record their peak spacing (i.e., the offset of the symmetry axis ) and the difference in symmetry. The difference in symmetry here refers to the absolute value of the difference between the symmetry parameters of two peaks (e.g., the peak asymmetry factor calculated through peak shape parameters). Assume that the symmetry parameter of peak 1 obtained through peak shape analysis is (1 represents perfect symmetry), and the symmetry parameter of peak 2 is , then the difference in symmetry is . In this way, over time, two sequences can be constructed: the peak spacing sequence (unit: nanometer): (corresponding to the peak pairs between ) and the difference in symmetry sequence (unitless): These two sequences together form the peak cluster alignment sequence for subsequent dynamic time warping.

[0062] Based on the peak cluster alignment sequence, the dynamic warping sub-module uses the dynamic time warping algorithm to match the peak cluster spacing sequences with different time dimensions, calculates the cumulative path cost of the peak height difference sequence and the symmetry sequence, and adjusts the cumulative amount of morphological differences through the bending path weight coefficient to generate a warping matching matrix; The peak cluster alignment sequence is the output of the previous sub-module, which contains the peak spacing sequence and the difference in symmetry sequence in the time dimension. Taking a peak spacing sequence of length 5 nanometers and a difference in symmetry sequence as an example, these two sequences represent the morphological evolution of a pair of tracked peak clusters in five consecutive time intervals. The purpose of dynamic time warping (DTW) is to find an optimal matching path to compare these two (or more) sequences of characteristics that change over time, even if their rates of change on the time axis are different. In this application, DTW is used to quantify the similarity or difference between the peak cluster evolution patterns in two different time periods (or two different experimental conditions).

[0063] Suppose we want to compare the above sequence (D, S) obtained in the current observation period with a similar sequence (D_ref, S_ref) obtained from a reference template (or the previous observation period). For simplicity, here we describe how to calculate the cumulative cost within a sequence or align it with another sequence. The core of DTW is to construct a cost matrix. Given two sequences and . The element of the cost matrix represents the local cost between points and , which is usually a measure of the difference between them, such as . In this module, the local cost function combines the peak spacing and the difference in symmetry: where is the difference in peak spacing and symmetry of one sequence (or time point), is the corresponding value of another sequence (or another time point). is the warping path weight coefficient, which is used to adjust the contribution of the symmetry difference to the total cost. This coefficient is set to 0.6. This value is determined through a series of simulation experiments: generating synthetic peak cluster evolution sequences with different degrees of time scale distortion and morphological changes, performing DTW alignment using different values (such as from 0.1 to 2.0, with a step size of 0.1), and comparing with the known true correspondences, and selecting the value that gives the highest alignment accuracy (or the smallest alignment error). Experiments show that when is in the range of 0.5 to 0.7, the alignment effect is more balanced in terms of sensitivity to changes in peak spacing and symmetry, so 0.6 is taken.

[0064] The DTW algorithm finds a path from to through dynamic programming, such that the sum of all local costs on the path is minimized. Each step of the path allows moving right, up, or right-up (diagonal). The cumulative cost is calculated as follows: ; Taking the comparison of the peak spacing sequence and another similar reference sequence as an example (here the symmetry is ignored to simplify the display of the DTW path calculation, and the actual calculation will include the symmetry term). The cost matrix : ; The cumulative cost matrix : ... And so on, until the entire matrix is calculated. The final bottom-right This is the DTW distance or the minimum cumulative path cost between these two sequences. The optimal path is found by backtracking.

[0065] Applying the DTW algorithm to the peak cluster alignment sequences, the calculated minimum cumulative path cost and the path itself are organized into a regularized matching matrix. If only the path is to be recorded, this matrix can be a binary matrix of the same size as the cost matrix, with elements on the path being 1 and others being 0. If the cost needs to be recorded, it can be the cumulative cost matrix itself. Here, the regularized matching matrix is mainly used to extract path-related statistics subsequently.

[0066] The drift coefficient generation sub-module calls the regularized matching matrix, extracts the product factor of the path bending cost and the cumulative amount of morphological differences, combines it with the feature-component weight vector for range normalization processing, calculates the periodic drift influence weights corresponding to multiple components, and generates component-related drift coefficients.

[0067] The regularized matching matrix is generated by the dynamic regularization sub-module, and the core is the minimum cumulative path cost and the optimal alignment path calculated by the DTW algorithm. Suppose the minimum cumulative path cost obtained by the previous sub-module (which comprehensively considers the peak spacing and symmetry differences and takes into account the bending path weight) is 0.06 (as in the simplified example ). The degree of path bending is measured by counting the number of non-diagonal steps (i.e., pure horizontal or pure vertical movements) in the optimal path. Suppose in the alignment, the optimal path contains 1 non-diagonal step. Then the path bending cost = the number of non-diagonal steps × the preset unit bending penalty (this penalty can be related to the bending path weight coefficient , or set independently, set to 0.1). The path bending cost = . The cumulative amount of morphological differences is the minimum cumulative path cost calculated by DTW, which is 0.06 here. Multiplying these two quantities gives a comprehensive product factor: Product factor = path bending cost × cumulative amount of morphological differences = .

[0068] Next, this product factor is combined with the feature-component weight vector (from the convolutional weight fusion sub-module, for example , and other values) for range normalization processing. The purpose of range normalization is to map parameters from different sources and different scales to a unified [0,1] interval. Suppose the value of the product factor being processed currently is . It is necessary to determine a reasonable variation range for this type of product factor, which can be obtained by performing DTW analysis on a large amount of historical data or simulated data, for example . Then the normalized value If the calculation result exceeds the range of [0, 1], it is truncated to 0 or 1.

[0069] Multiply this standardized product factor by the elements in the corresponding feature-component weight vector (for example, the component weight most relevant to the currently analyzed peak cluster , set to 0.20678) to obtain the weight of the periodic drift effect for this component: Weight of periodic drift effect = Generalize this process to all relevant components and all analysis time periods to generate a set of component-related drift coefficients. For example, for three main components, a sequence of drift coefficients may be obtained: .

[0070] Please refer to Figure 5 , the component quantification output module includes: The drift compensation sub-module calls the component-related drift coefficients, performs a Hadamard product operation on the initial content values and the drift coefficients, compensates for the concentration offset between adjacent time windows using piecewise linear interpolation, and adjusts the interpolation weights according to the compensation intensity coefficient of floating-point numbers in the range of 0.5 - 1.2 to generate a drift compensation sequence; The initial content value matrix is the concentration values of each component initially calculated based on the spectral data without drift correction through a conventional chemometric model (such as partial least squares regression). Suppose at a certain time point, the initial content values obtained by analyzing three main components are respectively: ppm (parts per million). The sequence of component-related drift coefficients generated by the previous sub-module is . These coefficients represent the degree to which each component is affected by periodic drift.

[0071] Perform a Hadamard product (element-wise product) operation on the initial content values and the corresponding drift coefficients to obtain the absolute drift amount for each component: Drift amount ppm. These drift amounts are the cumulative effects based on an observation period.

[0072] Next, use piecewise linear interpolation to compensate for the concentration offset between adjacent time windows. Here, "adjacent time windows" may refer to a finer time scale, or it may mean distributing the total drift amount to different parts within the measurement period. More directly, use the calculated drift amounts to correct the initial content values. If the drift coefficient indicates a relative change rate, the compensation method may be or , depending on the definition of the drift coefficient. Assume that the drift coefficient here represents the relative proportion to be deducted, and the concentration should not be negative. The compensated concentration . ppm ppm The concentration sequence after drift compensation in ppm: ppm (rounded to four decimal places).

[0073] The "interpolation weight adjustment of the compensation intensity coefficient for floating-point numbers in the range of 0.5 - 1.2" mentioned in the original text may refer to further amplitude adjustment of the calculated drift amount or the compensation behavior itself. There is a compensation intensity coefficient whose value is between 0.5 and 1.2 and is set based on prior knowledge of the system drift characteristics or real-time feedback, for example . This coefficient is used to scale the calculated drift effect. If the drift coefficient is already a correction factor, it may be . Using this adjustment: ppm ppm The final drift compensation sequence in ppm is ppm.

[0074] Based on the drift compensation sequence, the parameter update sub-module constructs the Kalman filter state equation and observation equation, calculates the prior estimate covariance matrix through the prediction step, and calculates the Kalman gain coefficient by fusing the observation noise covariance matrix in the update step, and iteratively outputs the updated concentration parameter set; Based on the drift compensation sequence obtained from the previous sub-module (i.e., the corrected concentration observation value), construct a Kalman filter to iteratively update the equation to obtain a better concentration estimate. Let be the true (but unknown) component concentration at time point , be the observed (drift-compensated) concentration at time point .

[0075] The state equation describes the dynamic evolution of the concentration: The observation equation describes the relationship between the observed value and the true state: Where: is the state transition matrix (or scalar), which describes the change law of the concentration from one time point to the next. Here, let , indicating that the concentration has a certain persistence but will slowly decay or change. The setting of this value is based on long-term observation and model fitting of the stability of the target substance in the monitoring system. For example, if no new substance is added and there is a small amount of consumption, will be slightly less than 1. is the observation matrix (or scalar), usually 1, indicating directly observing the state itself. Let . is the process noise, which follows a Gaussian distribution with a mean of 0 and a covariance of . Represents the uncertainty of the state transition model. Let . The value of is estimated by analyzing the random fluctuation amplitude of the concentration change in historical data or based on the theoretical analysis of the system's physical and chemical processes. For example, the variance of the concentration change introduced by process disturbances. is the observation noise, which follows a Gaussian distribution with a mean of 0 and a covariance of . Represents the uncertainty of the measurement process. Let . The value of is determined according to the square of the standard deviation of the concentration of the standard sample measured repeatedly, reflecting the precision of the measuring instrument and the residual error after drift compensation.

[0076] Kalman filtering process: Initialization: It is necessary to (Initial state estimate) and (Initial estimate covariance). Let ppm (Take the first value of the drift compensation sequence as the initial estimate of component 1). Let (The initial uncertainty is large).

[0077] Prediction step (Taking component 1 as an example, the time is from to ): Prior state estimate: ppm Prior estimate covariance: ; Update step (At when a new observation value is obtained): Let the new observation value (from the next time point of the drift compensation sequence or the observation of different components at the same time point) be ppm. Kalman gain: Posterior state estimate (Updated concentration): ppm Posterior estimate covariance: ; Take and as the and for the next iteration. Repeat the prediction and update steps for each time point in the drift compensation sequence (or for consecutive observations of each component), and iteratively output the updated set of concentration parameters. For example, after several iterations, the updated set of concentration parameters for the three components is: ppm.

[0078] The quantitative analysis sub-module calls the updated set of concentration parameters, calculates the kurtosis coefficient and skewness coefficient of the concentration values of multiple components, constructs a concentration probability density function in combination with the impurity determination threshold determined through three parallel experiments, and outputs a quantitative analysis report of impurities.

[0079] Call the updated concentration parameter set output by the previous sub-module. For example, for a certain key impurity component, the updated concentration values obtained in 5 consecutive measurements are as follows: ppm. Calculate the statistical characteristics of this concentration sequence, including kurtosis coefficient and skewness coefficient, to evaluate the shape of the concentration distribution. Mean ppm Standard deviation ppm Skewness coefficient . For these data, the calculated skewness is . (Positive skewness, slightly longer right tail) Kurtosis coefficient . For these data, the calculated kurtosis is . (Platykurtic peak, lower than the peakedness of the normal distribution) Combine the impurity determination thresholds determined through three parallel experiments to construct a concentration probability density function. Parallel experiments refer to conducting multiple measurements on negative samples known to contain no target impurities and positive samples known to contain target impurities with concentrations close to the critical concentration under the same nominal conditions.

[0080] Table 3 Experimental data table for determining impurity determination thresholds:

[0081] As shown in Table 3, by analyzing the data of three batches of parallel experiments, for example, using the receiver operating characteristic (ROC) curve analysis method, determine the concentration threshold for best distinguishing negative and positive samples in each batch. Combining the results of the three experiments, take the average value as the final impurity determination threshold, that is ppm. Based on this threshold and combined with the historical data or the impurity concentration distribution characteristics under specific process requirements, construct a concentration probability density function (PDF). For example, assume that the impurity concentration approximately follows a distribution with a certain lower mean and standard deviation when it is below the threshold in qualified products, and follows a distribution with a certain higher mean and standard deviation when it exceeds the standard. Or, more simply, the determination threshold can be regarded as a critical point.

[0082] For the currently obtained updated concentration value, for example, the latest ppm (taking the aforementioned mean), compare it with the determination threshold ppm. Since , it is initially judged that the impurity exceeds the standard. For a more refined quantitative analysis, the probability that the current concentration value falls into the "exceeding the standard" interval can be calculated, or its probability density value under a certain preset "impurity concentration distribution model". If it is assumed that the concentration exceeding the standard follows a distribution with as the lower limit, mean of ppm, and standard deviation of For a truncated normal distribution in ppm, calculations can be performed or The output impurity quantification analysis report will include: updated concentration values of current components, comparison with historical data, statistical characteristics (such as mean, standard deviation, skewness, kurtosis), whether the preset impurity determination threshold is exceeded, and risk assessment or confidence analysis based on the probability model.

[0083] The above are only the preferred embodiments of the present invention, and do not limit the present invention in other forms. Any person skilled in the art may use the disclosed technical content to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the technical solution content of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A spectral intelligent analysis system for impurities in chemical fiber raw materials, characterized in that, The system includes: A spectral feature deconstruction module, which is used to scan the full-band spectrum through a moving window, extract peak width and peak height parameters by using the Savitzky-Golay filtering algorithm, divide characteristic sections based on extreme value density, generate an impurity feature marking matrix, and transfer the impurity feature marking matrix to the impurity feature modeling module; An impurity feature modeling module, which is used to receive the impurity feature marking matrix, perform a convolution operation on the peak width change rate and the baseline offset of the characteristic section to generate a feature-component weight vector, and transfer the feature-component weight vector to the periodic drift analysis module; A periodic drift analysis module, which is used to call the feature-component weight vector, perform dynamic time warping matching on the peak height difference and symmetry of the continuous peak cluster, construct a warping matching matrix through path bending cost calculation and morphological difference cumulative amount, generate a component-related drift coefficient, and transfer the component-related drift coefficient to the component quantization output module; A component quantization output module, which is used to perform drift compensation on the initial content based on the component-related drift coefficient, update the concentration parameter by using the Kalman filtering algorithm, and output an impurity quantization analysis report.

2. The intelligent spectral analysis system for impurities in chemical fiber raw materials according to claim 1, wherein The impurity feature marking matrix specifically includes peak width parameters, peak height parameters, and extreme value density partitioning. The feature-component weight vector includes peak width change rate weight, baseline offset weight, and convolution operation coefficient. The component-related drift coefficient specifically refers to peak spacing coefficient of variation, fusion degree attenuation parameter, and dynamic warping matching degree. The impurity quantization analysis report includes xylene derivative concentration, ester by-product concentration, and Kalman gain parameter.

3. The intelligent spectral analysis system for impurities in chemical fiber raw materials according to claim 2, wherein The combined feature extraction method of the Savitzky-Golay filtering algorithm and extreme value density partitioning enhances the peak shape recognition accuracy through second derivative transformation and improves the characteristic section division efficiency by combining extreme point density statistics; The Kalman filtering algorithm uses a state equation and an observation equation , where A is the state transition matrix with a value in the range of 0.95 - 1.05, H is the observation matrix taking the identity matrix, and the process noise covariance Q and the observation noise covariance R are determined through a spectrometer calibration experiment.

4. The intelligent spectral analysis system for impurities in chemical fiber raw materials according to claim 3, characterized in that, The spectral feature deconstruction module includes: A window scanning sub-module acquires full-band spectral data, uses a fixed step size to move a rectangular window to cover the spectral band, records the spectral intensity sequence in each window, and linearly interpolates and splices the overlapping regions of adjacent windows to generate a window spectral sequence covering the full band; A peak shape parameter extraction sub-module, based on the window spectral sequence, performs second derivative transformation on each window spectrum by using the Savitzky-Golay filtering algorithm, calculates the spectral curvature change rate within the window, locates the coordinates of local extreme points, extracts the left half-peak width, right half-peak width, and peak top height values of multiple extreme points, and generates a peak shape parameter set; A characteristic section division sub-module calls the peak shape parameter set, counts the number of extreme points in the unit wavelength interval as the density benchmark, divides the continuous band according to the density threshold optimized by the gradient descent method, merges the boundaries of adjacent high-density intervals, and calculates the peak width mean and peak height coefficient of variation in multiple independent sections to generate an impurity feature marking matrix.

5. The intelligent spectral analysis system for impurities in chemical fiber raw materials according to claim 4, wherein The half-peak width measurement uses a nanometer wavelength unit and is dimensionally normalized by a standard wavelength calibrator.

6. The intelligent spectral analysis system for impurities in chemical fiber raw materials according to claim 5, wherein The impurity feature modeling module includes: Based on the impurity feature marking matrix, the peak width dynamic analysis sub-module extracts the characteristic sections marked as impurities and their corresponding sampling points, establishes a measurement mechanism for the spacing between adjacent sampling points, calculates the ratio of the peak width change in the current time window to the change in the previous time window, and uses the sliding window method to perform median filtering on the ratios of five consecutive time windows to generate a dynamic gradient matrix; The baseline offset quantization sub-module calls the baseline trajectory data of the characteristic section, constructs a moving average model of the reference baseline value, calculates the absolute deviation between the baseline values of multiple sampling points and the corresponding reference values, and performs arithmetic average operations on the deviations of three consecutive sampling points using a three-point sliding window to establish an offset intensity tensor; The convolution weight fusion sub-module splits the dynamic gradient matrix by row vectors and performs dot product operations with the column vectors of the offset intensity tensor, using the formula: ； Normalize the dot product results, accumulate the multi-vector components along the time dimension, and generate a feature-component weight vector; Among them, represents the weight vector of the i-th characteristic section and the j-th component, is the dimensionless peak width change rate, representing the ratio of the peak width change amount at the k-th sampling point in the i-th characteristic section to the change amount in the previous time window. i represents the characteristic section serial number, and k represents the sampling point serial number. is the baseline offset intensity, representing the absolute deviation amount between the baseline value at the k-th sampling point in the j-th characteristic section and the moving average value. j represents the characteristic section serial number, n represents the total number of sampling points, max(ΔP) represents the maximum peak width gradient in the current characteristic section, and max(ΔB) represents the maximum baseline offset intensity in the current characteristic section.

7. The intelligent spectral analysis system for impurities in chemical fiber raw materials according to claim 6, wherein The cycle drift analysis module includes: The peak cluster alignment sub-module calls the feature-component weight vector, locates the horizontal and vertical coordinates of the peak tops of consecutive peak clusters, calculates the morphological differences between adjacent peak clusters using the peak height difference and the symmetry axis offset, screens out valid peak cluster pairs based on the peak height difference threshold determined by spectral standard sample tests, constructs a sequence of peak spacing and symmetry difference according to the time series, and generates a peak cluster alignment sequence; The dynamic programming sub-module, based on the peak cluster alignment sequence, uses the dynamic time warping algorithm to match the peak cluster spacing sequences with different time dimensions, calculates the cumulative path cost of the peak height difference sequence and the symmetry sequence, and adjusts the cumulative amount of morphological differences through the bending path weight coefficient to generate a programming matching matrix; The drift coefficient generation sub-module calls the programming matching matrix, extracts the product factor of the path bending cost and the cumulative amount of morphological differences, performs range normalization processing in combination with the feature-component weight vector, calculates the cycle drift influence weights corresponding to multiple components, and generates component-related drift coefficients.

8. The intelligent spectral analysis system for impurities in chemical fiber raw materials according to claim 7, wherein The component quantization output module includes: The drift compensation sub-module calls the component-related drift coefficients, performs Hadamard product operations on the initial content values and the drift coefficients, compensates for the concentration offset between adjacent time windows using the piecewise linear interpolation method, and adjusts the interpolation weights according to the compensation intensity coefficient of floating-point numbers in the 0.5 - 1.2 interval to generate a drift compensation sequence; The parameter update sub-module constructs a Kalman filter state equation and an observation equation based on the drift compensation sequence, calculates the prior estimation covariance matrix through the prediction step, and fuses the observation noise covariance matrix in the update step to calculate the Kalman gain coefficient, and iteratively outputs an updated concentration parameter set; The quantization analysis sub-module calls the updated concentration parameter set, calculates the kurtosis coefficient and skewness coefficient of the concentration values of multiple components, constructs a concentration probability density function in combination with the impurity determination threshold determined by three parallel experiments, and outputs a quantitative analysis report on impurities.

Images