A method for multi-variable monitoring of a czochralski silicon single crystal growth process

Abstract

The application provides a multi-variable monitoring method for a Czochralski silicon single crystal growth process, and belongs to the technical field of the Czochralski silicon single crystal growth process. The method comprises the following steps: synchronously collecting and preprocessing time series data of various process variables in the Czochralski silicon single crystal growth process to obtain a multi-variable time series matrix; a multi-scale sliding window strategy is used to extract slow features from the multi-variable time series matrix, and static slow feature matrices and dynamic slow feature matrices are sequentially constructed; a static slow feature vector and a dynamic slow feature vector at each time step are spliced into a fused slow feature vector, all fused slow feature vectors form a fused slow feature matrix, and the fused slow feature matrix is divided into a fused training set and a fused test set; an autoencoder network model is constructed, and the autoencoder network model is trained and tested to obtain a multi-variable monitoring result. The application can realize accurate real-time monitoring of the Czochralski silicon single crystal growth process in a complex environment.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of Czochralski silicon single crystal growth, and in particular to a method for monitoring multiple variables in the process of Czochralski silicon single crystal growth. BACKGROUND

[0002] As a core basic material for integrated circuit manufacturing, the quality of silicon single crystal material directly determines the electrical performance, integration level and long-term reliability of chip devices. With the continuous development of semiconductor manufacturing processes to 7 nanometers and finer nodes, more stringent requirements are placed on the quality of silicon single crystal, particularly in terms of the suppression of internal defects (such as dislocations, vacancy clusters, etc.) in the silicon crystal, the uniformity of impurity (oxygen, carbon, etc.) concentration distribution, and the precise control of the morphology of the silicon crystal. For example, parameters such as the dislocation density and oxygen-carbon content in the silicon crystal can severely affect the breakdown characteristics, resistivity uniformity and service life of chip devices; for another example, geometric abnormalities such as diameter fluctuations and unstable growth interfaces of the silicon crystal can cause problems such as warping and cracking of the silicon wafer during processing, which severely restricts the yield and consistency of chip device manufacturing. Therefore, achieving precise monitoring and dynamic adjustment of key process variables during the growth of silicon single crystal is a core technical challenge to improve the overall performance of silicon crystals and meet the requirements of manufacturing processes.

[0003] In the process of Czochralski silicon single crystal growth, there is a complex coupling relationship between multiple process variables such as pulling speed, heater power, thermal field temperature, silicon crystal diameter and its derived parameters (such as V / G value, i.e. the ratio of growth rate to temperature gradient), which collectively affect the quality of the silicon crystal. For example, changes in pulling speed will change the cooling rate at the front of the growth interface, thereby affecting the generation and proliferation of silicon lattice defects; for example, heater power and thermal field temperature together determine the temperature gradient between the silicon melt and the crystal region, directly affecting the stability of the growth interface and impurity segregation; for another example, dynamic changes in the diameter of the silicon crystal will change the heat dissipation conditions of the exposed part of the silicon crystal, causing redistribution of thermal stress. The interrelation and dynamic interaction between these process variables make the silicon single crystal growth process exhibit strong nonlinearity, time variation and multivariable coupling characteristics, thus placing extremely high requirements on real-time perception and abnormality recognition of the process state.

[0004] As can be seen, achieving precise monitoring and abnormality recognition of the evolution process of process variables is a key prerequisite for ensuring the stability of silicon crystal quality and the reliability of manufacturing processes. However, current silicon single crystal growth process monitoring technologies have the following problems:

[0005] First, in the actual production of industry, the monitoring of the Czochralski silicon single crystal growth process highly depends on artificial experience and single variable threshold judgment. The operator judges whether the diameter of the silicon crystal is stable according to the change of the crystal image, or sets the process parameter range according to experience, so as to judge whether the system is in a normal state. This experience-dependent static judgment mode has obvious defects: first, the response is lagging, and it is difficult to capture rapid changes; second, it lacks multi-variable comprehensive perception ability and cannot effectively analyze the coupling effect between multiple process variables; third, it has very low sensitivity to small fluctuations or potential failure signs of the system, and often cannot be discovered until the problem is obvious or even causes loss. These defects seriously restrict the improvement of the stability and consistency of the quality of silicon crystals;

[0006] Secondly, for the monitoring needs of the Czochralski silicon single crystal growth process, multi-variable statistical process monitoring and data-driven methods are introduced in existing researches, for example, the method based on adaptive canonical variate analysis can realize global process state monitoring, but the sensitivity to local anomalies and small disturbances is low; for example, the method based on deep belief network has strong feature extraction ability, but has problems such as high false alarm rate and weak fault attribution ability. Therefore, it can be seen that the existing technology is generally difficult to achieve precise anomaly monitoring with high sensitivity and low false alarm rate in complex dynamic conditions, and lacks reliable identification ability for small anomalies and gradual drift under multi-variable coupling effect.

[0007] Therefore, it is necessary to propose a scheme to improve one or more problems in the above related technical solutions.

[0008] It should be noted that the information disclosed in the above background section is only used to strengthen the understanding of the background of the present application, and therefore can include information that does not constitute prior art known to those skilled in the art. SUMMARY

[0009] The present application provides a multi-variable monitoring method for the Czochralski silicon single crystal growth process, which comprises the following steps:

[0010] Synchronously collecting time series data of multiple process variables in the Czochralski silicon single crystal growth process, and respectively pre-processing all the time series data of process variables to obtain a multi-variable time series matrix;

[0011] Using a multi-scale sliding window strategy to extract slow features from the multi-variable time series matrix, and sequentially constructing a static slow feature matrix and a dynamic slow feature matrix;

[0012] Respectively splicing the static slow feature vector and the dynamic slow feature vector of each time step into a fusion slow feature vector, all fusion slow feature vectors form a fusion slow feature matrix, and the fusion slow feature matrix is divided into a fusion training set and a fusion test set;

[0013] The auto-encoder network model is constructed, the auto-encoder network model is trained by using the fusion training set, and the trained auto-encoder network model is tested by using the fusion test set;

[0014] All test results are used as multivariate monitoring results of the Czochralski silicon single crystal growth process;

[0015] The static slow feature matrix includes a plurality of static slow feature vectors, and the dynamic slow feature matrix includes a plurality of dynamic slow feature vectors.

[0016] All training samples in the fusion training set are fusion slow feature vectors in normal working conditions; a part of test samples in the fusion test set are fusion slow feature vectors in normal working conditions, and another part of test samples are fusion slow feature vectors in abnormal working conditions.

[0017] Further, the types of process variable time series data include silicon crystal diameter time series data, Value time series data, heater power time series data, silicon crystal pulling speed time series data, single crystal furnace pressure time series data and thermal field temperature time series data;

[0018] The preprocessing includes: respectively standardizing each process variable time series data to obtain the corresponding standardized time series data of each process variable time series data;

[0019] The expression of the standardized time series data is:

[0020] (1)

[0021] Wherein, The standardized time series data is represented by The process variable time series data is represented by The mean of the process variable time series data is represented by The standard deviation of the process variable time series data is represented by

[0022] All standardized time series data are combined according to different time steps to obtain a multivariate time series matrix.

[0023] Further, the slow feature extraction is performed on the multivariate time series matrix by using a multi-scale sliding window strategy, and the steps of sequentially constructing the static slow feature matrix and the dynamic slow feature matrix include:

[0024] Based on the total length of the time series in the multivariate time series matrix and the types of all process variable time series data, various window standard lengths are set, and under each window standard length, the time series is divided into multiple sliding windows connected end to end in sequence. The number of all sliding windows under each window standard length and the sample matrix in each sliding window are obtained respectively.

[0025] The sample matrix is ​​represented as follows: ,in, Indicates the first Under the standard window length, the first The sample matrix within a sliding window express 3D real vector space, Indicates the first Standard window length This indicates the type of all process variable time series data; the number of all sliding windows under the standard window length is represented as: ,in, Indicates the first The number of all sliding windows under a standard window length. Indicates the total length of the time series;

[0026] For each standard window length, a projection matrix is ​​solved using the sample matrix within each sliding window of that standard window length. Then, the sample vector of the last time step of the corresponding sliding window is projected using each projection matrix to obtain the unaligned static slow feature vector corresponding to each sliding window of that standard window length.

[0027] The projection matrix is ​​represented as: ,in, Indicates the first The first type of window standard length The projection matrix corresponding to each sliding window Indicates the first The number of all unaligned static slow features corresponding to a standard window length. express A 3D real vector space; the unaligned static slow eigenvector is represented as: , ,in, Indicates the first The first type of window standard length The unaligned static slow feature vectors corresponding to each sliding window Indicates the first The first type of window standard length The sample vector of the last time step within a sliding window. express a real number vector space;

[0028] define a set of all window standard lengths as , and the maximum value in the set is denoted as , and the minimum value in the set is denoted as ;

[0029] According to the maximum value and the minimum value , determine the start index and the end index of the common interval, the start index of the common interval is denoted as , the end index of the common interval is denoted as , and the length of the common interval is denoted as ; ;

[0030] According to the unaligned static slow feature vectors corresponding to all sliding windows under each window standard length, respectively obtain the end time step of the sliding window corresponding to each unaligned static slow feature vector;

[0031] Screen out all unaligned static slow feature vectors corresponding to the end time step falling within the common interval, and align and sort these unaligned static slow feature vectors according to the respective end time steps, to respectively obtain the aligned static slow feature vector corresponding to each time step under each window standard length;

[0032] Vertically stack all aligned static slow feature vectors under each window standard length according to the time dimension, to respectively obtain the static slow feature matrix corresponding to each window standard length;

[0033] The expression of the static slow feature matrix is as follows:

[0034] (2)

[0035] wherein, denotes the static slow feature matrix corresponding to the i-th window standard length, denotes the number of all aligned static slow features in the static slow feature matrix corresponding to the i-th window standard length, , denotes the number of all window standard lengths, , denotes the number of all window standard lengths, denotes a real number vector space;

[0036] ​​​For each window standard length, the first-order difference of each aligned static slow feature vector is calculated respectively, and the dynamic slow feature vector corresponding to each aligned static slow feature vector is obtained respectively;

[0037] The expression of the dynamic slow feature vector is:

[0038] (3)

[0039] wherein, represents the dynamic slow feature vector of the i-th time step under the j-th window standard length, represents the aligned static slow feature vector of the i-th time step under the j-th window standard length, represents the aligned static slow feature vector of the i-th time step under the j-th window standard length.

[0040] All dynamic slow feature vectors under each window standard length are vertically stacked in the time dimension respectively, and the dynamic slow feature matrix corresponding to each window standard length is obtained respectively.

[0041] The expression of the dynamic slow feature matrix is:

[0042] (4)

[0043] wherein, represents the dynamic slow feature matrix corresponding to the j-th window standard length, represents the number of all dynamic slow features in the dynamic slow feature matrix corresponding to the j-th window standard length, represents an n-dimensional real vector space.

[0044] Further, the expression of the fusion slow feature vector is:

[0045] z i, fused t =[ z i, static t , z i, dynamic ( t )] (5)

[0046] wherein, represents the dynamic slow feature vector of the i-th time step under the j-th window standard length, represents the aligned static slow feature vector of the i-th time step under the j-th window standard length, ​​​​​​​​​The fused slow feature vectors at each time step;

[0047] The expression for the fused slow feature matrix is:

[0048] (6)

[0049] in, Indicates the first The fusion slow feature matrix corresponding to the standard window length. This represents the total length of the time series. express A real vector space.

[0050] Furthermore, the autoencoder network model includes an encoder, a hybrid attention module, and a decoder connected in sequence;

[0051] The encoder includes multiple multi-scale temporal convolutional modules stacked in series. Each multi-scale temporal convolutional module includes a corresponding convolutional kernel and dilation rate. Each multi-scale temporal convolutional module has multiple independent feature channels. Each feature channel outputs a corresponding channel temporal feature vector. All channel temporal feature vectors form the corresponding channel temporal feature matrix.

[0052] The hybrid attention module includes interconnected channel attention submodules and temporal self-attention submodules;

[0053] Among them, the channel attention submodule performs compression and excitation operations in sequence;

[0054] The compression operation includes: performing global flat pooling on all feature channels in each multi-scale temporal convolution module, compressing each channel temporal feature vector into a corresponding channel description vector, and forming a corresponding channel description matrix from all channel description vectors;

[0055] Incentive operations include:

[0056] Each channel description vector is input into the first fully connected layer to obtain the initial channel weight vector corresponding to each channel description vector, and the negative values ​​of each initial channel weight vector are transformed into non-negative values ​​by the first activation function.

[0057] Each initial channel weight vector with non-negative values ​​is input into the second fully connected layer to obtain the final channel weight vector corresponding to each initial channel weight vector. Then, the second activation function maps each final channel weight vector to... Within the interval, the channel weight coefficients corresponding to each final channel weight vector are obtained;

[0058] Each final channel weight vector and its corresponding channel weight coefficient are weighted separately to obtain the corresponding channel weighted feature vector;

[0059] All channel weighting feature vectors constitute a channel weighting feature matrix;

[0060] The time sequence self-attention sub-module projects the channel weighting feature matrix into a query matrix, a key matrix and a value matrix respectively;

[0061] The attention coefficient matrix is calculated using the query matrix and the key matrix;

[0062] The attention coefficient matrix and the value matrix are weighted to generate a time sequence weighting feature matrix;

[0063] The decoder includes a plurality of multi-scale time sequence convolution modules, the order of the concatenation and stacking of all multi-scale time sequence convolution modules in the decoder is opposite to that of all multi-scale time sequence convolution modules in the encoder, and the decoder performs the following operations:

[0064] The global average pooling is performed on all time sequence weighting feature vectors in the time sequence weighting feature matrix in the time dimension to obtain a global feature vector, and the global feature vector is converted into a latent feature vector through a nonlinear transformation;

[0065] The latent feature vector is repeatedly expanded to form a time dimension tensor;

[0066] The time dimension tensor is sequentially subjected to multi-scale convolution processing and mapping processing to generate a reconstructed time vector.

[0067] Further, the expression of the channel time sequence feature vector is:

[0068] (7)

[0069] wherein, denotes the channel time sequence feature vector corresponding to the i-th feature channel of the j-th multi-scale time sequence convolution module, denotes a ReLU activation function, denotes a one-dimensional convolution operation of the i-th feature channel of the j-th multi-scale time sequence convolution module, denotes the channel time sequence feature vector corresponding to the i-th feature channel of the j-th multi-scale time sequence convolution module, denotes the number of all feature channels in the multi-scale time sequence convolution module; The channel time sequence feature matrix is expressed as: , denotes the channel time sequence feature matrix, denotes the channel time sequence feature matrix,

[0070] The channel time sequence feature matrix is expressed as: , denotes the channel time sequence feature matrix, ​​N denotes the number of all multi-scale temporal convolution modules in the encoder, N denotes the number of all multi-scale temporal convolution modules in the encoder, N denotes the number of all multi-scale temporal convolution modules in the encoder,

[0071] The expression of the channel description vector is:

[0072] (8)

[0073] wherein, denotes the channel description vector corresponding to the th feature channel, denotes the channel temporal feature of the th time step, the th feature channel, denotes the arbitrary th feature channel;

[0074] The expression of the channel description matrix is: , denotes the channel description matrix, denotes a -dimensional real vector space;

[0075]

[0076] (9)

[0077] wherein, denotes the channel weighted feature vector of the th time step, the th feature channel, denotes the channel weight coefficient of the th feature channel, s c = α [ ω 2 ∙ ReLU ( ω 1 ∙ z c )] , denotes an activation function, denotes an initial channel weight vector, , denotes a final channel weight vector, , denotes a compression rate, denotes a -dimensional real vector space;

[0078] The channel-weighted feature matrix is denoted as: , The channel-weighted feature matrix is denoted as denotes a d-dimensional real vector space.

[0079] Further, the query matrix is denoted as: , The query matrix is denoted as The channel-weighted feature matrix is denoted as The query projection weight matrix is denoted as , The attention embedding dimension is denoted as denotes a d-dimensional real vector space, denotes the number of all feature channels in the multi-scale temporal convolution module;

[0080] The key matrix is denoted as: , The key matrix is denoted as The key projection weight matrix is denoted as ;

[0081] The value matrix is denoted as: , The value matrix is denoted as The value projection weight matrix is denoted as ;

[0082] The expression of the attention coefficient matrix is:

[0083] (10)

[0084] wherein The attention coefficient matrix is denoted as denotes a normalization exponential function, denotes a transpose;

[0085] The expression of the temporal-weighted feature matrix is:

[0086] (11)

[0087] wherein The temporal-weighted feature matrix is denoted as , denotes the total length of the time series, denotes a d-dimensional real vector space.

[0088] Further, the expression of the latent feature vector is:

[0089] (12)

[0090] wherein, represents a latent feature vector, represents a global feature vector, represents a time step embedding, a time step embedding, represents the total length of the time series;

[0091] The expression of the time dimension tensor is:

[0092] B repect =[ B , … , B ]ϵ R L × d (13)

[0093] wherein, represents a time dimension tensor, represents an attention embedding dimension, represents a real vector space;

[0094] The expression of the reconstructed time vector is:

[0095] (14)

[0096] wherein, represents a reconstructed time vector, represents a mapping function of the decoder, represents the total dimension of the fusion slow feature matrix, represents a real vector space.

[0097] Further, the process of training the autoencoder network model by using the fusion training set and testing the trained autoencoder network model by using the fusion test set comprises:

[0098] In the training process, a loss function is set, the loss function is continuously monitored, the loss function includes a weighted mean square error loss term and an energy constraint term, and the autoencoder network model when the weighted mean square error loss term reaches a minimum value and the energy constraint term also reaches a minimum value within a preset training round is saved as the trained autoencoder network model;

[0099] Each training sample is re-input into the trained autoencoder network model respectively, and the reconstructed static slow feature vector, the reconstructed dynamic slow feature vector and the reconstructed fusion slow feature vector corresponding to each training sample are obtained respectively.

[0100] respectively calculate a static reconstruction error between the reconstructed static slow feature vector and the static slow feature vector corresponding to each training sample; the static reconstruction error is a root mean square error between the reconstructed static slow feature vector and the corresponding static slow feature vector;

[0101] respectively calculate a dynamic reconstruction error between the reconstructed dynamic slow feature vector and the dynamic slow feature vector corresponding to each training sample; the dynamic reconstruction error is a root mean square error between the reconstructed dynamic slow feature vector and the corresponding dynamic slow feature vector;

[0102] respectively calculate a fusion reconstruction error between the reconstructed fusion slow feature vector and the fusion slow feature vector corresponding to each training sample; the fusion reconstruction error is a root mean square error between the reconstructed fusion slow feature vector and the corresponding fusion slow feature vector;

[0103] respectively count the distribution of the global anomaly scores of all static reconstruction errors, all dynamic reconstruction errors and all fusion reconstruction errors, and set the upper threshold and lower threshold of the static anomaly score, the upper threshold and lower threshold of the dynamic anomaly score, and the upper threshold and lower threshold of the fusion anomaly score by using the percentile method respectively;

[0104] According to the distribution of the global anomaly scores of all static reconstruction errors, take the 98th percentile as the upper threshold of the static anomaly score, and take the 2nd percentile as the lower threshold of the static anomaly score;

[0105] According to the distribution of the global anomaly scores of all dynamic reconstruction errors, take the 98th percentile as the upper threshold of the dynamic anomaly score, and take the 2nd percentile as the lower threshold of the dynamic anomaly score;

[0106] According to the distribution of the global anomaly scores of all fusion reconstruction errors, take the 98th percentile as the upper threshold of the fusion anomaly score, and take the 2nd percentile as the lower threshold of the fusion anomaly score;

[0107] respectively weight the fusion anomaly score corresponding to each training sample by the sum variable weight vector to obtain the weighted fusion anomaly score corresponding to each training sample;

[0108] According to the distribution of the global anomaly scores of all weighted fusion anomaly scores, take the 98th percentile as the upper threshold of the weighted fusion anomaly score, and take the 2nd percentile as the lower threshold of the weighted fusion anomaly score;

[0109] respectively calculate the static reconstruction error, the dynamic reconstruction error, the fusion anomaly score and the weighted fusion anomaly score corresponding to each test sample;

[0110] The static reconstruction error corresponding to each test sample is compared with the upper and lower bound thresholds of the static anomaly score corresponding to the training sample at the corresponding time step.

[0111] The dynamic reconstruction error corresponding to each test sample is compared with the upper and lower bound thresholds of the dynamic anomaly score corresponding to the training sample at the corresponding time step.

[0112] The fusion anomaly score corresponding to each test sample is compared with the upper and lower bound thresholds of the fusion anomaly score corresponding to the training sample at the corresponding time step.

[0113] The weighted fusion anomaly score corresponding to each test sample is compared with the upper and lower bound thresholds of the weighted fusion anomaly score corresponding to the training sample at the corresponding time step.

[0114] If the comparison result exceeds any upper or lower bound threshold, the test sample is marked as an anomaly.

[0115] The frequency of anomalies in the static dimension, dynamic dimension, fusion dimension, and weighted fusion dimension of each abnormal test sample were counted separately.

[0116] Based on the anomaly frequencies of all static dimensions, all dynamic dimensions, all fusion dimensions, and all weighted fusion dimensions, the alarm percentage of each anomaly in the test samples is calculated, and the total alarm percentage is the test result.

[0117] Furthermore, the expression for the loss function is:

[0118] (15)

[0119] in, Represents the loss function. This represents the weighted mean square error loss term. , This indicates the number of all training samples in the fused training set. This represents the total length of the time series. This represents the total dimension of the fused slow feature matrix. Indicates the first The importance weights of each dimension Indicates the first The training sample at the th ... The time step, the first Input feature values ​​of each variable dimension In the autoencoder network model, the first... The training sample at the th ... The time step, the first The reconstructed output feature values ​​of each dimension a weight coefficient representing an energy constraint term, represents an energy constraint term, , represents the input feature vector of the i-th training sample at the j-th time step, represents the input feature vector of the i-th training sample at the j-th time step, represents the reconstructed output feature vector of the i-th training sample at the j-th time step in the autoencoder network model, represents the reconstructed output feature vector of the i-th training sample at the j-th time step in the autoencoder network model, represents the reconstructed output feature vector of the i-th training sample at the j-th time step in the autoencoder network model, represents the reconstructed output feature vector of the i-th training sample at the j-th time step in the autoencoder network model, represents the two-norm;

[0120] The reconstruction error includes a static reconstruction error, a dynamic reconstruction error, and a fusion reconstruction error;

[0121] The expression of the reconstruction error is:

[0122] (16)

[0123] wherein, represents the reconstruction error of the i-th sample, the j-th time step, and the k-th dimension, represents the static reconstruction error, represents the static reconstruction error, represents the dynamic reconstruction error, represents the fusion reconstruction error; The expression of the weighted fusion anomaly score is:

[0124]

[0125] (17)

[0126] wherein, represents the fusion training anomaly score of the i-th training sample, represents the i-th variable weight vector. The present application provides a Czochralski silicon single crystal growth process multivariate monitoring method, at least has the following beneficial effects:

[0127] (1) The present application extracts slow features from the multivariate time series matrix by using a multi-scale sliding window strategy, and sequentially constructs a static slow feature matrix and a dynamic slow feature matrix, which can systematically represent the long-term evolution law and local dynamic change of multivariate data;

[0128] (2) The present application can realize the complex nonlinear relationship between multivariate data and the joint modeling of different time scale features by constructing an autoencoder network model containing a multi-scale parallel time series convolution network and a hybrid attention mechanism;

[0129] (2) The present application can realize the complex nonlinear relationship between multivariate data and the joint modeling of different time scale features by constructing an autoencoder network model containing a multi-scale parallel time series convolution network and a hybrid attention mechanism; ​​​

[0130] (3) In the training process of the autoencoder network model, the application adopts a residual reconstruction mechanism to prevent the autoencoder network model from degenerating into a simple input identity mapping, and introduces an energy constraint term in the loss function, thereby improving the monitoring ability for micro anomalies, gradual anomalies and high-frequency anomalies by constraining the difference between high-frequency disturbances and global energy distribution;

[0131] (4) In the training stage, the application obtains the reconstruction error result by using the autoencoder network model, and introduces a variable weight vector to determine the upper threshold and lower threshold of the static anomaly score, dynamic anomaly score, fusion anomaly score and weighted fusion anomaly score of each training sample. In the test stage, the static reconstruction error, dynamic reconstruction error, fusion anomaly score and weighted fusion anomaly score of each test sample are calculated in real time, and the judgment is made based on the percentile method, realizing the multivariate global anomaly monitoring. At the same time, through the anomaly score attribution analysis, the key variables causing the anomaly are identified and located, and finally the precise anomaly alarm with pertinence is generated, providing intelligent decision support for subsequent process control and optimization. BRIEF DESCRIPTION OF DRAWINGS

[0132] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the application and, together with the description, serve to explain the principles of the application. It is apparent that the accompanying drawings in the following description are only some embodiments of the application, and other drawings can be obtained from these drawings without creative labor for those skilled in the art.

[0133] Figure 1 A step schematic diagram of a multivariate monitoring method for a Czochralski silicon single crystal growth process in an exemplary embodiment of the application is shown;

[0134] Figure 2 A flowchart of a multivariate monitoring method for a Czochralski silicon single crystal growth process in an exemplary embodiment of the application is shown;

[0135] Figure 3 A curve graph of the silicon crystal pulling speed time series data after standardization processing in an exemplary embodiment of the application is shown;

[0136] Figure 4 A curve graph of the static slow feature vector corresponding to the silicon crystal pulling speed time series data in an exemplary embodiment of the application is shown;

[0137] Figure 5 A curve graph of the dynamic slow feature vector corresponding to the silicon crystal pulling speed time series data in an exemplary embodiment of the application is shown;

[0138] Figure 6A structural schematic diagram of a self-encoder network model in an example embodiment of the present application is shown.

[0139] Figure 7 A curve diagram of values of a loss function in an example embodiment of the present application is shown.

[0140] Figure 8 A curve diagram of static abnormality scores corresponding to each process variable time series data in an example embodiment of the present application is shown.

[0141] Figure 9 A curve diagram of dynamic abnormality scores corresponding to each process variable time series data in an example embodiment of the present application is shown.

[0142] Figure 10 A curve diagram of fusion abnormality scores corresponding to each process variable time series data in an example embodiment of the present application is shown.

[0143] Figure 11 A bar chart of alarm times of test samples in a test set corresponding to each process variable time series data in an example embodiment of the present application is shown.

[0144] Figure 12 A schematic diagram of weighted fusion abnormality scores and alarm points of all test samples in a test set in an example embodiment of the present application is shown. DETAILED DESCRIPTION

[0145] Example implementations will now be described more fully with reference to the accompanying drawings. Example implementations may, however, be implemented in many different forms and should not be construed as limited to the implementations set forth herein; rather, these implementations are provided so that this disclosure will be thorough and complete, and will fully convey the inventive aspects to those skilled in the art. Features described in the description, structures, or characteristics may be combined in any suitable manner in one or more implementations.

[0146] In addition, the drawings are to be considered in all respects as illustrative and not restrictive; identical reference numerals have been used, where possible, to denote identical or similar features, and thus repetition of the description thereof will be omitted. Some of the blocks in the drawings are functional blocks that do not necessarily have a corresponding physical or logical entity in an actual device. These functional blocks may be implemented in software, or in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices.

[0147] In the following, a multi-variable monitoring method for a Czochralski silicon single crystal growth process in the present example embodiment will be described in more detail.

[0148] The example embodiment provides a method for multivariate monitoring of a Czochralski silicon single crystal growth process, which can include the following steps as shown in Figure 1 and Figure 2

[0149] The embodiment step S101: synchronously collect the time series data of multiple process variables in the Czochralski silicon single crystal growth process, and respectively pre-process all the time series data of process variables to obtain a multivariate time series matrix.

[0150] Further, the types of time series data of process variables include the time series data of silicon single crystal diameter, the time series data of temperature gradient, representing the growth rate, representing the temperature gradient, the ratio of the growth rate to the temperature gradient, the time series data of heater power, the time series data of silicon crystal pulling speed, the time series data of pressure in the single crystal furnace, and the time series data of thermal field temperature.

[0151] Further, each time series data of process variables is pre-processed, and no outlier rejection operation is performed in the pre-processing process, so as to retain potential abnormal condition information for subsequent learning and monitoring of the autoencoder network model. In the embodiment, all the time series data of process variables are shown in Table 1 as follows:

[0152] Table 1 Time series data of multiple process variables in the Czochralski silicon single crystal growth process

[0153]

[0154] Further, Z-score standardization method is used in the pre-processing to standardize the time series data of each process variable, so that these data have comparability, and the standardized time series data corresponding to each process variable is obtained.

[0155] The expression of the standardization processing is:

[0156] (1)

[0157] wherein, represents the standardized time series data, represents the time series data of process variables, represents the mean of the time series data of process variables, represents the standard deviation of the time series data of process variables.

[0158] As shown in Figure 3 ​As shown, the embodiment gives the curve change of the standardization processing of the time series data of the silicon single crystal pulling speed. All process variable time series data can be standardized by formula (1).

[0159] Then, all standardized time series data are combined according to different time steps to obtain a multivariate time series matrix.

[0160] The embodiment step S102: as shown in Figure 4 and Figure 5 , a multi-scale sliding window strategy is used to extract slow features from the multivariate time series matrix, and static slow feature matrix and dynamic slow feature matrix are constructed in turn.

[0161] The multi-scale sliding window strategy can fully utilize the dynamic evolution characteristics of all process variables, so as to extract the long-term trend and change rate of all process variables, and improve the perception ability of the autoencoder network model to complex patterns such as gradual change type anomaly and drift type anomaly. The embodiment step S102 can include the following sub-steps:

[0162] Sub-step S1021: according to the total length of the time series in the multivariate time series matrix and the type of all process variable time series data, a plurality of window standard lengths are set respectively, and the time series is divided into a plurality of sliding windows connected in turn under each window standard length, respectively obtaining the number of all sliding windows under each window standard length and the sample matrix in each sliding window.

[0163] In the embodiment, the window standard length is set to 30 time steps, 60 time steps and 120 time steps respectively.

[0164] Further, the sample matrix in the sliding window is represented as: , wherein, represents the sample matrix in the i-th sliding window under the j-th window standard length, represents a d-dimensional real vector space, represents the j-th window standard length, represents the type of all process variable time series data; the number of all sliding windows under the window standard length is represented as: , wherein, represents the number of all sliding windows under the j-th window standard length, represents the total length of the time series.

[0165] ​​​​In sub-step S1022, for each window standard length, a projection matrix is solved using the sample matrix in each sliding window of the window standard length, respectively, and each sample vector of the last time step of the corresponding sliding window is projected using each projection matrix, respectively, to obtain the unaligned static slow feature vector corresponding to each sliding window of the window standard length, respectively.

[0166] Further, the projection matrix is represented as: wherein, represents the projection matrix corresponding to the i-th sliding window of the j-th window standard length, represents the number of all unaligned static slow features corresponding to the j-th window standard length, represents the j-th window standard length, represents the i-th sliding window of the j-th window standard length, represents the number of all unaligned static slow features corresponding to the j-th window standard length, represents the j-th window standard length, represents the i-th sliding window of the j-th window standard length, represents the sample vector of the last time step in the i-th sliding window of the j-th window standard length, represents the j-th window standard length, represents the i-th sliding window of the j-th window standard length, represents the sample vector of the last time step in the i-th sliding window of the j-th window standard length, represents the j-th window standard length, represents the i-th sliding window of the j-th window standard length, represents the sample vector of the last time step in the i-th sliding window of the j-th window standard length, represents the j-th window standard length, represents the i-th sliding window of the j-th window standard length, represents the j-th window standard length,

[0167] In sub-step S1023, a set of all window standard lengths is defined as , the maximum value in the set is denoted as , and the minimum value in the set is denoted as .

[0168] According to the maximum value and the minimum value , the start index and the end index of the common interval are determined, the start index of the common interval is denoted as , the end index of the common interval is denoted as , and the length of the common interval is denoted as .

[0169] In sub-step S1024, according to the unaligned static slow feature vectors corresponding to all sliding windows of each window standard length, the end time step of the sliding window corresponding to each unaligned static slow feature vector is obtained, respectively.​​​

[0170] Sub-step S1025: Filter out all unaligned static slow feature vectors corresponding to the end time step falling within the common interval, and align and sort these unaligned static slow feature vectors according to their respective end time steps to obtain the aligned static slow feature vectors corresponding to each time step under each standard window length.

[0171] Sub-step S1026: Vertically stack all aligned static slow feature vectors for each standard window length along the time dimension to obtain the static slow feature matrix corresponding to each standard window length.

[0172] Furthermore, the expression for the static slow feature matrix is:

[0173] (2)

[0174] in, Indicates the first The static slow feature matrix corresponding to the standard window length. Indicates the first The number of all aligned static slow features in the static slow feature matrix corresponding to the standard window length. , This indicates the types of standard window lengths. express A real vector space.

[0175] Extracting the static slow feature matrix can characterize the long-term evolution trend of process variable time series data in the time dimension.

[0176] Sub-step S1026: In order to further enhance the ability to characterize the rate of change of variables, for each standard window length, calculate the first difference of each aligned static slow feature vector to obtain the dynamic slow feature vector corresponding to each static slow feature vector.

[0177] Furthermore, the expression for the dynamic slow feature vector is:

[0178] (3)

[0179] in, Indicates the first The first type of window standard length The dynamic slow feature vector at each time step Indicates the first The first type of window standard length Aligned static slow feature vectors at each time step Indicates the first The first type of window standard length aligned static slow feature vector of a time step.

[0180] Sub-step S1027: vertically stack all dynamic slow feature vectors of each window standard length along the time dimension, respectively, to obtain a dynamic slow feature matrix corresponding to each window standard length.

[0181] Further, the expression of the dynamic slow feature matrix is:

[0182] (4)

[0183] wherein, represents the dynamic slow feature matrix corresponding to the th window standard length, represents the number of all dynamic slow features in the dynamic slow feature matrix corresponding to the th window standard length, represents a dimensional real vector space.

[0184] The dynamic slow feature matrix is extracted to depict the local change rate and fluctuation characteristics of the process variable time series data.

[0185] Respectively extracting the static slow feature matrix and the dynamic slow feature can comprehensively retain the global evolution mode and local dynamic details of the process variable time series data.

[0186] In this embodiment, step S103: respectively concatenate the static slow feature vector and the dynamic slow feature vector of each time step into a fusion slow feature vector, all fusion slow feature vectors form a fusion slow feature matrix, and the fusion slow feature matrix is divided into a fusion training set and a fusion test set.

[0187] Further, the expression of the fusion slow feature vector is:

[0188] z i, fused t =[ z i, static t , z i, dynamic ( t )] (5)

[0189] wherein, represents the fusion slow feature vector of the th time step under the th window standard length.

[0190] Further, the expression of the fusion slow feature matrix is:

[0191] (6)

[0192] wherein, represents the fusion slow feature matrix corresponding to the th window standard length, represents a real vector space.

[0193] The fusion slow feature matrix is divided into a fusion training set and a fusion test set. In this embodiment, the proportion of the fusion training set and the fusion test set is 7:3. All training samples in the fusion training set are fusion slow feature vectors in a normal working condition, which are used for training and parameter adjustment of the autoencoder network model. Part of the test samples in the fusion test set are fusion slow feature vectors in a normal working condition, and the other part of the test samples are fusion slow feature vectors in an abnormal working condition, which are used for evaluating the final monitoring effect, so as to ensure that the autoencoder network model covers the normal and abnormal states in the entire process cycle.

[0194] In step S104 of this embodiment, as shown in Figure 6 , an autoencoder network model is constructed, the fusion training set is used to train the autoencoder network model, and the fusion test set is used to test the trained autoencoder network model.

[0195] Further, in this embodiment, a residual-free autoencoder network model integrating a multi-scale parallel time series convolution network and a channel-time hybrid attention mechanism is constructed, which is used to learn the time dependence structure between process variables and the importance features of variable channels, so as to realize abnormal monitoring of all process variables without labels. The training target is to minimize the error between the input sequence and the reconstructed sequence.

[0196] The autoencoder network model includes an encoder, a hybrid attention module and a decoder connected in sequence.

[0197] Firstly, the encoder includes a plurality of multi-scale time series convolution modules connected in series, each multi-scale time series convolution module includes corresponding convolution kernels and dilation rates, each multi-scale time series convolution module has a plurality of independent feature channels, each feature channel outputs a corresponding channel time series feature vector, and all channel time series feature vectors form a corresponding channel time series feature matrix. The encoder is used to extract multi-scale time dependence features in the input fusion slow feature vector.

[0198] Each multi-scale time series convolution module uses different convolution kernels and dilation rates to capture local fluctuations and global trend features in different time dimensions at the same time. All convolution outputs of each multi-scale time series convolution module are aggregated in the same feature channel to form multi-scale fusion feature expression.

[0199] Secondly, the hybrid attention module includes interconnected channel attention submodules and temporal self-attention submodules. This hybrid attention module integrates a channel-temporal hybrid attention mechanism, inserted between the encoder and decoder, to adaptively learn the channel importance and temporal dynamic weight distribution among variables, further enhancing the ability to identify anomalous features of process variables.

[0200] The channel attention submodule executes compression and excitation operations sequentially.

[0201] The compression operation includes: performing global flat pooling on all feature channels in each multi-scale temporal convolution module, compressing each channel temporal feature vector into a corresponding channel description vector, and forming a corresponding channel description matrix from all channel description vectors.

[0202] Incentive operations include:

[0203] Each channel description vector is input into the first fully connected layer to obtain the initial channel weight vector corresponding to each channel description vector, and the negative values ​​of each initial channel weight vector are transformed into non-negative values ​​by the first activation function.

[0204] Each initial channel weight vector with non-negative values ​​is input into the second fully connected layer to obtain the final channel weight vector corresponding to each initial channel weight vector. Then, the second activation function maps each final channel weight vector to... Within the interval, the channel weight coefficients corresponding to each final channel weight vector are obtained;

[0205] Each final channel weight vector and its corresponding channel weight coefficient are weighted separately to obtain the corresponding channel weighted feature vector;

[0206] The weighted eigenvectors of all channels form the channel-weighted eigenma matrix;

[0207] In this embodiment, the expression for the channel timing feature vector is:

[0208] (7)

[0209] in, Indicates the first The first multi-scale temporal convolution module The channel time-series feature vectors corresponding to each feature channel Represents the ReLU activation function. Indicates the first The first multi-scale temporal convolution module One-dimensional convolution operation for each feature channel Indicates the first The first multi-scale temporal convolution module The channel time-series feature vectors corresponding to each feature channel This indicates the number of all feature channels in the multi-scale temporal convolution module.

[0210] The channel time series feature matrix is ​​represented as follows: , Represents the channel time-series feature matrix. This indicates the number of multi-scale temporal convolutional modules in the encoder. This represents the total number of all input samples. This represents the total length of the time series. The input samples here can be either training samples or test samples.

[0211] In this embodiment, the expression for the channel description vector is:

[0212] (8)

[0213] in, Indicates the first The channel description vector corresponding to each feature channel. Indicates the first The time step, the first Channel temporal features of each feature channel Represents any number of Each feature channel.

[0214] The channel description matrix is ​​represented as follows: , This represents the channel description matrix. express A dimensional real vector space.

[0215] In this embodiment, the expression for the channel-weighted feature vector is:

[0216] (9)

[0217] in, Indicates the first The time step, the first Channel-weighted feature vectors of each feature channel Indicates the first Channel weight coefficients for each feature channel. s c = α [ ω 2 ∙ ReLU ( ω 1 ∙ z c )] denotes activation function, denotes an initial channel weight vector, denotes a final channel weight vector, denotes compression rate, denotes dimensional real vector space.

[0218] The channel-weighted feature matrix is denoted as: denotes channel-weighted feature matrix, denotes dimensional real vector space.

[0219] The temporal self-attention sub-module projects the channel-weighted feature matrix into a query matrix, a key matrix and a value matrix, respectively.

[0220] In this embodiment, the query matrix is denoted as: denotes query matrix, denotes channel-weighted feature matrix, denotes query projection weight matrix, denotes attention embedding dimension, denotes dimensional real vector space.

[0221] The key matrix is denoted as: denotes key matrix, denotes key projection weight matrix,

[0222] The value matrix is denoted as: denotes value matrix, denotes value projection weight matrix,

[0223] Using the query matrix and the key matrix, an attention coefficient matrix is calculated.

[0224] In this embodiment, the expression of the attention coefficient matrix is:

[0225] (10)

[0226] wherein, denotes attention coefficient matrix, denotes a normalization exponential function which can convert any real number into a number within the range of 0, 1.​​​​​​​​​​ This indicates transpose.

[0227] The attention coefficient matrix and the value matrix are weighted to generate a time-weighted feature matrix.

[0228] In this embodiment, the expression for the time-weighted feature matrix is:

[0229] (11)

[0230] in, Represents the time-weighted feature matrix. , This represents the total length of the time series. express A dimensional real vector space.

[0231] Finally, in this embodiment, the decoder structure is symmetrically designed with respect to the encoder structure for efficient fidelity preservation and reconstruction of the input sequence. The decoder includes multiple multi-scale temporal convolutional modules, and the cascaded stacking order of all multi-scale temporal convolutional modules in the decoder is the reverse of that in the encoder. The decoder performs the following operations:

[0232] Global average pooling is performed on all time-weighted eigenvectors in the time-weighted feature matrix along the time dimension to obtain global eigenvectors. Then, nonlinear transformation is used to convert the global eigenvectors into latent eigenvectors.

[0233] In this embodiment, the expression for the latent feature vector is:

[0234] (12)

[0235] in, Represents the latent feature vector. Represents the global feature vector. Indicates the first The time-weighted feature vector embedded at each time step This indicates the total length of the time series.

[0236] The latent feature vectors are repeatedly expanded to form a time-dimensional tensor.

[0237] In this embodiment, the expression for the time dimension tensor is:

[0238] B repect =[ B ,…, B ]ϵ R L x d (13)

[0239] in, Represents a time-dimension tensor. Indicates the attention embedding dimension. express A dimensional real vector space.

[0240] The time dimension tensor is sequentially subjected to multi-scale convolution and mapping processing to generate a reconstructed time vector.

[0241] In this embodiment, the expression for reconstructing the time vector is:

[0242] (14)

[0243] in, This represents the reconstructed time vector. The mapping function of the decoder. This represents the total dimension of the fused slow feature matrix. express A real vector space. Here, the total dimension of the fused slow feature matrix refers to the sum of the dimensions of all dynamic slow feature vectors and the dimensions of all static slow feature vectors.

[0244] The process of training an autoencoder network model using a fused training set and testing the trained autoencoder network model using a fused test set includes:

[0245] like Figure 7 As shown, a loss function is set during training and continuously monitored. The loss function includes a weighted mean square error loss term and an energy constraint term. The autoencoder network model that reaches the minimum value of both the weighted mean square error loss term and the energy constraint term within a preset training round is saved as the trained autoencoder network model.

[0246] Furthermore, the expression for the loss function is:

[0247] (15)

[0248] in, Represents the loss function. This represents the weighted mean square error loss term. , This indicates the number of all training samples in the fused training set. This represents the total length of the time series. This represents the total dimension of the fused slow feature matrix. Indicates the first The importance weights of each dimension Indicates the first The training sample at the th ... The time step, the first Input feature values ​​in each dimension In the autoencoder network model, the first... The training sample at the th ... The time step, the first The reconstructed output feature values ​​of each dimension This represents the weighting coefficient of the energy constraint term. Represents the energy constraint term. , Indicates the first The training sample at the th ... The input feature vector at each time step, In the autoencoder network model, the first... The training sample at the th ... The reconstructed output feature vector at each time step This represents the L2 norm.

[0249] Each training sample is re-input into the trained autoencoder network model to obtain the reconstructed static slow feature vector, reconstructed dynamic slow feature vector, and reconstructed fused slow feature vector for each training sample.

[0250] Calculate the static reconstruction error between the reconstructed static slow feature vector and the corresponding static slow feature vector for each training sample; the static reconstruction error is the root mean square error between the reconstructed static slow feature vector and the corresponding static slow feature vector.

[0251] Calculate the dynamic reconstruction error between the reconstructed dynamic slow feature vector and the corresponding dynamic slow feature vector for each training sample; the dynamic reconstruction error is the root mean square error between the reconstructed dynamic slow feature vector and the corresponding dynamic slow feature vector.

[0252] Calculate the fusion reconstruction error between the fused static slow feature vector and the fused slow feature vector for each training sample. The fusion reconstruction error is the root mean square error between the reconstructed fused slow feature vector and the corresponding fused slow feature vector.

[0253] Furthermore, the reconstruction error includes static reconstruction error, dynamic reconstruction error, and fusion reconstruction error. The expression for the reconstruction error is:

[0254] (16)

[0255] in, Indicates the first The sample, the first The time step, the first Reconstruction error in each dimension Indicates static reconstruction error. Indicates dynamic reconstruction error. denotes the fusion reconstruction error. The samples here can include training samples and test samples.

[0256] Further, the expression of the weighted fusion anomaly score is:

[0257] (17)

[0258] wherein, denotes the weighted fusion training anomaly score of the i-th training sample, denotes the i-th variable weight vector. As shown in ,

[0259] , Figure 8 , Figure 9 , Figure 10 , the distributions of the global anomaly scores of all static reconstruction errors, all dynamic reconstruction errors and all fusion reconstruction errors corresponding to the fusion training set are respectively counted, and the upper threshold and the lower threshold of the static anomaly score, the upper threshold and the lower threshold of the dynamic anomaly score, and the upper threshold and the lower threshold of the fusion anomaly score are respectively set by using the percentile method.

[0260] Figure 8 shows the curve change of the static anomaly score corresponding to each process variable time series data. Figure 9 shows the curve change of the dynamic anomaly score corresponding to each process variable time series data. Figure 10 shows the curve change of the fusion anomaly score corresponding to each process variable time series data. Among them, Figure 8 , Figure 9 and Figure 10 , a corresponds to the silicon crystal diameter time series data; Figure 8 , Figure 9 and Figure 10 , b corresponds to the value time series data; Figure 8 , Figure 9 and Figure 10 , c corresponds to the heater power time series data; Figure 8 , Figure 9 and Figure 10 , d corresponds to the silicon crystal pulling speed time series data; Figure 8 , Figure 9 and Figure 10 , e corresponds to the pressure in the single crystal furnace time series data; Figure 8 , Figure 9 and Figure 10 , f corresponds to the hot field temperature time series data. In Figure 8 , Figure 9 andFigure 10 In the curve change of the fusion abnormal score corresponding to each process variable time series data in the figure, the upper threshold and the lower threshold are respectively marked, and in particular Figure 8 The alarm points are respectively marked by black dots in the figure, which can effectively present the possible abnormal situation of each sample in each process variable time series data. According to the distribution of the global abnormal score of all static reconstruction errors, the 98% quantile is taken as the upper threshold of the static abnormal score, and the 2% quantile is taken as the lower threshold of the static abnormal score.

[0261] According to the distribution of the global abnormal score of all dynamic reconstruction errors, the 98% quantile is taken as the upper threshold of the dynamic abnormal score, and the 2% quantile is taken as the lower threshold of the dynamic abnormal score.

[0262] According to the distribution of the global abnormal score of all fusion reconstruction errors, the 98% quantile is taken as the upper threshold of the fusion abnormal score, and the 2% quantile is taken as the lower threshold of the fusion abnormal score.

[0263] The fusion abnormal score corresponding to each training sample and the variable weight vector are respectively weighted to obtain the weighted fusion abnormal score corresponding to each training sample.

[0264] According to the distribution of the global abnormal score of all weighted fusion abnormal scores, the 98% quantile is taken as the upper threshold of the weighted fusion abnormal score, and the 2% quantile is taken as the lower threshold of the weighted fusion abnormal score.

[0265] The static reconstruction error, the dynamic reconstruction error, the fusion abnormal score and the weighted fusion abnormal score corresponding to each test sample are respectively calculated.

[0266] The static reconstruction error corresponding to each test sample is compared with the upper threshold and the lower threshold of the static abnormal score corresponding to the training sample at the corresponding time step.

[0267] The dynamic reconstruction error corresponding to each test sample is compared with the upper threshold and the lower threshold of the dynamic abnormal score corresponding to the training sample at the corresponding time step.

[0268] The fusion abnormal score corresponding to each test sample is compared with the upper threshold and the lower threshold of the fusion abnormal score corresponding to the training sample at the corresponding time step.

[0269] The weighted fusion abnormal score corresponding to each test sample is compared with the upper threshold and the lower threshold of the weighted fusion abnormal score corresponding to the training sample at the corresponding time step.

[0270] If the comparison result exceeds any upper threshold or lower threshold, the test sample is marked as abnormal.

[0271] The static dimension abnormality frequency, dynamic dimension abnormality frequency, fusion dimension abnormality frequency and weighted fusion dimension abnormality frequency of each abnormality of the test sample are counted respectively.

[0272] As shown in Figure 9 , according to all static dimension abnormality frequencies, all dynamic dimension abnormality frequencies, all fusion dimension abnormality frequencies and all weighted fusion dimension abnormality frequencies, the alarm proportion of each abnormality of the test sample is calculated respectively, and all alarm proportions are the test results.

[0273] Figure 10 The frequency histogram of the time series data of different process variables being judged as abnormal is counted, which reveals the activity degree of each process variable time series data in the whole process.

[0274] The step S105 of the embodiment is as shown in Figure 10 , all test results are taken as the multivariate monitoring results of the Cz-Si single crystal growth process.

[0275] Figure 11 The time sequence curve of the weighted fusion abnormal score and the quantile threshold are shown, and the circle points represent the final alarm results. The dotted line in the figure represents the upper threshold of the 98% quantile, and the dotted line represents the lower threshold of the 2% quantile, which is used to determine the alarm threshold, indicating that the autoencoder network model adopts the quantile threshold strategy based on the statistical characteristics of the training set, which can reduce the influence of manually setting the threshold, so as to obtain more stable and objective abnormal monitoring results.

[0276] In addition, Figure 11 the fluctuation of the weighted fusion abnormal score of different samples in the time evolution process is clearly shown. When the weighted fusion abnormal score of a certain sample exceeds the upper threshold or is lower than the lower threshold, it will be judged as abnormal by the autoencoder network model, and is marked with a circle point in the figure, which indicates that the time step or sample has abnormal behavior in the whole monitoring feature and needs further attention, analysis and attribution.

[0277] The mean, standard deviation, maximum, skewness, kurtosis, and abnormal point proportion of the weighted fusion abnormal score are also given in Table 2 below, which can help researchers or engineers to fully understand the overall distribution of the global abnormal score output by the autoencoder network model. When the statistical indicators show abnormal or unreasonable characteristics, for example, when the abnormal point proportion is too high or the kurtosis is too large, the parameters, feature extraction, threshold setting and other processes of the autoencoder network model need to be further adjusted and optimized.

[0278] Table 2: Statistics of weighted fusion abnormal score

[0279]

[0280] Figure 8 Figure 12 Figure 12 Figure 12 And Table 2 fully shows the monitoring performance of the multi-variable monitoring method for the Czochralski silicon single crystal growth process by comprehensively showing the local and global, single variable and multi-variable, etc., which can provide strong support for multi-variable anomaly identification and variable attribution of the Czochralski silicon single crystal growth process.

[0281] In addition, the terms "first", "second", "third", etc. are used only for descriptive purposes and are not to be construed as indicating or implying relative importance or an ordered ranking of the indicated technical features. Thus, a feature defined with "first", "second", etc. can explicitly or implicitly include one or more of the features. In the description of the embodiments of the present application, the meaning of "a plurality of" is two or more, unless otherwise explicitly and specifically limited.

[0282] In the description of the present application, the description of the terms "one embodiment", "some embodiments", "an example", "a specific example" or "some examples" means that the specific features, structures, materials or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present application. In the present specification, the illustrative description of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art can combine and combine different embodiments or examples described in the present specification.

[0283] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto, and any skilled person in the art can easily think of various equivalent modifications or replacements within the technical scope disclosed by the present application, and these modifications or replacements should be covered within the protection scope of the present application.

[0284] Other embodiments of the present application will be readily apparent to those skilled in the art upon considering the disclosure herein, the principles and applications of which are intended to be used in the broadest sense. The present application is intended to cover any variations, uses or adaptive changes of the present application following the general principles thereof and including those art-known or customary technical means not disclosed in the present application.

Claims

1. A method of multivariate monitoring of a Czochralski silicon single crystal growth process, characterized by, The method comprises the following steps: Synchronously collecting time series data of multiple process variables in the growth process of a Czochralski silicon single crystal, and respectively preprocessing all the time series data of the process variables to obtain a multivariate time series matrix; Using a multi-scale sliding window strategy, slow feature extraction is performed on the multivariate time series matrix to sequentially construct a static slow feature matrix and a dynamic slow feature matrix, comprising: According to the total length of the time series in the multivariate time series matrix and the types of all the time series data of the process variables, a plurality of window standard lengths are respectively set, and under each window standard length, the time series is divided into a plurality of sequentially connected sliding windows to respectively obtain the number of all the sliding windows under each window standard length and the sample matrix in each sliding window; The sample matrix is ​​represented as follows: ,in, Indicates the first Under the standard window length, the first The sample matrix within a sliding window express 3D real vector space, Indicates the first Standard window length This indicates the type of all process variable time series data; the number of all sliding windows under the standard window length is represented as: ,in, Indicates the first The number of all sliding windows under a standard window length. Indicates the total length of the time series; For each window standard length, a projection matrix is respectively solved using the sample matrix in each sliding window under the window standard length, and a sample vector at the last time step of the corresponding sliding window is respectively projected using each projection matrix to respectively obtain an unaligned static slow feature vector corresponding to each sliding window under the window standard length; The projection matrix is ​​represented as: ,in, Indicates the first The first type of window standard length The projection matrix corresponding to each sliding window Indicates the first The number of all unaligned static slow features corresponding to a standard window length. express A real vector space; the unaligned static slow feature vector is represented as: , ,in, Indicates the first The first type of window standard length The unaligned static slow feature vectors corresponding to each sliding window Indicates the first The first type of window standard length The sample vector of the last time step within a sliding window. express 3D real vector space; define a set of all the window standard lengths as the maximum value in the set is denoted as the minimum value in the set is denoted as ; According to the maximum value and the minimum value , a start index and an end index of a common interval are determined, the start index of the common interval is represented by , , the end index of the common interval is represented by , , and a length of the common interval is represented by , ; According to the unaligned static slow feature vectors corresponding to all the sliding windows under each window standard length, the end time step of the sliding window corresponding to each unaligned static slow feature vector is respectively obtained; All unaligned static slow feature vectors corresponding to the end time steps falling within the common interval are screened out, and these unaligned static slow feature vectors are aligned and sorted according to their respective end time steps to respectively obtain an aligned static slow feature vector corresponding to each time step under each window standard length; All the aligned static slow feature vectors under each window standard length are respectively vertically stacked in the time dimension to respectively obtain a static slow feature matrix corresponding to each window standard length; The expression of the static slow feature matrix is: (2) wherein, represents a static slow feature matrix corresponding to a window standard length of the i-th window, represents a number of all aligned static slow features in the static slow feature matrix corresponding to the window standard length of the i-th window, , represents a kind of all window standard lengths, represents a d-dimensional real vector space;​​ For each window standard length, a first-order difference of each aligned static slow feature vector is respectively calculated to respectively obtain a dynamic slow feature vector corresponding to each aligned static slow feature vector; The expression of the dynamic slow feature vector is: (3) in, Indicates the first The first type of window standard length The dynamic slow feature vector at each time step Indicates the first The first type of window standard length Aligned static slow feature vectors at each time step Indicates the first The first type of window standard length Aligned static slow feature vectors at each time step; All the dynamic slow feature vectors under each window standard length are respectively vertically stacked in the time dimension to respectively obtain a dynamic slow feature matrix corresponding to each window standard length; The expression of the dynamic slow feature matrix is: (4) wherein, represents a dynamic slow feature matrix corresponding to a window standard length, represents a dynamic slow feature matrix corresponding to a window standard length, represents a number of all dynamic slow features in a dynamic slow feature matrix corresponding to a window standard length, represents a d-dimensional real vector space;​ The static slow feature vector and the dynamic slow feature vector of each time step are respectively spliced into a fusion slow feature vector, all the fusion slow feature vectors form a fusion slow feature matrix, and the fusion slow feature matrix is divided into a fusion training set and a fusion test set; An autoencoder network model is constructed, the fusion training set is used to train the autoencoder network model, and the fusion test set is used to test the trained autoencoder network model. all the test results are taken as multivariate monitoring results of the Czochralski silicon single crystal growth process; The static slow feature matrix comprises a plurality of static slow feature vectors, and the dynamic slow feature matrix comprises a plurality of dynamic slow feature vectors. All training samples in the fusion training set are the fusion slow feature vectors in the normal working condition; and part of the test samples in the fusion test set are the fusion slow feature vectors in the normal working condition, and the other part of the test samples are the fusion slow feature vectors in the abnormal working condition.

2. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 1, characterized in that, types of the process variable time series data include silicon crystal diameter time series data, value time series data, heater power time series data, silicon crystal pulling speed time series data, pressure in a single crystal furnace time series data, and thermal field temperature time series data; The preprocessing comprises: respectively performing standardization processing on each process variable time series data to obtain standardized time series data corresponding to each process variable time series data; The expression of the standardized time series data is: (1) wherein, denotes the standardized time series data, denotes the process variable time series data, denotes the mean of the process variable time series data, denotes the standard deviation of the process variable time series data; All the standardized time series data are combined according to different time steps to obtain the multivariate time series matrix.

3. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 1, characterized in that, The expression of the fusion slow feature vector is: (5) wherein, represents the fusion slow feature vector at the i-th time step under the j-th window criterion length; represents the fusion slow feature vector at the i-th time step under the j-th window criterion length; represents the fusion slow feature vector at the i-th time step under the j-th window criterion length; The expression of the fusion slow feature matrix is: (6) wherein, represents a fusion slow feature matrix corresponding to a window standard length, represents a d-dimensional real vector space.​ 4. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 1, characterized in that, The autoencoder network model comprises an encoder, a hybrid attention module and a decoder connected in sequence. The encoder comprises a plurality of multi-scale time series convolution modules connected in series and stacked, each multi-scale time series convolution module comprises a corresponding convolution kernel and a dilation rate, each multi-scale time series convolution module has a plurality of independent feature channels, each feature channel outputs a corresponding channel time series feature vector, and all channel time series feature vectors form a corresponding channel time series feature matrix. The hybrid attention module comprises a channel attention submodule and a time series self-attention submodule connected with each other. The channel attention submodule sequentially performs a compression operation and an excitation operation. The compression operation comprises: respectively performing global average pooling on all feature channels in each multi-scale time series convolution module, and respectively compressing each channel time series feature vector into a corresponding channel description vector, and all channel description vectors form a corresponding channel description matrix. The excitation operation comprises: Each channel description vector is input into a first full connection layer to obtain an initial channel weight vector corresponding to each channel description vector, and a negative value of each initial channel weight vector is converted into a non-negative value by a first activation function. respectively inputting all the initial channel weight vectors with non-negative values into a second fully connected layer to obtain a final channel weight vector corresponding to each of the initial channel weight vectors, and respectively mapping each of the final channel weight vectors to obtaining a channel weight coefficient corresponding to each of the final channel weight vectors within an interval. Each final channel weight vector and a corresponding channel weight coefficient are weighted to obtain a corresponding channel weighted feature vector. All channel weighted feature vectors form a channel weighted feature matrix. The time series self-attention submodule projects the channel weighted feature matrix into a query matrix, a key matrix and a value matrix, respectively. An attention coefficient matrix is calculated by using the query matrix and the key matrix. The attention coefficient matrix and the value matrix are weighted to generate a time series weighted feature matrix. The decoder comprises a plurality of multi-scale time series convolution modules, the serial and stacked order of all multi-scale time series convolution modules in the decoder is opposite to that of all multi-scale time series convolution modules in the encoder, and the decoder performs the following operations: performing global average pooling on all time-series weighted feature vectors in the time-series weighted feature matrix in a time dimension to obtain a global feature vector, and converting the global feature vector into a latent feature vector through a nonlinear transformation; repeatedly expanding the latent feature vector to form a time dimension tensor; performing multi-scale convolution processing and mapping processing on the time dimension tensor in sequence to generate a reconstructed time vector.

5. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 4, characterized in that, An expression of the channel time-series feature vector is: (7) wherein, represents a channel time sequence feature vector corresponding to the i-th feature channel of the j-th multi-scale time sequence convolution module, represents a ReLU activation function, represents a one-dimensional convolution operation of the i-th feature channel of the j-th multi-scale time sequence convolution module, represents a channel time sequence feature vector corresponding to the i-th feature channel of the j-th multi-scale time sequence convolution module, represents the number of all feature channels in the multi-scale time sequence convolution module.​​​​​​ The channel timing feature matrix is represented as: , represents a channel timing feature matrix, represents the number of all multi-scale timing convolution modules in the encoder, represents the number of all input samples, represents the total length of the time series; An expression of the channel description vector is: (8) wherein, represents a channel description vector corresponding to the th feature channel, represents a channel temporal feature of the th time step, the th feature channel, represents an arbitrary th feature channel; The channel description matrix is represented as: , The channel description matrix is represented as The channel description matrix is represented as a d-dimensional real vector space; An expression of the channel weighted feature vector is: (9) wherein, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step, represents the channel-weighted feature vector of the i-th feature channel at the t-th time step,​​​​​​​ The channel-weighted feature matrix is represented as: , represents the channel-weighted feature matrix, represents a d-dimensional real vector space.

6. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 4, characterized in that, The query matrix is represented as: , represents the query matrix, represents the channel-weighted feature matrix, represents the query projection weight matrix, , represents the attention embedding dimension, represents a real-valued vector space, represents the number of all feature channels in the multi-scale time convolution module; The key matrix is represented as: , The key matrix is represented as, The key projection weight matrix is represented as, ; The value matrix is represented as: , The value matrix is represented as, The value projection weight matrix is represented as, ; An expression of the attention coefficient matrix is: (10) wherein, denotes the attention coefficient matrix, denotes the normalized exponential function, denotes the transpose; An expression of the time-series weighted feature matrix is: (11) wherein, denotes a time-sequential weighting feature matrix, , denotes the total length of the time series, denotes dimensional real vector space.

7. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 4, characterized in that, An expression of the latent feature vector is: (12) in, Represents the latent feature vector. Represents the global feature vector. Indicates the first The time-weighted feature vector embedded at each time step Indicates the total length of the time series; An expression of the time dimension tensor is: (13) wherein, denotes a time dimension tensor, denotes an attention embedding dimension, denotes dimensional real vector space; An expression of the reconstructed time vector is: (14) wherein, denotes a reconstruction time vector, denotes a mapping function of the decoder, denotes the total dimension of the fused slow feature matrix, denotes dimensional real vector space.

8. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 1, characterized in that, The process of training the autoencoder network model using the fusion training set and testing the trained autoencoder network model using the fusion test set includes: setting a loss function during the training process, continuously monitoring the loss function, the loss function including a weighted mean square error loss term and an energy constraint term, and saving the autoencoder network model when the weighted mean square error loss term reaches a minimum value and the energy constraint term also reaches a minimum value within a preset training round as the trained autoencoder network model; re-inputting each training sample into the trained autoencoder network model to obtain a reconstructed static slow feature vector, a reconstructed dynamic slow feature vector, and a reconstructed fusion slow feature vector corresponding to each training sample, respectively; calculating a static reconstruction error between the reconstructed static slow feature vector and the static slow feature vector corresponding to each training sample, respectively; the static reconstruction error is a root mean square error between the reconstructed static slow feature vector and the corresponding static slow feature vector; calculating a dynamic reconstruction error between the reconstructed dynamic slow feature vector and the dynamic slow feature vector corresponding to each training sample, respectively; the dynamic reconstruction error is a root mean square error between the reconstructed dynamic slow feature vector and the corresponding dynamic slow feature vector; calculating a fusion reconstruction error between the reconstructed fusion slow feature vector and the fusion slow feature vector corresponding to each training sample, respectively; the fusion reconstruction error is a root mean square error between the reconstructed fusion slow feature vector and the corresponding fusion slow feature vector; statistically analyzing the distribution of global anomaly scores of all static reconstruction errors, all dynamic reconstruction errors, and all fusion reconstruction errors, and setting upper and lower threshold values of static anomaly scores, upper and lower threshold values of dynamic anomaly scores, and upper and lower threshold values of fusion anomaly scores using the percentile method, respectively; taking the 98th percentile as the upper threshold value of the static anomaly score and the 2nd percentile as the lower threshold value of the static anomaly score according to the distribution of global anomaly scores of all static reconstruction errors; According to the distribution of the global anomaly scores of all the dynamic reconstruction errors, the 98th percentile is taken as the upper threshold of the dynamic anomaly score, and the 2nd percentile is taken as the lower threshold of the dynamic anomaly score; According to the distribution of the global anomaly scores of all the fusion reconstruction errors, the 98th percentile is taken as the upper threshold of the fusion anomaly score, and the 2nd percentile is taken as the lower threshold of the fusion anomaly score; The fusion anomaly score and the variable weight vector corresponding to each training sample are weighted respectively to obtain a weighted fusion anomaly score corresponding to each training sample; According to the distribution of the global anomaly scores of all the weighted fusion anomaly scores, the 98th percentile is taken as the upper threshold of the weighted fusion anomaly score, and the 2nd percentile is taken as the lower threshold of the weighted fusion anomaly score; The static reconstruction error, the dynamic reconstruction error, the fusion anomaly score and the weighted fusion anomaly score corresponding to each test sample are calculated respectively; The static reconstruction error corresponding to each test sample is compared with the upper threshold and the lower threshold of the static anomaly score corresponding to the training sample corresponding to the time step; The dynamic reconstruction error corresponding to each test sample is compared with the upper threshold and the lower threshold of the dynamic anomaly score corresponding to the training sample corresponding to the time step; The fusion anomaly score corresponding to each test sample is compared with the upper threshold and the lower threshold of the fusion anomaly score corresponding to the training sample corresponding to the time step; The weighted fusion anomaly score corresponding to each test sample is compared with the upper threshold and the lower threshold of the weighted fusion anomaly score corresponding to the training sample corresponding to the time step; If the comparison result exceeds any upper threshold or lower threshold, the test sample is marked as abnormal; The static dimension abnormal frequency, the dynamic dimension abnormal frequency, the fusion dimension abnormal frequency and the weighted fusion dimension abnormal frequency of the test sample of each anomaly are counted respectively; According to all the static dimension abnormal frequencies, all the dynamic dimension abnormal frequencies, all the fusion dimension abnormal frequencies and all the weighted fusion dimension abnormal frequencies, the alarm proportion of the test sample of each anomaly is counted respectively, and all the alarm proportions are the test results.

9. The multivariate monitoring method for the Czochralski silicon single crystal growth process according to claim 8, characterized in that, The expression of the loss function is: (15) wherein, represents a loss function, represents a weighted mean squared error loss term, , represents the number of all training samples in the fusion training set, represents the total length of the time series, represents the total dimension of the fusion slow feature matrix, represents the importance weight of the th dimension, represents the input feature value of the th training sample at the th time step, the th dimension, represents the reconstructed output feature value of the th training sample in the autoencoder network model at the th time step, the th dimension, represents the weight coefficient of the energy constraint term, represents the energy constraint term, , represents the input feature vector of the th training sample at the th time step, represents the reconstructed output feature vector of the th training sample in the autoencoder network model at the th time step, represents a two-norm; The reconstruction error includes the static reconstruction error, the dynamic reconstruction error and the fusion reconstruction error; The expression of the reconstruction error is: (16) wherein, represents the reconstruction error of the th sample, the th time step, the th dimension, represents the static reconstruction error, represents the dynamic reconstruction error, represents the fused reconstruction error; The expression of the weighted fusion anomaly score is: (17) in, Indicates the first The weighted fusion anomaly score of each training sample. Indicates the first A vector of variable weights.

Images