Control method of composite copper foil thickness control equipment

Abstract

The invention provides a control method of composite copper foil thickness control equipment, which comprises the following steps: calculating the deposition thickness of each abnormal concentration group at a future moment according to a deposition rate change rule diagram to obtain the copper layer thickness distribution of an abnormal region, and comparing the copper layer thickness distribution of the abnormal region with the thickness distribution of a normal region to obtain the thickness distribution of the abnormal region. Calculating a thickness deviation value of the abnormal area; according to the thickness deviation value, a deviation compensation algorithm is adopted, electroplating process parameters of the abnormal area are adjusted according to the size and the sign of the deviation value, and the electroplating process parameters comprise the current density and the electrolyte flow speed; and according to the optimized electroplating process parameters, performing parameter setting and adjustment on a power supply system and an electrolyte circulation system of the composite copper foil thickness control equipment, so that the composite copper foil thickness control equipment performs production according to the optimized electroplating process parameters.

Description

Technical Field

[0001] The invention relates to the field of information technology, and in particular to a control method for composite copper foil thickness control equipment. Background Art

[0002] There is a key technical problem in the production process of composite copper foil electroplating: how to accurately control the uniformity of copper layer thickness. This problem stems from the complex coupling relationship between the three physical fields of electric field distribution, electrolyte flow and copper ion concentration distribution. Specifically, the inhomogeneity of electric field distribution will lead to electric field distortion in local areas, which will affect the deposition rate of copper ions; at the same time, the flow disturbance of electrolyte will also cause the inhomogeneity of copper ion concentration distribution, which will cause local deviation of deposition thickness. These factors influence and restrict each other, making it extremely difficult to accurately control the thickness of the copper layer. More deeply, the morphology and thickness parameters of the interface layer will directly affect the electric field distribution and the flow state of the electrolyte, and then affect the mass transfer process of copper ions. Interface layers of different shapes will produce different electric field distribution and flow field characteristics, resulting in local anomalies in copper ion concentration. The deposition rate of these abnormal areas is significantly different from that of normal areas, which ultimately causes the inhomogeneity of copper layer thickness. In addition, the time evolution law of copper ion concentration is also an important issue. Since the electroplating process is a dynamically changing process, the concentration change trend in the abnormal area will directly affect the future deposition thickness. How to accurately predict this dynamic change and adjust the electroplating process parameters accordingly is the key to achieving precise thickness control. In general, a series of problems in the production process of composite copper foil electroplating, such as multi-physical field coupling effects, interface layer influence, abnormal area identification, dynamic concentration prediction, etc., together constitute the core technical difficulty of controlling the uniformity of copper layer thickness. Summary of the invention

[0003] The present invention provides a control method for composite copper foil thickness control equipment, which mainly includes:

[0004] According to the mutual coupling relationship between the electric field, flow field and mass transfer process of the electroplating process, a multi-physics field coupling model describing the composite copper foil production process is constructed, and the interface layer morphology and thickness parameters of the composite copper foil are integrated into the multi-physics field coupling model. The electric field distribution, electrolyte flow state and copper ion concentration distribution data on the interface layer with different shapes are obtained through the multi-physics field coupling model simulation.

[0005] Based on the electric field distribution data obtained by simulation, determine whether there is a local electric field distortion phenomenon. If there is a local electric field distortion, determine the position coordinates of the distortion area, and mark the electric field data in the distortion area as abnormal data. At the same time, based on the electrolyte flow state data, determine whether there is a flow disturbance phenomenon. If there is a disturbance, determine the position coordinates of the disturbance area, and mark the flow field data in the disturbance area as abnormal data.

[0006] The copper ion concentration is grouped according to the copper ion concentration distribution data to obtain copper ion groups of different concentration levels, and the central coordinates of the copper ion groups of each concentration level are calculated. According to the position coordinates of the distortion area and the disturbance area, it is determined whether the central coordinates of the copper ion groups of each concentration level are located in the abnormal area. If they are located in the abnormal area, the copper ion group of this concentration level is marked as an abnormal concentration group;

[0007] For abnormal concentration groups, a time series prediction algorithm is used to predict the concentration change at future moments according to the copper ion concentration change trend at historical moments, and the concentration change curve of each abnormal concentration group is obtained. At the same time, according to the concentration change curve, a curve fitting algorithm is used to obtain the deposition rate change law diagram of each abnormal concentration group;

[0008] According to the deposition rate variation law diagram, the deposition thickness of each abnormal concentration group at a future time is calculated to obtain the copper layer thickness distribution in the abnormal area, and the copper layer thickness distribution in the abnormal area is compared with the thickness distribution in the normal area to calculate the thickness deviation value in the abnormal area;

[0009] A deviation compensation algorithm is used for the thickness deviation value, and the electroplating process parameters of the abnormal area are adjusted according to the size and sign of the deviation value, wherein the electroplating process parameters include current density and electrolyte flow rate;

[0010] According to the optimized electroplating process parameters, the parameters of the power supply system and the electrolyte circulation system of the composite copper foil thickness control equipment are set and adjusted so that the composite copper foil thickness control equipment can be produced according to the optimized electroplating process parameters.

[0011] The technical solution provided by the embodiment of the present invention may have the following beneficial effects:

[0012] The invention discloses a control method for a composite copper foil thickness control device. The method constructs a multi-physical field coupling model to simulate the interaction of an electric field, a flow field and a mass transfer process, and integrates the interface layer morphology and thickness parameters into the model. The electric field distribution, electrolyte flow and copper ion concentration data are obtained by simulation to identify abnormal areas. For the copper ion concentration group in the abnormal area, a time series prediction and a curve fitting algorithm are used to predict the future deposition rate and thickness distribution. The thickness of the abnormal area is compared with that of the normal area, and a deviation value is calculated. According to the deviation value, a compensation algorithm is used to adjust the electroplating process parameters, including current density and electrolyte flow rate. Finally, the composite copper foil thickness control device is parameterized and adjusted, and the anode plate position angle and electrolyte flow rate pressure are automatically adjusted to make the current density distribution uniform and the electrolyte flow stable. The invention realizes precise control of the thickness in the composite copper foil production process through multi-physical field coupling analysis and intelligent algorithm optimization, and improves product quality and production efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 The present invention is a flow chart of a control method of a composite copper foil thickness control device. DETAILED DESCRIPTION

[0014] The technical solution of the present invention will be clearly and completely described below in conjunction with the embodiments. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0015] like Figure 1 The control method of the composite copper foil thickness control device of this embodiment may specifically include:

[0016] Step S101, based on the mutual coupling relationship between the electric field, flow field and mass transfer process of the electroplating process, a multi-physical field coupling model describing the composite copper foil production process is constructed, the interface layer morphology and thickness parameters of the composite copper foil are integrated into the multi-physical field coupling model, and the electric field distribution, electrolyte flow state and copper ion concentration distribution data on the interface layers of different shapes are obtained through the multi-physical field coupling model simulation.

[0017] The surface of the composite copper foil substrate is divided into quadrilateral grid units by a finite element numerical calculation method, and the interface layer thickness distribution function is calculated according to the three-dimensional coordinates of the boundary points of the grid units; a three-dimensional geometric structure digital model of the composite copper foil is established according to the interface layer thickness distribution function, and the Navier-Stokes equation is discretely solved by the finite volume method to obtain the electrolyte velocity field and pressure field distribution data near the interface layer; based on the electrolyte velocity field and pressure field distribution data, the Poisson equation is solved by the finite element method to obtain the potential distribution, and the current density at the interface layer is calculated in combination with the electrolyte conductivity distribution function; a copper ion concentration control equation is established for the current density at the interface layer, and the mass transfer equation is solved by the finite difference method in combination with the velocity field transport term and the diffusion coefficient diffusion term to obtain the ion concentration distribution data.

[0018] Specifically, the surface of the composite copper foil substrate is divided into quadrilateral grid units by finite element numerical calculation, the interface layer thickness distribution function is calculated according to the three-dimensional coordinates of the grid unit boundary points, and the three-dimensional geometric structure digital model of the composite copper foil is constructed by the substrate surface roughness parameters and the interface layer thickness data. The dynamic parameters such as electrolyte density, viscosity, and temperature are imported, and the Navier-Stokes equation is discretely solved by the finite volume method. The three-dimensional geometric structure digital model of the composite copper foil is used as the boundary condition to obtain the electrolyte velocity field and pressure field distribution data near the interface layer. Based on the three-dimensional geometric structure of the composite copper foil and the electrolyte flow field data, the finite element method is used to solve the Poisson equation to calculate the potential distribution, and the electric field intensity distribution is calculated in combination with the electrolyte conductivity distribution function, and the current density calculation expression at the interface layer is established. The current density distribution at the interface layer is used as the source term of the electrochemical reaction, and the copper ion concentration control equation is established. The mass transfer equation is solved by the finite difference method in combination with the velocity field transport term and the diffusion coefficient diffusion term. According to the electric field intensity, fluid velocity and ion concentration distribution data, the Newton iteration method is used to solve the multi-physics field coupling equations, and the steady-state solution is obtained as the multi-physics field coupling model of the composite copper foil electroplating process. The multi-physics field coupling model is substituted into different interface layer morphology parameters, and the coupled equations are solved by time advancement using the Runge-Kutta method to obtain the evolution data of each physical field over time. In the production process of composite copper foil, the key factors affecting the electroplating quality are the interface layer structure and the multi-physics field distribution. Taking the interface layer thickness distribution function h(x) as an example, when the surface of the composite copper foil substrate presents a periodic rough structure, the Fourier series expansion can be used to describe its height distribution: h(x)=h0+Σ(Ancos(2πnx / L)+Bnsin(2πnx / L)), where h0 is the average thickness, An and Bn are Fourier coefficients, L is the period length, and x is the position parameter of different interface layer morphologies. By performing least squares fitting on the measured profile data, the first five-order Fourier coefficients are obtained, and an accurate three-dimensional geometric model is established. After obtaining the geometric model, the calculation of the electrolyte flow field is particularly important. For the typical laminar flow state in the electroplating tank, the Reynolds number Re = ρvD / μ, where ρ is the electrolyte density, v is the characteristic velocity, D is the hydraulic diameter, and μ is the dynamic viscosity. The Navier-Stokes equations are discretized by the finite volume method, and the grid size is taken as 1 / 10 of the thickness of the interface layer, that is, 1 micron. The SIMPLE algorithm is used to solve the pressure-velocity coupling equation to obtain the three-dimensional velocity field distribution near the interface layer. The calculation of the electric field distribution is based on the Poisson equation of the electric potential: The electrolyte conductivity κ changes with the ion concentration c: κ = κ0 (1 + α (c-c0)), where κ0 is the standard state conductivity 5.8S / m, c0 is the standard ion concentration 0.6mol / L, and α is the concentration coefficient 0.25L / mol. A potential of -0.3V is applied to the cathode surface and 0V to the anode surface. The potential distribution φ is obtained by the finite element method, and then the electric field strength is calculated. The mass transfer process is governed by the convection-diffusion equation: The diffusion coefficient D = 6.5 × 10 -10 m 2 / s. At the interface layer, the current density and ion flux satisfy the relationship: j = nFN, n is the valence number 2, F is the Faraday constant 96485C / mol, and N is the ion flux mol / (m 2 ·s). The staggered grid finite difference format is used, the time step is 0.001s, and the space step is consistent with the flow field calculation to obtain the concentration field evolution. When solving the multi-physics field coupling equations, the relative error tolerance is selected as 10 -6 , Newton iteration method is used to solve. For a given interface layer morphology, such as a periodic groove structure with a 30-degree inclination angle and a spacing of 75 microns, calculations show that the current density at the bottom of the groove decreases by 25%, while the current density at the edge increases by 35%, which is consistent with the experimentally observed coating thickness distribution trend. The dynamic evolution characteristics of each physical field can be obtained by time advancement using the fourth-order Runge-Kutta method with a time step of 0.01s and a calculation period of 10s.

[0019] Step S102, for the electric field distribution data obtained by simulation, determine whether there is local electric field distortion. If there is local electric field distortion, determine the position coordinates of the distortion area, and mark the electric field data of the distortion area as abnormal data. At the same time, determine whether there is flow disturbance based on the electrolyte flow state data. If there is disturbance, determine the position coordinates of the disturbance area, and mark the flow field data of the disturbance area as abnormal data.

[0020] The electric field intensity gradient value is calculated according to the electric field distribution data of the grid points. If the gradient value exceeds three times the standard deviation of the local average value, a three-dimensional space clustering method is used to obtain the coordinate set of the boundary points of the distortion area. The velocity gradient tensor and the curl field are calculated from the grid point velocity vector field. If the curl value is greater than twice the laminar reference value, a boundary tracking algorithm is used to obtain the coordinate set of the boundary points of the disturbance area. Grid point matching is performed on the coordinate set of the boundary points of the distortion area and the coordinate set of the boundary points of the disturbance area to obtain the volume fraction of the overlapping area. If the volume fraction is greater than a critical value, the area is determined to be a strongly coupled anomaly area. A thermocouple array is used to obtain the temperature field distribution of the strongly coupled anomaly area, and the local conductivity is obtained according to the temperature field distribution and the conductivity temperature coefficient.

[0021] Specifically, the electric field intensity gradient value g(x,y,z) of each grid point is calculated based on the electric field distribution data. If the gradient value g(x,y,z) exceeds three times the standard deviation of the local average value, the three-dimensional spatial clustering method is used to determine the coordinate set P(x,y,z) of the boundary points of the distortion region, and the distortion region position marking matrix M1 is generated. The velocity gradient tensor and the curl field ω(x,y,z) are calculated from the velocity vector field v(x,y,z) of the fluid grid points. If the curl value is greater than twice the laminar reference value or the velocity gradient exceeds the Reynolds stress threshold, the boundary tracking algorithm is used to determine the coordinate set Q(x,y,z) of the boundary points of the disturbance region, and the disturbance region position marking matrix M2 is generated. The electric field distortion region marking matrix M1 and the flow field disturbance region marking matrix M2 are matched by grid points, and the overlapping region volume fraction α(x,y,z) is calculated. If the volume fraction α(x,y,z) is greater than the critical value, the region is marked as a strong coupling abnormal region. The temperature field distribution T(x, y, z) in the abnormal area is obtained by using a thermocouple array, and the local conductivity κ(x, y, z) is calculated in combination with the conductivity temperature coefficient to generate the temperature-compensated electric field intensity distribution E(x, y, z). According to the distorted area position marker matrix M1 and the disturbed area position marker matrix M2, the electric field intensity E(x, y, z), fluid velocity v(x, y, z) and temperature T(x, y, z) data of the abnormal area are extracted, and the Hilbert transform is used to calculate the instantaneous frequency characteristics of the abnormal data. The abnormal area data is decomposed by wavelet, the abnormal characteristic coefficients at different scales are extracted, the multi-scale abnormal characteristic vectors are generated, and the abnormal area feature database is established. In the process of electric field distortion detection, the electric field intensity gradient is calculated using the central difference format:

[0022] g(x,y,z)=[(E(x+h,y,z)-E(xh,y,z)) / 2h,(E(x,y+h,z)-E(x,yh,z)) / 2h,(E(x,y,z+h)-E(x,y,zh)) / 2h], where h is the grid spacing value of 0.1mm. For a typical electroplating tank, when the gradient value of the local area exceeds the mean plus 3 times the standard deviation, that is, |g(x,y,z)|>μg+3σg, the point is marked as a distortion point, where μg is the mean value of the gradient field 0.5kV / m 2 , σg is the standard deviation 0.15kV / m 2 In flow field disturbance judgment, the velocity gradient tensor The calculation of the curl field also uses the central difference method. Characterize the local vortex intensity. Under laminar conditions, the reference curl value is 10s-1. When the actual curl exceeds 20s-1, it indicates that there is a significant disturbance. The Reynolds stress threshold is set to 0.1Pa, corresponding to a turbulence intensity of 5%. The regional matching process uses grid overlap analysis, and the volume fraction α(x, y, z) represents the proportion of electric field distortion and flow field disturbance in a unit volume. When the α value exceeds 0.3, it indicates that there is a strong coupling anomaly in the area. In the actual electroplating process, the size of a typical strong coupling anomaly area is about 5mm×5mm×2mm. The influence of temperature field distribution on conductivity is described by a linear relationship: κ(T)=κ 0 [1+β(TT 0 )], where κ 0 is the reference temperature T 0 = The conductivity is 5.8S / m at 298K, and β is the temperature coefficient of 0.02K-1. When the temperature rises to 308K, the conductivity increases by about 20%, resulting in a corresponding decrease in the local electric field strength. The Hilbert transform is used to extract the instantaneous frequency characteristics of the abnormal signal. After transforming the electric field strength data E(t), the analytical signal z(t)=E(t)+jH[E(t)] is obtained, where H[E(t)] is the Hilbert transform result. The instantaneous frequency f(t)=(1 / 2π)·d(argz(t)) / dt reflects the frequency characteristics of the electric field fluctuations in the abnormal area, and the typical value is in the range of 0.1-10Hz. The wavelet decomposition uses the db4 wavelet basis to decompose the data in the abnormal area into 4 layers to obtain coefficients of different scales. The first layer of coefficients reflects high-frequency disturbances (>5Hz), and the fourth layer of coefficients reflects low-frequency changes (<0.5Hz). By analyzing the energy distribution of each scale coefficient, the dominant frequency component of the disturbance can be identified. When the high-frequency component energy accounts for more than 50%, it indicates that there is a violent local fluctuation. This multi-scale feature extraction method can comprehensively characterize the dynamic characteristics of the abnormal area and provide a basis for subsequent process parameter optimization.

[0023] Step S103, grouping the copper ion concentration according to the copper ion concentration distribution data to obtain copper ion groups of different concentration levels, and calculating the central coordinates of the copper ion groups at each concentration level, and judging whether the central coordinates of the copper ion groups at each concentration level are located in the abnormal area according to the position coordinates of the distortion area and the disturbance area. If located in the abnormal area, the copper ion group at that concentration level is marked as an abnormal concentration group.

[0024] A three-dimensional concentration field is constructed according to the copper ion concentration distribution data, and the three-dimensional concentration field obtains three concentration level areas of high, medium and low through a density clustering algorithm; the centroid coordinates and point density are obtained for the concentration level areas, and the group spatial distribution function and the group characteristic radius are obtained through the Gaussian kernel function of the centroid coordinates and point density; the abnormal area spatial envelope surface is constructed using the boundary point set of the electric field distortion area and the fluid disturbance area, and the shortest distance value is obtained between the abnormal area spatial envelope surface and the centroid of the concentration group; whether the ratio of the shortest distance value to the group characteristic radius is less than a critical value is judged, and if the ratio is less than the critical value, the concentration group is marked as an abnormal group and a spatial distribution characteristic matrix of the abnormal group is generated.

[0025] Specifically, a three-dimensional concentration field is constructed based on the copper ion concentration distribution data c(x, y, z), and a density clustering algorithm based on Euclidean distance is used to spatially cluster the concentration field. The minimum point number threshold n is used to divide the three concentration level areas C into high, medium and low 1 , C 2 , C 3 . Calculate the centroid coordinates r = (x, y, z) for each concentration level area C, count the point density ρ(r) in the area, and use the Gaussian kernel function to calculate the population spatial distribution function φ(r) and the population characteristic radius R. Based on the electric field distortion area D 1 and fluid disturbance area D 2 The boundary point set is used to construct the spatial envelope S(x, y, z) of the abnormal area, and the shortest distance d from the centroid of the concentration group to the envelope is calculated. The ratio λ=d / R of the shortest distance d from the centroid of the concentration group to the abnormal area and the characteristic radius R of the group is used as the judgment parameter. If λ is less than the critical value λc, the concentration group is marked as an abnormal group. For concentration groups with overlapping spatial positions, the overlapping volume ratio γ is calculated. j , if γ j If the overlap is greater than the overlap threshold γc, the overlapping areas are merged to form a new concentration group and the group characteristic parameters are updated. A coding mapping is established for the marked abnormal concentration group to generate the abnormal group spatial distribution characteristic matrix M (x, y, z). The matrix elements contain characteristic quantities such as concentration value, centroid position, and group radius. In the copper ion concentration distribution analysis, the density clustering algorithm distinguishes areas with different concentration levels by setting the concentration gradient threshold and the minimum number of points. Taking the actual electroplating tank as an example, when the concentration gradient threshold is set When the minimum point threshold n = 100, the concentration field can be divided into high concentration area (c>0.8mol / L), medium concentration area (0.4-0.8mol / L) and low concentration area (c<0.4mol / L). For the divided concentration groups, the centroid coordinates are calculated using the weighted average method: r = Σ(c j r j ) / Σc j , where cj is the concentration value at point j, r j is its position vector. The population spatial distribution function adopts the Gaussian kernel form: φ(r)=exp[-(rr) 2 / 2σ 2 ], σ is the characteristic width and takes the value of 2mm. By integrating ∫φ(r)d 3 r is calculated to obtain the characteristic radius R of the group, and the typical value is in the range of 3-5mm. In the judgment of abnormal areas, if the shortest distance d between the center of mass of a concentration group and the boundary of the abnormal area is 2mm, and its characteristic radius R is 4mm, the judgment parameter λ=0.5. The critical value λc=0.8 is set, and the concentration group is marked as an abnormal group. This judgment method based on distance ratio takes into account the spatial scale of the concentration group and avoids the deviation that may be caused by a simple distance threshold. For the case of spatial overlap, the overlapping volume ratio γj is calculated by Monte Carlo integration: points are randomly scattered in the overlapping area, and the proportion of points that are simultaneously in the two concentration groups is calculated. When γ j When >0.3, it indicates that there is significant overlap and the group needs to be merged. The characteristic parameters of the new concentration group after the merger are obtained by weighted averaging, and the weight is proportional to the number of points in the original group. The abnormal group characteristic matrix M(x,y,z) adopts a multidimensional data structure, and each matrix element contains: concentration value c(x,y,z), center of mass position r=(x,y,z), group radius R, point density ρ and other information. Through this matrix, the evolution process of the abnormal concentration group can be tracked. It is observed that in the electric field distortion area, the ion concentration often shows a local accumulation phenomenon, and the concentration gradient can reach 0.2mol / (L·mm); while in the fluid disturbance area, the ion distribution tends to be uniform, and the concentration gradient is reduced to below 0.02mol / (L·mm). This multi-scale feature analysis method can comprehensively characterize the concentration distribution characteristics of the abnormal area and provide an important basis for the optimization of process parameters.

[0026] Step S104, for the abnormal concentration group, a time series prediction algorithm is used to predict the concentration change at future moments according to the copper ion concentration change trend at historical moments, and a concentration change curve of each abnormal concentration group is obtained. At the same time, a curve fitting algorithm is used according to the concentration change curve to obtain a deposition rate change law diagram of each abnormal concentration group.

[0027] According to the historical data of abnormal concentration groups, resampling is performed at a fixed sampling interval to obtain equally spaced time series data and a concentration change rate sequence; the concentration change rate sequence is used to train a long short-term memory network model, and the concentration prediction value at the next moment is output through the network model; the concentration difference at adjacent moments is calculated for the concentration prediction value, and a theoretical deposition thickness increment sequence is obtained according to the concentration difference and Faraday's law; the theoretical deposition thickness increment sequence is subjected to Fourier transform to extract frequency components, and a polynomial fitting is performed on the frequency components to obtain a deposition rate change law function, and if the root mean square error of the polynomial fitting is greater than a preset threshold, the polynomial order is increased and the fitting process is repeated.

[0028] Specifically, the time series {t 1 ,c 1}, resample the concentration data with a fixed sampling interval δt to generate an equally spaced time series data {t,c}, and obtain the first-order concentration change rate sequence {t,Δc / Δt} through differential operation. Set the time window width τ, and use the long short-term memory network to model the concentration change rate sequence. The network input is the concentration change rate value at τ consecutive moments, and the output is the concentration prediction value at the next moment. The network parameters are optimized by the back propagation algorithm. Calculate the concentration difference at adjacent moments for the predicted concentration sequence, calculate the theoretical deposition thickness increment Δh at each moment based on Faraday's law and current density distribution, and generate the cumulative thickness sequence {t,h} over time. Perform Fourier transform on the cumulative thickness sequence, extract the main frequency component ω and its corresponding amplitude A in the amplitude spectrum, and reconstruct the deposition rate time function v(t) through inverse transformation. Use the n-order polynomial function to perform least squares fitting on the deposition rate v(t), obtain the polynomial coefficient {a}, and calculate the root mean square error σ between the fitting curve and the measured data. According to the comparison result of the root mean square error σ and the preset threshold, the polynomial order n is adjusted and the fitting process is repeated until the error meets the requirements to obtain the final deposition rate change law function f(t). In the prediction of copper ion concentration, the construction of time series is crucial to the prediction accuracy. Taking the actual electroplating process as an example, when the sampling interval δt is set to 0.5 seconds, the original concentration data points {0s, 0.6mol / L; 0.47s, 0.58mol / L; 0.91s, 0.55mol / L} are resampled to obtain an equally spaced sequence {0s, 0.6mol / L; 0.5s, 0.59mol / L; 1.0s, 0.55mol / L}, and the concentration change rates are -0.02mol / (L·s) and -0.08mol / (L·s) through differential calculation. The time window width τ of the long short-term memory network is set to 10 sampling points (i.e. 5 seconds), the input layer contains 10 neurons corresponding to the concentration change rate value in the window, the hidden layer is set to 64 neurons, and the output layer has 1 neuron to predict the concentration value at the next moment. After 500 iterative trainings, the prediction error dropped from the initial 0.15 mol / L to 0.01 mol / L. Based on Faraday's law, the calculation formula for the incremental thickness of copper ion deposition is: Δh = (M / nF·ρ)·j·Δt, where M is the molar mass of copper 63.5 g / mol, n is the valence number 2, F is the Faraday constant 96485 C / mol, and ρ is the copper density 8.9 g / cm 3 , j is the local current density. When the current density is 50mA / cm 2The theoretical deposition thickness increment within 1 second is 0.017μm. The accumulated thickness sequence is subjected to fast Fourier transform with 1024 sampling points, and the main fluctuation periods in the spectrum are 20 seconds and 60 seconds, with corresponding amplitudes of 0.05μm and 0.03μm, respectively. This indicates that there are two main components in the deposition process: short-period fluctuations and long-period drifts. The deposition rate function reconstructed by inverse transform is v(t)=0.017+0.005sin(0.314t)+0.003sin(0.105t)μm / s. The fourth-order polynomial was used for least square fitting, and the initial coefficient values ​​were {0.017, -0.001, 0.0002, -0.00001, 0.0000002}. After least square iterative optimization, the correction coefficients were {0.0172, -0.00095, 0.00018, -0.000012, 0.00000018}. The root mean square error between the fitting curve and the measured data was 0.002μm / s. When the order of the polynomial was increased to 5, the error was only reduced to 0.0019μm / s, and the improvement was not significant. Therefore, the fourth-order polynomial was determined as the final deposition rate variation function. The function reveals the nonlinear variation characteristics of the deposition rate over time: the rate drops rapidly in the initial stage (0-20 seconds), stabilizes in the middle stage (20-60 seconds), and shows a slow upward trend in the later stage (>60 seconds), which is closely related to the evolution of the electrode surface state and the dynamic changes of the ion concentration distribution. Through this multi-time scale analysis method, we can capture both the rapid fluctuation characteristics and the long-term change trends.

[0029] Step S105, calculating the deposition thickness of each abnormal concentration group at a future time according to the deposition rate variation law diagram, obtaining the copper layer thickness distribution in the abnormal area, comparing the copper layer thickness distribution in the abnormal area with the thickness distribution in the normal area, and calculating the thickness deviation value of the abnormal area.

[0030] A fourth-order Runge-Kutta integral equation is constructed according to the deposition rate variation law function, and the cumulative deposition thickness distribution function of the abnormal concentration group is obtained by the weighted average method; the radial basis function is used to interpolate and fit the copper layer thickness in the normal area, and the benchmark thickness distribution function is obtained by three-dimensional grid division; the difference between the copper layer thickness distribution function in the abnormal area and the benchmark thickness distribution function is calculated to obtain the thickness deviation function; the Gaussian kernel function is used to estimate the probability density of the thickness deviation function, and the deviation distribution map is obtained according to the kernel function weight calculation.

[0031] Specifically, the fourth-order Runge-Kutta integral equation is constructed according to the sedimentation rate variation function v(t), and the time step h is set to calculate the state quantity k 1 , k 2 , k 3 , k 4The weighted average method is used to calculate the cumulative deposition thickness H(x,y,z,t) of each abnormal concentration group. The coordinate set P(x,y,z) of the discrete points in the normal area is obtained by three-dimensional grid division, and the radial basis function is used to calculate the copper layer thickness H in the normal area. 0 (x,y,z) is interpolated and fitted to generate a continuous benchmark thickness distribution function B(x,y,z). The mean μ and standard deviation σ are calculated based on the benchmark thickness distribution function B(x,y,z) in the normal area, and the thickness fluctuation range [μ-3σ,μ+3σ] in the normal area is established using the triple standard deviation principle as the benchmark standard. The difference between the copper layer thickness distribution H(x,y,z,t) in the abnormal area and the benchmark thickness distribution function B(x,y,z) is calculated to generate a thickness deviation function D(x,y,z,t)=H(x,y,z,t)-B(x,y,z). The Gaussian kernel function is used to estimate the probability density of the thickness deviation value, and the bandwidth parameter w is selected to calculate the kernel function weight to generate the deviation value probability density function P(D). The deviation value interval is divided according to the probability density function P(D), and the spatial distribution characteristics of the deviation points in each interval are statistically analyzed. The contour method is used to generate the deviation distribution map G(x,y,z). In the process of calculating the copper layer deposition thickness, the fourth-order Runge-Kutta method is used to numerically integrate the deposition rate function. Taking an actual case as an example, given the deposition rate function v(t) = 0.015 + 0.002sin(0.1t) μm / s, the time step h = 0.1s, the state quantity calculation formula is: k 1 = h·v(t), k 2 =h·v(t+h / 2), k 3 =h·v(t+h / 2), k 4 =h·v(t+h), the final integral result H=(k 1 +2k 2 +2k 3 +k 4 ) / 6. After calculation, the cumulative thickness reaches 0.152μm at t=10s. The normal area benchmark thickness distribution is interpolated using a multi-quadratic radial basis function, and the function form is B(r)=(1+(εr) 2 )^(1 / 2), where r is the spatial distance and ε is the shape parameter, which is 0.1. For the typical normal area measurement data, the mean μ=0.15μm and the standard deviation σ=0.005μm are calculated, and the reference range [0.135μm, 0.165μm] is established. In the thickness deviation calculation, local thickness mutations occurred in the abnormal area. Taking an abnormal concentration group as an example, the thickness at the center reached 0.18μm, which is 0.03μm different from the reference value, far beyond the normal fluctuation range. Through the Gaussian kernel function K(x)=exp(-x 2 / 2h2) is used for probability density estimation, and the kernel function bandwidth parameter h is determined to be 0.0075μm using the Silverman criterion. The probability density distribution of the deviation value shows an obvious bimodal feature, with the main peak at 0.02μm and the secondary peak at -0.015μm, indicating that thickening and thinning occur simultaneously in the abnormal area. In terms of spatial distribution, the thickened area is mainly concentrated at the edge of the electric field distortion area, with the maximum deviation reaching 0.035μm; while the thinned area appears at the center of the fluid disturbance, with the maximum deviation being -0.025μm. The deviation distribution map is drawn by the contour line method, and the contour line spacing is set to 0.005μm. Projection on the xy plane shows that the thickened area is elliptical in shape, with a major axis of about 5mm and a minor axis of about 3mm; the thinned area is irregularly circular with a diameter of about 4mm. There is a transition zone between the two areas, with a width of about 1mm, and the deviation value gradually transitions from -0.01μm to 0.01μm. The z-direction profile shows that the deviation value decays exponentially with depth, tends to be stable at 50μm from the surface, and reaches the reference thickness level. This three-dimensional distribution feature reflects the superposition of electric field effect and fluid action, providing a quantitative basis for subsequent process optimization.

[0032] Step S106, using a deviation compensation algorithm for the thickness deviation value, and adjusting the electroplating process parameters of the abnormal area according to the size and sign of the deviation value, the electroplating process parameters including current density and electrolyte flow rate.

[0033] A proportional-integral compensation equation is established according to the thickness deviation value, and the current density compensation coefficient and the flow velocity compensation coefficient are obtained by integrating the deviation between the target value and the actual value; a piecewise linear feedback function is established for the current density compensation coefficient, and if the absolute value of the compensation coefficient is less than a preset threshold, linear compensation is used to obtain a corrected current density distribution; an exponential nonlinear mapping function is used to transform the flow velocity compensation coefficient, and the corrected flow velocity distribution is obtained by calculating the preset mapping coefficient and the original flow velocity value; a compensation parameter evaluation function including the squared thickness deviation term, the current density gradient term and the flow velocity gradient term is constructed, and a deep neural network is used to optimize the compensation parameters to obtain the optimized compensation coefficient.

[0034] Specifically, a proportional-integral compensation equation is established based on the thickness deviation value Δh(x,y,z): u(t)=Kp·e(t)+Ki∫e(t)dt, where u(t) is the output of the controller, Kp and Ki are compensation coefficients, e(t) is the deviation between the target value and the actual value, and t is time. The current density compensation coefficient α(x,y,z) and the flow rate compensation coefficient β(x,y,z) are calculated. A piecewise linear feedback function is established for the current density compensation coefficient α(x,y,z). If |α| is less than the threshold α 0 , linear compensation is used; if |α| is greater than α 0, saturation compensation is used to generate the corrected current density distribution j′(x,y,z). The velocity compensation coefficient β(x,y,z) is converted using an exponential nonlinear mapping function: v′=v·exp(λβ), where λ is the mapping coefficient, and the corrected velocity distribution v′(x,y,z) is calculated. Construct the compensation parameter evaluation function:

[0035] Where J(α,β) is the optimization target value of the compensation parameter, Δh is the thickness deviation, j is the current density distribution, v is the flow velocity distribution, and w 1 、w 2 、w 3 is the weight coefficient, and the optimization target value of the compensation parameter is calculated. A five-layer deep neural network is used to optimize the compensation parameters. The input layer is the deviation value and the original parameter, the hidden layer uses the ReLU activation function, and the output layer is the optimized compensation coefficient. The process parameters are recalculated according to the optimized compensation coefficient, and the compensation effect is verified by numerical simulation to generate the compensation parameter verification index R(x, y, z). In the process of electroplating process parameter compensation, the deviation compensation is calculated using the proportional integral algorithm. When the thickness deviation Δh=5μm at a certain point is measured, the proportional coefficient Kp=0.2 and the integral coefficient Ki=0.05 are set, and the corresponding compensation output is u(t)=0.2×5+0.05×∫5dtμm. The integral term accumulates over time, and the compensation output gradually increases from 1μm to 2.5μm within 10 seconds. The current density compensation is processed by a piecewise function, and the threshold α is set. 0 =0.3, when the compensation coefficient α=0.2, it is in the linear interval, and the corrected current density j′=j×(1+0.2); when α=0.4, it exceeds the threshold and enters the saturation zone, and the saturation value j′=j×1.3 is taken. This segmented processing method avoids over-compensation. The flow rate compensation adopts an exponential mapping relationship, and the mapping coefficient λ=0.5. When β=0.4, the corrected flow rate v′=v×exp(0.5×0.4)=1.22v, showing a nonlinear growth characteristic. This exponential relationship responds smoothly in the area of ​​small compensation coefficients and responds quickly in the area of ​​large compensation coefficients. The weight coefficient in the evaluation function reflects the importance of different optimization objectives, and the typical value is w 1 =0.5, w 2 =0.3, w 3=0.2. Taking a certain abnormal area as an example, the evaluation value before compensation is J=2.5, of which the thickness deviation term contributes 1.5, the current gradient term contributes 0.6, and the flow velocity gradient term contributes 0.4. The deep neural network adopts a five-layer structure. The 8 neurons in the input layer receive the deviation value and the original parameters. The three hidden layers contain 32, 16, and 8 neurons respectively. The 2 neurons in the output layer output the optimized compensation coefficient. The ReLU activation function f(x)=max(0,x) is used to process the hidden layer output to avoid the gradient disappearance problem. In the verification stage, the effect of the compensation parameters is evaluated by the indicator R(x,y,z). The pre-compensation indicator R=0.85 of a certain test point was increased to R=0.92 after a round of optimization. Analysis of the optimization process found that current density compensation dominated the improvement of thickness uniformity, while flow velocity compensation significantly reduced local disturbances. The synergistic effect of the two compensation mechanisms is reflected in multiple components of the evaluation function: the sum of squares of thickness deviations is reduced by 45%, the sum of squares of current gradients is reduced by 35%, and the sum of squares of flow velocity gradients is reduced by 30%. This multi-objective optimization strategy not only improves the deposition uniformity but also maintains a smooth distribution of process parameters.

[0036] Step S107, according to the optimized electroplating process parameters, the power supply system and the electrolyte circulation system of the composite copper foil thickness control device are parameterized and adjusted so that the composite copper foil thickness control device can be produced according to the optimized electroplating process parameters.

[0037] The theoretical position coordinates of the anode plate in the target area are calculated according to the current density distribution function, the actual position of the anode plate is measured by the spatial distance sensor array, and the position deviation compensation amount is obtained by the proportional integral controller; the distance change rate matrix between the anode plate and the cathode plate is obtained by the laser displacement sensor array, the anode plate inclination is calculated according to the distance change rate matrix, and the angle adjustment signal is output through the angle closed-loop controller; if the current density uniformity index is less than the preset threshold, the anode plate position compensation mechanism is driven to adjust until the uniformity index is greater than the preset threshold; the electrolyte flow rate set value is calculated according to the flow field distribution function, the actual flow rate value is measured by the flow sensor, the valve opening is adjusted by the pressure closed-loop controller, and the flow ratio of the inlet pipeline and the return pipeline is adjusted according to the valve opening.

[0038] Specifically, the theoretical position coordinates (x 0 ,y 0 ,z 0), measure the actual position (x, y, z) of the anode plate through a spatial distance sensor array, and use a proportional-integral controller to calculate the position deviation compensation Δr. Use a laser displacement sensor array to measure the distance change rate between the anode plate and the cathode plate in real time, calculate the optimal inclination angle θ(x, y) of the anode plate based on the distance change rate matrix D(x, y, z), and output the adjustment signal through an angle closed-loop controller. Calculate the anode plate fine-tuning parameters based on the current density uniformity index. If the uniformity index is less than the preset threshold η0, drive the anode plate position compensation mechanism for fine adjustment until the uniformity index is greater than the threshold η0. Use the optimized flow field distribution function v(x, y, z) to calculate the electrolyte flow rate setting value v 0 , the actual flow velocity v is measured by the flow sensor, and according to the deviation Δv=v 0 -v adjusts the speed of the circulating pump n. The electrolyte pipeline pressure p(t) is collected in real time, and the pressure signal is filtered by a third-order Butterworth low-pass filter. The valve opening compensation φ is calculated by a pressure closed-loop controller. The flow ratio γ between the inlet pipeline and the return pipeline is adjusted according to the valve opening compensation φ. The liquid level h(t) of the electroplating tank is monitored by a liquid level sensor to establish a liquid level stability control loop. In the process of adjusting the parameters of the electroplating equipment, the anode plate position control adopts a multi-level closed-loop feedback mechanism. When the theoretically calculated optimal position coordinates of the anode plate are (150mm, 200mm, 80mm), the actual position measured by the spatial distance sensor is (152mm, 198mm, 82mm), and the position deviation vector Δr=(2,-2,2)mm. A proportional-integral controller is used, and the proportional coefficient Kp=0.8 and the integral coefficient Ki=0.2 are set. The position compensation is calculated and the servo motor is driven for fine-tuning. During the adjustment of the anode plate inclination angle, the laser displacement sensor array is arranged around the cathode plate, and the sampling frequency is 100Hz. When the measured distance change rate matrix shows that the right distance increase rate is 0.5mm / s and the left distance decrease rate is -0.4mm / s, it indicates that the anode plate has a tendency to tilt to the right. The optimal inclination angle θ = 2.5 degrees is obtained through matrix calculation, and the angle controller outputs the corresponding stepper motor pulse signal. The current density uniformity index η is calculated using the ratio of the standard deviation to the average value, and the preset threshold η 0 =0.05. The measured data shows that the uniformity index of a certain area is η=0.08, which exceeds the threshold range. At this time, the precision position compensation mechanism is activated and fine-tuned with a step size of 0.1mm. After 3 iterations, the uniformity index drops to 0.045. In terms of flow field control, the optimized flow velocity distribution function gives the set value v 0=0.15m / s, the flow sensor measures the actual flow velocity v=0.12m / s. The PI controller calculates the circulating pump speed increment Δn=300rpm, and the frequency converter adjusts the output frequency accordingly. At the same time, the pressure sensor monitors the pipeline pressure fluctuation, which contains 50Hz power frequency interference. The pressure signal is processed by a third-order Butterworth filter with a cutoff frequency of 10Hz, and a stable pressure reading p=0.2MPa is obtained after filtering out the high-frequency noise. The liquid level control adopts a dual-pipeline balancing strategy, and the flow ratio γ between the inlet pipeline and the return pipeline is achieved by adjusting the valve opening. When the inlet valve opening is 60%, the return valve opening should be adjusted to 55% to keep the liquid level at the target height h 0 =500mm. If the liquid level fluctuates by ±10mm, the level controller will automatically adjust the valve opening, and the compensation amount φ is determined by the deviation value and the rate of change. This multi-parameter coupling control scheme achieves a dynamic balance of process parameters, and each control loop cooperates with each other to maintain the stable operation of the electroplating process.

[0039] The preferred embodiments of the present invention disclosed above are only used to help explain the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to only specific implementation methods. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. A control method for a composite copper foil thickness control device, characterized in that: The method comprises: According to the mutual coupling relationship between the electric field, flow field and mass transfer process of the electroplating process, a multi-physics field coupling model describing the composite copper foil production process is constructed, and the interface layer morphology and thickness parameters of the composite copper foil are integrated into the multi-physics field coupling model. The electric field distribution, electrolyte flow state and copper ion concentration distribution data on the interface layer with different shapes are obtained through the multi-physics field coupling model simulation. Based on the electric field distribution data obtained by simulation, determine whether there is a local electric field distortion phenomenon. If there is a local electric field distortion, determine the position coordinates of the distortion area, and mark the electric field data in the distortion area as abnormal data. At the same time, based on the electrolyte flow state data, determine whether there is a flow disturbance phenomenon. If there is a disturbance, determine the position coordinates of the disturbance area, and mark the flow field data in the disturbance area as abnormal data. The copper ion concentration is grouped according to the copper ion concentration distribution data to obtain copper ion groups of different concentration levels, and the central coordinates of the copper ion groups of each concentration level are calculated. According to the position coordinates of the distortion area and the disturbance area, it is determined whether the central coordinates of the copper ion groups of each concentration level are located in the abnormal area. If they are located in the abnormal area, the copper ion group of this concentration level is marked as an abnormal concentration group; For abnormal concentration groups, a time series prediction algorithm is used to predict the concentration change at future moments according to the copper ion concentration change trend at historical moments, and the concentration change curve of each abnormal concentration group is obtained. At the same time, according to the concentration change curve, a curve fitting algorithm is used to obtain the deposition rate change law diagram of each abnormal concentration group; According to the deposition rate variation law diagram, the deposition thickness of each abnormal concentration group at a future time is calculated to obtain the copper layer thickness distribution in the abnormal area, and the copper layer thickness distribution in the abnormal area is compared with the thickness distribution in the normal area to calculate the thickness deviation value in the abnormal area; A deviation compensation algorithm is used for the thickness deviation value, and the electroplating process parameters of the abnormal area are adjusted according to the size and sign of the deviation value, wherein the electroplating process parameters include current density and electrolyte flow rate; According to the optimized electroplating process parameters, the parameters of the power supply system and the electrolyte circulation system of the composite copper foil thickness control equipment are set and adjusted so that the composite copper foil thickness control equipment can be produced according to the optimized electroplating process parameters.

2. The method according to claim 1, characterized in that According to the mutual coupling relationship between the electric field, flow field and mass transfer process of the electroplating process, a multi-physics field coupling model describing the composite copper foil production process is constructed, the interface layer morphology and thickness parameters of the composite copper foil are integrated into the multi-physics field coupling model, and the electric field distribution, electrolyte flow state and copper ion concentration distribution data on the interface layer of different shapes are obtained through the multi-physics field coupling model simulation, including: The finite element numerical calculation method is used to divide the surface of the composite copper foil substrate into quadrilateral grid units, and the interface layer thickness distribution function is calculated according to the three-dimensional coordinates of the grid unit boundary points; A three-dimensional geometric structure digital model of the composite copper foil is established according to the interface layer thickness distribution function, and the Navier-Stokes equation is discretely solved by the finite volume method to obtain the distribution data of the electrolyte velocity field and pressure field near the interface layer; According to the velocity field and pressure field distribution data of the electrolyte, the finite element method is used to solve the Poisson equation to obtain the potential distribution, and the current density at the interface layer is calculated in combination with the electrolyte conductivity distribution function; According to the current density at the interface layer, a copper ion concentration control equation is established, and the mass transfer equation is solved by the finite difference method in combination with the velocity field transport term and the diffusion coefficient diffusion term to obtain the ion concentration distribution data.

3. The method according to claim 1, characterized in that The method comprises: judging whether there is a local electric field distortion phenomenon for the electric field distribution data obtained by simulation, and if there is a local electric field distortion, determining the position coordinates of the distortion area, and marking the electric field data of the distortion area as abnormal data; and judging whether there is a flow disturbance phenomenon according to the electrolyte flow state data, and if there is a disturbance, determining the position coordinates of the disturbance area, and marking the flow field data of the disturbance area as abnormal data, including: Calculate the electric field intensity gradient value according to the electric field distribution data of the grid points, and if the gradient value exceeds three times the standard deviation of the local average value, use a three-dimensional spatial clustering method to obtain a set of coordinates of the boundary points of the distortion area; The velocity gradient tensor and curl field are calculated from the velocity vector field of the grid points. If the curl value is greater than twice the laminar reference value, the boundary tracking algorithm is used to obtain the coordinate set of the boundary points of the disturbance area. Grid point matching is performed on the coordinate set of the boundary points of the distortion region and the coordinate set of the boundary points of the disturbance region to obtain a volume fraction of the overlapping region. If the volume fraction is greater than a critical value, the region is determined to be a strongly coupled abnormal region. The temperature field distribution in the strongly coupled abnormal region is obtained by using a thermocouple array, and the local conductivity is obtained according to the temperature field distribution and the conductivity temperature coefficient.

4. The method according to claim 1, characterized in that: The method comprises: grouping the copper ion concentration according to the copper ion concentration distribution data to obtain copper ion groups of different concentration levels, calculating the central coordinates of the copper ion groups at each concentration level, and judging whether the central coordinates of the copper ion groups at each concentration level are located in an abnormal area according to the position coordinates of the distortion area and the disturbance area. If the central coordinates of the copper ion groups at each concentration level are located in the abnormal area, marking the copper ion group at the concentration level as an abnormal concentration group comprises: Constructing a three-dimensional concentration field according to the copper ion concentration distribution data, wherein the three-dimensional concentration field obtains three concentration level areas of high, medium and low through a density clustering algorithm; Obtaining centroid coordinates and point density for the concentration level area, and obtaining a population spatial distribution function and a population characteristic radius through a Gaussian kernel function using the centroid coordinates and point density; The boundary point set of the electric field distortion region and the fluid disturbance region is used to construct the abnormal region space envelope surface, and the abnormal region space envelope surface and the concentration group centroid obtain the shortest distance value; It is determined whether the ratio of the shortest distance value to the group characteristic radius is less than a critical value. If the ratio is less than the critical value, the concentration group is marked as an abnormal group and a spatial distribution characteristic matrix of the abnormal group is generated.

5. The method according to claim 1, characterized in that For the abnormal concentration group, a time series prediction algorithm is used to predict the concentration change at a future moment according to the copper ion concentration change trend at a historical moment, and a concentration change curve of each abnormal concentration group is obtained. At the same time, according to the concentration change curve, a curve fitting algorithm is used to obtain a deposition rate change law diagram of each abnormal concentration group, including: According to the historical data of abnormal concentration groups, resampling is performed with a fixed sampling interval to obtain equal-interval time series data and concentration change rate series; The concentration change rate sequence is used to train a long short-term memory network model, and a predicted concentration value at the next moment is output through the network model; Calculating the concentration difference at adjacent moments according to the concentration prediction value, and obtaining a theoretical deposition thickness increment sequence according to the concentration difference and Faraday's law; The theoretical deposition thickness increment sequence is subjected to Fourier transform to extract frequency components, and a deposition rate variation law function is obtained by performing polynomial fitting on the frequency components. If the root mean square error of the polynomial fitting is greater than a preset threshold, the polynomial order is increased and the fitting process is repeated.

6. The method according to claim 1, characterized in that The method of calculating the deposition thickness of each abnormal concentration group at a future time according to the deposition rate variation law diagram, obtaining the copper layer thickness distribution in the abnormal area, comparing the copper layer thickness distribution in the abnormal area with the thickness distribution in the normal area, and calculating the thickness deviation value in the abnormal area includes: The fourth-order Runge-Kutta integral equation is constructed according to the sedimentation rate variation law function, and the cumulative sedimentation thickness distribution function of the abnormal concentration group is obtained by the weighted average method. The radial basis function is used to interpolate and fit the copper layer thickness in the normal area, and the reference thickness distribution function is obtained through three-dimensional grid division; Performing difference calculation on the copper layer thickness distribution function in the abnormal area and the reference thickness distribution function to obtain a thickness deviation function; A Gaussian kernel function is used to perform probability density estimation on the thickness deviation function, and a deviation distribution map is obtained according to kernel function weight calculation.

7. The method according to claim 1, characterized in that The deviation compensation algorithm is used for the thickness deviation value, and the electroplating process parameters of the abnormal area are adjusted according to the size and sign of the deviation value. The electroplating process parameters include current density and electrolyte flow rate, including: A proportional-integral compensation equation is established according to the thickness deviation value, and a current density compensation coefficient and a flow velocity compensation coefficient are obtained by integrating the deviation between the target value and the actual value; A piecewise linear feedback function is established for the current density compensation coefficient, and if the absolute value of the compensation coefficient is less than a preset threshold, linear compensation is used to obtain a corrected current density distribution; The flow velocity compensation coefficient is transformed by using an exponential nonlinear mapping function, and a modified flow velocity distribution is obtained by calculating the preset mapping coefficient and the original flow velocity value; A compensation parameter evaluation function including the squared thickness deviation term, the current density gradient term and the flow velocity gradient term is constructed, and a deep neural network is used to optimize the compensation parameters to obtain the optimized compensation coefficient.

8. The method according to claim 7, characterized in that The proportional-integral compensation equation is: u(t)=Kp·e(t)+Ki∫e(t)dt Where u(t) is the output of the controller, Kp and Ki are compensation coefficients, e(t) is the deviation between the target value and the actual value, and t is time.

9. The method according to claim 7, characterized in that: The compensation parameter evaluation function is: Where J(α, β) is the optimization target value of the compensation parameter, Δh is the thickness deviation, j is the current density distribution, v is the flow velocity distribution, and w1, w2, and w3 are weight coefficients.

10. The method according to claim 1, characterized in that According to the optimized electroplating process parameters, the power supply system and the electrolyte circulation system of the composite copper foil thickness control device are parameterized and adjusted so that the composite copper foil thickness control device is produced according to the optimized electroplating process parameters, including: The theoretical position coordinates of the anode plate in the target area are calculated according to the current density distribution function, the actual position of the anode plate is measured by the spatial distance sensor array, and the position deviation compensation is obtained by using a proportional integral controller; A laser displacement sensor array is used to obtain a distance change rate matrix between the anode plate and the cathode plate, the anode plate inclination angle is calculated according to the distance change rate matrix, and an angle adjustment signal is output through an angle closed-loop controller; If the current density uniformity index is less than a preset threshold, the anode plate position compensation mechanism is driven to adjust until the uniformity index is greater than the preset threshold; The electrolyte flow rate setting value is calculated according to the flow field distribution function, the actual flow rate value is measured by a flow sensor, the valve opening is adjusted by a pressure closed-loop controller, and the flow ratio of the liquid inlet pipeline and the liquid return pipeline is adjusted according to the valve opening.

Images