IGBT module life prediction method based on multi-condition load spectrum

Abstract

The application discloses an IGBT module life prediction method based on a multi-working condition load spectrum, and relates to the technical field of IGBT module life prediction. The method comprises the following steps: collecting the transient junction temperature sequence and the collector current waveform of the IGBT module under multiple working conditions in real time through an embedded sensor network, and constructing an original multi-working condition load spectrum; performing phase space reconstruction on the original multi-working condition load spectrum based on chaos theory, extracting load attractor geometric features, and using a stochastic resonance model to enhance weak fatigue feature components to generate an enhanced load feature set; inputting the enhanced load feature set into a sparse representation classifier, separating out damage sparse coefficients associated with the bond wire failure mechanism through dictionary learning, and determining cumulative damage primitives; adopting a graph neural network to map the cumulative damage primitives to a stress space-time evolution graph of each material layer interface, and outputting a fatigue crack propagation rate field; aligning the fatigue crack propagation rate field to a reference life database, constructing a residual life distribution probability graph, and outputting a predicted residual effective life.

Description

Technical Field

[0001] This invention relates to the field of IGBT module life prediction technology, specifically to an IGBT module life prediction method based on multi-condition load spectrum. Background Technology

[0002] IGBT modules in power electronic systems are subjected to alternating loads under multiple operating conditions for extended periods. Bond wire failure and solder layer fatigue are the main mechanisms for terminating IGBT module lifespan. Existing IGBT module lifespan prediction methods largely rely on empirical lifespan models under a single operating condition or cumulative damage methods based on cycle counting. Conventional methods process junction temperature fluctuation signals at the time-domain statistical or frequency-domain analysis level, treating junction temperature changes as linear or quasi-steady-state processes, ignoring the inherent nonlinear chaotic dynamics of the junction temperature sequence during multi-condition switching. This approach fails to extract the attractor geometry formed by load fluctuations in high-dimensional phase space, resulting in insensitivity to subtle fatigue characteristics and significantly reduced prediction accuracy under multi-condition coupling.

[0003] Existing technologies lack effective utilization of the structural correspondence between load characteristics and physical failure mechanisms in the life prediction process. Conventional methods typically input extracted load characteristics directly into data-driven models for life mapping, failing to clearly separate the damage contributions of different failure modes. In life mapping, existing solutions either perform statistical regression in an abstract feature space, detaching from the actual physical process of stress evolution in each material layer within the IGBT module; or rely on simplified finite element models that are ill-suited to the complex temporal evolution of multi-condition loads; or directly apply models trained on historical data in cross-condition scenarios without considering data distribution shifts. These combined problems make it difficult for existing methods to accurately capture the entire process from load input to internal damage evolution under complex multi-condition conditions, thus limiting the reliability and generalization ability of IGBT module remaining life prediction. Summary of the Invention

[0004] This invention provides an IGBT module life prediction method based on multi-condition load spectrum. It extracts nonlinear fatigue characteristics from the transient junction temperature sequence under multiple conditions and establishes a physical mapping from load characteristics to internal stress evolution, aiming to solve the defects of existing IGBT module life prediction methods.

[0005] To achieve the above objectives, this invention provides the following technical solution: This invention discloses a method for predicting the lifetime of IGBT modules based on multi-condition load spectra. This method uses an embedded sensor network to collect transient junction temperature sequences and collector current waveforms of IGBT modules under multiple operating conditions in real time, constructing an original multi-condition load spectrum. For the constructed original multi-condition load spectrum, phase space reconstruction is performed based on chaos theory, extracting the geometric features of the load attractor. A stochastic resonance model is then used to enhance the weak fatigue characteristic components in these geometric features, generating an enhanced load feature set. Preferably, in this process, a time-delay embedding operation is first performed on the transient junction temperature sequence, calculating its correlation dimension and Lyapunov exponent, thereby determining the embedding dimension and time delay parameters required for phase space reconstruction. Then, the transient junction temperature sequence is projected into the corresponding high-dimensional phase space to form a load phase trajectory. By calculating the local neighborhood radius and phase point density of each phase point in the phase trajectory, the geometric features of the load attractor, including the attractor dimension, Kolmogorov entropy, and the maximum Lyapunov exponent, are extracted. After obtaining these features, they are input into a bistable stochastic resonance model, which exhibits stochastic resonance characteristics determined by the barrier height and noise intensity parameters. By adjusting these two parameters, the model is made to operate in a stochastic resonance state, thereby significantly enhancing the weak fatigue feature components hidden in the geometric features of the load attractor through signal resonance. The obtained resonant enhanced signal is then subjected to wavelet packet decomposition, and the energy proportion and peak frequency of each decomposition node are calculated. Only node signals with energy proportions exceeding the global mean and peak frequencies falling within the fatigue-sensitive frequency band are selected to form the enhanced load feature set.

[0006] After obtaining the enhanced load feature set, it is input into a sparse representation classifier. Through dictionary learning, damage sparse coefficients associated with the bond line failure mechanism of the IGBT module are separated, and the cumulative damage primitives are determined based on the energy distribution of these coefficients. As a specific technical solution of this invention, a portion of the enhanced load feature set is allocated as a training dictionary, and the remainder as test samples. The training dictionary is iteratively optimized using K-singular value decomposition to learn an overcomplete dictionary, where each atom corresponds to a specific load fluctuation mode. The test samples are input into the sparse representation classifier, and the orthogonal matching pursuit algorithm is used to sparsely decompose the overcomplete dictionary to obtain a sparse coefficient vector. For the bond line failure mechanisms of the IGBT module, namely bond line peeling, aluminum reconstruction, and solder fatigue, corresponding failure mechanism labels are pre-assigned to each atom in the overcomplete dictionary, and corresponding atom index intervals are set, thereby accurately selecting damage sparse coefficients belonging to these three mechanisms from the sparse coefficient vector. The Lp norm energy distribution of the above damage sparsity coefficients is calculated. The atomic index and coefficient magnitude corresponding to the energy peak are extracted by local peak detection, and the product sequence of the two is used as the cumulative damage primitive.

[0007] Subsequently, a graph neural network is used to map the accumulated damage primitives to the spatiotemporal evolution map of interface stress in each material layer within the IGBT module, thereby outputting the fatigue crack propagation rate field for each material layer. Preferably, a topological graph structure of each material layer within the IGBT module is first constructed, defining discrete units on the interfaces of different material layers as nodes, and the thermal coupling relationship between adjacent discrete units as edges. Node attributes include the temperature gradient and equivalent stress of the corresponding unit. The accumulated damage primitives are fed into a graph convolutional network as temporal input features. This network performs layer-by-layer convolution operations on the topological graph structure, aggregating the attribute features of neighboring nodes in each convolution to update the hidden state vector of the node. After each graph convolution operation, a graph-gated recurrent unit is applied to capture the memory and forgetting characteristics of the node state in the time dimension, thereby outputting the interface stress prediction value of each node in continuous time steps. The interface stress prediction values ​​of all nodes are arranged according to the material layer order and spatial coordinates, thus forming the spatiotemporal evolution map of interface stress. By performing edge detection and region growth segmentation on this map, the fatigue crack front position and propagation direction of each material layer are accurately extracted, ultimately generating the fatigue crack propagation rate field for each material layer.

[0008] Based on this, transfer learning techniques are used to align the generated fatigue crack propagation rate field to a pre-defined baseline lifetime database. A probability map of the remaining lifetime distribution is constructed by analyzing the difference between the rate fields before and after alignment, and the predicted remaining effective lifetime of the IGBT module is ultimately output based on this. As a preferred implementation of this step, the baseline lifetime database stores a large number of historical IGBT module fatigue crack propagation rate field sequences under different operating conditions and their corresponding actual lifetime labels. A domain-adaptive transfer learning method is used to calculate the maximum mean difference distance between the current IGBT module's fatigue crack propagation rate field and each historical rate field sequence in the baseline database, and the marginal distribution of the current rate field is aligned to the source domain distribution by minimizing this distance. In the aligned common feature space, the difference between the current rate field and the corresponding rate field in the baseline database at each spatial location is calculated, forming a rate field difference tensor. This tensor is decomposed to obtain the main difference pattern and the residual difference pattern. These two patterns are concatenated and input into a Bayesian deep regression network. The weight parameters of this network are assigned a Gaussian prior distribution, and Monte Carlo dropout layers are inserted between fully connected layers as Bayesian approximate inference modules. By randomly discarding some neurons during multiple forward propagation processes, multiple output activation values ​​are obtained. The arithmetic mean of these activation values ​​is calculated as the predicted mean of the remaining lifetime, and its standard deviation is calculated as the predicted variance. Using the predicted mean as the center and the predicted variance as the scale, a probability density function of the remaining lifetime on the time axis is constructed under the Gaussian distribution assumption, and this function is discretized to obtain a probability map of the remaining lifetime distribution. Finally, the lifetime time coordinate corresponding to the maximum probability density value in this probability map is retrieved and output as the predicted remaining effective lifetime of the IGBT module. Preferably, the aforementioned embedded sensing network for synchronously acquiring transient junction temperature sequences and collector current waveforms consists of multiple thermocouple sensors and Rogowski coil sensors integrated inside the IGBT module.

[0009] The technical effects and advantages provided by the present invention in the above technical solution are as follows: Phase space reconstruction of multi-condition load spectra using chaos theory is employed to project transient junction temperature sequences onto a high-dimensional phase space, forming load phase trajectories. Geometric features of the load attractors, including attractor dimension, Kolmogorov entropy, and the maximum Lyapunov exponent, are extracted. A bistable stochastic resonance model is used to enhance weak fatigue characteristic components within these geometric features. By adjusting the barrier height and noise intensity parameters, the system enters a stochastic resonance state, transferring noise energy to the weak characteristic components and significantly amplifying fatigue features previously submerged in the noise background. Compared to existing techniques that directly perform statistical or spectral analysis on time-domain signals, this method fully utilizes the deterministic chaotic structure inherent in the load signal, transforming the system's nonlinear dynamic behavior into a feature space usable for fatigue state identification. It can effectively capture weak signal components closely related to fatigue damage in complex, noisy environments under multiple conditions, overcoming the insensitivity of conventional methods to weak fatigue features.

[0010] By combining a sparse representation classifier with dictionary learning techniques, the enhanced load feature set is sparsely decomposed on an overcomplete dictionary. Each dictionary atom corresponds to a load fluctuation mode and is pre-labeled with failure mechanisms such as bond line peeling, aluminum reconstruction, and solder fatigue. The damage sparse coefficients corresponding to each failure mechanism are directly separated from the sparse coefficient vector according to the atom index interval, and the cumulative damage primitives are determined based on the energy distribution of these coefficients. Furthermore, a graph neural network is used to build the topological graph structure of each material layer inside the IGBT module. The cumulative damage primitives are used as time-series inputs and mapped to the spatiotemporal evolution diagram of interface stress, outputting the fatigue crack propagation rate field of each material layer. In the transfer learning stage, the edge distribution of the crack propagation rate field is aligned to the source domain distribution of the baseline lifetime database through a domain adaptation method. In the common feature space of the aligned distribution, a Bayesian deep regression network is used to calculate the predicted mean and predicted variance of the remaining lifetime. This combined approach connects the complete causal link from load signal to internal physical damage and then to cross-condition life inference. While maintaining the separability of damage contributions of each failure mode, it uses transfer learning to eliminate the domain bias problem caused by the difference in data distribution under different conditions, so that life prediction can still maintain physical consistency when the conditions change, and output a life estimate with uncertainty measurement. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0012] Figure 1 This is a flowchart of an IGBT module lifetime prediction method based on multi-condition load spectrum; Figure 2 This is a flowchart of the load feature enhancement process based on phase space reconstruction and bistable random resonance; Figure 3 This is a flowchart of the method for generating fatigue crack propagation rate field in IGBT modules. Detailed Implementation

[0013] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0014] See Figure 1 This invention provides a method for predicting the lifetime of an IGBT module based on a multi-condition load spectrum, comprising the following steps: Real-time acquisition of transient junction temperature sequences and collector current waveforms of the IGBT module under multiple conditions using an embedded sensor network to construct an original multi-condition load spectrum; Phase space reconstruction of the original multi-condition load spectrum based on chaos theory, extraction of load attractor geometric features, and enhancement of weak fatigue feature components in the load attractor geometric features using a stochastic resonance model to generate an enhanced load feature set; Inputting the enhanced load feature set into a sparse representation classifier, separating damage sparse coefficients associated with the bond line failure mechanism of the IGBT module through dictionary learning, and determining cumulative damage primitives based on the energy distribution of the damage sparse coefficients; Mapping the cumulative damage primitives to the spatiotemporal evolution diagram of interface stress in each material layer inside the IGBT module using a graph neural network, outputting the fatigue crack propagation rate field of each material layer; Aligning the fatigue crack propagation rate field to a preset baseline lifetime database using transfer learning, constructing a remaining lifetime distribution probability map by comparing the difference in rate fields before and after alignment, and outputting the predicted remaining effective lifetime of the IGBT module accordingly.

[0015] In specific implementation, please refer to Figure 2The transient junction temperature sequence in the original multi-condition load spectrum is subjected to a time-delay embedding operation. The correlation dimension and maximum Lyapunov exponent of the transient junction temperature sequence are calculated. Based on the calculated correlation dimension and maximum Lyapunov exponent, the embedding dimension and time delay parameters required for phase space reconstruction are determined. More specifically, the transient junction temperature sequence is zero-mean normalized to obtain a normalized transient junction temperature sequence. The mutual information function of the normalized transient junction temperature sequence is calculated, and the time delay corresponding to the first local minimum of the mutual information function is used as the time delay parameter for phase space reconstruction. The pseudo-nearest neighbor method is used to determine the embedding dimension parameters. Specifically, starting with a small initial embedding dimension, the normalized transient junction temperature sequence is reconstructed in phase space with a fixed time delay parameter. The proportion of pseudo-nearest neighbors in each embedding dimension is calculated. The embedding dimension corresponding to the first drop in the proportion of pseudo-nearest neighbors below a preset saturation threshold is selected as the embedding dimension parameter for phase space reconstruction. The preset saturation threshold is set according to the noise level of the transient junction temperature sequence. When the signal-to-noise ratio (SNR) of the transient junction temperature sequence is higher than 20 dB, the preset saturation threshold is set to 0.05; when the SNR is between 10 dB and 20 dB, the preset saturation threshold is set to 0.1. For the reconstructed phase space matrix obtained using the same time delay and embedding dimension parameters, the correlation dimension is calculated using the GP algorithm: a series of increasing neighborhood radius values ​​are selected in the reconstructed phase space matrix, the number of phase point pairs contained within each neighborhood radius is counted, and a double logarithmic curve of the correlation integral and neighborhood radius is plotted. The slope of the linear scaling interval in the double logarithmic curve is taken as the correlation dimension of the transient junction temperature sequence. For each phase point in the reconstructed phase space matrix, the maximum Lyapunov exponent is calculated using a small data method: the evolution trajectory of each phase point and its nearest neighbor in the phase space is traced, the distance growth rate between phase points at each time step is calculated, and the slope of the average distance growth rate over time is calculated after taking the natural logarithm. This slope is the maximum Lyapunov exponent of the transient junction temperature sequence. After obtaining the correlation dimension and the maximum Lyapunov exponent, the correlation dimension is used to constrain the lower limit of the integer value of the minimum embedding dimension, i.e., the embedding dimension parameter is not less than the smallest integer greater than the correlation dimension. Simultaneously, the maximum Lyapunov exponent is retained as one of the characteristic quantities characterizing the chaotic properties of the load attractor.

[0016] Based on the determined time delay and embedding dimension parameters, the transient junction temperature sequence is projected into a high-dimensional phase space, forming a load phase trajectory within the high-dimensional phase space. For each phase point in the load phase trajectory, the nearest several phase points are searched to form a local neighborhood. The geometric radius of the local neighborhood is calculated as the local neighborhood radius of that phase point, and the number of phase points contained within a unit volume of phase space is counted as the phase point density of that phase point. Based on the distribution of local neighborhood radii and phase point densities of all phase points, the geometric features of the load attractor are extracted. Specifically, the geometric features of the load attractor include the attractor dimension, Kolmogorov entropy, and the maximum Lyapunov exponent. The attractor dimension is calculated based on the scaling law of the phase point density in the phase space; the Kolmogorov entropy is obtained by calculating the statistical average of the divergence rates of adjacent phase point trajectories in the phase trajectory, reflecting the rate of information loss within the load attractor; the maximum Lyapunov exponent uses the value calculated in the previous step, reflecting the predictability of the long-term evolution of the phase trajectory.

[0017] The geometric features of the load attractor are input into the bistable stochastic resonance model. The processing procedure of the bistable stochastic resonance model is as follows: A bistable system described by the nonlinear Langevin equation is constructed. The potential function of the bistable system contains two steady-state wells and a potential barrier. The barrier height is used as an adjustment parameter of the model. The system is simultaneously injected with a useful signal composed of the geometric features of the load attractor and environmental noise. The intensity of the environmental noise is used as another adjustment parameter of the model. The value of the attractor dimension is used as a reference for determining the initial value of the barrier height, and the value of the Kolmogorov entropy is used as a reference for determining the initial value of the noise intensity parameter. Near the initial values ​​of the barrier height and the initial value of the noise intensity parameter, a two-dimensional parameter scan is performed on the barrier height and the noise intensity according to a preset adjustment step size. Under various parameter combinations, the bistable stochastic resonance model is driven to generate an output signal, and the signal-to-noise ratio and peak frequency of the power spectral density of the output signal are calculated. When a certain parameter combination causes the signal-to-noise ratio of the output signal to reach a local maximum, and the deviation between the peak frequency of the power spectral density of the output signal and the preset fatigue characteristic frequency is less than a preset frequency deviation tolerance, the bistable stochastic resonance model is determined to be operating in a stochastic resonance state. The fatigue characteristic frequency was extracted from the characteristic signal of the bond wire failure process by performing a standard power cycle accelerated aging test on the IGBT module. The preset frequency deviation tolerance was set to ±1% of the fatigue characteristic frequency. The barrier height parameter and noise intensity parameter values ​​that put the bistable stochastic resonance model into a stochastic resonance state were locked. The maximum Lyapunov exponent in the geometric features of the load attractor was used as the input signal and re-inputted into the parameter-locked bistable stochastic resonance model. The output signal of the bistable stochastic resonance model is the resonance enhancement signal of the weak fatigue characteristic component.

[0018] The resonant enhancement signal is processed by wavelet packet decomposition. A 5-level wavelet packet decomposition using the db6 wavelet basis function yields 32 decomposition nodes, each corresponding to a sub-band signal. The energy proportion of each decomposition node is calculated, defined as the ratio of the sub-band signal energy of that node to the total energy of the resonant enhancement signal. The peak frequency of each decomposition node is calculated, defined as the frequency corresponding to the maximum power spectral density of its sub-band signal. The global mean of the energy proportions of all decomposition nodes is calculated; this global mean is the arithmetic mean of the energy proportions of all decomposition nodes. Based on the pre-obtained fatigue-sensitive frequency band range of the IGBT module bonding wire material, the fatigue-sensitive frequency band is calibrated using the crack propagation characteristic frequency band obtained from bonding wire material fatigue tests. This band typically falls within the range of 100Hz to 2kHz and depends on specific material parameters; in this embodiment, the fatigue-sensitive frequency band is taken as 200Hz to 800Hz. From the 32 decomposition nodes, select the decomposition nodes that simultaneously satisfy the condition that the energy proportion exceeds the global average and the peak frequency falls within the fatigue-sensitive frequency band of 200Hz to 800Hz. Combine the sub-frequency band signals corresponding to these selected decomposition nodes to form an enhanced load feature set.

[0019] In the specific implementation, the enhanced load feature set is divided into a training dictionary and test samples. A predetermined proportion of 80% of the samples are randomly selected from the enhanced load feature set as the training dictionary, and the remaining 20% ​​are used as test samples. Each sample in the training dictionary is column-vectorized, and all the column-vectorized training dictionary samples are arranged column-wise to form a training dictionary matrix. The number of columns in the overcomplete dictionary matrix is ​​set to be greater than the number of rows, and the number of rows in the overcomplete dictionary matrix is ​​equal to the dimension of a single training dictionary sample after column vectorization. The number of columns in the overcomplete dictionary matrix, i.e., the number of atoms, is set to 512 based on the diversity requirements of load fluctuation modes. An initial overcomplete dictionary matrix is ​​generated using a random initialization method. Each column vector of the initial overcomplete dictionary matrix corresponds to one atom, and each column vector is normalized using the L2 norm. The K-Singular Value Decomposition (KSVD) algorithm is used to perform overcomplete dictionary learning on the training dictionary matrix. The optimization objective of KSVD is to minimize the reconstruction error between the training dictionary matrix and the overcomplete dictionary matrix multiplied by the sparse coefficient matrix, while constraining the sparsity of each column of the sparse coefficient matrix to not exceed a preset sparsity upper limit, which is set to one-tenth of the dimension value of the training dictionary samples. The execution process of KSVD consists of alternating iterative sparse encoding and dictionary update phases. In the sparse encoding phase, the overcomplete dictionary matrix obtained in the previous iteration is fixed. For each column vector in the training dictionary matrix, the Orthogonal Matching Pursuit (ORP) algorithm is used to solve for the corresponding sparse coefficient vector. The convergence condition of ORP is that the L2 norm of the current residual decreases to less than one-hundredth of the L2 norm of the initial residual, or the number of selected atoms reaches the preset sparsity upper limit. In the dictionary update phase, the sparse coefficient matrix is ​​fixed, and each atom in the overcomplete dictionary matrix is ​​updated column by column. When updating the j-th atom in the overcomplete dictionary matrix, the residual matrix after removing the contribution of the j-th atom is calculated. A local residual matrix is ​​formed by selecting the residual columns corresponding to the training dictionary samples that used the j-th atom. Singular value decomposition (SVD) is performed on the local residual matrix to obtain a left singular matrix, a singular value diagonal matrix, and a right singular matrix. The first column of the left singular matrix is ​​used as the updated j-th atom vector. The product of the first singular value of the singular value diagonal matrix and the transpose of the first row of the right singular matrix is ​​used as the coefficient value at the non-zero position in the j-th row of the updated sparse coefficient matrix. The sparse coding stage and the dictionary update stage are repeated alternately until the Frobenius norm of the changes in all atoms in the overcomplete dictionary matrix is ​​less than a preset convergence threshold, which is set to 1×10⁻⁶. -6At this point, the final overcomplete dictionary is output. During the execution of the K-singular value decomposition algorithm, a bond line failure type label is pre-assigned to each sample in the training dictionary. These labels include bond line peeling, aluminum reconstruction, and solder fatigue. In the sparse coding stage, for training dictionary samples with specific bond line failure type labels, their sparse decomposition is restricted to selecting atoms only within a preset atomic index range. Atomic index ranges 1 to 200 are pre-assigned to the bond line peeling label, 201 to 400 to the aluminum reconstruction label, and 401 to 512 to the solder fatigue label. Through these constraints, a definite correspondence is formed between the atomic index ranges and the bond line failure mechanisms of the IGBT module in the final learned overcomplete dictionary. The load fluctuation pattern corresponding to each atom in the overcomplete dictionary is attributed to an associated bond line failure mechanism.

[0020] Each column vector in the test sample is then input into a sparse representation classifier, which incorporates the overcomplete dictionary learned earlier. The orthogonal matching pursuit algorithm is used to perform sparse decomposition on the test sample column vectors within the overcomplete dictionary. During the solution process, the sparsity is also limited to a preset upper limit, resulting in a sparse coefficient vector corresponding to the test sample. The positions of the non-zero coefficients in the sparse coefficient vector represent which atoms in the overcomplete dictionary match the load fluctuation pattern of the current test sample.

[0021] Damage sparse coefficients associated with the bond line failure mechanism of the IGBT module are selected from the sparse coefficient vector. Each non-zero coefficient in the sparse coefficient vector is traversed, and the corresponding atomic index value is obtained. The preset index interval is determined. When the atomic index value is in the interval 1 to 200, the associated failure mechanism is bond line peeling; when the atomic index value is in the interval 201 to 400, the associated failure mechanism is aluminum reconstruction; and when the atomic index value is in the interval 401 to 512, the associated failure mechanism is solder fatigue. All non-zero coefficient amplitudes within these three intervals are extracted to form a damage sparse coefficient set. Each coefficient amplitude in the damage sparse coefficient set is recorded with its corresponding atomic index value and failure mechanism label.

[0022] The Lp norm energy distribution of the impaired sparse coefficients is calculated, with p=2. The reason for setting p to 2 is that the L2 norm energy can directly represent the power contribution intensity of the signal components during the reconstruction process. The Lp norm energy distribution of the impaired sparse coefficients consists of the energy value of each impaired sparse coefficient in the set. The formula for calculating the energy value of the m-th impaired sparse coefficient is: ; in, This represents the magnitude of the m-th damage sparse coefficient in the set of damage sparse coefficients. This represents the energy value corresponding to the m-th damage sparsity coefficient. For the energy value sequence... Local peak detection is performed, where Q is the total number of coefficients in the sparse coefficient set of the damage. The method for local peak detection is as follows: the energy value sequence is arranged in index order, and for index m, the following condition is met: For each position, if simultaneously satisfying and Then determine the energy value at the m-th position. This is a local peak. For m equals 1, if... Then it is considered a local peak; for m equal to Q, if These are considered local peaks. The atomic indices and coefficient magnitudes corresponding to all local peaks are extracted. The extracted atomic index value at each local peak is multiplied by the coefficient magnitude to obtain a product term. All product terms are arranged in ascending order of atomic index values ​​to form a cumulative damage primitive sequence.

[0023] In specific implementation, please refer to Figure 3 A topology graph of the material layers within the IGBT module is constructed. From top to bottom, the material layers within the IGBT module include the chip layer, chip solder layer, upper copper layer, ceramic layer, lower copper layer, and system solder layer. The interface regions between each material layer and its adjacent material layers are discretized, dividing the interface region into quadrilateral discrete units of the same size. Each discrete unit serves as a node in the topology graph. All discrete units on the material layer interfaces together constitute the node set of the topology graph, and each node is assigned a unique node index. Edges in the topology graph represent the thermal coupling relationship between adjacent discrete units. The rules for establishing edges are as follows: if two discrete units share the same unit boundary in space and are located on the same material layer interface, an undirected edge is established between the corresponding two nodes; if two discrete units belong to two adjacent material layer interfaces (upper and lower) and their projections overlap in the normal direction, an undirected edge is established between the corresponding two nodes. Each edge is assigned a weight, which is the product of the contact area and thermal conductivity between the two discrete units divided by the Euclidean distance between the center points of the two discrete units, representing the strength of the thermal coupling. The node attributes of the topology graph structure are the temperature gradient and equivalent stress of the corresponding discrete element. The temperature gradient is obtained by measuring the temperature value through thermocouple sensors placed at the center of the discrete element and then calculating the difference along the spatial coordinate direction. The equivalent stress is calculated using the strain data measured by strain gauges near the discrete element through the Von Mises stress formula.

[0024] The cumulative damage primitives are used as temporal input features. The cumulative damage primitive sequence consists of cumulative damage primitives from multiple consecutive time steps. Each cumulative damage primitive at a time step is a row vector, and the vector length is equal to the number of product terms extracted at that time step. A graph neural network model consisting of two graph convolutional layers and one graph-gated recurrent unit layer is used to map the cumulative damage primitives to the spatiotemporal evolution diagram of interface stress in each material layer within the IGBT module. The complete forward propagation flow of the graph neural network model is as follows: Within a single time step, the cumulative damage primitive vector for the current time step and the node attribute vectors of all nodes in the topological graph are input into the graph neural network model. The cumulative damage primitive vector is concatenated with the node attribute vector of each node. The node attribute vector of each node has a dimension of 2 and includes temperature gradient values ​​and equivalent stress values. The concatenation results in the initial hidden state vector of each node, with a dimension equal to 2 plus the length of the cumulative damage primitive vector. The first layer of graph convolution is performed on the topological graph structure. The first layer of graph convolution aggregates and transforms the set of neighboring nodes and the initial hidden state vector of each node using the mean aggregation function. After aggregation, a linear transformation and the ReLU activation function are applied to generate the first hidden state vector. The first hidden state vector output by the first layer of graph convolution is fed into the second layer of graph convolution. The structure of the second layer of graph convolution is the same as that of the first layer of graph convolution, but it uses independent weight matrices and bias vectors. The aggregation method is still mean aggregation, and the second hidden state vector of each node is output.

[0025] In the first layer graph convolution operation, the first hidden state vector of node v is calculated as follows: ; in, This represents the initial hidden state vector of node u at time step t; MEAN represents the set of neighboring nodes of node v that are directly connected to node v by an edge in the topological graph structure; MEAN represents the element-wise average of all vectors within the curly braces. It is the weight matrix for the first layer graph convolution operation. It has 64 rows and the number of columns is equal to the dimension of the initial hidden state vector. The initial values ​​of all elements in the weight matrix are randomly sampled from a Gaussian distribution with a mean of 0 and a standard deviation of 0.01. It is the bias vector of the first layer graph convolution operation, with a dimension of 64, and all elements are initialized to 0; This represents the modified linear unit activation function. This is the first hidden state vector of node v. The second-layer graph convolution operation uses the same aggregation and transformation to obtain the second hidden state vector of node v. The output dimension of the second-layer graph convolution operation is set to 32, and the initialization method of the weight matrix and bias vector is the same as that of the first-layer graph convolution operation.

[0026] After the second layer of graph convolution is completed, a graph-gated recurrent unit (GRU) layer is applied to the second hidden state vector of each node. The GRU layer receives the second hidden state vector of node v at the current time step t. And the hidden state vector of node v output from the previous time-step graph-gated recurrent unit layer. Calculate the hidden state vector of the current time step graph gated recurrent unit layer. The calculation process for the gated cyclic cell layer includes the following steps: First, calculate the reset gate vector. and updating the gate vector : ; ; Then calculate the candidate hidden state vector. : ; Finally, calculate and output the hidden state vector: ; In the above calculations, The sigmoid activation function is represented by sigmoid; the hyperbolic tangent activation function is represented by tanh. This represents element-wise multiplication of vectors. , , For the second hidden state vector respectively The weight matrices are all 32×32 in shape; , , These are the hidden state vectors from the previous time step. The weight matrices are all 32×32 in shape; , , The bias vectors are all 32 in dimension. The initial values ​​of the weight matrices of all graph-gated recurrent unit layers are randomly sampled from a Gaussian distribution with a mean of 0 and a standard deviation of 0.01, and the initial values ​​of the bias vectors are all set to zero. This is the hidden state vector of node v at time step t, which serves as one of the inputs to the graph-gated cyclic unit layer for the next time step, and also as the interface stress encoding representation of node v at the current time step.

[0027] The graph neural network model is trained on historical operating data before use. A training sample set is constructed by extracting temperature gradient equivalent stress, cumulative damage primitives, and the actual interface stress values ​​of corresponding discrete units calibrated by ultrasonic scanning or destructive testing from several time steps throughout the entire lifespan of multiple IGBT modules from a baseline lifespan database. During training, the sequence of all node attribute vectors and cumulative damage primitive vectors from multiple consecutive time steps is input into the graph neural network model. The hidden state vector output from the graph-gated recurrent unit layer is passed through a linear mapping layer to generate scalar interface stress prediction values. The weight vector of the linear mapping layer has a dimension of 32, a scalar bias, and is initialized to 0. The mean squared error loss between the predicted interface stress value and the actual interface stress value is calculated. The Adam optimizer with a learning rate of 0.001 is used to update all parameters of the graph neural network model through gradient backpropagation, iterating until the relative decrease in the loss value over 10 consecutive training cycles is less than 1%.

[0028] After training, the multi-time-step cumulative damage primitive sequence of the IGBT module to be predicted and the real-time acquired node attribute data are input into the graph neural network model. The graph neural network model outputs the hidden state vector of each node at consecutive time steps, which is transformed into the interface stress prediction value through a linear mapping layer. The interface stress prediction values ​​of all nodes are arranged according to the material layer order and spatial coordinates: the nodes of each material layer interface are organized into a two-dimensional matrix according to the spatial row and column coordinates of their discrete units. The two-dimensional matrices of all material layer interfaces are stacked along the depth direction to form a data cube with three dimensions: material layer index, spatial row coordinates, and spatial column coordinates. Each time step corresponds to a data cube, and the time step order constitutes the time dimension. The overall data is the spatiotemporal evolution diagram of interface stress.

[0029] Edge detection and region growth segmentation are performed on the spatiotemporal evolution map of interface stress. Edge detection uses the Sobel operator to convolve the two-dimensional stress distribution of each material layer interface at each time step, calculates the stress gradient magnitude, and extracts pixels whose gradient magnitude exceeds 0.3 times the mean stress of the material layer interface to form an edge point set. Region growth segmentation uses the location identified as an edge point in three consecutive time steps as the seed point. The growth criterion is that the absolute value of the difference between the stress value of an adjacent pixel and the stress value of the seed point is less than 0.5 times the standard deviation of the stress at the interface where the seed point is located. The growth region expands continuously until no new pixels meet the criterion; the resulting growth region is the fatigue crack region. Comparing the fatigue crack regions in adjacent time steps, the boundary of the newly added crack region is extracted as the fatigue crack front position. The fatigue crack propagation direction is defined as the normal vector direction from the inside of the crack region to the outside at each point on the fatigue crack front contour. The propagation rate of each point is calculated based on the time interval between adjacent time steps and the displacement of the crack front position. The propagation rate values ​​at each discrete unit are arranged according to spatial coordinates to generate the fatigue crack propagation rate field for each material layer.

[0030] In practical implementation, a pre-defined baseline lifespan database is obtained. The construction process of the baseline lifespan database is as follows: At least 50 historical IGBT modules that have completed their full lifespan operation and whose failure mode has been confirmed as bond wire failure through disassembly and analysis are collected, covering different power levels, different packaging forms, and different operating conditions. For each historical IGBT module, during its entire lifespan operation, the fatigue crack propagation rate field corresponding to each time step is recorded at a fixed sampling interval. The fatigue crack propagation rate fields of all time steps are arranged in chronological order to form a fatigue crack propagation rate field sequence. Simultaneously, the total operating time experienced by the historical IGBT module from the start of operation to the occurrence of bond wire functional failure is recorded as an actual lifespan tag. The fatigue crack propagation rate field sequences of all historical IGBT modules are paired and stored with their corresponding actual lifespan tags to form the baseline lifespan database.

[0031] A domain-adaptive transfer learning method is employed to align the fatigue crack propagation rate field of the IGBT module to be predicted to the source domain distribution of the baseline life database. The fatigue crack propagation rate field of the IGBT module to be predicted is considered as the target domain sample, and the fatigue crack propagation rate field sequences of all historical IGBT modules in the baseline life database are flattened into several source domain samples. The maximum mean difference distance between the target domain samples and all source domain samples is calculated. The maximum mean difference distance is calculated using a Gaussian kernel function, with the bandwidth parameter of the Gaussian kernel function set to the median of the pairwise Euclidean distances between all source domain samples. The formula for calculating the maximum mean difference distance is: ; in, This represents the set of feature mapping functions, specifically all functions within the unit sphere in the reproducing kernel Hilbert space; The source domain sample set is composed of flattened vectors of all fatigue crack propagation rate fields in the baseline life database. The target domain sample set is composed of the flattened vector of the fatigue crack propagation rate field of the IGBT module to be predicted. This represents the total number of samples in the source domain sample set; This represents the total number of samples in the target domain sample set; This represents the i-th sample vector in the source domain sample set; This represents the j-th sample vector in the target domain sample set; The feature mapping function represents the feature space from the original feature space to the reproducing kernel Hilbert space; Represents the regenerated nucleus Hilbert space The norm in.

[0032] By minimizing the maximum mean difference distance, the edge distribution of the fatigue crack propagation rate field of the IGBT module to be predicted is aligned to the source domain distribution of the benchmark lifetime database. Specifically, a correction perturbation is applied to the flattened vector of the fatigue crack propagation rate field of the IGBT module to be predicted. This perturbation is opposite in direction to the gradient of the maximum mean difference distance relative to the flattened vector. The perturbation step size is set to 0.01. This perturbation is iteratively applied until the decrease in the maximum mean difference distance is less than one-thousandth of its current value, resulting in the aligned fatigue crack propagation rate field. After alignment, the fatigue crack propagation rate field of the IGBT module to be predicted and the source domain samples in the benchmark lifetime database are in a common feature space.

[0033] In the common feature space, the difference between each spatial location of the aligned fatigue crack propagation rate field and the mean rate field value at the corresponding spatial location of all source domain samples in the baseline life database is calculated, forming a rate field difference tensor. The dimension of the rate field difference tensor is the same as the spatial dimension of the fatigue crack propagation rate field, and each element in the tensor is a scalar difference value. Tensor decomposition is performed on the rate field difference tensor using the Tucker decomposition method, decomposing it into a kernel tensor and the sum of the products of factor matrices along each mode. The principal component cutoff ranks for each mode are set as follows: material layer mode cutoff rank is 3, row space mode cutoff rank is 5, and column space mode cutoff rank is 5. The reconstructed low-rank approximate tensor is used as the principal difference mode, and the difference between the original rate field difference tensor and the principal difference mode is used as the residual difference mode.

[0034] A Bayesian deep regression network is constructed, consisting of four fully connected layers and one output layer. The input dimension of the first fully connected layer is equal to the total dimension of the concatenated main difference pattern and residual difference pattern. The output dimension of the first fully connected layer is 128; the output dimension of the second fully connected layer is 64; the output dimension of the third fully connected layer is 32; the output dimension of the fourth fully connected layer is 16; and the output layer has a dimension of 2. The two output neurons are used to generate the predicted mean and predicted variance of the remaining lifetime, respectively. The weight parameters of each fully connected layer follow a Gaussian prior distribution with a mean of 0 and a standard deviation of 0.1. Dropout layers are inserted as Monte Carlo dropout modules between the first and second fully connected layers, between the second and third fully connected layers, and between the third and fourth fully connected layers. The dropout probability of the three dropout layers is uniformly set to 0.2 during the training and forward inference phases. The dropout operation randomly sets the activation values ​​of the input neurons to zero according to a Bernoulli distribution and scales the activation values ​​of neurons with a ratio of 0.8.

[0035] The principal difference pattern vector and the residual difference pattern vector are concatenated end-to-end to generate a joint feature vector. The dimension of the joint feature vector is equal to the sum of the number of principal difference pattern elements and the number of residual difference pattern elements. This joint feature vector is input into a Bayesian deep regression network, and the output layer activation values ​​are calculated through forward propagation. During forward propagation, Monte Carlo dropout inference is performed 50 times. In each Monte Carlo dropout inference, all three dropout layers independently and randomly drop a corresponding proportion of neurons with a dropout probability of 0.2, resulting in 50 different sets of output activation values. The arithmetic mean of the first output neuron among the 50 sets of output activation values ​​is calculated as the predicted mean of remaining lifetime, and the standard deviation of the first output neuron among the 50 sets of output activation values ​​is calculated as the predicted variance of remaining lifetime. The logarithm of the predicted variance of the second output neuron is taken as the exponent, and then weighted and averaged with the standard deviation of the first output neuron, with a weighting coefficient of 0.5, to obtain the final predicted variance.

[0036] Using the predicted mean as the center and the predicted variance as the scale, a probability density function of the remaining lifetime on the time axis is constructed under the Gaussian distribution assumption. The probability density function is discretized and sampled at equal intervals within three standard deviations of the predicted mean, with the step size set to one-tenth of the final predicted variance, resulting in a remaining lifetime probability distribution plot. The horizontal axis of the remaining lifetime probability distribution plot represents the remaining lifetime time value, and the vertical axis represents the probability density value. The time coordinate corresponding to the maximum probability density value in the remaining lifetime probability distribution plot is retrieved, and this time coordinate is output as the predicted remaining effective lifetime of the IGBT module.

[0037] In practical implementation, the transient junction temperature sequence in the original multi-condition load spectrum is normalized to zero mean. The arithmetic mean of all sampled values ​​in the transient junction temperature sequence is calculated. This arithmetic mean is then subtracted from each sampled value in the transient junction temperature sequence, and the result is divided by the standard deviation of the transient junction temperature sequence to obtain the normalized transient junction temperature sequence. The autocorrelation function of the normalized transient junction temperature sequence is then calculated. The autocorrelation function is calculated as follows: for the delay... Take the normalized transient junction temperature sequence and its own translation. The Pearson correlation coefficient of the sequence after sampling points, iterating through... An autocorrelation function sequence is generated using integer values ​​ranging from 1 to one-third of the length of the normalized transient junction temperature sequence. A preset threshold is set to 1 / e of the maximum value of the autocorrelation function sequence, where e is the natural constant, approximately 2.718. The autocorrelation function sequence is scanned, and the delay corresponding to the first time the autocorrelation function value drops below the preset threshold from a positive value is recorded. This delay is used as a candidate time delay. The mutual information function of the normalized transient junction temperature sequence is calculated as follows: for the delay... The normalized transient junction temperature sequence is compared with its own delay. The sequence after each sampling point forms a two-dimensional scatter plot. The plane containing the scatter plot is uniformly divided into 256 grid cells. The frequency of the scatter points in each grid cell is counted as the joint probability density estimate. Simultaneously, the marginal probability density estimate of the normalized transient junction temperature sequence is calculated. The mutual information value under this delay is calculated according to the mutual information definition formula. This process is iterated through... A mutual information function sequence is generated by taking integer values ​​from 1 to one-third of the length of the normalized transient junction temperature sequence. The mutual information function sequence is scanned, and the delay corresponding to the first occurrence of a local minimum value in the mutual information function is recorded. This delay is then used as another candidate time delay. The criterion for determining a local minimum is: the mutual information function value at that point is less than the mutual information function values ​​of the two points immediately preceding and following it.

[0038] The pseudo-neighborhood method is used to determine the final embedding dimension. Starting with a candidate embedding dimension d=2, the value is incremented sequentially to 15. For each candidate embedding dimension d, the smaller of the two candidate time delays is selected as the current computational time delay. The normalized transient junction temperature sequence is then reconstructed in phase space according to the current computational time delay and the candidate embedding dimension d. After reconstruction, in the d-dimensional phase space, for each phase point, the phase point with the nearest Euclidean distance is searched, and the distance between this phase point and its nearest neighbor is calculated and recorded. Subsequently, the embedding dimension is increased to d+1, and the nearest neighbor phase point with the same phase point index is searched again in the d+1-dimensional phase space. The nearest neighbor distance at this time is calculated and recorded as follows. Define the pseudo-neighborhood criterion ratio. If P is greater than the pseudo-neighborhood threshold, then the phase point is determined to be a pseudo-nearest neighbor in the embedding dimension d. The pseudo-neighborhood threshold is set to 10. The proportion of pseudo-nearest neighbors among all phase points in the embedding dimension d is calculated to obtain the pseudo-neighborhood proportion. When the pseudo-neighborhood proportion first decreases below 0.01 with the increase of the embedding dimension and remains stable without significant fluctuation, the corresponding embedding dimension d is taken as the final embedding dimension. The criterion for determining that the pseudo-neighborhood proportion is stable without significant fluctuation is that the change in the pseudo-neighborhood proportion in three consecutive embedding dimensions is less than 0.001.

[0039] The final time delay is determined from two candidate time delays using the average displacement method. The normalized transient junction temperature sequence is reconstructed into a phase space spanned by the final embedding dimension and the candidate time delays. For all phase points in the reconstructed phase space, the Euclidean distance between each phase point and its temporally adjacent phase points is calculated. The arithmetic mean of the corresponding Euclidean distances for all phase points is taken as the average displacement. The values ​​of the candidate time delays are varied, and the average displacement is calculated for each time delay by traversing the parameter interval formed by the two candidate time delays with a step size of 1, resulting in an average displacement variation curve. The slope of the average displacement variation curve at each time delay point is calculated. When the slope at a certain time delay point is first less than a preset slope threshold of 0.01, that time delay point is selected as the final time delay.

[0040] The normalized transient junction temperature sequence is reconstructed in phase space based on the final embedding dimension and the final time delay. The reconstruction method is as follows: the normalized transient junction temperature sequence is truncated sequentially according to the final time delay. Each truncated segment contains sampling points of the final embedding dimension, and the starting sampling points of adjacent segments differ by one sampling interval. All truncated segments are arranged in rows to form the reconstructed phase space matrix. Singular value decomposition is performed on the reconstructed phase space matrix to obtain a sequence of singular values ​​arranged in descending order. A logarithmic curve is plotted with the singular value ordinal number as the x-axis and the natural logarithm of the singular values ​​as the y-axis. The slope of the linear segment in the logarithmic curve, where the singular values ​​decrease rapidly and then flatten out, is calculated, and the absolute value of this slope is used as the correlation dimension. The linear segment is selected by removing the first singular value and one-fifth of the total number of the last singular values ​​from the singular value sequence, and using the remaining singular values ​​to form the fitting interval for the linear segment.

[0041] For each phase point in the reconstructed phase space matrix, the evolution trajectory of its nearest neighbor is traced. The number of tracing steps, Q, is set to one-twentieth of the total number of rows in the reconstructed phase space matrix. For the u-th phase point in the reconstructed phase space matrix, the nearest neighbor phase point with the smallest Euclidean distance to the u-th phase point and whose index interval in the original normalized transient junction temperature sequence is greater than the average orbital period is searched. The average orbital period is obtained by taking the period corresponding to the main peak of the spectrum after performing a Fourier transform on the autocorrelation function of the transient junction temperature sequence. The initial distance between the u-th phase point and its nearest neighbor is denoted as . Within the range of tracking steps k from 1 to Q, calculate the distance between two phase points after advancing k steps in their respective time delay directions. Calculate the exponential separation rate. ,in The sampling time interval is denoted as . The maximum Lyapunov exponent is obtained by averaging the exponential separation rates over all phase points and all tracking steps.

[0042] In the tuning part of the stochastic resonance model, a bistable stochastic resonance model described by the Langevin equation is constructed, whose potential function is... It contains two stable equilibrium points and one unstable equilibrium point, with a potential barrier height. The potential difference between two stable equilibrium points in the potential function is determined by the potential energy difference, where α is the coefficient determining the width of the potential well, and β is the coefficient determining the steepness of the potential well; both α and β are positive real numbers. The value of the correlation dimension is directly used as the initial reference value for determining the α coefficient, and the value of the Kolmogorov entropy is directly used as the initial reference value for determining the noise intensity parameter D of the Gaussian white noise injected into the bistable stochastic resonance model. The noise intensity parameter D is the root mean square amplitude of the Gaussian white noise. The β coefficient in the bistable stochastic resonance model is fixed at 1, and the barrier height is... The relationship between the coefficient and α satisfies Following a preset adjustment step size of 0.05, the coefficients of α are iterated within a ±30% neighborhood of the initial reference value α, and the noise intensity parameter D is iterated within a ±30% neighborhood of the initial reference value D. For each (α, D) parameter combination, the Langevin equation is numerically solved using the fourth-order Runge-Kutta method to drive the bistable stochastic resonance model to generate the output signal. The discrete time step of the Langevin equation is set to one-fiftieth of the sampling time interval of the transient junction temperature sequence. The power spectral density of the output signal of the bistable stochastic resonance model is calculated, and the peak value of the power spectral density curve is taken as the peak power spectral density. The frequency corresponding to the peak power spectral density is extracted as the peak frequency of the power spectral density. Simultaneously, the signal-to-noise ratio (SNR) of the output signal of the bistable stochastic resonance model is calculated. The SNR is defined as the ratio of the signal power at the peak frequency of the power spectral density to the average noise power within a ±50Hz frequency band near the peak frequency of the power spectral density. When the signal-to-noise ratio reaches a local maximum under a certain parameter combination and the deviation between the peak frequency of the power spectral density and the preset fatigue characteristic frequency of 600Hz is less than 6Hz, the bistable stochastic resonance model is determined to be operating in a stochastic resonance state under this parameter combination. The values ​​of the α coefficient and the noise intensity parameter D in this parameter combination are locked. The value of the maximum Lyapunov exponent is used as the input signal to the parameter-locked bistable stochastic resonance model, driving it to generate an output signal. This output signal is the resonance enhancement signal for the weak fatigue characteristic component.

[0043] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for predicting the lifetime of an IGBT module based on multi-condition load spectrum, characterized in that, Includes the following steps: The transient junction temperature sequence and collector current waveform of the IGBT module under multiple operating conditions are collected in real time by an embedded sensor network to construct the original multi-condition load spectrum. Based on chaos theory, the phase space of the original multi-condition load spectrum is reconstructed, the geometric features of the load attractor are extracted, and the weak fatigue feature components in the geometric features of the load attractor are enhanced by using a stochastic resonance model to generate an enhanced load feature set. The enhanced load feature set is input into the sparse representation classifier, and the damage sparse coefficients associated with the bond line failure mechanism of the IGBT module are separated through dictionary learning. The cumulative damage primitive is determined based on the energy distribution of the damage sparse coefficients. A graph neural network is used to map the cumulative damage primitives to the spatiotemporal evolution diagram of the interface stress of each material layer inside the IGBT module, and outputs the fatigue crack propagation rate field of each material layer. Transfer learning is used to align the fatigue crack propagation rate field to a preset baseline lifetime database. By comparing the difference in rate fields before and after alignment, a probability map of remaining lifetime distribution is constructed, and the predicted remaining effective lifetime of the IGBT module is output accordingly.

2. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 1, characterized in that, The phase space reconstruction of the original multi-condition load spectrum based on chaos theory is performed to extract the geometric features of the load attractor. Then, a stochastic resonance model is used to enhance the weak fatigue feature components in the geometric features of the load attractor, generating an enhanced load feature set. Specifically: A time-delay embedding operation is performed on the transient junction temperature sequence in the original multi-condition load spectrum. The correlation dimension and Lyapunov exponent of the transient junction temperature sequence are calculated. The embedding dimension and time delay parameters of the phase space reconstruction are determined based on the correlation dimension and Lyapunov exponent. The transient junction temperature sequence is projected into a high-dimensional phase space according to the embedding dimension and time delay parameters to form the load phase trajectory. The local neighborhood radius and phase point density of each phase point in the load phase trajectory are calculated, and the attractor dimension, Kolmogorov entropy and maximum Lyapunov exponent are extracted from the geometric features of the load attractor. The geometric features of the load attractor are input into the bistable stochastic resonance model. The barrier height and noise intensity parameters of the bistable stochastic resonance model are adjusted to make the bistable stochastic resonance model run in a stochastic resonance state and output the resonance enhancement signal of the weak fatigue characteristic component. Wavelet packet decomposition is performed on the resonant enhancement signal. The energy proportion and peak frequency of each decomposition node are calculated. Decomposition nodes with energy proportions exceeding the global mean and peak frequencies falling within the fatigue-sensitive frequency band are selected as the enhanced load feature set.

3. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 2, characterized in that, The enhanced load feature set is input into a sparse representation classifier, and damage sparse coefficients associated with the IGBT module bond line failure mechanism are separated through dictionary learning. The cumulative damage primitive is determined based on the energy distribution of the damage sparse coefficients. Specifically: The enhanced load feature set is divided into a training dictionary and test samples. The training dictionary is iteratively optimized by K-singular value decomposition to learn an overcomplete dictionary, where each atom of the overcomplete dictionary corresponds to a load fluctuation mode. The test samples are input into the sparse representation classifier, and the orthogonal matching pursuit algorithm is used to perform sparse decomposition on the overcomplete dictionary to obtain the sparse coefficient vector. Damage sparse coefficients associated with the bond line failure mechanism of the IGBT module are selected from the sparse coefficient vector. The bond line failure mechanism includes bond line peeling, aluminum reconstruction and solder fatigue, which correspond to the preset atomic index intervals in the sparse coefficient vector. The Lp norm energy distribution of the damage sparsity coefficients is calculated, local peak detection is performed on the energy distribution, the atomic index and coefficient amplitude corresponding to the energy peak are extracted, and the product sequence of the atomic index and coefficient amplitude is used as the cumulative damage primitive.

4. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 3, characterized in that, The method employs a graph neural network to map the accumulated damage primitives to the spatiotemporal evolution diagram of the interface stress of each material layer within the IGBT module, outputting the fatigue crack propagation rate field of each material layer. Specifically: Construct a topology graph of each material layer inside the IGBT module. The nodes of the topology graph represent discrete units on the interface of different material layers, the edges of the topology graph represent the thermal coupling relationship between adjacent discrete units, and the node attributes are the temperature gradient and equivalent stress of the corresponding discrete unit. The cumulative damage primitive is used as the temporal input feature. A graph convolutional network is used to perform layer-by-layer convolution operations on the topological graph structure. In each convolution operation, the temperature gradient and equivalent stress features of the neighboring nodes are aggregated to update the hidden state vector of each node. After each graph convolution operation, a graph-gated recurrent unit is applied to the updated hidden state vector to capture the memory and forgetting characteristics of the node state in the time dimension, and output the interface stress prediction value of each node in consecutive time steps. The predicted interfacial stress values ​​of all nodes are arranged according to the material layer order and spatial coordinates to form a spatiotemporal evolution map of interfacial stress. Edge detection and region growth segmentation are performed on the spatiotemporal evolution map of interfacial stress to extract the fatigue crack front position and propagation direction of each material layer, thereby generating a fatigue crack propagation rate field.

5. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 4, characterized in that, The process involves aligning the fatigue crack propagation rate field to a preset baseline lifetime database using transfer learning, constructing a remaining lifetime probability map by comparing the differences in the rate fields before and after alignment, and outputting the predicted remaining effective lifetime of the IGBT module accordingly. Obtain a preset baseline life database, which contains fatigue crack propagation rate field sequences of multiple historical IGBT modules under different operating conditions and corresponding actual life tags. A domain-adaptive transfer learning method is adopted to calculate the maximum mean difference distance between the fatigue crack propagation rate field and each rate field sequence in the benchmark life database. By minimizing the maximum mean difference distance, the edge distribution of the fatigue crack propagation rate field is aligned to the source domain distribution of the benchmark life database. In the aligned common feature space, the velocity field difference tensor between each spatial location of the fatigue crack propagation rate field and the corresponding velocity field in the benchmark life database is calculated. Tensor decomposition is performed on the velocity field difference tensor to obtain the main difference mode and the residual difference mode. The main difference pattern and residual difference pattern are input into the Bayesian deep regression network. The Bayesian deep regression network outputs the predicted mean and predicted variance of the remaining lifetime. Based on the predicted mean and predicted variance, a probability distribution map of the remaining lifetime is constructed on the lifetime time axis. The lifetime value corresponding to the maximum probability density in the probability distribution map of the remaining lifetime is taken as the predicted remaining effective lifetime.

6. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 2, characterized in that, The process involves performing a time-delay embedding operation on the transient junction temperature sequence in the original multi-condition load spectrum, calculating the correlation dimension and Lyapunov exponent of the transient junction temperature sequence, and determining the embedding dimension and time delay parameters for phase space reconstruction based on the correlation dimension and Lyapunov exponent. Specifically: The transient junction temperature sequence in the original multi-condition load spectrum is normalized to zero mean. The autocorrelation function and mutual information function of the normalized transient junction temperature sequence are calculated. The time delay when the autocorrelation function first drops to a preset threshold is taken as the candidate time delay. The time delay when the mutual information function first appears at a local minimum value is taken as the candidate time delay. The pseudo-neighborhood method is used to calculate the pseudo-neighborhood ratio under each candidate embedding dimension. The candidate embedding dimension when the pseudo-neighborhood ratio drops to saturation is selected as the final embedding dimension. The average displacement method is used to calculate the average displacement under each candidate time delay. The candidate time delay corresponding to the first time when the slope of the average displacement change curve is less than the preset slope threshold is selected as the final time delay. The normalized transient junction temperature sequence is reconstructed into a phase space based on the final embedding dimension and the final time delay to obtain the reconstructed phase space matrix. The reconstructed phase space matrix is ​​then subjected to singular value decomposition, and the correlation dimension is calculated by taking the logarithmic slope of the non-zero singular values. For each phase point in the reconstructed phase space matrix, trace the evolution trajectory of its nearest neighbor, calculate the exponential separation rate between adjacent phase points, and take the average of all exponential separation rates to obtain the maximum Lyapunov exponent.

7. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 6, characterized in that, The process involves inputting the geometric features of the load attractor into a bistable stochastic resonance model, adjusting the barrier height and noise intensity parameters of the bistable stochastic resonance model to operate in a stochastic resonance state, and outputting a resonance enhancement signal for weak fatigue characteristic components. Specifically: Construct a bistable stochastic resonance model described by the Langevin equation, wherein the potential function of the bistable stochastic resonance model contains two stable equilibrium points and one unstable equilibrium point, and the height of the potential barrier is determined by the potential energy difference between the two stable equilibrium points in the potential function. The attractor dimension in the geometric features of the load attractor is used as the initial reference value for determining the barrier height, and the Kolmogorov entropy is used as the initial reference value for determining the noise intensity parameter. The barrier height and noise intensity parameter are scanned in the neighborhood of the initial reference value according to the preset adjustment step size. The bistable stochastic resonance model is driven to run under each scan combination. The signal-to-noise ratio and peak frequency of the power spectral density output by the bistable stochastic resonance model are calculated. When the signal-to-noise ratio reaches a local maximum and the peak frequency of the power spectral density matches the preset fatigue characteristic frequency, the bistable stochastic resonance model is determined to be in a stochastic resonance state. The barrier height and noise intensity parameters that keep the bistable stochastic resonance model in a stochastic resonance state are locked. The maximum Lyapunov exponent in the geometric features of the load attractor is used as the input signal to the locked bistable stochastic resonance model, and the resonance enhancement signal of the weak fatigue characteristic component is output.

8. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 3, characterized in that, The enhanced load feature set is divided into a training dictionary and test samples. The training dictionary is iteratively optimized using K-singular value decomposition to learn an overcomplete dictionary. Each atom in the overcomplete dictionary corresponds to a load fluctuation mode, specifically: A predetermined proportion of samples are randomly selected from the enhanced payload feature set as training dictionaries, and the remaining samples are used as test samples. Each sample in the training dictionary is processed into a column vector to form a training dictionary matrix. The overcomplete dictionary matrix is ​​initialized as a random matrix, where the number of columns in the overcomplete dictionary matrix is ​​greater than the number of rows, and each column vector corresponds to an atom. The K-singular value decomposition algorithm is used to perform alternating iterations of sparse encoding stage and dictionary update stage on the training dictionary matrix. In the sparse coding stage, the overcomplete dictionary matrix is ​​fixed, and the orthogonal matching pursuit algorithm is used to calculate the sparse coefficient matrix of the training dictionary matrix on the overcomplete dictionary matrix. In the dictionary update stage, the sparse coefficient matrix is ​​fixed, and each atom of the overcomplete dictionary matrix is ​​updated column by column. When updating a column, the residual matrix corresponding to that column is calculated. Singular value decomposition is performed on the residual matrix, and the left singular vector corresponding to the largest singular value is taken as the updated atom. The sparse encoding stage and dictionary update stage are iterated alternately until the change in the column vectors of the overcomplete dictionary matrix is ​​less than the preset convergence threshold, and the final overcomplete dictionary is output.

9. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 8, characterized in that, The process involves selecting damage sparse coefficients from the sparse coefficient vector that are associated with the bond wire failure mechanism of the IGBT module. This bond wire failure mechanism includes bond wire peeling, aluminum reconstruction, and solder fatigue, each corresponding to a preset atomic index interval in the sparse coefficient vector. Specifically: Each atom in the overcomplete dictionary is pre-assigned a failure mechanism label, which includes bond line stripping label, aluminum reconstruction label and solder fatigue label. Each atom is assigned to the corresponding failure mechanism label based on the physical characteristics of the load fluctuation mode used during its training. After the sparse coefficient vector of the test sample is generated, each sparse coefficient in the sparse coefficient vector is traversed in the order of atomic index, and the failure mechanism label of the atom corresponding to each sparse coefficient is queried. The sparse coefficients of failure mechanisms labeled as bond line stripping are extracted into a subset of bond line stripping sparse coefficients, the sparse coefficients of failure mechanisms labeled as aluminum reconstruction are extracted into a subset of aluminum reconstruction sparse coefficients, and the sparse coefficients of failure mechanisms labeled as solder fatigue are extracted into a subset of solder fatigue sparse coefficients. The Lp norm normalization process was performed on the subsets of bond line stripping sparse coefficients, aluminum reconstruction sparse coefficients, and solder fatigue sparse coefficients, respectively. Sparse coefficients whose normalized coefficient values ​​exceeded a preset threshold were retained as damage sparse coefficients, and the atomic index and failure mechanism label corresponding to each damage sparse coefficient were recorded.

10. The IGBT module lifetime prediction method based on multi-condition load spectrum according to claim 1, characterized in that, The embedded sensing network includes multiple thermocouple sensors and Rogowski coil sensors integrated inside the IGBT module, used to synchronously acquire transient junction temperature sequences and collector current waveforms.

Images