An automatic peak searching method for XRD patterns

Abstract

This invention discloses an automatic peak-finding method for X-ray diffraction (XRD) spectra, relating to the field of phase identification based on XRD spectra. The method mainly includes the following steps: S1: Constructing wide and narrow window symmetrical zero-area transformation functions to transform the spectral data and form a transformation result group; S2: Comparison method and peak flattening processing; S3: Background subtraction, intensity screening, and width margin calculation; S4: Data fusion of the transformation results, using the "AND" fusion criterion to screen out channels not within the peak position characteristic range; S5: Cluster analysis to further divide peak positions into peak and non-peak positions, obtaining the final peak position determination result. Using the automatic peak-finding method proposed in this invention for XRD spectrum peak finding eliminates the need for manual parameter setting, accurately identifying spectral peak positions. This solves the problem of traditional peak-finding methods requiring frequent parameter adjustments based on input spectral data characteristics, thus limiting their application and bringing convenience to phase identification.

Description

Technical Field

[0001] This invention belongs to the field of phase identification based on XRD patterns, and specifically relates to an automatic peak finding method for XRD patterns. Background Technology

[0002] X-ray diffraction (XRD) is a common phase analysis method used for mineral detection. Diffraction patterns are related to the phase composition; different phases produce different diffraction patterns. Comparing the diffraction pattern with a standard powder diffraction card (PDF card) can determine the phase composition of a sample. Locating spectral peaks and comparing them with the PDF card is a crucial step in determining the presence of a specific phase. Common peak-finding methods often suffer from parameter setting issues in practical applications. For example, the high-performance symmetric zero-area transform method requires manually setting three parameters: full width at half maximum (FWHM), half-window width, and threshold. These parameter settings need to consider the spectral data characteristics and be continuously adjusted according to actual requirements, limiting the method's applicability in different scenarios. Therefore, this invention proposes an automatic peak-finding method for XRD patterns that can identify peak positions in spectral data without parameter setting, overcoming the limitations of existing peak-finding methods and paving the way for further phase identification and automatic detection. Summary of the Invention

[0003] The purpose of this invention is to provide an automatic peak finding method for XRD patterns. This method filters out channels within a limited feature range by examining the results of various transformations, and further divides peak positions and non-peak positions using improved K-means clustering analysis. Ultimately, given an input pattern, peak position information can be obtained without setting any parameters.

[0004] The technical solution of the present invention is: an automatic peak finding method for XRD patterns, the method comprising the following steps:

[0005] Step 1: Symmetrical zero-area transformation between wide and narrow windows

[0006] Gaussian linear function Constructing transformation functions The transformation function C simultaneously satisfies zero area and symmetry. C j =C -j The transformation function is determined by the parameter H (full width and height). G Determined by the half-window width m, set H G =5, m are 2, 3, 4 and 5 respectively, construct transformation functions covering different window widths, comprehensively utilize wide window functions to suppress noise, and narrow window functions to improve peak recognition.

[0007] The spectral data are transformed sequentially using a sliding window function, and the transformation results are grouped together. According to the characteristics of the window function, the transformation result at the peak is a large positive number, while the transformation result at the non-peak position is small, or even negative;

[0008] Step 2: Comparison method and peak / substrate processing

[0009] A comparison method is introduced to ensure that the peak position satisfies the peak position characteristics within the left and right cmpN addresses. cmpN = 2 is selected, and the intensity magnitudes of the two channels to the left and right of channel i are compared to establish a relationship vector cmpA(i) = [a i1 a i2 a i3 a i4 ],

[0010] a i1 The intensity comparison results for track address i and track address i-1 are given by a. i2 The intensity comparison results for channel address i and channel address i+1 are given by a. i3 The strength comparison results for track address i-1 and track address i-2 are shown in a. i4 The intensity comparison results for channel addresses i+1 and i+2 are given by sgn(), which is the sign function. If there is an element -1 in cmpA(i), then channel address i does not meet the peak position requirement; if all elements in cmpA(i) are 1, then channel address i is a standard peak position, and the maximum value exists only on one channel address; when the peak value exists on more than one channel, there is an element 0 in cmpA(i), and flat-top peak identification is required.

[0011] Establish a search matrix If cmpA(i) is a row element in srchM, then the address i is a flat peak;

[0012] When the identification result is a standard peak or a flat-topped peak, add element 1 to the end of the result group; otherwise, add element -1.

[0013] Step 3: Baseline deduction, strength screening, and width margin calculation

[0014] Step 3.1 Press The raw spectral data is processed to reduce the large fluctuations in the counts of each channel. y(i) is the count of channel address i, and V(i) is the transformation result.

[0015] Step 3.2 According to Iteratively obtain the full spectrum background, where p is the iteration number and V is the background. p (i) is the result of the p-th iteration, V0(i) = V(i), and the total number of iterations is 10;

[0016] Step 3.3 Press The processed data is inversely transformed back to its original form, y SNIP The final result of background subtraction using the improved Simpson-SNIP algorithm;

[0017] Step 3.4 Set an intensity threshold based on the statistical index of spectral data for intensity screening. If the intensity of the channel address after background subtraction is greater than the threshold, the result is 1; otherwise, it is 0 or -1. Add the result to the result group and take the intensity threshold as 0.8 times the mean intensity of the spectral data after background subtraction.

[0018] Step 3.5 Calculate the width margin of the road address after screening, that is, the number of road addresses that meet the strength screening on the left and right. Compare the width margin with the lower limit of the margin. If the width margin is greater than the lower limit of the margin, the result is 1; otherwise, it is 0 or -1. Add the result to the result group and take the lower limit of the margin as 4.

[0019] Step 4: Data fusion, screening and filtering out road addresses outside the feature range.

[0020] Based on the transformation process of the first three steps, the peak position characteristic range is when all elements in the channel address result group are positive. Using the "AND" fusion criterion, if there is a non-positive number in a certain channel address result group, it is judged to be a non-peak position channel address.

[0021] Step 5: Cluster analysis to further distinguish between peak positions and non-peak positions, and obtain the final peak position determination result.

[0022] Step 5.1 Select the intensity value after background subtraction and the average value of the SZAC results as clustering features, and use the Z-score method to standardize the data to eliminate the influence of units. * = (x-μ) / σ, where μ is the sample mean and σ is the sample variance;

[0023] Step 5.2 Calculate the distance between each pair of sample points and the sum of the distances from each object to the other objects. Select the point corresponding to the maximum value as the first isolated point. Delete the isolated point and repeat step 5.2 until I isolated points are removed. Let I = 4.

[0024] Step 5.3 Set the number of clusters to 2, and find the two sample points that are farthest apart after removing outliers as the initial cluster centers;

[0025] Step 5.4: Using the conventional K-means algorithm, traverse the samples, find the cluster center closest to the sample point, assign the sample to that cluster, and update the cluster center until the maximum number of iterations is reached or the cluster center no longer changes. Save the final cluster center information.

[0026] Step 5.5 Calculate the inter-class distance D between the two classes. b The sum of the average distances within the class D w According to D b / D w Determine the number of clusters. If the result is greater than 1.25, then the optimal number of clusters is determined to be 2; otherwise, the number of clusters is determined to be 1.

[0027] Step 5.6 If the number of clusters is 2, the isolated points are classified according to the nearest distance priority principle; if the number of clusters is 1, all sample points are classified into the same class.

[0028] Step 5.7 If the number of clusters is 2, the sample with the largest average value of the SZAC transformation result is labeled as the peak class, and the other class is labeled as the non-peak class. If the number of clusters is 1, all samples are labeled as peak classes.

[0029] Compared with the closest prior art, the beneficial effects of the present invention are as follows:

[0030] 1. By introducing the comparison method and transformation function group, using the wide window function to suppress noise and the narrow window function to improve the peak recognition ability, the problems of misjudgment and missed judgment are solved, and the accuracy of recognition results is improved;

[0031] 2. Processing spectral data based on a set of transformation functions reduces the impact of source spectrum randomness, solves the problem of frequent parameter adjustments for different spectral signals, enhances the applicability of the method, and facilitates automated detection.

[0032] 3. By introducing a width margin while performing intensity screening, the peak positions within the feature range are guaranteed to have a certain width, effectively reducing the misjudgment of peak positions caused by abrupt noise.

[0033] 4. Introducing cluster analysis further divides peak positions and non-peak positions, reducing the possibility of misclassifying non-significant peak positions among the remaining positions after limited screening as peak positions. Under this premise, the feature range can be expanded to a certain extent, retaining more possible positions during the limited screening process and reducing the peak position omission rate. Attached Figure Description

[0034] Figure 1 This is a peak finding result diagram of the XRD pattern of a sample with known phase;

[0035] Figure 2 This is a flowchart of the automatic peak finding algorithm;

[0036] Figure 3 This is a graph showing the results of background subtraction;

[0037] Figure 4 This is a diagram illustrating the elements contained in the transformation result group and their meanings;

[0038] Figure 5 This is a schematic diagram of the data transformation and fusion process;

[0039] Figure 6 This is a clustering result diagram; Detailed Implementation

[0040] The technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings.

[0041] This invention provides an automatic peak finding method for XRD patterns. Figure 1 The peak-finding results of the method of the present invention for XRD patterns are shown. The sample crystal composition includes kaolinite and lithium. The X-ray diffractometer model is D8, with copper target, voltage 35KV, voltage 30mA, and step size 0.02. Comparison with the standard card of crystal composition shows that the peak positions of the two are basically consistent, and the automatic peak-finding effect is good.

[0042] Figure 2 The implementation process of the method described in this embodiment is illustrated, including the following steps:

[0043] Step 1: Symmetrical zero-area transformation between wide and narrow windows

[0044] Gaussian linear function Constructing transformation functions The transformation function C simultaneously satisfies zero area and symmetry. C j =C -j The transformation function is determined by the parameter H (full width and height). G Determined by the half-window width m, set H G =5, m are 2, 3, 4 and 5 respectively, construct transformation functions covering different window widths, comprehensively utilize wide window functions to suppress noise, and narrow window functions to improve peak recognition.

[0045] The spectral data are transformed sequentially using a sliding window function, and the transformation results are grouped together. According to the characteristics of the window function, the transformation result at the peak is a large positive number, while the transformation result at the non-peak position is small, or even negative;

[0046] Step 2: Comparison method and peak / substrate processing

[0047] A comparison method is introduced to ensure that the peak position satisfies the peak position characteristics within the left and right cmpN addresses. cmpN = 2 is selected, and the intensity magnitudes of the two channels to the left and right of channel i are compared to establish a relationship vector cmpA(i) = [a i1 a i2 a i3 a i4 ],

[0048] a i1 The intensity comparison results for track address i and track address i-1 are given by a. i2 The intensity comparison results for channel address i and channel address i+1 are given by a. i3 The strength comparison results for track address i-1 and track address i-2 are shown in a. i4The intensity comparison results for channel addresses i+1 and i+2 are given by sgn(), which is the sign function. If there is an element -1 in cmpA(i), then channel address i does not meet the peak position requirement; if all elements in cmpA(i) are 1, then channel address i is a standard peak position, and the maximum value exists only on one channel address; when the peak value exists on more than one channel, there is an element 0 in cmpA(i), and flat-top peak identification is required.

[0049] Establish a search matrix If cmpA(i) is a row element in srchM, then the address i is a flat peak;

[0050] When the identification result is a standard peak or a flat-topped peak, add element 1 to the end of the result group; otherwise, add element -1.

[0051] Step 3: Baseline deduction, strength screening, and width margin calculation

[0052] Step 3.1 Press The raw spectral data is processed to reduce the large fluctuations in the counts of each channel. y(i) is the count of channel address i, and V(i) is the transformation result.

[0053] Step 3.2 According to Iteratively obtain the full spectrum background, where p is the iteration number and V is the background. p (i) is the result of the p-th iteration, V0(i) = V(i), and the total number of iterations is 10;

[0054] Step 3.3 Press The processed data is inversely transformed back to its original form, y SNIP The final result of background subtraction using the improved Simpson-SNIP algorithm; Figure 3 The background subtraction effect in this embodiment is shown. Subfigure (a) is the original XRD pattern, and subfigure (b) is the background subtraction result.

[0055] Step 3.4 Set an intensity threshold based on the statistical index of spectral data for intensity screening. If the intensity of the channel address after background subtraction is greater than the threshold, the result is 1; otherwise, it is 0 or -1. Add the result to the result group and take the intensity threshold as 0.8 times the mean intensity of the spectral data after background subtraction.

[0056] Step 3.5 Calculate the width margin of the road address after screening, that is, the number of road addresses that meet the strength screening on the left and right. Compare the width margin with the lower limit of the margin. If the width margin is greater than the lower limit of the margin, the result is 1; otherwise, it is 0 or -1. Add the result to the result group and take the lower limit of the margin as 4.

[0057] Step 4: Data fusion, screening and filtering out road addresses outside the feature range.

[0058] Based on the transformation process in the first three steps, the peak position characteristic range is defined as a channel address result group where all elements are positive. Using the "AND" fusion criterion, if a channel address result group contains a non-positive number, it is determined to be a non-peak position channel address. Figure 4 The transformation result group shows the elements and meanings contained in it. Subgraph (a) is the peak position channel address result group in this example, all of which are positive. Subgraph (b) is the non-peak position channel address result group, at least one of which is a non-positive number.

[0059] Figure 5 The data transformation and fusion process was demonstrated;

[0060] Step 5: Cluster analysis to further distinguish between peak positions and non-peak positions, and obtain the final peak position determination result.

[0061] Step 5.1 Select the intensity value after background subtraction and the average value of the SZAC results as clustering features, and use the Z-score method to standardize the data to eliminate the influence of units. * = (x-μ) / σ, where μ is the sample mean and σ is the sample variance;

[0062] Step 5.2 Calculate the distance between each pair of sample points and the sum of the distances from each object to the other objects. Select the point corresponding to the maximum value as the first isolated point. Delete the isolated point and repeat step 5.2 until I isolated points are removed. Let I = 4.

[0063] Step 5.3 Set the number of clusters to 2, and find the two sample points that are farthest apart after removing outliers as the initial cluster centers;

[0064] Step 5.4: Using the conventional K-means algorithm, traverse the samples, find the cluster center closest to the sample point, assign the sample to that cluster, and update the cluster center until the maximum number of iterations is reached or the cluster center no longer changes. Save the final cluster center information.

[0065] Step 5.5 Calculate the inter-class distance D between the two classes. b The sum of the average distances within the class D w According to D b / D w Determine the number of clusters. If the result is greater than 1.25, then the optimal number of clusters is determined to be 2; otherwise, the number of clusters is determined to be 1.

[0066] Step 5.6 If the number of clusters is 2, the isolated points are classified according to the nearest distance priority principle; if the number of clusters is 1, all sample points are classified into the same class.

[0067] Step 5.7 If the number of clusters is 2, the sample with the largest average value of the SZAC transformation result is labeled as the peak class, and the other class is labeled as the non-peak class. If the number of clusters is 1, all samples are labeled as peak classes.

[0068] Figure 6 The clustering results for this example are shown in the diagram. In this embodiment, the optimal number of clusters is 2.

[0069] The program in this example runs on MATLAB R2017a.

[0070] Finally, it should be noted that the described embodiments are only some, not all, of the embodiments in this application. All other embodiments obtained by those skilled in the art based on the embodiments in this application without inventive effort are within the scope of protection of this application.

Claims

1. An automatic peak searching method of XRD pattern, characterized in that, The method comprises the following steps: First step: wide and narrow window symmetric zero area conversion (SZAC) A group of peak-like symmetrical zero-area transformation functions is constructed by using Gaussian linear function, and the full width at half maximum H G =5, and the half window width m is 2, 3, 4, and 5 respectively, the spectrum data is transformed, and a group of transformation results is formed. Second step: comparison method and flat top peak processing The comparison method is introduced to ensure that the peak position meets the peak position characteristics within the left and right cmpN trace addresses, cmpN=2 is selected, the intensity of the left and right trace addresses i is compared, and a relationship vector cmpA(i) is established, cmpA(i) = [a i1 a i2 a i3 a i4 ], a i1 is the intensity comparison result of track address i and track address i-1, a i2 is the intensity comparison result of track address i and track address i+1, a i3 is the intensity comparison result of track address i-1 and track address i-2, a i4 is the intensity comparison result of track address i+1 and track address i+2, sgn() is the sign function, if there is an element -1 in cmpA(i), the track address i does not meet the peak position requirement; if all elements in cmpA(i) are 1, the track address i is a standard peak position; if there is an element 0 in cmpA(i) and the peak width is greater than 1, flat top peak identification is required; Establish search matrix If cmpA(i) is a row element in srchM, then trace i is a flat peak. When the identification result is a standard peak or a flat top peak, an element 1 is added at the end of the result group, otherwise an element -1 is added; Third step: background subtraction, intensity screening and width margin calculation Step 3.1 is as follows The raw spectral data is processed to reduce large fluctuations in the counts of each channel, y(i) is the count at channel i, and V(i) is the result of the transformation; Step 3.2 according to The full spectrum background is iteratively obtained, p is the iteration number, V p (i) is the result of the pth iteration, V0(i) = V(i), and the total iteration number is 10; Step 3.3 is performed as in The processed data is inverse transformed to the original form, y SNIP Final result of the background subtraction using the improved Simpson-SNIP algorithm; Step 3.4 sets an intensity threshold based on the spectral data statistical index for intensity screening. The trace address with an intensity greater than the threshold after background subtraction is processed as 1, otherwise as 0 or -1. The result is added to the result group. The intensity threshold is 0.8 times the average intensity of the spectral data after background subtraction; Step 3.5 calculates the width margin of the trace address after screening, that is, the number of left and right trace addresses that meet the intensity screening. The width margin is compared with the margin lower limit. If the width margin is greater than the margin lower limit, the processing result is 1, otherwise 0 or -1. The result is added to the result group. The margin lower limit is 4; Fourth step: data fusion, investigation and screening of trace addresses not in the characteristic range According to the transformation processing of the previous three steps, the peak position trace address characteristic range is that the elements in the result group are all positive. The "and" fusion criterion is adopted. If there is a non-positive number in the result group of a trace address, it is judged to be a non-peak position trace address; Fifth step: cluster analysis, further division of peak positions and non-peak positions, and final peak position determination result Step 5.1 selects the average value of the intensity after background subtraction and the SZAC result as the clustering feature, and uses the Z-score method for data standardization to eliminate the influence of dimension; Step 5.2 calculates the distance between each pair of sample points, and calculates the sum of the distances from each object to the remaining objects. The point corresponding to the maximum value is selected as the first isolated point. After deleting the isolated point, step 5.2 is repeated until I isolated points are removed. I=4 is selected; Step 5.3 sets the number of clusters to 2, and finds the two sample points with the farthest distance after removing the isolated points as the initial cluster centers; Step 5.4 conventional K-means algorithm, traverse the sample, find the nearest cluster center of the sample point, classify the sample into this class, and update the cluster center until the maximum iteration number is reached or the cluster center no longer changes. The final cluster center information is saved; Step 5.5 Calculate the inter-class distance D of two classes b and the sum of intra-class average distance D w , according to D b / D w Determine the number of classes, if the result is greater than 1.25, determine the optimal number of clusters as 2, otherwise determine the number of clusters as 1; Step 5.6 if the number of clusters is 2, the isolated points are classified according to the nearest distance priority principle. If the number of clusters is 1, all sample points are classified into the same class; Step 5.7 if the number of clusters is 2, the class with the maximum average value of the SZAC transformation result is marked as the peak position class, and the other class is marked as the non-peak position class. If the number of clusters is 1, all samples are marked as the peak position class.

Images