Explanatable data-driven semiconductor silicon single crystal growth diameter prediction method

Abstract

The invention discloses an interpretable data-driven semiconductor silicon single crystal growth diameter prediction method, which comprises the following steps of: firstly, preprocessing experimental data, and calculating a maximum mutual information coefficient (MIC) to obtain an auxiliary characteristic variable related to a crystal diameter variable; estimating a time delay order between the auxiliary characteristic variable and the crystal diameter, and performing time sequence alignment of the data sequence; and establishing a crystal diameter deep learning prediction model TCN-LSTM and the like based on the time convolution network TCN and the long and short term memory network LSTM, and finally obtaining a crystal diameter prediction result. The problems that in a crystal diameter prediction method in the prior art, the number of mechanism model assumed conditions is large, the prediction result is inaccurate, and a data driving model does not have interpretability are solved.

Description

Technical Field

[0001] The invention belongs to the technical field of silicon single crystal Czochralski method, and in particular relates to an interpretable data-driven semiconductor silicon single crystal growth diameter prediction method. Background Art

[0002] Silicon single crystal is the most important semiconductor material at present. More than 90% of integrated circuit products are made of silicon single crystal as the base material. As a leading industry, the quality of semiconductor silicon single crystal material preparation industry determines the performance of integrated circuit chips processed and manufactured later. With the rapid development of Internet of Things, artificial intelligence, and 5G communication technology, in order to improve the integration and intelligence of various smart terminal devices, mobile wearable devices, and electronic communication equipment, high-quality silicon single crystal preparation is required. To this end, how to improve the process control technology of growing silicon single crystals by Czochralski method is particularly important. The optimization control of crystal diameter during the growth of Czochralski silicon single crystals has a direct and close relationship with the final crystal quality. However, with the increase of silicon single crystal diameter, the size of the thermal field structure of the single crystal furnace also increases, and the heat transfer process in crystal growth becomes more complex and changeable, and presents process characteristics such as nonlinearity, large hysteresis, and uncertainty. These characteristics make the traditional method of estimating crystal diameter based on experience increasingly inaccurate, making it increasingly difficult to achieve accurate crystal diameter control. Therefore, in order to solve the problem of crystal diameter detection and control, it is necessary to establish a more accurate crystal growth process diameter prediction model to achieve high-precision diameter control goals.

[0003] Although there are many simplified mathematical mechanism models based on various physical assumptions in theory to describe the change behavior of crystal diameter, due to the high complexity of the growth process of CZ silicon single crystals, such as multi-field (thermal field, magnetic field, flow field) coupling, multi-phase (solid phase, liquid phase, gas phase) action, dynamic time-varying thermal field, high uncertainty of the model and large time lag, the prediction accuracy of the mechanism model is low, and the established mechanism model is difficult to apply in practice due to too many assumptions. However, compared with the mechanism modeling method, data-driven modeling shows its superiority. The data-driven silicon single crystal diameter modeling method does not require detailed and accurate prior knowledge of crystal growth. By collecting data during the crystal growth process, the offline modeling and online prediction functions of the crystal diameter parameters can be realized. At present, with the widespread application of deep learning technology in industrial processes, it provides a good data-driven modeling strategy for solving modeling problems such as nonlinearity and time-varying in the crystal growth process. In addition, considering that the data-driven modeling method based on deep learning is a black box model and does not have the function of explaining the relationship between input variables and predicted variables, it is critical to develop a deep learning crystal diameter prediction model with interpretable functions. Therefore, the present invention designs an interpretable data-driven CZ silicon single crystal diameter prediction method, which can not only solve the shortcomings of the traditional mechanism model that cannot be predicted, but also explain the degree of effect of different operating variables on the crystal diameter prediction results. In short, this inventive method can effectively guide the actual CZ silicon single crystal growth site and assist engineering technicians and crystal growth control systems to make accurate decisions. Summary of the invention

[0004] The purpose of the present invention is to provide an interpretable data-driven semiconductor silicon single crystal growth diameter prediction method, which solves the problems in the prior art of crystal diameter prediction methods that the mechanism model assumes many conditions and the prediction results are inaccurate, and the data-driven model is not interpretable.

[0005] The technical solution adopted by the present invention is an interpretable data-driven semiconductor silicon single crystal growth diameter prediction method, which is specifically implemented according to the following steps:

[0006] Step 1, preprocessing the experimental data, and obtaining the auxiliary characteristic variables related to the crystal diameter variable by calculating the maximum mutual information coefficient MIC;

[0007] Step 2, estimating the time delay order between the auxiliary characteristic variable and the crystal diameter, and performing time alignment of the data sequence;

[0008] Step 3: Establish a crystal diameter deep learning prediction model TCN-LSTM based on the temporal convolutional network TCN and the long short-term memory network LSTM;

[0009] Step 4: Train the TCN-LSTM network model and use the Grey Wolf Optimization Algorithm (GWO) to optimize the network model hyperparameters.

[0010] Step 5: Obtain the crystal diameter prediction result based on the trained crystal diameter deep learning prediction model TCN-LSTM.

[0011] The present invention is also characterized in that:

[0012] Step 1 is implemented according to the following steps:

[0013] Step 1.1: Use the triple standard deviation principle to eliminate outliers:

[0014] The three-times standard deviation principle is used to eliminate outliers in the crystal diameter observation data. First, the mean μ and standard deviation σ of the crystal diameter data y are calculated according to the following formula:

[0015]

[0016]

[0017] In the formula, y i is the i-th crystal diameter data sample, and n is the total number of sample data;

[0018] Determine the threshold of outliers. The threshold range of outliers is determined by adding or subtracting three times the standard deviation from the mean, that is, the lower limit is μ-3σ and the upper limit is μ+3σ. Check whether there are values ​​in the data set that are less than the lower threshold or greater than the upper threshold. Values ​​that are less than the lower threshold or greater than the upper threshold are considered outliers and are deleted from the data set.

[0019] Step 1.2: Use the minimum-maximum scaling method to normalize the collected auxiliary variables. When normalizing the crystal diameter data, the specific steps are as follows:

[0020] Find the maximum value y in the crystal diameter data max Minimum value y min , linearly scale the original crystal diameter data to between 0 and 1, as follows:

[0021]

[0022] In the formula, y nor represents the normalized data;

[0023] The normalized data is used as training data for the crystal diameter prediction model;

[0024] Step 1.3: Use the maximum mutual information coefficient (MIC) to screen auxiliary variables with a greater correlation with crystal diameter. The specific steps are as follows:

[0025] The mutual information used in the MIC principle is expressed by the following equation:

[0026]

[0027] In the formula, x represents the operating variable related to the silicon single crystal, y represents the actual crystal diameter variable, p(x,y) represents the joint probability density distribution function of x and y, and p(x) and p(y) distributions represent the marginal probability density functions of x and y. The calculation formula of MIC is as follows:

[0028]

[0029] In the formula, a represents the number of grid divisions in the x direction, b represents the number of grid divisions in the y direction, and B is a set variable. B is the 0.6th power of the data capacity of the crystal diameter variable involved in calculating MIC. The range of MIC is [0, 1]. The larger the value of MIC, the stronger the correlation between the two variables.

[0030] Step 2 is implemented as follows:

[0031] N features are screened out through MIC, and the original time series data matrix X is obtained by sampling M times on the time axis:

[0032]

[0033] The subscript i of the element in the original data matrix represents the process variable, and the superscript t represents the sampling time, where x0 represents the time series data of the zero-delay reference variable, that is, the target variable, and x i represents the time series data of the i-th process variable;

[0034] Assume that the multiple correlation delay parameters τ between the target variable and each process variable are represented by the delay sequence Γ, that is:

[0035] Γ=[τ1,τ2,...,τ i ,...,τ N ]

[0036] The sampling period of each process variable in the original time series data matrix is ​​T, then the delay parameter of the i-th process variable relative to the target variable needs to satisfy:

[0037] Γ=DT=[d1T,d2T,...,d i T,...,d N T]

[0038] Where D is the time base sequence, d i is the time base of the ith process variable, d iis a dimensionless integer, which ensures that the multiple correlation delay parameters of the process variables are integer multiples of the sampling period, ensuring the feasibility of the process variables in the time series reconstruction;

[0039] Reconstruct the time series correlation matrix composed of process variables and target variables under different time base sequences in the space-time dimension. The reconstruction process is as follows: First, start from the zero-delay reference variable Select F data with sampling frequency T from the tth time to form a time series

[0040]

[0041] FT≥max(τ1,τ2,…,τ i ,…,τ N )

[0042] In the formula, F must be greater than the delay of each auxiliary variable;

[0043] Relative to the target variable, the process variable x i The delay coefficient relative to the zero-delay reference variable is τ i , so take the process variable x i From t-τ i Generate a new time series with F data starting at time That is, from td i F consecutive sampling data starting from T:

[0044]

[0045] The temporal correlation matrix reconstructed in time and space according to the time delay parameters of each process variable

[0046]

[0047] Let the target variable in the reconstructed time series correlation matrix be the parent sequence, and the multivariate process variable whose delay parameter needs to be estimated be the subsequence. First, calculate the difference between the corresponding elements of the parent sequence and the subsequence, and express it with matrix A:

[0048]

[0049] In the formula,

[0050] Let a=min{a ji}, b = max{a ji}, then the grey correlation coefficient matrix is ​​expressed by γ:

[0051]

[0052] In the formula, ρ is the resolution coefficient. The smaller ρ is, the greater the resolution is. When ρ is 0.5, the grey correlation coefficient matrix can better reflect the actual correlation between multivariate variables.

[0053] The elements in the grey correlation coefficient matrix are averaged by column, that is, the grey correlation degree of the subsequence relative to the parent sequence is defined as:

[0054] r=[r1,r2,...,r i ,...,r N ]

[0055] In the formula,

[0056] The sum of grey correlation R is used to describe the overall correlation between multivariate process variables and target variables, where R is expressed as:

[0057]

[0058] Step 3 is implemented as follows:

[0059] Based on the historical data of silicon single crystal growth, a crystal diameter prediction model TCN-LSTM is established by combining the temporal convolutional network TCN and the long short-term memory network LSTM. The hyperparameters are optimized using the GWO optimization algorithm. The TCN-LSTM model is used to learn the nonlinear relationship between the auxiliary variables and the crystal diameter. The nonlinear relationship is described as follows:

[0060] y(t+1)=f(y(t),y(t-1),x(t-τ),x(t-τ-1))

[0061] Where f(·) is the nonlinear mapping function to be learned, τ represents the system lag, x(t) represents the auxiliary variable input at time t, and y(t+1) represents the predicted value of the crystal diameter at time t+1. Specifically, the current and previous moments of each auxiliary variable, plus the current and previous moments of the crystal diameter output, are used to predict the next moment of the crystal diameter.

[0062] Step 3.1, TCN consists of three parts: residual connection, causal convolution and dilated convolution;

[0063] Causal convolution: The value before time t of the previous layer determines the value at time t of the next layer. Causal convolution ensures that information at future moments will not be leaked. The convolution kernel of causal convolution is based on the current and past time steps in the time series. Specifically, causal convolution ensures that the output of each time step depends only on the current and previous input data;

[0064] Dilated Convolution: The dilated convolution operation is defined as follows:

[0065]

[0066] In the formula, * represents the convolution operation, k represents the convolution kernel size, and d represents the dilation factor;

[0067] Convolutional networks achieve a larger perceptual range with fewer layers, and dilated convolution increases the receptive field by inserting holes between the elements of the convolution kernel, thereby capturing dependencies over longer time spans.

[0068] Residual connection: In each residual block, the input signal passes through one or more convolution layers and is added to the original input signal to form a residual. This residual connection helps alleviate the gradient vanishing problem and makes deeper network training more efficient and stable.

[0069] The calculation formula is as follows:

[0070] o=Activation(F(x)+x)

[0071] In the formula, F(x) represents the output of the convolutional layer, x represents the input, and Activation represents the activation function;

[0072] Step 3.2, the memory unit consists of a forget gate, an input gate, and an output gate. Input x t The flow through the memory unit is as follows:

[0073] Step 3.2.1. Input information passes through the forget gate. The forget gate determines how much information the memory unit forgets from the previous state. The calculation formula is:

[0074] f t =σ[W f (h t-1 ,x t )+b f ]

[0075] Where W f is the weight matrix of the forget gate, b f is the offset vector of the forget gate, σ is the sigmoid activation function, f t After the sigmoid function operation, it is constrained between (0, 1), 1 means "keep everything" and 0 means "forget everything";

[0076] Step 3.2.2, then pass through the input gate, which determines which information will be retained. The calculation formula is:

[0077] i t =σ[W f (h t-1 ,x t )+b i ]

[0078]

[0079] Where W f is the weight matrix of the input gate; b i is the offset vector of the input gate; is the state variable for acquiring new information; tanh is the hyperbolic tangent curve;

[0080] Step 3.2.3, after determining the information to be forgotten and remembered, update the memory unit state variable. The calculation formula is:

[0081]

[0082] In the formula, C t is the updated memory unit state variable; ⊙ is the dot multiplication operation;

[0083] Step 3.2.4: The processed data finally passes through the output gate O t , the calculation of the output gate determines the hidden layer state variable h t , the calculation formula is:

[0084] O t =σ[W O (h t-1 ,x t )+b O ]

[0085] h t =O t ⊙tanh(C t )

[0086] Where W O is the weight matrix of the output gate; b O is the output gate offset vector;

[0087] Step 3.2.5, finally pass through a fully connected layer to get the final prediction result. The calculation formula of the fully connected layer is as follows:

[0088] y t =σ(W y [h t ,h t ]+b y )

[0089] Where W y is the output state weight, b y is the output offset vector, y t represents the predicted crystal diameter;

[0090] Step 3.3, TCN and LSTM are combined to obtain TCN-LSTM. The specific structure of TCN-LSTM includes: (1) TCN layer: First, the input data sequence is processed by a series of TCN layers. The TCN layer uses causal convolution and dilated convolution to capture the long-term dependencies in the sequence; (2) LSTM layer: The output of the TCN layer is passed to the LSTM layer for feature fusion operation. The LSTM layer further processes the feature sequence data output by the TCN layer and uses memory units to capture complex time dependencies and sequence patterns; (3) Output layer: The output of the LSTM layer passes through a fully connected layer or other types of output layers to generate the final crystal diameter prediction result.

[0091] Step 4 is as follows:

[0092] Use the GWO algorithm to optimize the hyperparameters of TCN-LSTM:

[0093] The hunting process of the gray wolf is as follows:

[0094] (1) Surrounding the prey

[0095] In the GWO algorithm, the gray wolf uses the following position update formula to surround the prey during the hunting process:

[0096]

[0097]

[0098] In the formula, is the distance between the gray wolf and its prey, Update the formula for the gray wolf's position, and are the position vector of the prey and the position vector of the gray wolf, respectively, and t is the current iteration number;

[0099] and To determine the coefficient, the calculation formula is as follows:

[0100]

[0101]

[0102] In the formula, and are two one-dimensional random number vectors whose components are between [0, 1]. Used to simulate the attack behavior of gray wolves on prey. The value is subject to The influence of convergence factor The value of decreases linearly from 2 to 0 as the number of iterations increases;

[0103] (2) Chasing prey

[0104] In the iterative process, α, β and δ are used to guide the movement of ω, thereby achieving global optimization, using the positions of α, β and δ Update the positions of all gray wolves using the following equation:

[0105]

[0106]

[0107]

[0108] In the formula, Respectively represent the distances of the individual ω gray wolf from the α, β and δ wolf packs:

[0109]

[0110]

[0111]

[0112]

[0113] In the formula, They represent the position that the individual ω gray wolf needs to adjust under the influence of α, β and δ wolf packs, and ω takes the average value

[0114] (3) Attacking prey

[0115] In the following formula, t represents the current number of iterations, and T represents the maximum number of iterations. When the value of a decreases from 2 to 0, the corresponding value of A also changes in the interval [-a, a]. The larger the value of a, the farther away the gray wolf will be from the prey, hoping to find a more suitable prey, thus prompting the wolf pack to conduct a global search. If the value of a is smaller, the gray wolf will be closer to the prey, prompting the wolf pack to conduct a local search.

[0116]

[0117] Step 5 is implemented according to the following steps:

[0118] The SHAP algorithm is used to interpret the crystal diameter prediction results. The Shapley value of each feature is calculated through samples to measure the contribution of each feature to the model prediction results. The calculation of the Shapley value is as follows:

[0119]

[0120] In the formula, φ i(f) is the Shapley value of the i-th sample; N is the set of participants, including i; S is a subset of N, excluding i, and |S| is the size of the set S; represents the feature subset that does not include feature i; f(S) represents the model prediction result containing feature subset S.

[0121] The beneficial effect of the present invention is that the interpretable data-driven semiconductor silicon single crystal growth diameter prediction method overcomes the shortcomings of traditional crystal diameter mechanism modeling, inaccurate model and high modeling cost. At the same time, compared with the existing data-driven crystal diameter prediction method, it can take into account the mutual influence between multiple auxiliary characteristic variables, ensure that the diameter prediction is more accurate, and can explain the contribution of different auxiliary characteristic variables to the crystal diameter prediction results. BRIEF DESCRIPTION OF THE DRAWINGS

[0122] Figure 1 It is a flow chart of a method for predicting the diameter of a CZ silicon single crystal with interpretable data-driven method of the present invention;

[0123] Figure 2 It is a MIC auxiliary characteristic variable analysis result diagram of an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0124] Figure 3 It is a time delay identification result diagram between each auxiliary characteristic variable and the crystal diameter in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0125] FIG4( a ) is a structural diagram of causal convolution in a TCN network in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0126] FIG4( b ) is a structural diagram of the dilated convolution in the TCN network in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0127] FIG4( c ) is a structural diagram of a residual block in TCN in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0128] FIG4( d ) is a structural diagram of an LSTM unit in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0129] FIG4(e) is a schematic diagram of the structure of a TCN-LSTM model in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0130] Figure 5 It is a comparison diagram of the predictive effect of an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention and other models;

[0131] Figure 6 It is a characteristic importance diagram for explaining each auxiliary variable based on the SHAP algorithm in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention;

[0132] Figure 7 It is a force-directed graph based on the SHAP algorithm to explain various auxiliary variables in an interpretable data-driven CZ silicon single crystal diameter prediction method of the present invention. DETAILED DESCRIPTION

[0133] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0134] The present invention can interpret the data driven semiconductor silicon single crystal growth diameter prediction method, the flow chart is as follows Figure 1 As shown, the specific implementation steps are as follows:

[0135] Step 1, preprocessing the experimental data, and obtaining the auxiliary characteristic variables related to the crystal diameter variable by calculating the maximum mutual information coefficient MIC (Maximal Information Coefficient);

[0136] Step 1 is implemented according to the following steps:

[0137] Step 1.1: Use the triple standard deviation principle to eliminate outliers:

[0138] The three-times standard deviation principle is used to eliminate outliers in the crystal diameter observation data. First, the mean μ and standard deviation σ of the crystal diameter data y are calculated according to the following formula:

[0139]

[0140]

[0141] In the formula, y i is the i-th crystal diameter data sample, and n is the total number of sample data;

[0142] Determine the threshold of outliers. The threshold range of outliers is determined by adding or subtracting three times the standard deviation from the mean, that is, the lower limit is μ-3σ and the upper limit is μ+3σ. Check whether there are values ​​in the data set that are less than the lower threshold or greater than the upper threshold. Values ​​that are less than the lower threshold or greater than the upper threshold are considered outliers and are deleted from the data set.

[0143] Step 1.2: Use the minimum-maximum scaling method to normalize the collected auxiliary variables. When normalizing the crystal diameter data, the specific steps are as follows:

[0144] Find the maximum value y in the crystal diameter datamax Minimum value y min , linearly scale the original crystal diameter data to between 0 and 1, as follows:

[0145]

[0146] In the formula, y nor represents the normalized data;

[0147] The normalized data is used as training data for the crystal diameter prediction model;

[0148] Step 1.3: Use the maximum mutual information coefficient (MIC) to screen auxiliary variables that are highly correlated with crystal diameter. The specific steps are as follows:

[0149] The mutual information used in the MIC principle is expressed by the following equation:

[0150]

[0151] Where x represents the operating variables related to silicon single crystal, such as power, crucible rise, crystal rise and thermal field temperature; y represents the actual crystal diameter variable, p(x, y) represents the joint probability density distribution function of x and y, and p(x) and p(y) distribution represent the marginal probability density functions of x and y. The calculation formula of MIC is as follows:

[0152]

[0153] In the formula, a represents the number of grid divisions in the x direction, b represents the number of grid divisions in the y direction, and B is a set variable. B is the 0.6th power of the data capacity of the crystal diameter variable involved in calculating MIC. The range of MIC is [0, 1]. The larger the value of MIC, the stronger the correlation between the two variables.

[0154] Step 2: Estimate the time delay order between the auxiliary characteristic variable and the crystal diameter based on the multiple correlation time delay estimation method of the time correlation matrix, and perform time alignment of the data sequence;

[0155] Step 2 is implemented as follows:

[0156] N features are screened out through MIC, and the original time series data matrix X is obtained by sampling M times on the time axis:

[0157]

[0158] The subscript i of the element in the original data matrix represents the process variable, and the superscript t represents the sampling time, where x0 represents the time series data of the zero-delay reference variable, that is, the target variable, and x i represents the time series data of the i-th process variable;

[0159] Assume that the multiple correlation delay parameters τ between the target variable and each process variable are represented by the delay sequence Γ, that is:

[0160] Γ=[τ1,τ2,...,τ i ,...,τ N ]

[0161] The sampling period of each process variable in the original time series data matrix is ​​T, then the delay parameter of the i-th process variable relative to the target variable needs to satisfy:

[0162] Γ=DT=[d1T,d2T,...,d i T,...,d N T]

[0163] Where D is the time base sequence, d i is the time base of the ith process variable, d i is a dimensionless integer, which ensures that the multiple correlation delay parameters of the process variables are integer multiples of the sampling period, ensuring the feasibility of the process variables in the time series reconstruction;

[0164] In order to estimate the multiple correlation delay parameters of different process variables, the time series correlation matrix composed of process variables and target variables under different time base sequences is reconstructed in the time and space dimensions, and the feedback information with different delays is identified by using the time base. The reconstruction process is as follows: First, from the zero delay reference variable Select F data with sampling frequency T from the tth time to form a time series

[0165]

[0166] FT≥max(τ1,τ2,…,τ i ,…,τ N )

[0167] In the formula, F must be greater than the delay of each auxiliary variable;

[0168] According to the multi-correlated time delay sequence defined above, relative to the target variable, the process variable x i The delay coefficient relative to the zero-delay reference variable is τ i , so take the process variable x i From t-τ i Generate a new time series with F data starting at time That is, from td i F consecutive sampling data starting from T:

[0169]

[0170] The temporal correlation matrix reconstructed in time and space according to the time delay parameters of each process variable

[0171]

[0172] Let the target variable in the reconstructed time series correlation matrix be the parent sequence, and the multivariate process variable whose delay parameter needs to be estimated be the subsequence. First, calculate the difference between the corresponding elements of the parent sequence and the subsequence, and express it with matrix A:

[0173]

[0174] In the formula,

[0175] Let a=min{a ji}, b = max{a ji}, then the grey correlation coefficient matrix is ​​expressed by γ:

[0176]

[0177] In the formula, ρ is the resolution coefficient. The smaller ρ is, the greater the resolution is. When ρ is 0.5, the grey correlation coefficient matrix can better reflect the actual correlation between multivariate variables.

[0178] In order to describe the overall correlation between the subsequence and the parent sequence, the elements in the grey correlation coefficient matrix are averaged by column, that is, the grey correlation degree of the subsequence relative to the parent sequence is defined as:

[0179] r=[r1,r2,...,r i ,...,r N ]

[0180] In the formula,

[0181] In order to quantitatively describe the multiple correlation degrees between each process variable and the target variable in the reconstructed temporal correlation matrix, the sum of the grey correlation degrees R is used to describe the overall correlation between the multivariate process variables and the target variables, where R is expressed as:

[0182]

[0183] The sum of the grey correlation degrees R can measure the correlation between the multivariate time series in the time series correlation matrix to a certain extent. The larger the value of R, the stronger the correlation of the multivariate time series. Finally, the grey wolf optimization algorithm is used to optimize the R value. When the R value reaches the maximum, it is the optimal time delay.

[0184] Step 3, establishing a crystal diameter deep learning prediction model TCN-LSTM based on the temporal convolutional network TCN (Temporal Convolutional Network) and the long short-term memory network LSTM (Long Short-Term Memory);

[0185] Step 3 is implemented as follows:

[0186] Since the growth process of CZ silicon single crystal is a nonlinear and large-hysteresis dynamic process, based on the historical data of silicon single crystal growth, a crystal diameter prediction model TCN-LSTM is established by combining the time convolution network TCN and the long short-term memory network LSTM. The GWO optimization algorithm is used to optimize its hyperparameters. The TCN-LSTM model is used to learn the nonlinear relationship between the auxiliary variables and the crystal diameter. The nonlinear relationship is described as follows:

[0187] y(t+1)=f(y(t),y(t-1),x(t-τ),x(t-τ-1))

[0188] Where f(·) is the nonlinear mapping function to be learned, τ represents the system lag, x(t) represents the auxiliary variable input at time t, and y(t+1) represents the predicted value of the crystal diameter at time t+1. Specifically, the current and previous moments of each auxiliary variable, plus the current and previous moments of the crystal diameter output, are used to predict the next moment of the crystal diameter.

[0189] Step 3.1, TCN is a deep learning network framework specifically used to extract features of time series data. It consists of three parts: residual connection, causal convolution and dilated convolution;

[0190] Causal convolution: Causal convolution is a one-way structure that cannot see future data. The value of the previous layer before time t determines the value of the next layer at time t. Causal convolution ensures that information at future times will not be leaked. The causal convolution structure is shown in Figure 4(a). The convolution kernel of causal convolution is based on the current and past time steps in the time series to prevent future information from leaking to the current time step. Specifically, causal convolution ensures that the output of each time step depends only on the current and previous input data, and does not involve future data.

[0191] Dilated convolution: Since the size of the convolution kernel limits the feature extraction of causal convolution, which requires linear stacking of more layers to capture longer dependencies, dilated convolution solves this problem well. The definition of the dilated convolution operation is as follows:

[0192]

[0193] In the formula, * represents the convolution operation, k represents the convolution kernel size, and d represents the dilation factor;

[0194] In this way, the convolutional network obtains a larger perceptual range through fewer layers. The dilated convolution structure is shown in Figure 4(b). Dilated convolution increases the receptive field by inserting holes between the elements of the convolution kernel, thereby capturing dependencies over a longer time span. Dilated convolution allows the network to cover a larger input range without increasing the amount of computation, thereby more effectively processing long-term dependencies. In the application of the convolution kernel, the dilation factor determines the size of the interval between the elements of the convolution kernel.

[0195] Residual connection: Residual connection helps to avoid the gradient vanishing or gradient exploding problem when training the network. The residual block structure is shown in Figure 4(c). Residual block is a method to simplify network training and optimization by introducing short-circuit connections. In each residual block, the input signal passes through one or more convolution layers and is added to the original input signal to form a residual. This residual connection helps alleviate the gradient vanishing problem and makes deeper network training more efficient and stable.

[0196] The calculation formula is as follows:

[0197] o=Activation(F(x)+x)

[0198] In the formula, F(x) represents the output of the convolutional layer, x represents the input, and Activation represents the activation function;

[0199] Step 3.2, LSTM is a network model improved on the basis of recurrent neural network. LSTM introduces a gating mechanism to replace the unit structure in RNN, thereby improving the problem of gradient vanishing or gradient exploding in recurrent neural network. The structure of the memory unit in LSTM is shown in Figure 4(d). The memory unit consists of a forget gate, an input gate, and an output gate. The input x t The flow through the memory unit is as follows:

[0200] Step 3.2.1. Input information passes through the forget gate. The forget gate determines how much information the memory unit forgets from the previous state. The calculation formula is:

[0201] f t =σ[W f (h t-1 ,x t )+b f ]

[0202] Where W f is the weight matrix of the forget gate, b f is the offset vector of the forget gate, σ is the sigmoid activation function, f tAfter the sigmoid function operation, it is constrained between (0, 1), 1 means "keep everything" and 0 means "forget everything";

[0203] Step 3.2.2, then pass through the input gate, which determines which information will be retained. The calculation formula is:

[0204] i t =σ[W f (h t-1 ,x t )+b i ]

[0205]

[0206] Where W f is the weight matrix of the input gate; b i is the offset vector of the input gate; is the state variable for acquiring new information; tanh is the hyperbolic tangent curve;

[0207] Step 3.2.3, after determining the information to be forgotten and remembered, update the memory unit state variable. The calculation formula is:

[0208]

[0209] In the formula, C t is the updated memory unit state variable; ⊙ is the dot multiplication operation;

[0210] Step 3.2.4: The processed data finally passes through the output gate O t , the calculation of the output gate determines the hidden layer state variable h t , the calculation formula is:

[0211] O t =σ[W O (h t-1 ,x t )+b O ]

[0212] h t =O t ⊙tanh(C t )

[0213] Where W O is the weight matrix of the output gate; b O is the output gate offset vector;

[0214] Step 3.2.5, finally pass through a fully connected layer to get the final prediction result. The calculation formula of the fully connected layer is as follows:

[0215] y t=σ(W y [h t ,h t ]+b y )

[0216] Where W y is the output state weight, b y is the output offset vector, y t represents the predicted crystal diameter;

[0217] Step 3.3, TCN takes into account both long-term and short-term dependencies between data and simultaneously utilizes the parallel computing advantages of TCN and the sequence data processing capabilities of LSTM. TCN and LSTM are combined to obtain TCN-LSTM, and its structure diagram is shown in Figure 4(e); the TCN-LSTM structure contains two network architectures, TCN and LSTM, aiming to utilize the advantages of both. The specific structure includes: (1) TCN layer: First, the input data sequence (including operational variable data such as crystal rise, thermal field temperature, heating power, etc.) is processed through a series of TCN layers. The TCN layer uses causal convolution and dilated convolution to capture the long-term dependencies in the sequence; (2) LSTM layer: The output of the TCN layer is passed to the LSTM layer for feature fusion operation. The LSTM layer further processes the feature sequence data output by the TCN layer and uses memory units to capture complex time dependencies and sequence patterns; (3) Output layer: The output of the LSTM layer passes through a fully connected layer or other types of output layers to generate the final crystal diameter prediction result. This structure combines the powerful convolution characteristics of TCN and the time series modeling capabilities of LSTM to improve the effect of sequence prediction and processing.

[0218] Step 4: Train the TCN-LSTM network model to learn the nonlinear relationship between input variables and output variables, and use the Grey Wolf Optimizer (GWO) algorithm to optimize the network model hyperparameters.

[0219] Step 4 is as follows:

[0220] The settings of some hyperparameters of the TCN-LSTM network generally directly affect the prediction results. In order to achieve the best prediction results, the GWO algorithm is used to optimize the hyperparameters of TCN-LSTM. The principle of GWO is as follows:

[0221] The hunting process of the gray wolf is as follows:

[0222] (1) Surrounding the prey

[0223] In the GWO algorithm, the gray wolf uses the following position update formula to surround the prey during the hunting process:

[0224]

[0225]

[0226] In the formula, is the distance between the gray wolf and its prey, Update the formula for the gray wolf's position, and

[0227] are the position vector of the prey and the position vector of the gray wolf, respectively, and t is the current iteration number;

[0228] and

[0229] To determine the coefficient, the calculation formula is as follows:

[0230]

[0231]

[0232] In the formula, and

[0233] are two one-dimensional random number vectors whose components are between [0, 1]. Used to simulate the attack behavior of gray wolves on prey. The value is subject to The influence of convergence factor The value of decreases linearly from 2 to 0 as the number of iterations increases;

[0234] (2) Chasing prey

[0235] In the iterative process, α, β and δ are used to guide the movement of ω, thereby achieving global optimization, using the positions of α, β and δ Update the positions of all gray wolves using the following equation:

[0236]

[0237]

[0238]

[0239] In the formula, Respectively represent the distances of the individual ω gray wolf from the α, β and δ wolf packs:

[0240]

[0241]

[0242]

[0243]

[0244] In the formula, They represent the position that the individual ω gray wolf needs to adjust under the influence of α, β and δ wolf packs, and ω takes the average value

[0245] (3) Attacking prey

[0246] In the following formula, t represents the current number of iterations, T represents the maximum number of iterations set, and when the value of a decreases from 2 to 0, the corresponding value of A also changes in the interval [-a, a]: the larger the value of a, the farther the gray wolf will move from the prey, hoping to find a more suitable prey, thus prompting the wolf pack to conduct a global search (|A|>1); if the value of a is smaller, the gray wolf will move closer to the prey, prompting the wolf pack to conduct a local search (|A|<1);

[0247]

[0248] Step 5: Obtain the crystal diameter prediction result based on the trained crystal diameter deep learning prediction model TCN-LSTM. Then, use the SHAP (SHapley Additive exPlanations) algorithm to interpret the crystal diameter prediction result, so as to obtain the contribution of different auxiliary feature variables to the crystal diameter prediction result.

[0249] Step 5 is implemented according to the following steps:

[0250] Since the deep learning network model is a black box model, the present invention uses the SHAP algorithm to explain the crystal diameter prediction results. SHAP is a post hoc interpretation method for explaining the prediction results of the deep learning model. The idea is to introduce the Shapley value in cooperative game theory into the deep learning model, and calculate the Shapley value of each feature through samples to measure the contribution of each feature to the model prediction results. The calculation of the Shapley value is as follows:

[0251]

[0252] In the formula, φ i (f) is the Shapley value of the i-th sample; N is the set of participants, including i; S is a subset of N, excluding i, and |S| is the size of the set S; represents the feature subset that does not include feature i; f(S) represents the model prediction result containing feature subset S.

[0253] Example 1

[0254] The present invention predicts crystal diameter by establishing a data-driven deep learning model. The flow chart is as follows: Figure 1 As shown, the specific implementation steps are as follows:

[0255] Step 1: Collect process variable data during silicon single crystal growth through sensors, including 15 process variables such as power, crucible lifting, crucible rotation, argon flow, crystal lifting, crystal rotation, thermal field temperature, main chamber low vacuum, auxiliary chamber vacuum, main heater voltage, main heater current, main heater resistance, crystal weight, crystal length, average growth rate and crystal diameter. Perform outlier processing and normalization respectively. The specific steps are as follows:

[0256] Step 1.1: Remove outliers using the three standard deviation principle

[0257] To remove outliers from crystal diameter using the three standard deviation principle, we must first calculate the mean (μ) and standard deviation (σ) of the crystal diameter data (y) according to the following formula:

[0258]

[0259]

[0260] In the formula, y i is the ith sample data, n is the total number of samples, and the present invention takes n=6000.

[0261] Step 1.2: Use the minimum-maximum scaling method to normalize the collected auxiliary variables, such as normalizing the crystal diameter data. The specific steps are as follows:

[0262] Find the maximum value (y max )Minimum value (y min ), and scale it linearly to between 0 and 1, as follows:

[0263]

[0264] In the formula, y nor Represents normalized data.

[0265] The normalization process of other process variables is the same. The normalized data is used as the input and output of the model.

[0266] Step 2: Use historical silicon single crystal growth data to screen several auxiliary variables with the greatest correlation with crystal diameter through the maximum mutual information coefficient (MIC).

[0267] The main idea of ​​MIC is: if there is a certain correlation between two variables, then after some grid division is performed on the scatter plot of the two variables, the mutual information of the two variables can be calculated according to the approximate probability density distribution of the two variables in the grid. After regularization, this value can be used to measure the correlation between the two variables. The MIC principle uses the concept of mutual information, which can be expressed by the following equation:

[0268]

[0269] In the formula, x and y represent variables, p(x, y) represents the joint probability density distribution function of x and y, and p(x) and p(y) represent the marginal probability density functions of x and y. The calculation formula of MIC is as follows:

[0270]

[0271] In the formula, a and b represent the number of grid divisions in the x and y directions, B is a variable that is considered to be set, which is generally 0.6 power of the data size. The range of MIC is [0, 1]. The larger the value of MIC, the stronger the correlation between the two variables.

[0272] The principle of using MIC to screen the main features is: (1) the correlation between the selected features and the crystal diameter should be as large as possible, that is, the MIC should be as large as possible; (2) the correlation between the selected features should be as small as possible, that is, the MIC should be as small as possible. The auxiliary variables with the greatest correlation with the crystal diameter are screened out based on the MIC.

[0273] Based on the 16 process variables collected above, including power, crucible lifting, crucible rotation, argon flow rate, crystal lifting, crystal rotation, thermal field temperature, main chamber low vacuum, auxiliary chamber vacuum, main heater voltage, main heater current, main heater resistance, crystal weight, crystal length, average growth rate and crystal diameter, the maximum mutual information coefficients between the first 15 and the diameter and the maximum mutual information coefficients between them were calculated by Python program design. The MIC analysis results are shown in the figure. Figure 2 ,According to the two principles of MIC screening features mentioned above, 5 ,features are screened out, namely thermal field temperature, crystal ,transformation, crystal rise, power, and crucible rise.

[0274] Step 3: Use the multiple correlation delay estimation method based on the timing correlation matrix to identify the lag order of each auxiliary variable relative to the crystal diameter.

[0275] Step 4: Determine the model structure and establish a TCN-LSTM network model to learn the nonlinear mapping relationship between input and output, and use the Grey Wolf Optimizer (GWO) algorithm to optimize some hyperparameters of the network.

[0276] Step 5: Use the SHAP algorithm to interpret the crystal diameter prediction results and obtain the influence of different process variables on the crystal diameter.

[0277] Example 2

[0278] The present invention predicts the crystal diameter by establishing a data-driven deep learning model. The flow chart is as follows: Figure 1 As shown, the specific steps are as follows:

[0279] Step 1: Preprocess the collected crystal growth process data and divide it into training set and test set for training the TCN-LSTM model.

[0280] Step 2: Based on the historically accumulated industrial process data, MIC is used to screen out auxiliary variables related to crystal diameter.

[0281] Step 2 is implemented as follows:

[0282] The maximum mutual information coefficient between each process variable and the crystal diameter is calculated respectively, and the maximum mutual information coefficient between each variable is calculated respectively. According to the principles that (1) the correlation between each variable and the crystal diameter is as large as possible and (2) the correlation between each process variable is as small as possible, several auxiliary variables are selected as the input of the crystal diameter prediction model.

[0283] According to the principle of MIC screening characteristics, 5 auxiliary variables with the greatest correlation with crystal diameter were screened out through MIC, namely thermal field temperature, crystal rotation, crystal rise, power, and crucible rise. The MIC analysis results are shown in the figure below: Figure 2 , these five process variables are used to predict the crystal diameter.

[0284] Step 3: Use the multiple correlation delay estimation method based on the time series correlation matrix to estimate the time lag order of the model's input and output.

[0285] Step 3 is implemented as follows:

[0286] From step 2, we know that 5 features are screened out through MIC, so N=5 in the following formula, and the original time series data matrix X is obtained by sampling M times on the time axis. In the present invention, M=6000, and the original time series data matrix X is as follows:

[0287]

[0288] In the formula, N = 5, indicating 5 process variables, M = 6000, indicating 6000 sampling moments, and x0 represents the time series data of the zero-delay reference variable, i.e., the target variable. i Represents the time series data of the ith process variable.

[0289] Assume that the multiple associated delay parameters between the target variable and each process variable are represented by a delay sequence, that is:

[0290] Γ=[τ1,τ2,...,τ i ,...,τ N ]

[0291] Wherein, N=5.

[0292] Each process variable in the original time series data matrix collects data every 2 seconds, that is, T = 2, then the time delay parameter of each i-th process variable relative to the target variable needs to satisfy:

[0293] Γ=DT=[d1T,d2T,...,d i T,...,d N T]

[0294] Where, T = 2, D is the time base sequence, d i is the time base of each process variable of the ith dimension, which is a dimensionless integer. This ensures that the multiple correlation delay parameters of the process variables are integer multiples of the sampling period, which is to ensure the feasibility of the process variables in the timing reconstruction.

[0295] In order to estimate the multiple correlation delay parameters of different process variables, the time series correlation matrix composed of process variables and target variables under different time base sequences is reconstructed in the time and space dimensions, and the feedback information with different delays is identified by using the time base. The reconstruction process is as follows: First, from the zero delay reference variable Select F data with sampling frequency T from the tth time to form a time series

[0296]

[0297] FT≥max(τ1,τ2,…,τ i ,…,τ N )

[0298] The data contained in the sequence must contain at least one complete cycle, and the F data selected from the sampling time t must be guaranteed to contain a transition state of the zero-delay reference variable.

[0299] According to the multiple correlation time delay sequences defined above, relative to the target variable, the process variable x i The delay coefficient relative to the zero-delay reference variable is τ i , so take the process variable x i From t-τ i Generate a new time series with F data starting at time That is, from td iF consecutive sampling data starting from T:

[0300]

[0301] Since the new time series is from td i F consecutive sampling data starting from T, and the data contained must contain at least one complete delay cycle, that is, the F data must ensure a transition state of the zero-delay reference variable (target variable), so F in the present invention is taken as 600.

[0302] The temporal correlation matrix reconstructed in time and space according to the time delay parameters of each process variable

[0303]

[0304] Let the target variable in the reconstructed time series correlation matrix be the parent sequence, and the multivariate process variable whose delay parameter needs to be estimated be the subsequence. First, calculate the difference between the corresponding elements of the parent sequence and the subsequence, and express it with matrix A:

[0305]

[0306] In the formula,

[0307] Let a=min{a ji}, b = max{a ji}, then the grey correlation coefficient matrix can be expressed by γ:

[0308]

[0309] In the formula, ρ is the resolution coefficient. The smaller ρ is, the greater the resolution is. When ρ is 0.5, the grey correlation coefficient matrix can better reflect the actual correlation between multivariate variables. Therefore, the present invention selects ρ=0.5.

[0310] In order to describe the overall correlation between the subsequence and the parent sequence, the elements in the grey correlation coefficient matrix are averaged by column, that is, the grey correlation degree of the subsequence relative to the parent sequence is defined as:

[0311] r=[r1,r2,...,r i ,...,r N ]

[0312] In the formula,

[0313] In order to quantitatively describe the multiple correlation degrees between each process variable and the target variable in the reconstructed temporal correlation matrix, the sum of the grey correlation degrees R is used to describe the overall correlation between the multivariate process variables and the target variable, where R is expressed as:

[0314]

[0315] The sum of the grey correlation R can measure the correlation between the multivariate time series in the time series correlation matrix to a certain extent. The larger the value, the stronger the correlation of the multivariate time series. The present invention uses the GWO optimization algorithm to optimize the R value. When the R value reaches the maximum, the corresponding multi-correlation time delay is the real multi-correlation time delay of the CZ silicon single crystal growth process. The delay result is as follows: Figure 3 .

[0316] Step 4: Determine the model structure and establish a TCN-LSTM network model to learn the nonlinearity between input and output, and use the Grey Wolf Optimizer (GWO) algorithm to optimize some hyperparameters of the network.

[0317] Step 5: Use the SHAP algorithm to interpret the prediction results and obtain the contribution model of different process variables to the crystal diameter prediction results, thereby helping on-site engineering and technical personnel make correct control decisions.

[0318] Embodiment 3:

[0319] The present invention predicts the crystal diameter of the CZ silicon single crystal growth process by establishing a data-driven deep learning model. The modeling flow chart is as follows: Figure 1 As shown, the specific implementation steps are as follows:

[0320] Step 1: Collect process variable data at the same time during the growth of silicon single crystals, first use the three standard deviation principle to process the data for outliers, and then use the minimum-maximum scaling method to normalize the data.

[0321] Step 2: Using the normalized data, MIC was used to select the auxiliary variable with the greatest correlation with the crystal diameter from 15 process variables, including power, crucible lifting, crucible rotation, argon flow, crystal lifting, crystal rotation, thermal field temperature, main chamber low vacuum, auxiliary chamber vacuum, main heater voltage, main heater current, main heater resistance, crystal weight, crystal length, and average growth rate. The MIC analysis results are shown in the figure. Figure 2 Finally, five auxiliary variables were selected, namely, thermal field temperature, crystal rotation, crystal rise, power, and crucible rise.

[0322] Step 3: Through the large amount of data accumulated during the crystal pulling process, a time-series-based correlation matrix is ​​constructed, and then the time-series correlation matrix-based multiple correlation delay estimation method is used to estimate the time lag order between the auxiliary variable and the crystal diameter. The results of the time lag order identification are shown in Figure 3. Figure 3 , the figure shows the time lag order between the five auxiliary variables and the crystal diameter.

[0323] Step 4: Determine the model structure and establish the nonlinearity between the TCN-LSTM neural network learning input and output, and use the Grey Wolf Optimizer (GWO) algorithm to optimize some hyperparameters of the network.

[0324] The present invention reconstructs data through the above-screened process variables and the identified delay coefficients between each process variable and the crystal diameter, and finally obtains 5000 groups of samples. The first 4000 groups of samples are used as the training set of the model to train the model, and the last 1000 groups of sample data are used as the test set of the model. In addition, the GWO algorithm is used to optimize the hyperparameters.

[0325] The present invention uses the GWO optimization algorithm to optimize the four hyperparameters of minimum batch size, learning rate, TCN convolution kernel size, and LSTM hidden layer neuron number. The optimized hyperparameter results are shown in Table 1:

[0326] Table 1 Hyperparameters after GWO optimization

[0327]

[0328] The prediction effect of the method of the present invention is compared with other data-driven models (such as TCN model and LSTM model). Figure 5 As shown in Table 2, the evaluation indicators are shown in Table 2. It can be seen that the prediction effect of the method of the present invention is better than that of other methods.

[0329] Table 2 Comparison of prediction effect evaluation indicators of each model

[0330]

[0331] Step 5: Use the SHAP algorithm to interpret the crystal diameter prediction results and obtain the influence of different process variables on the crystal diameter. The results are as follows: Figure 6 and Figure 7 . Figure 6 The size of the median value indicates the influence of the auxiliary feature variable on the prediction result. The crystal transformation variable has the greatest influence on the prediction result, and the power variable has the smallest influence. Figure 7 It shows whether the auxiliary characteristic variable has a positive or negative effect on the crystal diameter prediction result. In the figure, red indicates a positive effect on the prediction result, and blue indicates a negative effect on the prediction result.

Claims

1. An interpretable data-driven semiconductor silicon single crystal growth diameter prediction method, characterized in that: Follow the steps below to implement it: Step 1, preprocessing the experimental data, and obtaining the auxiliary characteristic variables related to the crystal diameter variable by calculating the maximum mutual information coefficient MIC; Step 2, estimating the time delay order between the auxiliary characteristic variable and the crystal diameter, and performing time alignment of the data sequence; Step 3: Establish a crystal diameter deep learning prediction model TCN-LSTM based on the temporal convolutional network TCN and the long short-term memory network LSTM; Step 4: Train the TCN-LSTM network model and use the Grey Wolf Optimization Algorithm (GWO) to optimize the network model hyperparameters. Step 5: Obtain the crystal diameter prediction result based on the trained crystal diameter deep learning prediction model TCN-LSTM.

2. The interpretable data-driven semiconductor silicon single crystal growth diameter prediction method according to claim 1, characterized in that: The step 1 is specifically implemented according to the following steps: Step 1.1: Use the triple standard deviation principle to eliminate outliers: The three-times standard deviation principle is used to eliminate outliers in the crystal diameter observation data. First, the mean μ and standard deviation σ of the crystal diameter data y are calculated according to the following formula: In the formula, y i is the i-th crystal diameter data sample, and n is the total number of sample data; Determine the threshold of outliers. The threshold range of outliers is determined by adding or subtracting three times the standard deviation from the mean, that is, the lower limit is μ-3σ and the upper limit is μ+3σ. Check whether there are values ​​in the data set that are less than the lower threshold or greater than the upper threshold. Values ​​that are less than the lower threshold or greater than the upper threshold are considered outliers and are deleted from the data set. Step 1.2: Use the minimum-maximum scaling method to normalize the collected auxiliary variables. When normalizing the crystal diameter data, the specific steps are as follows: Find the maximum value y in the crystal diameter data max Minimum value y min , linearly scale the original crystal diameter data to between 0 and 1, as follows: In the formula, y nor represents the normalized data; The normalized data is used as training data for the crystal diameter prediction model; Step 1.3: Use the maximum mutual information coefficient MIC to screen auxiliary variables with a greater correlation with crystal diameter. The specific steps are as follows: The mutual information used in the MIC principle is expressed by the following equation: In the formula, x represents the operating variable related to the silicon single crystal, y represents the actual crystal diameter variable, p(x,y) represents the joint probability density distribution function of x and y, and p(x) and p(y) distributions represent the marginal probability density functions of x and y. The calculation formula of MIC is as follows: In the formula, a represents the number of grid divisions in the x direction, b represents the number of grid divisions in the y direction, and B is a set variable. B is the 0.6th power of the data capacity of the crystal diameter variable involved in calculating MIC. The range of MIC is [0, 1]. The larger the value of MIC, the stronger the correlation between the two variables.

3. The interpretable data-driven semiconductor silicon single crystal growth diameter prediction method according to claim 2, characterized in that: The step 2 is specifically implemented according to the following steps: N features are screened out through MIC, and the original time series data matrix X is obtained by sampling M times on the time axis: The subscript i of the element in the original data matrix represents the process variable, and the superscript t represents the sampling time, where x0 represents the time series data of the zero-delay reference variable, that is, the target variable, and x i represents the time series data of the i-th process variable; Assume that the multiple correlation delay parameters τ between the target variable and each process variable are represented by the delay sequence Γ, that is: C=[τ1,τ2,...,τ i ,...,t N ] The sampling period of each process variable in the original time series data matrix is ​​T, then the delay parameter of the i-th process variable relative to the target variable needs to satisfy: Γ=DT=[d1T,d2T,...,d i T,...,d N T] Where D is the time base sequence, d i is the time base of the ith process variable, d i is a dimensionless integer, which ensures that the multiple correlation delay parameters of the process variables are integer multiples of the sampling period, ensuring the feasibility of the process variables in the time series reconstruction; Reconstruct the time series correlation matrix composed of process variables and target variables under different time base sequences in the space-time dimension. The reconstruction process is as follows: First, start from the zero-delay reference variable Select F data with sampling frequency T from the tth time to form a time series FT≥max(τ1,τ2,…,τ i ,,t N ) In the formula, F must be greater than the delay of each auxiliary variable; Relative to the target variable, the process variable x i The delay coefficient relative to the zero-delay reference variable is τ i , so take the process variable x i From t-τ i Generate a new time series with F data starting at time That is, from td i F consecutive sampling data starting from T: The temporal correlation matrix reconstructed in time and space according to the time delay parameters of each process variable Let the target variable in the reconstructed time series correlation matrix be the parent sequence, and the multivariate process variable whose delay parameter needs to be estimated be the subsequence. First, calculate the difference between the corresponding elements of the parent sequence and the subsequence, and express it with matrix A: In the formula, Let a=min{a ji }, b = max{a ji }, then the grey correlation coefficient matrix is ​​expressed by γ: In the formula, ρ is the resolution coefficient. The smaller ρ is, the greater the resolution is. When ρ is 0.5, the grey correlation coefficient matrix can better reflect the actual correlation between multivariate variables. The elements in the grey correlation coefficient matrix are averaged by column, that is, the grey correlation degree of the subsequence relative to the parent sequence is defined as: r＝[r1,r2,...,r i ,...,r N ] In the formula, The sum of grey correlation R is used to describe the overall correlation between multivariate process variables and target variables, where R is expressed as:

4. The interpretable data-driven semiconductor silicon single crystal growth diameter prediction method according to claim 3, characterized in that: The step 3 is specifically implemented according to the following steps: Based on the historical data of silicon single crystal growth, a crystal diameter prediction model TCN-LSTM is established by combining the temporal convolutional network TCN and the long short-term memory network LSTM. The hyperparameters are optimized using the GWO optimization algorithm. The TCN-LSTM model is used to learn the nonlinear relationship between the auxiliary variables and the crystal diameter. The nonlinear relationship is described as follows: y(t+1)=f(y(t),y(t-1),x(t-τ),x(t-τ-1)) In the formula, f(·) is the nonlinear mapping function to be learned, τ represents the system lag, x(t) represents the auxiliary variable input at time t, and y(t+1) represents the predicted value of the crystal diameter at time t+1. Specifically, the current and previous moments of each auxiliary variable, plus the current and previous moments of the crystal diameter output, are used to predict the next moment of the crystal diameter; Step 3.1, TCN consists of three parts: residual connection, causal convolution and dilated convolution; Causal convolution: The value before time t of the previous layer determines the value at time t of the next layer. Causal convolution ensures that information at future moments will not be leaked. The convolution kernel of causal convolution is based on the current and past time steps in the time series. Specifically, causal convolution ensures that the output of each time step depends only on the current and previous input data; Dilated Convolution: The dilated convolution operation is defined as follows: In the formula, * represents the convolution operation, k represents the convolution kernel size, and d represents the dilation factor; Convolutional networks achieve a larger perceptual range with fewer layers, and dilated convolution increases the receptive field by inserting holes between the elements of the convolution kernel, thereby capturing dependencies over longer time spans. Residual connection: In each residual block, the input signal passes through one or more convolution layers and is added to the original input signal to form a residual. This residual connection helps alleviate the gradient vanishing problem and makes deeper network training more efficient and stable. The calculation formula is as follows: o=Activation(F(x)+x) In the formula, F(x) represents the output of the convolutional layer, x represents the input, and Activation represents the activation function; Step 3.2, the memory unit consists of a forget gate, an input gate, and an output gate. Input x t The flow through the memory unit is as follows: Step 3.2.

1. Input information passes through the forget gate. The forget gate determines how much information the memory unit forgets from the previous state. The calculation formula is: f t =σ[W f (h t-1 ,x t )+b f ] Where W f is the weight matrix of the forget gate, b f is the offset vector of the forget gate, σ is the sigmoid activation function, f t After the sigmoid function operation, it is constrained between (0, 1), 1 means "keep everything" and 0 means "forget everything"; Step 3.2.2, then pass through the input gate, which determines which information will be retained. The calculation formula is: i t =σ[W f (h t-1 ,x t )+b i ] Where W f is the weight matrix of the input gate; b i is the offset vector of the input gate; is the state variable for acquiring new information; tanh is the hyperbolic tangent curve; Step 3.2.3, after determining the information to be forgotten and remembered, update the memory unit state variable. The calculation formula is: In the formula, C t is the updated memory unit state variable; ⊙ is the dot multiplication operation; Step 3.2.4: The processed data finally passes through the output gate O t , the calculation of the output gate determines the hidden layer state variable h t , the calculation formula is: The t =σ[W O (h t-1 ,x t )+b O ] h t =O t ⊙tanh(C t ) Where W O is the weight matrix of the output gate; b O is the output gate offset vector; Step 3.2.5, finally pass through a fully connected layer to get the final prediction result. The calculation formula of the fully connected layer is as follows: y t =σ(W y [h t ,h t ]+b y ) Where W y is the output state weight, b y is the output offset vector, y t represents the predicted crystal diameter; Step 3.3, TCN and LSTM are combined to obtain TCN-LSTM. The specific structure of TCN-LSTM includes: (1) TCN layer: First, the input data sequence is processed by a series of TCN layers. The TCN layer uses causal convolution and dilated convolution to capture the long-term dependencies in the sequence; (2) LSTM layer: The output of the TCN layer is passed to the LSTM layer for feature fusion operation. The LSTM layer further processes the feature sequence data output by the TCN layer and uses memory units to capture complex time dependencies and sequence patterns; (3) Output layer: The output of the LSTM layer passes through a fully connected layer or other types of output layers to generate the final crystal diameter prediction result.

5. The interpretable data-driven semiconductor silicon single crystal growth diameter prediction method according to claim 4, characterized in that: The step 4 is specifically as follows: Use the GWO algorithm to optimize the hyperparameters of TCN-LSTM: The hunting process of the gray wolf is as follows: (1) Surrounding the prey In the GWO algorithm, the gray wolf uses the following position update formula to surround the prey during the hunting process: In the formula, is the distance between the gray wolf and its prey, Update the formula for the gray wolf's position, and are the position vector of the prey and the position vector of the gray wolf, respectively, and t is the current iteration number; and To determine the coefficient, the calculation formula is as follows: In the formula, and are two one-dimensional random number vectors whose components are between [0, 1]. Used to simulate the attack behavior of gray wolves on prey. The value is subject to The value of the convergence factor a decreases linearly from 2 to 0 as the number of iterations increases. (2) Chasing prey In the iterative process, α, β and δ are used to guide the movement of ω, thereby achieving global optimization, using the positions of α, β and δ Update the positions of all gray wolves using the following equation: In the formula, Respectively represent the distances of the individual ω gray wolf from the α, β and δ wolf packs: In the formula, They represent the position that the individual ω gray wolf needs to adjust under the influence of α, β and δ wolf packs, and ω takes the average value (3) Attacking prey In the following formula, t represents the current number of iterations, and T represents the maximum number of iterations. When the value of a decreases from 2 to 0, the corresponding value of A also changes in the interval [-a, a]. The larger the value of a, the farther away the gray wolf will be from the prey, hoping to find a more suitable prey, thus prompting the wolf pack to conduct a global search. If the value of a is smaller, the gray wolf will be closer to the prey, prompting the wolf pack to conduct a local search.

6. The interpretable data-driven semiconductor silicon single crystal growth diameter prediction method according to claim 5, characterized in that: The step 5 is specifically implemented according to the following steps: The SHAP algorithm is used to interpret the crystal diameter prediction results. The Shapley value of each feature is calculated through samples to measure the contribution of each feature to the model prediction results. The calculation of the Shapley value is as follows: In the formula, φ i (f) is the Shapley value of the i-th sample; N is the set of participants, including i; S is a subset of N, excluding i, and |S| is the size of the set S; represents the feature subset that does not include feature i; f(S) represents the model prediction result containing feature subset S.

Images