A similarity visual analysis method for a momordica grosvenori raw material fingerprint
By using curvature energy constraint-based map alignment and 3D topographic map analysis, the quality control challenge in the analysis of large-scale data of monk fruit raw materials was solved, enabling rapid identification of abnormal batches and deviations in process parameters, and improving the intelligence and traceability of quality control.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- GUILIN SANLENG BIOTECH CO LTD
- Filing Date
- 2026-04-02
- Publication Date
- 2026-07-24
Smart Images

Figure CN122449047A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of visualization analysis technology, specifically a method for visual analysis of the similarity of fingerprint spectra of monk fruit raw materials. Background Technology
[0002] Currently, those skilled in the art widely use high performance liquid chromatography or ultra-high performance liquid chromatography to construct chemical fingerprints of monk fruit. They then evaluate the consistency of quality by calculating the similarity values (such as the cosine of the angle and the correlation coefficient) between the sample chromatogram and the standard fingerprint chromatogram, or by using the "overlay method" that directly superimposes multiple chromatograms.
[0003] However, existing technologies have significant limitations when processing large-scale data analysis of monk fruit raw materials. First, traditional similarity evaluation only outputs a macroscopic statistical value, often masking the microscopic differences in specific chromatographic peaks. This makes it impossible for technicians to identify which specific component (such as mogroside V or specific impurities) caused the decrease in similarity. Second, when analyzing dozens or even hundreds of batches of raw materials, the traditional overlay method results in densely interwoven and overlapping chromatographic curves, easily causing visual confusion. This makes it difficult for those skilled in the art to quickly identify abnormal batches or subtle evolutionary patterns between batches with the naked eye. Therefore, existing technologies lack an analytical method that can both preserve the rich information of fingerprint spectra and intuitively reduce the dimensionality of multidimensional data differences, failing to meet the needs for precise and efficient quality control of monk fruit raw materials.
[0004] To address this, a visualization analysis method for the similarity of fingerprint spectra of monk fruit raw materials is proposed. Summary of the Invention
[0005] The purpose of this invention is to provide a similarity visualization analysis method for the fingerprint spectrum of monk fruit raw materials. By introducing rigid weight constraints based on curvature energy, the fingerprint spectrum is elastically aligned, and the alignment result is mapped to the spatial deformation of a ring elastic grid. Multiple batches of three-dimensional topographic maps are constructed, and combined with Laplace field analysis, the quality difference is visualized and the process fault is diagnosed.
[0006] To achieve the above objectives, the present invention provides the following technical solution: A method for visual analysis of the similarity of fingerprint spectra of monk fruit raw materials, comprising: Smooth filtering and second derivative calculation are performed on the standard monk fruit fingerprint spectrum to obtain the curvature energy value at each time point; the curvature energy value is then mapped to a rigid weight vector. Using the rigid weight vector as the path penalty weight of the dynamic time warping algorithm, the test map and the standard map are elastically aligned; the aligned test map is projected onto the standard map and the projection coefficient is calculated; the difference between the test map and the projection component is calculated to obtain the residual vector orthogonal to the standard map. A ring-shaped elastic mesh model is constructed, and the projection coefficients are transformed into centripetal contraction loads acting on all nodes. The magnitude of the residual vector is transformed into the normal pushing load of the corresponding node. The elastic equilibrium equation is solved to obtain the node displacements and generate the deformed topological surface. Multiple batches of topological surfaces are stacked and interpolated along the batch axis to form a three-dimensional terrain map. Calculate the Laplace field of the 3D topographic map, extract the time coordinates and batch range of the Laplace extreme points, perform pattern matching with the preset process fault knowledge base, and output process diagnosis conclusions.
[0007] Preferably, the process of obtaining the rigid weight vector includes: smoothing the absorbance time series data of the standard monk fruit fingerprint using a multinomial convolution filtering method, and simultaneously calculating the second derivative value at each time point; squaring the second derivative value to obtain the curvature energy value at each time point, wherein the curvature energy value at the characteristic peak position is higher than the curvature energy value of the baseline region; normalizing the curvature energy value, and converting the normalized curvature energy value into a rigid weight value through linear mapping to generate a rigid weight vector corresponding one-to-one with the time axis.
[0008] Preferably, the elastic alignment process includes: constructing a distance matrix between the spectrum to be tested and the standard spectrum, wherein each element in the distance matrix is the product of the square of the absorbance difference at the corresponding time point and the rigidity weight value at that time point; searching for the optimal curved path from the starting point to the ending point in the distance matrix, wherein the optimal curved path minimizes the cumulative distance and is subject to path continuity constraints and boundary constraints; resampling the spectrum to be tested along the time axis along the optimal curved path to obtain spectrum data aligned with the time axis of the standard spectrum; wherein, the rigidity weight value serves as a path penalty weight, which increases the path curvature cost in the characteristic peak region with high rigidity weight, thereby limiting the temporal distortion amplitude in that region, while the baseline region with low rigidity weight allows for temporal distortion amplitude to absorb the overall drift.
[0009] Preferably, the annular elastic mesh model includes a node layer, a connection layer, and a load layer, specifically: The node layer is used to generate multiple nodes in a ring distribution within a virtual plane. The angular position of each node corresponds to the retention time of the fingerprint map, and the initial radial position of each node is equal. The connection layer is used to establish virtual elastic connections between adjacent nodes, and the stiffness coefficient of each connection is taken from the stiffness weight value at the corresponding time point in the stiffness weight vector; The load layer is used to calculate the centripetal contraction load based on the projection coefficient and apply it to all nodes. The magnitude of the centripetal contraction load is proportional to the difference between the projection coefficient and a preset threshold. The load layer is also used to calculate the normal pushing load based on the residual vector and apply it to nodes with non-zero residuals. The magnitude of the normal pushing load is proportional to the magnitude of the residual vector at the corresponding node and its direction is a radial direction away from the center of the circle. The node layer provides the mesh topology, the stiffness distribution of the connection layer determines the local resistance to deformation of the mesh, and the load layer drives the mesh to deform.
[0010] Preferably, the process of forming the three-dimensional topographic map includes: assembling the overall load vector of the grid system based on the centripetal contraction load and normal jacking load applied by the load layer; assembling the overall stiffness matrix of the grid system based on the stiffness coefficient of the connection layer; solving the linear equation system composed of the overall stiffness matrix and the overall load vector to obtain the displacement of each node in the radial direction; superimposing the initial coordinates of each node with the displacement to obtain the final coordinates of each node after deformation, the final coordinates constituting a single batch of topological surfaces; arranging the topological surfaces of multiple batches in batch order along the direction perpendicular to the surface, and using an interpolation method to generate a continuous transition surface between corresponding nodes of adjacent batches, forming a three-dimensional topographic map with batch as one coordinate axis, retention time as another coordinate axis, and radial displacement as the height value.
[0011] Preferably, the pattern matching process includes: calculating the second-order partial derivatives of the height field data of the three-dimensional topographic map along the batch direction and the time direction respectively, and adding the second-order partial derivatives in the two directions to obtain the Laplace value of each grid point; detecting local extreme points in the Laplace field, and recording the Laplace value sign, time coordinate, and batch coordinate of each extreme point; comparing the time coordinate of each extreme point with the feature time window in the preset process fault knowledge base, wherein the feature time window corresponds to the peak time range of chemical components; associating the Laplace value sign of the extreme point with the defect type in the knowledge base, wherein a positive Laplace extreme value corresponds to a concave defect and a negative Laplace extreme value corresponds to a convex defect.
[0012] Preferably, the process of obtaining the process diagnostic conclusion includes: determining the chemical composition type and defect type corresponding to the extreme point based on the pattern matching result; retrieving fault cause entries that match the chemical composition type and defect type from the process fault knowledge base, wherein the fault cause entries include process parameter names and deviation directions; determining the batch range where the anomaly occurred based on the batch coordinates of the extreme point; and generating and outputting a process diagnostic conclusion that includes an abnormal batch identifier, defect location, cause, and adjustment suggestions by combining the chemical composition type, defect type, fault cause entries, and abnormal batch range.
[0013] In summary, the present invention adopts the above technical solution, and the present invention has the following technical effects: 1. This invention constructs a rigid weight vector based on curvature energy by smoothing and filtering the standard Luo Han Guo fingerprint spectrum and performing second-order derivative operations. This weight vector is then introduced into a dynamic time warping algorithm as a path penalty weight, enabling the elastic alignment process to have regional differential constraint capabilities. This invention can effectively suppress excessive stretching or compression of the time axis in the characteristic peak region, while allowing the baseline region to absorb overall drift errors. Therefore, even in the presence of retention time drift, injection differences, or system fluctuations, it can still maintain the structural consistency of the peak positions of key chemical components, significantly improving the physical rationality and repeatability of the fingerprint spectrum alignment results.
[0014] 2. This invention overcomes the limitations of traditional fingerprint map similarity methods, which only present similarity through numerical indicators or simple curve overlays. It introduces a ring-shaped elastic grid model, mapping the map projection coefficients and residual vectors to centripetal contraction loads and normal pushing loads, respectively. By solving the elastic equilibrium equations, the continuous deformation of the grid nodes is obtained, and a three-dimensional topographic map is further constructed. This method can intuitively transform the similarity differences, local deviations, and batch evolution trends between maps into topographic undulations and spatial structure changes. This allows quality fluctuations to move beyond abstract judgments of "similarity" and instead form a visualized result with spatial continuity and directionality, significantly improving the intuitiveness and engineering readability of fingerprint map analysis in quality control scenarios.
[0015] 3. Based on obtaining multiple batches of three-dimensional topographic maps, this invention further calculates the Laplace field of the height field and extracts extreme points. The time coordinates, batch ranges, and convex / concave morphology of these extreme points are then matched with a pre-set process fault knowledge base to automatically identify the peak times of abnormal chemical components and defect types. This invention can directly correlate structural differences in fingerprint maps to the deviation direction of specific process parameters and abnormal batch intervals, forming diagnostic conclusions that include defect location, cause analysis, and adjustment suggestions. This constructs a closed-loop mechanism from data analysis to process decision-making, significantly improving the intelligence and traceability of Luo Han Guo raw material quality control. Attached Figure Description
[0016] Figure 1 A flowchart of a method for visualizing the similarity of fingerprint spectra of monk fruit raw materials provided by the present invention; Figure 2 This is a schematic diagram of the topological surface acquisition process provided by the present invention; Figure 3 This is a schematic diagram of the pattern matching process provided by the present invention. Detailed Implementation
[0017] 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 preferred embodiments. However, it should be noted that many details listed in the specification are merely to provide the reader with a thorough understanding of one or more aspects of the present invention, and these aspects of the invention can be implemented even without these specific details.
[0018] Example 1 like Figure 1 As shown, the similarity visualization analysis method for the fingerprint spectrum of Luo Han Guo raw materials provided by the present invention includes the following steps: smoothing filtering and second derivative calculation of the standard Luo Han Guo fingerprint spectrum to obtain the curvature energy value at each time point; mapping the curvature energy value into a rigid weight vector; using the rigid weight vector as the path penalty weight of the dynamic time warping algorithm to elastically align the test spectrum and the standard spectrum; projecting the aligned test spectrum onto the standard spectrum and calculating the projection coefficient; calculating the difference between the test spectrum and the projection component to obtain the similarity with the standard spectrum. The residual vectors are orthogonal to the spectrum; a ring-shaped elastic mesh model is constructed, and the projection coefficients are transformed into centripetal contraction loads acting on all nodes, and the magnitude of the residual vectors is transformed into normal jacking loads for the corresponding nodes; the elastic equilibrium equations are solved to obtain the node displacements, and the deformed topological surface is generated; multiple batches of topological surfaces are stacked and interpolated along the batch axis to form a three-dimensional topographic map; the Laplace field of the three-dimensional topographic map is calculated, the time coordinates and batch ranges of the Laplace extreme points are extracted, and pattern matching is performed with a pre-set process fault knowledge base to output process diagnosis conclusions.
[0019] Furthermore, the process of obtaining the rigid weight vector includes: smoothing the absorbance time series data of the standard monk fruit fingerprint using a multinomial convolution filtering method, and simultaneously calculating the second derivative value at each time point; squaring the second derivative value to obtain the curvature energy value at each time point, wherein the curvature energy value at the characteristic peak position is higher than the curvature energy value of the baseline region; normalizing the curvature energy value, and converting the normalized curvature energy value into a rigid weight value through linear mapping to generate a rigid weight vector corresponding one-to-one with the time axis.
[0020] Specifically, the absorbance time series data of the standard monk fruit fingerprint spectrum is smoothed and filtered. The purpose of smoothing and filtering is to suppress the influence of high-frequency noise on the subsequent calculation of the second derivative and to avoid amplifying the noise.
[0021] The polynomial convolution filtering method is adopted. The basic principle of this method is: at each data point, a certain number of adjacent data points are taken to form a sliding window. Within the window, a polynomial is used to perform least squares fitting on the data, and the function value of the fitting polynomial at the center point of the window is used as the smoothed data value.
[0022] The selection of the filter window width needs to balance smoothing effect and signal fidelity. A window that is too narrow will result in insufficient noise suppression, while a window that is too wide may lead to peak broadening and reduced peak height. For routine HPLC data, a window width of 5 to 15 data points is recommended. Specifically, if the half-width at half-maximum (WHM) of a characteristic peak corresponds to n sampling points, then the window width should be the larger of n divided by 3 and 5. For example, if the WHM of a characteristic peak corresponds to 15 sampling points, the window width can be set to 5 points; if the WHM of a characteristic peak corresponds to 30 sampling points, the window width can be set to 10 points. This achieves a balance between noise suppression and peak shape preservation.
[0023] While performing smoothing filtering, the second derivative values at each time point are calculated simultaneously using the differential properties of the polynomial convolution filter. Specifically, the convolution coefficients used to obtain the second derivative are pre-calculated, and these coefficients are convolved with the original data within the sliding window to directly obtain the second derivative value at the window's center point.
[0024] The curvature energy values at each time point are obtained by squaring the obtained second derivative values. The calculation of the curvature energy values is based on the following: Mathematically, the curvature of a plane curve is related to its second derivative; the larger the absolute value of the second derivative, the more severe the curvature of the curve. For chromatographic peaks, the curvature is greatest when the curve changes from rising to falling at the peak apex, and the second derivative is a negative extreme value; in the region where the slope changes significantly on both sides of the peak, the absolute value of the second derivative is also relatively large; in the region where the baseline is flat, the curve is approximately a straight line, and the second derivative approaches zero.
[0025] The squaring operation eliminates the influence of the sign of the second derivative, transforming both negative extreme values at the peak and positive values elsewhere into positive curvature energy values. A larger curvature energy value indicates a more drastic signal change at that location, making it more likely to be a region containing a characteristic peak. To avoid noise amplification during the second derivative calculation, it is recommended to perform thresholding on the second derivative value before squaring. Specifically, calculate the standard deviation of the second derivative values in the flat baseline region, set values whose absolute value is less than twice this standard deviation to zero, and then perform the squaring operation. This effectively suppresses noise interference in the baseline region.
[0026] The curvature energy value is normalized to map its range to between zero and one. The normalization method is as follows: find the maximum and minimum values of curvature energy value at all time points, and for each time point, calculate its curvature energy value, subtract the minimum value, and divide by the difference between the maximum and minimum values to obtain the normalized curvature energy value.
[0027] The purpose of normalization is to eliminate the differences in absolute values between different batches of standard spectra, so that the subsequent mapping of rigid weights is consistent.
[0028] The normalized curvature energy values are converted into rigid weight values through linear mapping. The mapping method is as follows: set a lower limit and an upper limit for the rigid weights. For each time point, the rigid weight value is equal to the lower limit plus the normalized curvature energy value multiplied by the difference between the upper and lower limits. The lower limit ranges from 0.1 to 0.3, and the upper limit ranges from 1.0 to 2.0. The lower limit represents the weight of the baseline region. If the value is too small, the baseline region will be almost unconstrained, allowing excessive time distortion. If the value is too large, it will over-constrain the baseline and affect the ability to compensate for overall drift. The upper limit represents the weight of the characteristic peak region. If the value is too small, it will not be enough to protect the temporal consistency of the characteristic peak region. If the value is too large, it may lead to over-constraint of the peak region and be unable to adapt to reasonable retention time changes. For routine analysis with small peak shape changes, a lower limit of 0.2 and an upper limit of 1.5 can be selected. For analysis with large peak shape changes or high precision requirements, a lower limit of 0.1 and an upper limit of 2.0 can be selected.
[0029] The mapping result is as follows: feature peak regions with higher curvature energy values receive higher rigidity weight values, while baseline regions with lower curvature energy values receive lower rigidity weight values. Finally, a rigidity weight vector corresponding one-to-one with the time axis is generated, which will be used to constrain the time distortion amplitude in the subsequent elastic alignment step.
[0030] Furthermore, the elastic alignment process includes: constructing a distance matrix between the target spectrum and the standard spectrum, wherein each element in the distance matrix is the product of the square of the absorbance difference at the corresponding time point and the rigidity weight value at that time point; searching for the optimal curved path from the starting point to the ending point in the distance matrix, wherein the optimal curved path minimizes the cumulative distance and is subject to path continuity constraints and boundary constraints; resampling the target spectrum along the time axis along the optimal curved path to obtain target spectrum data aligned with the time axis of the standard spectrum; wherein, the rigidity weight value serves as a path penalty weight, which increases the path curvature cost in the characteristic peak region with high rigidity weight, thereby limiting the temporal distortion amplitude in that region, while the baseline region with low rigidity weight allows for temporal distortion amplitude to absorb the overall drift.
[0031] Specifically, a distance matrix is constructed between the spectrum to be tested and the standard spectrum. The spectrum to be tested contains absorbance data at M time points, and the standard spectrum contains absorbance data at N time points. The distance matrix is a two-dimensional array of M rows and N columns, where the element in the i-th row and j-th column represents the weighted distance between the i-th time point of the spectrum to be tested and the j-th time point of the standard spectrum.
[0032] The weighted distance is calculated as follows: First, the square of the absorbance difference between two time points is calculated, and then multiplied by the rigidity weight value corresponding to the j-th time point of the standard spectrum. The weight value of the standard spectrum is used here instead of the weight value of the spectrum under test because the standard spectrum represents an ideal reference state, and its weight vector is used to constrain the degree to which the spectrum under test undergoes elastic deformation towards the standard spectrum. The rigidity weight value acts as a multiplication factor, amplifying the distance in high-rigidity weight regions and reducing the distance in low-rigidity weight regions.
[0033] The optimal curved path from the starting point to the ending point is searched in the distance matrix; the starting point is defined as the lower left element of the distance matrix, corresponding to the pairing of the first time point of the test map and the first time point of the standard map; the ending point is defined as the upper right element of the distance matrix, corresponding to the pairing of the last time point of the test map and the last time point of the standard map.
[0034] Path search follows these constraints: a path must start from a starting point and end at a destination; the difference between the row index and the column index between any two adjacent paired points on the path does not exceed one, ensuring the continuity of the path; and the row index and column index on the path are monotonically increasing, ensuring the unidirectional nature of time.
[0035] The optimal path is defined as the path that minimizes the sum of distances between paired points among all paths that satisfy the constraints. Dynamic programming can be used to efficiently search for the optimal path: starting from the starting point, the minimum cumulative distance to each matrix element is calculated sequentially, and the optimal path is eventually obtained by backtracking.
[0036] The impact of rigid weights on path search is reflected in the following ways: Since the distance is amplified in regions with high rigid weights, the algorithm tends to avoid generating large time offsets in these regions during the search process, as offsets lead to large cumulative distances. Conversely, in regions with low rigid weights, even if there are large time offsets, their contribution to the cumulative distance is relatively small, and the algorithm can tolerate larger time distortions in these regions to absorb the overall drift.
[0037] The time axis of the spectrum under test is resampled along the optimal curved path; the optimal path effectively establishes a nonlinear mapping relationship between the time axis of the spectrum under test and the time axis of the standard spectrum. For each time point of the standard spectrum, a corresponding time point of the spectrum under test can be found based on the path. When the pairing is one-to-one, the absorbance value of the spectrum under test at that time point is directly taken; when the pairing is one-to-many, the absorbance values of multiple time points of the spectrum under test are averaged.
[0038] After resampling, the target spectrum data is obtained and aligned with the time axis of the standard spectrum. The aligned target spectrum has the same number of time points and time scale as the standard spectrum, eliminating the influence of retention time drift on subsequent analysis.
[0039] Furthermore, the ring-shaped elastic mesh model includes a node layer, a connection layer, and a load layer, specifically: The node layer is used to generate multiple nodes in a ring distribution within a virtual plane. The angular position of each node corresponds to the retention time of the fingerprint map, and the initial radial position of each node is equal. The connection layer is used to establish virtual elastic connections between adjacent nodes, and the stiffness coefficient of each connection is taken from the stiffness weight value at the corresponding time point in the stiffness weight vector; The load layer is used to calculate the centripetal contraction load based on the projection coefficient and apply it to all nodes. The magnitude of the centripetal contraction load is proportional to the difference between the projection coefficient and a preset threshold. The load layer is also used to calculate the normal pushing load based on the residual vector and apply it to nodes with non-zero residuals. The magnitude of the normal pushing load is proportional to the magnitude of the residual vector at the corresponding node and its direction is a radial direction away from the center of the circle. The node layer provides the mesh topology, the stiffness distribution of the connection layer determines the local resistance to deformation of the mesh, and the load layer drives the mesh to deform.
[0040] Specifically, multiple nodes are generated in a ring-shaped distribution within a virtual plane. The number of nodes corresponds to the number of time sampling points in the fingerprint spectrum and can be appropriately simplified according to the required computational accuracy. The virtual plane uses a Cartesian coordinate system, with a polar coordinate mapping relationship established with the origin as the center.
[0041] A linear mapping relationship is established between the angular position of each node and the retention time of the fingerprint spectrum. If the analysis time range of the fingerprint spectrum is zero to T, then the angle corresponding to the node with a retention time of t is t divided by T and then multiplied by 360 degrees.
[0042] The initial radial position of each node is set to an equal baseline radius, forming a regular circular distribution. The recommended range for the baseline radius is 50 to 200 units. The specific value can be set according to visualization requirements; for example, if the display window width is W pixels, the baseline radius can be set to W divided by 4. The value of the baseline radius does not affect the relative results of deformation calculations, only the scale of the visualization display. A baseline radius that is too small may lead to numerical accuracy issues, while a baseline radius that is too large may lead to calculation overflow.
[0043] Virtual elastic connections are established between adjacent nodes. Each node establishes an elastic connection with each of its two adjacent nodes in the circumferential direction, forming a closed loop structure that connects end to end.
[0044] The stiffness coefficient of each elastic connection is taken from the stiffness weight value at the corresponding time point in the stiffness weight vector. Specifically, the stiffness coefficient of the elastic connection connecting the i-th node and the (i+1)-th node is taken as the average of the stiffness weight values at the i-th time point and the (i+1)-th time point.
[0045] The physical meaning of stiffness coefficient is the ability of an elastic connection to resist expansion and contraction. The higher the stiffness coefficient, the smaller the elongation or shortening of the connection under the same load; the lower the stiffness coefficient, the more easily the connection deforms.
[0046] The stiffness distribution of the connecting layer makes the mesh in the region corresponding to the characteristic peak more holistic and less susceptible to damage by local loads; while the mesh in the region corresponding to the baseline is relatively soft and more likely to undergo local deformation in response to loads.
[0047] The load layer calculates and applies two types of loads based on the projection coefficient and the residual vector.
[0048] The first type is the global centripetal load. When the projection coefficient is less than a preset threshold (usually set to 1.0, indicating consistency with the standard content), a global centripetal load is calculated and applied to all nodes. The magnitude of the global centripetal load is equal to the degree of insufficient projection coefficient (1 minus the projection coefficient) multiplied by the load coefficient. The recommended value range for the load coefficient is 5 to 20, and the specific value should be adjusted according to the desired deformation amplitude. The load direction is the radial direction from each node to the center of the circle, and its component in the Cartesian coordinate system is: the negative value of the node's initial coordinate multiplied by the load magnitude and then divided by the distance from the node to the origin. The effect of the global centripetal load is to cause the entire mesh to tend to shrink towards the center, simulating the physical effect of insufficient overall content.
[0049] The second type is the local normal load. For time points where the magnitude of the residual vector is non-zero, a local normal load is applied at the corresponding node. The magnitude of the local normal load is equal to the magnitude of the residual vector at the corresponding time point multiplied by the load factor. The recommended value range for the load factor is 2 to 8; the larger the value, the more pronounced the deformation.
[0050] The load direction is radial, pointing from the center of the circle to the node, i.e., away from the center. Its component in the Cartesian coordinate system is: the initial coordinate of the node multiplied by the load magnitude and then divided by the distance from the node to the origin. The effect of the local normal load is to cause the corresponding node to bulge outward, simulating the physical effect of the presence of local impurities.
[0051] When a node is simultaneously subjected to a global centripetal load and a local normal load, the two types of loads are vectored together, and the final net load depends on the magnitude comparison of the two types of loads.
[0052] Furthermore, the formation process of the three-dimensional topographic map includes: assembling the overall load vector of the mesh system based on the centripetal contraction load and normal jacking load applied by the load layer; assembling the overall stiffness matrix of the mesh system based on the stiffness coefficient of the connection layer; solving the linear equation system composed of the overall stiffness matrix and the overall load vector to obtain the displacement of each node in the radial direction; superimposing the initial coordinates of each node with the displacement to obtain the final coordinates of each node after deformation, the final coordinates constituting a single batch of topological surfaces; arranging the topological surfaces of multiple batches in batch order along the direction perpendicular to the surface, and using an interpolation method to generate a continuous transition surface between corresponding nodes of adjacent batches, forming a three-dimensional topographic map with batch as one coordinate axis, retention time as another coordinate axis, and radial displacement as the height value.
[0053] Specifically, the overall load vector of the mesh system is assembled based on the centripetal and normal loads applied by the load layer; the dimension of the overall load vector is twice the number of nodes, because each node has two degrees of freedom in the plane, namely the displacement component along the horizontal direction and the displacement component along the vertical direction.
[0054] For each node, the net load acting on that node is decomposed into horizontal and vertical components, which are then filled into the corresponding positions of the overall load vector.
[0055] The overall stiffness matrix of the mesh system is assembled based on the stiffness coefficients of the connection layers; the dimension of the overall stiffness matrix is a square matrix with twice the number of nodes multiplied by twice the number of nodes. The assembly of the stiffness matrix follows the standard procedure of the finite element method: first, the element stiffness matrix of each elastic connection is calculated, and then the stiffness matrices of each element are superimposed onto the corresponding positions of the overall stiffness matrix according to the node numbers.
[0056] The element stiffness matrix of a single elastic connection is related to the stiffness coefficient, connection direction, and connection length of the connection. The larger the connection stiffness coefficient, the larger the corresponding stiffness matrix element values, indicating a stronger constraint effect of the connection on nodal displacements.
[0057] Solving the system of linear equations consisting of the global stiffness matrix and the global load vector yields the displacements at each node. The linear equations are in the form that the stiffness matrix multiplied by the displacement vector equals the load vector. This system of equations can be solved using either a direct method or an iterative method. For a small number of nodes, direct methods such as Gaussian elimination can obtain an exact solution; for a large number of nodes, iterative methods such as the conjugate gradient method can obtain a sufficiently accurate approximate solution.
[0058] The resulting displacement vector contains the horizontal and vertical displacement components of each node. Based on the initial angular positions of the nodes, the displacement components in the Cartesian coordinate system can be converted into radial displacement components in the polar coordinate system.
[0059] The initial coordinates of each node are superimposed with the displacement to obtain the final coordinates of each node after deformation. For each node, its final radial position is equal to the initial reference radius plus the radial displacement. When the radial displacement is negative, the node contracts towards the center; when the radial displacement is positive, the node bulges outward.
[0060] The coordinates of all nodes after deformation form a single batch of topological surfaces, as referenced. Figure 2 The deformed node coordinates are connected sequentially to form a closed curve, which visually demonstrates the differences between this batch of samples and the standard.
[0061] Multiple batches of topological surfaces are arranged in batch order along a direction perpendicular to the surface, and interpolation methods are used to generate continuous three-dimensional terrain maps.
[0062] Suppose there are K batches of samples to be analyzed, then K topological surfaces are generated. These surfaces are arranged at equal intervals along the vertical direction according to their batch numbers, from one to K, with the spacing between adjacent surfaces representing the batch interval.
[0063] Interpolation is used to generate a smooth transition between corresponding nodes in adjacent batches. The interpolation method can be linear interpolation or bilinear interpolation. The former is simpler to calculate but the surface is not smooth enough, while the latter produces a smooth surface but requires more computation.
[0064] After interpolation, a 3D topographic map is generated, with the batch as one coordinate axis, the angle of the retained time mapping as another coordinate axis, and the radial displacement as the height value. This topographic map allows for a visual observation of the mass distribution of each batch of samples and the trend of mass variation across batches.
[0065] Furthermore, the pattern matching process refers to Figure 3 The process includes: calculating the second-order partial derivatives of the height field data of the three-dimensional topographic map along the batch direction and the time direction respectively, and adding the second-order partial derivatives in the two directions to obtain the Laplace value of each grid point; detecting local extreme points in the Laplace field, and recording the Laplace value sign, time coordinate, and batch coordinate of each extreme point; comparing the time coordinate of each extreme point with the characteristic time window in the preset process fault knowledge base, wherein the characteristic time window corresponds to the peak time range of chemical components; and associating the Laplace value sign of the extreme point with the defect type in the knowledge base, wherein positive Laplace extreme values correspond to concave defects and negative Laplace extreme values correspond to convex defects.
[0066] Specifically, the Laplace field is calculated from the height field data of the 3D topographic map; the height values of the 3D topographic map are stored as a two-dimensional array, with row indices corresponding to batches and column indices corresponding to time sampling points. Second-order partial derivatives are calculated for this two-dimensional array along both row and column directions.
[0067] The second partial derivative along the row direction is calculated using the central difference scheme: for the element in the i-th row and j-th column, its second partial derivative along the row direction is equal to the sum of the elements in the (i+1)-th row and j-th column and the elements in the (i-1)-th row and j-th column, minus twice the sum of the elements in the i-th row and j-th column, and then divided by the square of the row spacing.
[0068] The second-order partial derivatives along the column direction are calculated using the same central difference scheme; the second-order partial derivatives in both directions are added to obtain the Laplace value for each grid point. The Laplace value characterizes the degree and direction of the difference between the height of that point and the average height of its neighborhood.
[0069] Local extrema are detected in a Laplace field. Each grid point in the Laplace field is traversed, and its Laplace value is compared with the Laplace values of its neighboring grid points. Neighboring grid points are defined as the eight adjacent points surrounding the current point. If a point's Laplace value is greater than all its neighbors and greater than a positive threshold, it is identified as a positive extrema; if a point's Laplace value is less than all its neighbors and less than a negative threshold, it is identified as a negative extrema. The thresholds are set as follows: calculate the mean and standard deviation of the Laplace values of all internal grid points excluding the boundary region; set the positive extrema threshold to twice the mean plus the standard deviation, and the negative extrema threshold to twice the mean minus the standard deviation. This adaptive thresholding method automatically adjusts the detection sensitivity based on the actual noise level of each batch of data, avoiding false alarms or abnormal omissions.
[0070] For each detected extreme point, its Laplacian value sign, time coordinate (i.e., the retention time corresponding to the column index), and batch coordinate (i.e., the batch number corresponding to the row index) are recorded. The time coordinates of each extreme point are compared with the characteristic time windows in the pre-set process fault knowledge base; the process fault knowledge base pre-stores the peak time ranges of each characteristic component in the mogro fingerprint spectrum, including but not limited to: the peak time window of mogroside V, the peak time window of symmenidine I, and the peak time window of mogroside alcohol, etc.
[0071] For each extreme point, determine which feature time window its time coordinate falls into. If it falls into the time window of a certain feature component, mark the extreme point as an anomaly related to that component; otherwise, mark it as an unknown anomaly. Associate the Laplacian value of the extreme point with the defect type in the knowledge base: a positive Laplacian extreme value indicates that the point is lower than the surrounding area, corresponding to a concave feature on the topological surface, associated with a concave defect, physically meaning that the content of that component is insufficient. A negative Laplacian extreme value indicates that the point is higher than the surrounding area, corresponding to a convex feature on the topological surface, associated with a convex defect, physically meaning that there is an abnormal component or excessive impurities at that location.
[0072] Furthermore, the process of obtaining the process diagnostic conclusion includes: determining the chemical composition type and defect type corresponding to the extreme point based on the pattern matching result; retrieving fault cause entries that match the chemical composition type and defect type from the process fault knowledge base, wherein the fault cause entries include the process parameter name and deviation direction; determining the batch range where the anomaly occurred based on the batch coordinates of the extreme point; and generating and outputting a process diagnostic conclusion that includes an abnormal batch identifier, defect location, cause, and adjustment suggestions by combining the chemical composition type, defect type, fault cause entries, and abnormal batch range.
[0073] Based on the pattern matching results, the chemical composition type and defect type corresponding to each extreme point are determined; by combining the time window comparison results and the defect type association results, complete label information is generated for each extreme point, including: chemical composition name or unknown, defect type is concave or convex.
[0074] Retrieve fault cause entries from the process fault knowledge base that match the chemical composition type and defect type; the fault cause entries in the process fault knowledge base are stored in a structured manner, and each entry includes: applicable chemical composition, applicable defect type, fault cause description, name of the process parameter involved, and direction of parameter deviation.
[0075] The search logic is as follows: using the chemical composition and defect type corresponding to the extreme point as search criteria, the system searches the knowledge base for completely matching fault cause entries. If multiple matching entries exist, all are returned for technical personnel to refer to; if no matching entries exist, a default general prompt is returned.
[0076] Example fault cause entries include: When mogroside V corresponds to a concave defect, the fault cause may be that the extraction temperature is too high, leading to glycoside degradation. The relevant parameter is the extraction temperature, and the deviation direction is too high. When the elution end corresponds to a convex defect, the fault cause may be that the macroporous resin is not fully regenerated, resulting in the residue of non-polar impurities. The relevant parameter is the amount of regeneration solution, and the deviation direction is too low.
[0077] The batch range where the anomaly occurred is determined based on the batch coordinates of the extreme point. For a single, isolated extreme point, the abnormal batch is the single batch containing that extreme point. For multiple extreme points occurring consecutively in the batch direction, they are merged into one abnormal event, and the abnormal batch range is from the batch containing the first extreme point to the batch containing the last extreme point. Consecutive anomalies may indicate a systematic process deviation and require close monitoring.
[0078] Based on the comprehensive chemical composition type, defect type, fault cause entries, and abnormal batch range, a complete process diagnostic conclusion is generated and output. The format of the process diagnostic conclusion includes the following fields: Abnormal Batch Identifier, listing all batch numbers or batch ranges where the abnormality occurred; Defect Location, indicating the time position corresponding to which chemical composition the abnormality occurred; Defect Type, indicating whether it is a content deficiency type or an impurity excess type; Possible Cause, referencing the fault cause description matched in the knowledge base; Adjustment Suggestions, providing corresponding process parameter adjustment directions based on the fault cause. The diagnostic conclusion can be output in text report form or in structured data form for downstream systems to use.
[0079] This invention first performs smoothing filtering and second-derivative calculations on the standard Luo Han Guo fingerprint spectrum, using curvature energy to characterize the degree of signal change, and constructs a rigid weight vector corresponding one-to-one with the time axis, enabling differentiated constraints between the characteristic peak region and the baseline region in subsequent analysis. Based on this rigid weight vector, constrained dynamic time warping is performed on the test spectrum and the standard spectrum to effectively correct for retained time drift while maintaining the structural stability of key characteristic peaks. Subsequently, the aligned test spectrum is decomposed into projection components and residual components, and a ring-shaped elastic mesh model is constructed, mapping the overall similarity difference and local anomalies to the centripetal contraction deformation and local convex deformation of the mesh, respectively, generating single-batch topological surfaces. By stacking and interpolating multiple batch topological surfaces, a three-dimensional topographic map reflecting the batch evolution characteristics is formed. Furthermore, the Laplace field is calculated on the three-dimensional topographic map and extreme points are extracted. The Laplace field is then matched with the process fault knowledge base to automatically output abnormal components, abnormal batches and corresponding process diagnostic conclusions. This achieves a comprehensive evaluation effect of the Luo Han Guo raw material fingerprint map from similarity analysis and visualization to process tracing.
[0080] Example 2 The overall method flow of this invention includes four stages: rigid weight benchmark modeling, orthogonal residual decoupling, finite element deformation visualization, and process fault tracing. Example 1 details the specific steps of each stage. This example enhances the three stages of extreme point detection, anomaly feature analysis, and diagnostic conclusion generation.
[0081] In extreme point detection, a fixed threshold is difficult to adapt to changes in noise levels across different batches of data. This embodiment employs an adaptive thresholding method based on background statistics. The adaptive threshold specifically includes: Statistical analysis is performed on all grid points in the Laplace field. After excluding grid points located within a preset boundary range, the mean and standard deviation of the Laplace values of the remaining grid points are calculated as an estimate of the background noise level. The positive extreme value detection threshold is set to a preset positive multiple of the mean plus the standard deviation, and the negative extreme value detection threshold is set to a preset negative multiple of the mean minus the standard deviation. The adaptive threshold is used instead of the fixed threshold for extreme point detection, so that the detection sensitivity is automatically adjusted according to the actual noise level of each batch of data.
[0082] Specifically, the grid points in the Laplace field are first statistically analyzed, and grid points located within a preset boundary range are excluded. The mean and standard deviation of the Laplace values of the remaining grid points are then calculated as an estimate of the background noise level. Boundary points are excluded because the central difference scheme lacks complete neighborhood information at the boundaries, which may lead to boundary effects.
[0083] Then, the positive extreme value detection threshold is set to a preset positive multiple of the mean plus the standard deviation, and the negative extreme value detection threshold is set to a preset negative multiple of the mean minus the standard deviation. The selection of the preset multiples is based on statistical principles to keep the probability of exceeding the threshold within the expected significance level.
[0084] An adaptive threshold is used instead of a fixed threshold for extreme point detection. When the background noise level is high, the threshold is automatically increased to avoid false alarms; when the background noise level is low, the threshold is automatically decreased to allow subtle anomalies to be detected.
[0085] Single-resolution analysis may miss both macroscopic trends and microscopic local anomalies. This embodiment introduces a multi-resolution pyramid analysis method. Multi-scale analysis steps: The three-dimensional terrain map is downsampled step by step. Each downsampling step reduces the resolution in both the batch direction and the time direction to half of the previous step, generating a multi-resolution pyramid that includes the original resolution layer and multiple low-resolution layers. The Laplace field is calculated and extreme point features are extracted at each level of the multi-resolution pyramid. The extreme point features extracted from each level are correlated across scales. Extreme points that appear only in high-resolution layers are marked as local micro anomalies, while extreme points that appear in multiple resolution layers are marked as macro trend anomalies. The labeled multi-scale anomaly features are then uniformly input into the subsequent pattern matching steps.
[0086] Specifically, the 3D topographic map is first downsampled step by step. Each downsampling step reduces the resolution in both the batch direction and the time direction to half that of the previous step, generating a multi-resolution pyramid containing the original resolution layer and multiple low-resolution layers. Low-pass filtering can be performed before downsampling to avoid aliasing.
[0087] Then, the Laplace field is calculated at each level and extreme point features are extracted. High-resolution layers preserve detailed information and are suitable for detecting microscopic anomalies; low-resolution layers highlight trend information and are suitable for detecting macroscopic anomalies.
[0088] Finally, the extreme point features extracted from each level are correlated across scales, and the positions of extreme points at each level are uniformly mapped to the original resolution coordinate system for comparison. Extreme points appearing only in high-resolution layers are marked as local micro-anomalies, while extreme points appearing in multiple resolution layers are marked as macro-trend anomalies. The marked multi-scale anomaly features are then uniformly input into the subsequent pattern matching steps.
[0089] Process diagnostic conclusions primarily focus on the spatial characteristics of anomalies, paying insufficient attention to the evolution patterns at the batch level. This embodiment introduces temporal correlation analysis to distinguish between systematic drift and sporadic failures. Batch-to-batch temporal correlation analysis: The height difference between adjacent batches is calculated along the batch direction of the three-dimensional topographic map to obtain the batch difference field; The cumulative difference at each time position is calculated for the batch differential field to characterize the changing trend of the component corresponding to that time position with the batch. The anomaly type is determined based on the change pattern of the differential cumulative amount: if the differential cumulative amount shows a monotonically increasing or monotonically decreasing trend, it is determined to be a gradual anomaly, indicating the existence of a systematic process drift; if the differential cumulative amount shows a pulse pattern of sudden change followed by recovery, it is determined to be a sudden anomaly, indicating the existence of an intermittent fault. The anomaly type determination results are added to the process diagnosis conclusion.
[0090] Specifically, the height difference between adjacent batches is first calculated along the batch direction on the 3D topographic map to obtain the batch difference field. A positive difference value indicates that the radial displacement of the current batch at that location is greater than that of the previous batch, and a negative difference value indicates that it is less than that of the previous batch.
[0091] Then, the cumulative difference for each time position in the batch difference field is calculated, which is to sum the differences of that position in all batch differences sequentially. The cumulative difference reflects the changing trend of the corresponding component at that position across batches.
[0092] The anomaly type is determined based on the changing pattern of the differential cumulative amount. If the differential cumulative amount shows a monotonically increasing or decreasing trend, it is identified as a gradual anomaly, indicating a systematic process drift. If the differential cumulative amount shows a pulse pattern of sudden change followed by recovery, it is identified as a sudden anomaly, indicating an intermittent fault. The anomaly type determination results are added to the process diagnostic conclusion, with an anomaly evolution type field added, and targeted adjustment suggestions provided.
[0093] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.
Claims
1. A method for visual analysis of the similarity of fingerprint spectra of monk fruit raw materials, characterized in that, include: Smoothing filtering and second derivative calculation were performed on the standard Luo Han Guo fingerprint spectrum to obtain the curvature energy value at each time point; Map the curvature energy value to a rigid weight vector; The rigid weight vector is used as the path penalty weight of the dynamic time warping algorithm to flexibly align the test map with the standard map. Project the aligned spectrum to be tested onto the standard spectrum direction and calculate the projection coefficients; Calculate the difference between the spectrum to be measured and the projected components, and obtain the residual vector orthogonal to the standard spectrum; A ring-shaped elastic mesh model is constructed, and the projection coefficients are transformed into centripetal contraction loads acting on all nodes. The magnitude of the residual vector is transformed into the normal pushing load of the corresponding node. The elastic equilibrium equation is solved to obtain the node displacements and generate the deformed topological surface. Multiple batches of topological surfaces are stacked and interpolated along the batch axis to form a three-dimensional terrain map. Calculate the Laplace field of the 3D topographic map, extract the time coordinates and batch range of the Laplace extreme points, perform pattern matching with the preset process fault knowledge base, and output process diagnosis conclusions.
2. The method for visualizing the similarity of fingerprint spectra of monk fruit raw materials according to claim 1, characterized in that, The process of obtaining the rigid weight vector includes: smoothing the absorbance time series data of the standard monk fruit fingerprint using a multinomial convolution filter method, and simultaneously calculating the second derivative value at each time point; squaring the second derivative value to obtain the curvature energy value at each time point, wherein the curvature energy value at the characteristic peak position is higher than the curvature energy value of the baseline region; normalizing the curvature energy value, and converting the normalized curvature energy value into a rigid weight value through linear mapping to generate a rigid weight vector corresponding one-to-one with the time axis.
3. The method for visualizing the similarity of fingerprint spectra of monk fruit raw materials according to claim 1, characterized in that, The elastic alignment process includes: constructing a distance matrix between the target spectrum and the standard spectrum, where each element of the distance matrix is the product of the square of the absorbance difference at the corresponding time point and the rigidity weight value at that time point; searching for the optimal curved path from the starting point to the ending point in the distance matrix, where the optimal curved path minimizes the cumulative distance and is subject to path continuity constraints and boundary constraints; resampling the target spectrum along the time axis along the optimal curved path to obtain target spectrum data aligned with the time axis of the standard spectrum; wherein, the rigidity weight value serves as a path penalty weight, increasing the path curvature cost in the characteristic peak region with high rigidity weight, thereby limiting the temporal distortion amplitude in that region, while the baseline region with low rigidity weight allows for temporal distortion amplitude to absorb the overall drift.
4. The method for visual analysis of the similarity of fingerprint spectra of monk fruit raw materials according to claim 1, characterized in that, The ring-shaped elastic mesh model includes a node layer, a connection layer, and a load layer, specifically: The node layer is used to generate multiple nodes in a ring distribution within a virtual plane. The angular position of each node corresponds to the retention time of the fingerprint map, and the initial radial position of each node is equal. The connection layer is used to establish virtual elastic connections between adjacent nodes, and the stiffness coefficient of each connection is taken from the stiffness weight value at the corresponding time point in the stiffness weight vector; The load layer is used to calculate the centripetal contraction load based on the projection coefficient and apply it to all nodes. The magnitude of the centripetal contraction load is proportional to the difference between the projection coefficient and a preset threshold. The load layer is also used to calculate the normal pushing load based on the residual vector and apply it to nodes with non-zero residuals. The magnitude of the normal pushing load is proportional to the magnitude of the residual vector at the corresponding node and its direction is a radial direction away from the center of the circle. The node layer provides the mesh topology, the stiffness distribution of the connection layer determines the local resistance to deformation of the mesh, and the load layer drives the mesh to deform.
5. The method for visualizing the similarity of fingerprint spectra of monk fruit raw materials according to claim 4, characterized in that, The formation process of the three-dimensional topographic map includes: assembling the overall load vector of the grid system based on the centripetal contraction load and normal jacking load applied by the load layer; assembling the overall stiffness matrix of the grid system based on the stiffness coefficient of the connection layer; solving the linear equation system composed of the overall stiffness matrix and the overall load vector to obtain the displacement of each node in the radial direction; superimposing the initial coordinates of each node with the displacement to obtain the final coordinates of each node after deformation, the final coordinates constituting a single batch of topological surfaces; arranging the topological surfaces of multiple batches in batch order along the direction perpendicular to the surface, and using an interpolation method to generate continuous transition surfaces between corresponding nodes of adjacent batches, forming a three-dimensional topographic map with batch as one coordinate axis, retention time as another coordinate axis, and radial displacement as the height value.
6. The method for visualizing the similarity of fingerprint spectra of monk fruit raw materials according to claim 5, characterized in that, The pattern matching process includes: calculating the second-order partial derivatives of the height field data of the three-dimensional topographic map along the batch direction and the time direction respectively, and adding the second-order partial derivatives in the two directions to obtain the Laplace value of each grid point; detecting local extreme points in the Laplace field, and recording the Laplace value sign, time coordinate, and batch coordinate of each extreme point; comparing the time coordinate of each extreme point with the feature time window in the preset process fault knowledge base, wherein the feature time window corresponds to the peak time range of chemical components; associating the Laplace value sign of the extreme point with the defect type in the knowledge base, wherein positive Laplace extreme values correspond to concave defects and negative Laplace extreme values correspond to convex defects.
7. The method for visualizing the similarity of fingerprint spectra of monk fruit raw materials according to claim 6, characterized in that, The process of obtaining the process diagnostic conclusion includes: determining the chemical composition type and defect type corresponding to the extreme point based on the pattern matching result; retrieving fault cause entries that match the chemical composition type and defect type from the process fault knowledge base, wherein the fault cause entries include the process parameter name and deviation direction; determining the batch range where the anomaly occurred based on the batch coordinates of the extreme point; and generating and outputting a process diagnostic conclusion that includes an abnormal batch identifier, defect location, cause, and adjustment suggestions by combining the chemical composition type, defect type, fault cause entries, and abnormal batch range.