Multi-Dimensional Uniform Initialization Gaussian Mixture Model for Spar Crack Quantification under Uncertainty

Guided Wave (GW)-based crack monitoring method as a promising method has been widely studied, as this method is sensitive to small cracks and can cover a wide monitoring range. Online crack quantification is difficult as the initiation and growth of crack are affected by various uncertainties. In addition, crack-sensitive GW features are influenced by time-varying conditions which further increase the difficulty in crack quantification. Considering these uncertainties, the Gaussian mixture model (GMM) is studied to model the probability distribution of GW features. To further improve the accuracy and stability of crack quantification under uncertainties, this paper proposes a multi-dimensional uniform initialization GMM. First, the multi-channel GW features are integrated to increase the accuracy of crack quantification, as GW features from different channels have different sensitivity to cracks. Then, the uniform initialization method is adopted to provide more stable initial parameters in the expectation-maximization algorithm. In addition, the relationship between the probability migration index of GMMs and crack length is calibrated with fatigue tests on prior specimens. Finally, the proposed method is applied for online crack quantification on the notched specimen of an aircraft spar with complex fan-shaped cracks under uncertainty.


Introduction
Structural health monitoring (SHM) technology has gradually developed from basic theoretical research in the laboratory to real aircraft applications recently [1][2][3][4]. Among them, the guided wave (GW) based SHM method has been widely studied because this method can cover a wide monitoring range and it is sensitive to small damage [5][6][7][8][9][10][11]. In addition, the GW-based SHM method can also be applied on-line both for damage monitoring and impact event monitoring.
Fatigue crack is one of the primary damage types in engineering structures. Accurate fatigue crack quantification is of great significance for ensuring the safety and reliability of aircraft structures. Moreover, it can realize condition-based maintenance and extend structure service life. However, there are some obstacles to crack monitoring in real engineering applications. First, initiation and growth of crack are affected by various uncertainties [12,13], such as inhomogeneity of material, scatter of mechanical properties of the material, the randomness of cracking process, and quality of manufacturing, etc. These uncertain factors lead to the scatter of both crack size and fatigue lifetime. Second, various time-varying conditions introduce uncertainty effects on GW features. These timevarying conditions include random dynamic load, varying environmental temperature, changing structural boundary conditions, noise, etc. GW features sensitive to the crack are also usually sensitive to time-varying conditions. Therefore, the uncertain and random variations of the GW features caused by the time-varying conditions lead to lower reliability and stability for crack monitoring in real engineering applications.
Considering these uncertainties, some probability statistical models are widely studied for fatigue crack diagnosis and prognosis recently. Probability statistical models can be introduced to construct a mathematical model of damage features to estimate the damage of structures. Yuan and Mei et al. [14] proposed the GW-Hidden Markov model based method to achieve a probabilistic evaluation of the propagation state of cracks. Uncertainties such as changing load and structural boundary conditions were considered in this research. For crack prognosis, Chen et al. [15] improved the particle filter based methods. In addition, some other probability statistical models, such as the stochastic global model [16] and Gaussian process [17] are also able to characterize uncertainties of GW signals. These probabilistic statistical models focus on addressing different aspects of crack diagnosis and prediction under uncertainty. However, more research is still needed for real engineering applications.
In recent years, the Gaussian mixture model (GMM) has been proven to be an effective approach to model the complicated probability distributions of GW features. GMM is a data-driven model without physical modeling. It can model arbitrary probability distributions by decomposing a non-Gaussian distribution into a combination of a finite number of Gaussian components (GCs) based on unsupervised learning [18]. Yuan and Qiu et al. [19] proposed the GW-GMM based method to model the probability characteristic of GW features under time-varying conditions. This method is validated in the full-scale aircraft fatigue test. Chakraborty et al. [20] proposed a Dirichlet process GMM and it was validated on an aluminum compact tension specimen for crack monitoring subjected to variable-amplitude cyclic load. Banerjee et al. [21] applied the GMM for crack propagation monitoring under the condition of sudden temperature changes. All these researches have proved the potential to apply the GMM to deal with the time-varying problem. However, deep research still needs to be performed concerning further improve the accuracy and stability of crack quantification under uncertainties.
Some structures of aircraft are very complex, and their cracking forms are also complex. For example, the crack form is complex fan-shaped at the notch-edge in an aircraft landing gear spar. There are cracks in multiple directions at the same crack initiation point. At this time, a single channel cannot meet the accuracy requirements of crack quantification. Some researchers reported that different piezoelectric (PZT) excitation-sensing channels can have different sensitivity to damage [22]. Therefore, it is necessary to fuse multi-channel GW features for more accurate crack quantification.
Furthermore, the GMM construction results are influenced by the initial parameters. The parameters of the GMM can be obtained by using the expectation-maximization (EM) algorithm, which is sensitive to the initial value [23]. If the initialization parameters are not set properly, it is easy to make the EM algorithm converge to a relatively poor local optimal solution. EM algorithm initialized by the k-means method is widely used to construct a GMM [18]. However, the centroids are chosen randomly in the k-means method, even with the same training data, k-means can lead to different results and thus provide different GMM initial parameters.
Thus, this paper proposes a multi-dimensional uniform initialization Gaussian mixture model (MdUI-GMM) to improve the accuracy and stability of on-line crack quantification under uncertainty. The multi-dimensional GMM is established by integrating GW features of multiple channels. In this paper, the damage index (DI) based on normalized correlation moment (NCM) is selected as the GW feature. Considering the different sensitivity of different PZT excitation-sensing channels, it is possible to increase the accuracy of crack quantification by integrating the GW features of each channel. Moreover, a uniform initialization method is adopted to improve the stability of multi-dimensional GMM. Compared with the traditional k-means initialization method, the stability of modeling results is improved. Finally, the proposed method is applied for online crack quantification on notched specimens of an aircraft spar with complex fan-shaped cracks under uncertainty. This paper is organized as follows. Section 2 proposes the MdUI-GMM based crack quantification method, including multi-channel fusion GW features extraction, the model- ing process of MdUI-GMM, and the migration measuring method. The implementation process of the MdUI-GMM based crack quantification process is summarized as well. In Section 3, the notched specimen of an aircraft spar fatigue test is introduced, and the multi-channel GW features acquired under dynamic load are given and discussed. Then, crack quantitative calibration results by the MdUI-GMM and crack monitoring results are shown. In addition, comparisons between single-channel and multi-channel fusion measurement are discussed. Finally, conclusions are given in Section 4.

MdUI-GMM Based Crack Quantification Method
In this section, the multi-channel fusion GW features extraction is described first. Then, the modeling process and the migration measuring method of MdUI-GMM are proposed. Finally, the implementation process of the MdUI-GMM based crack quantification process is summarized as well.

Multi-Channel Fusion GW Features Extraction
There are great uncertainties in the crack initiation point and crack propagation direction of each specimen because of the scatter of the manufacturing process and external fatigue load, etc. It is not reliable for crack quantification only based on a single excitation-sensing channel. Because different PZT excitation-sensing channels have different sensitivity to crack, it is possible to increase the accuracy of crack quantification by integrating multi-channel GW features.
The GW-based SHM is deemed as one of the most appealing methods due to its merits including sensitivity to small damage and the ability to monitor a region when forming a sensor network. Damage such as cracks within a structure can cause scattering and energy absorption of GW propagating through the structure. So the crack can be detected when it is located in the monitoring area which is defined by the sensor network. Generally speaking, when performing the GW based SHM technology, a lot of GW signals can be acquired during the online damage monitoring process and the corresponding GW features can be extracted using different kinds of signal processing methods. The GW features can be extracted in the time domain, frequency domain, and time-frequency domain [24][25][26]. The DI based on NCM [27] is used in this article. DI NCM can be calculated by Equation (1), which can be used to obtain the variation of GW signals phase and amplitude. where where b(t) and m(t) represent baseline signal and monitoring signal respectively, t 0 and t 1 are the start and stop time corresponding to the selected signal segment. In this paper, the baseline signal b(t) is the average of the first 10 dynamic load signals at the initial stage of fatigue loading when the structure is in a healthy state.

Modeling Process of MdUI-GMM
The GW signals acquired under aircraft structure service conditions can be considered as a mixture of uncertain changes. Consequently, a GW feature can be considered as a random variable. Let X = [X 1 , X 2 , . . . , X r , . . . , X R ] be a GW feature sample set composed of R independent features that are obtained from R GW signals. Xr denotes a d-dimensional GW feature in the sample set, where Xr = [DI 1 , DI 2 , . . . , DI d ] T and r = 1, 2, . . . , R.
Generally, the variation of GW features introduced by time-varying conditions is random and complex. Consequently, when using a single probability density function to describe the distribution of the GW feature sample set, it may yield unexpected and inaccurate results. Therefore, the mixture probability model is applied to address this problem. The joint probability density function (PDF) of the GW feature sample set can be approximately modeled by a GMM which is considered as a finite weighted sum of GCs. Similar to [18,28], the GMM is expressed as Equation (3): where C is the number of GCs, i = 1, 2, . . . , C. µ i and Σ i are the mean and the covariance matrix of the ith GC. w i is the mixture weight of the ith GC. The PDF of each GC Φ i is a d-dimensional Gaussian function which is expressed as Equation (4) [15] In this paper, the GW features of three channels are extracted to construct multidimensional GMM, where d = 3.
The value of w i , µ i , and Σ i can be obtained based on the sample set by using the EM algorithm. To ensure the stability of training, a uniform initialization process instead of the k-means initialization method is performed. The main modeling process of MdUI-GMM is as follows: Step 1: Before the EM algorithm begins, the uniform weights and posterior probabilities are given for each GC, as indicated in Equations (5) and (6) where C is the number of GCs, i = 1, 2, . . . , C. P ri,0 represents the probability that the initialized rth GW feature belongs to the ith GC and ∑ C i=1 P ri,0 = 1. w i,0 is the initialized mixture weight of the ith GC and ∑ C i=1 w i,0 = 1. Because the total number of GC is C. Therefore, the posterior probabilities P ri and mixture weight w i can be initialized uniformly, as shown in Equations (5) and (6). Based on the uniform initial values of the P ri,0 and w i,0 , the parameters of GMM can be calculated iteratively by the EM algorithm in step 2 and step 3. Since this initialization is only related to the number of GC, a stable initialization result can be obtained after the number of GC is determined.
Step 2: Based on the uniform initial values of the weights and posterior probabilities, the mean vectors and covariance matrices of UIGMM are calculated using Equations (7) and (8) where t is the number of iterations and R is the total number of GW features used to construct GMM. w i,t−1, and P ri,t−1 represent the weight and posterior probability of the last iteration, respectively. Step 3: The posterior probability and weight of GCs are updated according to Equations (9) and (10) where P ri,t is the probability that the rth GW feature belongs to the ith GC in the tth iteration.
Step 4: Repeat Step 2 and Step 3 until the log-likelihood value converges.
L t−1 is the value of the log-likelihood function in the previous iteration step, L t is the log-likelihood function value of the current step. The value of ε is set to be 1 × 10 −10 .

On-Line Migration Measuring Method of the MdUI-GMM
The migration measuring result between the on-line monitoring MdUI-GMM and the baseline MdUI-GMM is defined as the probability migration index (PMI). If the monitored structure keeps healthy, the result of PMI will be maintained at a low level. If a crack occurs in the structure and propagates continuously, a cumulative migration trend will appear in both the MdUI-GMM and the PMI. Based on the result of PMI, reliable damage monitoring can be realized. Considering that the MdUI-GMM is composed of many GCs, a Kullback-Leibler (KL) divergence [29] based method is used to calculate the PMI.
The PDF of the baseline GMM and the nth monitoring GMM are Φ 0 and Φ n , respectively. The KL divergence-based PMI can be calculated by Equation (13).
where d is the dimension of the GW feature vector. µ 0 and Σ 0 represent the weighted sum of the mean and covariance matrices in the baseline GMM, respectively. µ n and Σ n represent the weighted sum of the mean and covariance matrices in the nth monitoring GMM, respectively. tr is the matrix trace.

The MdUI-GMM Based Crack Quantification Process
The effect of time-varying conditions on GW features is random and uncertain. However, there is a cumulative migration trend in the effect of cracks on GW features. Considering this, the GMM is used to model the probability distribution of GW features and suppress the influence of uncertainty. If the structure is in a healthy state or if the cracks are not extended, each GC of the GMM can cover the random distribution of the GW features due to environmental changes. Time-varying conditions (e.g., random changes with dynamic load, or temperature, etc.) do not change the form of the GMM distribution. If the crack propagates, the form of the distribution of the GCs in the GMM undergoes a cumulative migration change. That is, under the influence of cracks, The PDF of monitoring GMM undergoes a cumulative migration change compared to the PDF of the baseline GMM.
The calibrating process of the MdUI-GMM based crack quantification method is shown schematically in Figure 1. The process of calibrating PMI with crack using prior specimens includes three parts: scope and scale lines on one surface of the specimen. After crack length measurement is completed, fatigue testing is resumed. By calculating KL divergence between the baseline GMM in part 1 and the monitoring GMM in part 2, PMI can be obtained. Finally, in combination with the measured crack length, a calibration between PMI and crack length can be performed.  Part 1 is the baseline MdUI-GMM. GW signal acquisition is performed under timevarying conditions when the structure is in a healthy state. GW features calculated by Equation (1) are extracted to construct a GW baseline sample set. After that, the parameters of MdUI-GMM, including weight, mean, and covariance matrix, are trained by an EM algorithm with a uniform initialization process. Finally, the PDF of the baseline GMM can be obtained.
Part 2 is online monitoring MdUI-GMM. Once a GW signal is acquired during the online damage monitoring process, new GW features can be obtained, and the GW feature sample set is updated first. The updating rule is called first in and first out which means that the new feature is added to be the last (newest) feature and the first (oldest) feature is removed. Hence, the number of GW features in the GW feature sample set is consistent. Then, the parameters of monitoring GMM are trained by an EM algorithm with a uniform initialization process. Finally, the PDF of monitoring GMM can be obtained.
Part 3 is crack quantification calibration. When a certain number of load cycles is reached, the fatigue test is paused to measure the crack length using an electron microscope and scale lines on one surface of the specimen. After crack length measurement is completed, fatigue testing is resumed. By calculating KL divergence between the baseline GMM in part 1 and the monitoring GMM in part 2, PMI can be obtained. Finally, in combination with the measured crack length, a calibration between PMI and crack length can be performed.
The process of online crack quantification of the new specimen is shown schematically in Figure 2. Under time-varying conditions, the GW features of the new specimen are obtained online to construct the online monitoring MdUI-GMM and then PMI can be calculated. Finally, according to the quantitative calibration relationship between PMI and crack length in Part 3, the real-time crack length can be calculated directly.
The process of online crack quantification of the new specimen is shown schema cally in Figure 2. Under time-varying conditions, the GW features of the new specim are obtained online to construct the online monitoring MdUI-GMM and then PMI can calculated. Finally, according to the quantitative calibration relationship between PMI a crack length in Part 3, the real-time crack length can be calculated directly.

Validation on the Notched Specimen of an Aircraft Spar under Fatigue Load
In this section, the fatigue test of notched specimens of an aircraft spar is introduce Then, the multi-channel GW features acquired under dynamic load are given and d cussed. After that, crack quantitative calibration results of prior specimens by the MdU GMM, and monitoring results are shown.

Experimental Setup
A typical load-transmitting component in aircraft structures is studied in this pap as illustrated in Figure 3. The thickness of the specimen is 10 mm. There are two obliq notches with a radius of 8 mm on both sides of the length direction. The oblique axis the notch on both sides is 102.72°. To increase the strength of the specimen, two bosses a in the thickness direction. The radius of the outer edge of the boss is around 20 mm, a the height is around 4 mm. Four specimens, with the same size and material paramete are tested and marked as T1-T4, respectively. The first three specimen T1-T3 are used prior specimens for quantitative calibration, and T4 is used for testing.

Validation on the Notched Specimen of an Aircraft Spar under Fatigue Load
In this section, the fatigue test of notched specimens of an aircraft spar is introduced. Then, the multi-channel GW features acquired under dynamic load are given and discussed. After that, crack quantitative calibration results of prior specimens by the MdUI-GMM, and monitoring results are shown.

Experimental Setup
A typical load-transmitting component in aircraft structures is studied in this paper as illustrated in Figure 3. The thickness of the specimen is 10 mm. There are two oblique notches with a radius of 8 mm on both sides of the length direction. The oblique axis of the notch on both sides is 102.72 • . To increase the strength of the specimen, two bosses are in the thickness direction. The radius of the outer edge of the boss is around 20 mm, and the height is around 4 mm. Four specimens, with the same size and material parameters, are tested and marked as T 1 -T 4 , respectively. The first three specimen T 1 -T 3 are used as prior specimens for quantitative calibration, and T 4 is used for testing.
The material of the specimen is aluminum alloy 7A04 T7351. Table 1 lists the mechanical properties of aluminum alloy 7A04. PZT layers are attached to each specimen and used as the actuators and the sensors, whose positions are illustrated in Figure 4. PZT2 and PZT4 are used as the actuators. PZT1 and PZT3 are used as sensors. The distance between Channel 2-1 and Channel 4-3 is 15mm. Because there are two notches in the specimen, cracks may be initiated at both notches, so another four PZT layers are also arranged at the other notch. When the crack starts from Notch 1, the GW features extracted from channel 2-1, channel 2-3, and channel 4-1 will be integrated to construct multi-dimensional GMM. When the crack starts from Notch 2, the GW features extracted from channel 6-5, channel 6-7, and channel 8-5 will be integrated to construct multi-dimensional GMM. the notch on both sides is 102.72°. To increase the strength of the specimen, two in the thickness direction. The radius of the outer edge of the boss is around 2 the height is around 4 mm. Four specimens, with the same size and material p are tested and marked as T1-T4, respectively. The first three specimen T1-T3 a prior specimens for quantitative calibration, and T4 is used for testing.  The material of the specimen is aluminum alloy 7A04 T7351. Table 1 lists t ical properties of aluminum alloy 7A04. PZT layers are attached to each specimen and used as the actuators and whose positions are illustrated in Figure 4. PZT2 and PZT4 are used as the actu and PZT3 are used as sensors. The distance between Channel 2-1 and Cha 15mm. Because there are two notches in the specimen, cracks may be initia notches, so another four PZT layers are also arranged at the other notch. Whe starts from Notch 1, the GW features extracted from channel 2-1, channel 2-3, a 4-1 will be integrated to construct multi-dimensional GMM. When the crack Notch 2, the GW features extracted from channel 6-5, channel 6-7, and channe integrated to construct multi-dimensional GMM.  The material of the specimen is aluminum alloy 7A04 T7351. Table 1 lists the mechanical properties of aluminum alloy 7A04. PZT layers are attached to each specimen and used as the actuators and the sensors, whose positions are illustrated in Figure 4. PZT2 and PZT4 are used as the actuators. PZT1 and PZT3 are used as sensors. The distance between Channel 2-1 and Channel 4-3 is 15mm. Because there are two notches in the specimen, cracks may be initiated at both notches, so another four PZT layers are also arranged at the other notch. When the crack starts from Notch 1, the GW features extracted from channel 2-1, channel 2-3, and channel 4-1 will be integrated to construct multi-dimensional GMM. When the crack starts from Notch 2, the GW features extracted from channel 6-5, channel 6-7, and channel 8-5 will be integrated to construct multi-dimensional GMM. Fan-shaped crack initiation location The experimental setup is illustrated in Figure 5. An MTS electro-hydraulic servo tensile machine is used to apply the variable amplitude fatigue load. The fatigue load spectrum is generated by using the in-service flight data to simulate different flight modes of the aircraft such as takeoff, level flight, overloading, landing, on-ground, etc. A digital microscope is applied to measure crack length. The multi-channel PZT array scanning system developed by the authors' group is employed to perform the active GW-based structural health monitoring. A three-cycle tone-burst signal with the center frequency of 110 kHz and ±70 V amplitude is used as the excitation signal. The sampling frequency is set to be 50M samples s −1 , where the high time resolution of the GW signals enables accurate damage feature extraction. The experimental setup is illustrated in Figure 5. An MTS electro-hydraulic servo tensile machine is used to apply the variable amplitude fatigue load. The fatigue load spectrum is generated by using the in-service flight data to simulate different flight modes of the aircraft such as takeoff, level flight, overloading, landing, on-ground, etc. A digital microscope is applied to measure crack length. The multi-channel PZT array scanning system developed by the authors' group is employed to perform the active GW-based structural health monitoring. A three-cycle tone-burst signal with the center frequency of 110 kHz and ±70 V amplitude is used as the excitation signal. The sampling frequency Based on the in-service experience, the edge fan-shaped cracks of oblique not the main damage forms, because the specimen is thick, and it is the primary loadstructure of the aircraft. There are front crack and side crack at the same crack in point. The front crack propagates along the width direction in the front of the sp The side crack extends along the thickness direction. The crack starts from a certa in the notch and gradually propagates to the inside of the specimen.
The front crack and the side crack have the same crack initiation point. Th between the front crack and the side crack is approximately 90°. Taking the crack in point as the center of the circle, the length of the front crack and the side crack is re as the semi-major axis and the semi-minor axis of the ellipse. The front crack len the side crack length are recorded as a and b, respectively. The damaged area A calculated by Equation (14).
Generally, with the propagation of the side cracks, the thickness direction is ally penetrated through, and then the fracture rapidly expands along the width di resulting in the fracture of the specimen. During the loading process of the fatig cracks of the T1-T4 are all initiated at Notch 1. Due to the difference in material pr and external load, the crack initiation position of each specimen is different.
To quantitatively monitor the crack growth, the crack length is measured by a microscope. Figure 6 shows the relationship between the front crack length, the sid length, and the damaged area with the number of load cycles. Based on the in-service experience, the edge fan-shaped cracks of oblique notches are the main damage forms, because the specimen is thick, and it is the primary load-bearing structure of the aircraft. There are front crack and side crack at the same crack initiation point. The front crack propagates along the width direction in the front of the specimen. The side crack extends along the thickness direction. The crack starts from a certain point in the notch and gradually propagates to the inside of the specimen.
The front crack and the side crack have the same crack initiation point. The angle between the front crack and the side crack is approximately 90 • . Taking the crack initiation point as the center of the circle, the length of the front crack and the side crack is regarded as the semi-major axis and the semi-minor axis of the ellipse. The front crack length and the side crack length are recorded as a and b, respectively. The damaged area A can be calculated by Equation (14).
Generally, with the propagation of the side cracks, the thickness direction is gradually penetrated through, and then the fracture rapidly expands along the width direction, resulting in the fracture of the specimen. During the loading process of the fatigue test, cracks of the T 1 -T 4 are all initiated at Notch 1. Due to the difference in material properties and external load, the crack initiation position of each specimen is different.
To quantitatively monitor the crack growth, the crack length is measured by a digital microscope. Figure 6 shows the relationship between the front crack length, the side crack length, and the damaged area with the number of load cycles. It can be seen that the growth rate of the front crack and side crack is relatively slow at first. With the increase of load cycle number and crack length, the crack growth rate increases. However, the uncertainty leads to the scatter of both crack size and fatigue lifetime. Figure 7a shows the typical S0 mode direct wave signal obtained from T2 channel 2-1 without time-varying conditions. When the excitation center frequency is 110 kHz, the complete S0 mode direct wave can be obtained. When there is no crack in the specimen, the peak amplitude of the S0 mode direct wave is 193.7 mV and the corresponding phase is 92.5 μs. When the front crack extends to 20 mm, the peak amplitude of the S0 mode direct wave is 109.4 mV and the corresponding phase is 93.1 μs. For every 2 mm crack propagation, the peak amplitude of the S0 mode direct wave decreases by about 20 mV. Therefore, without the influence of time-varying conditions, the crack can be monitored by the variation of GW signals phase and amplitude. As shown in Figure 7b, without environmental effects, there are differences in the GW features acquired by different channels when only the effects of cracks are considered. Considering that different channels have different sensitivity to cracks, it is possible to increase the accuracy of crack quantification by integrating multi-channel GW features. The GW feature extracted in this paper is DINCM by Equation (1). It can be seen that the growth rate of the front crack and side crack is relatively slow at first. With the increase of load cycle number and crack length, the crack growth rate increases. However, the uncertainty leads to the scatter of both crack size and fatigue lifetime. Figure 7a shows the typical S 0 mode direct wave signal obtained from T 2 channel 2-1 without time-varying conditions. When the excitation center frequency is 110 kHz, the complete S 0 mode direct wave can be obtained. When there is no crack in the specimen, the peak amplitude of the S 0 mode direct wave is 193.7 mV and the corresponding phase is 92.5 µs. When the front crack extends to 20 mm, the peak amplitude of the S 0 mode direct wave is 109.4 mV and the corresponding phase is 93.1 µs. For every 2 mm crack propagation, the peak amplitude of the S 0 mode direct wave decreases by about 20 mV. Therefore, without the influence of time-varying conditions, the crack can be monitored by the variation of GW signals phase and amplitude. As shown in Figure 7b, without environmental effects, there are differences in the GW features acquired by different channels when only the effects of cracks are considered. Considering that different channels have different sensitivity to cracks, it is possible to increase the accuracy of crack quantification by integrating multi-channel GW features. The GW feature extracted in this paper is DI NCM by Equation (1). However, the structure is rarely only affected by static load without time-varying conditions in real engineering applications. Therefore, crack monitoring under time-varying dynamic load will be discussed. The direct wave signal segment of S0 mode is intercepted to calculate the DI under dynamic load by the Equation (1). t0 and t1 of the signal segment are 0.8 × 10 −4 s and 1.05 × 10 −4 s, respectively. The baseline signal b(t) is the average of the first 10 dynamic load signals at the initial stage of fatigue loading when the structure is in a healthy state.

GW Features Extraction under Uncertainty
For T1, all DIs extracted from channel 2-1, channel 2-3, and channel 4-1 under timevarying dynamic load are illustrated in Figure 8a. The variation of the front crack length and the side crack length with the number of load cycles is illustrated in Figure 8b.    However, the structure is rarely only affected by static load without time-varying conditions in real engineering applications. Therefore, crack monitoring under timevarying dynamic load will be discussed. The direct wave signal segment of S 0 mode is intercepted to calculate the DI under dynamic load by the Equation (1). t 0 and t 1 of the signal segment are 0.8 × 10 −4 s and 1.05 × 10 −4 s, respectively. The baseline signal b(t) is the average of the first 10 dynamic load signals at the initial stage of fatigue loading when the structure is in a healthy state.
For T 1 , all DIs extracted from channel 2-1, channel 2-3, and channel 4-1 under timevarying dynamic load are illustrated in Figure 8a. The variation of the front crack length and the side crack length with the number of load cycles is illustrated in Figure 8b. However, the structure is rarely only affected by static load without time-varying conditions in real engineering applications. Therefore, crack monitoring under time-varying dynamic load will be discussed. The direct wave signal segment of S0 mode is intercepted to calculate the DI under dynamic load by the Equation (1). t0 and t1 of the signal segment are 0.8 × 10 −4 s and 1.05 × 10 −4 s, respectively. The baseline signal b(t) is the average of the first 10 dynamic load signals at the initial stage of fatigue loading when the structure is in a healthy state.
For T1, all DIs extracted from channel 2-1, channel 2-3, and channel 4-1 under timevarying dynamic load are illustrated in Figure 8a. The variation of the front crack length and the side crack length with the number of load cycles is illustrated in Figure 8b.    As shown in Figures 8-10, the effect of time-varying conditions on GW features is random and uncertain. Therefore, the traditional DI method has limitations in crack quantification. However, there is a cumulative migration trend in the effect of cracks on GW features.

Establishment and Stability of MdUI-GMM
Considering that the traditional DI method has limitations in crack quantification under time-varying conditions. In this section, MdUI-GMM is used to model the probability distribution of GW features and suppress the influence of uncertainty.  As shown in Figures 8-10, the effect of time-varying conditions on GW features is random and uncertain. Therefore, the traditional DI method has limitations in crack quantification. However, there is a cumulative migration trend in the effect of cracks on GW features.

Establishment and Stability of MdUI-GMM
Considering that the traditional DI method has limitations in crack quantification under time-varying conditions. In this section, MdUI-GMM is used to model the probability distribution of GW features and suppress the influence of uncertainty. As shown in Figures 8-10, the effect of time-varying conditions on GW features is random and uncertain. Therefore, the traditional DI method has limitations in crack quantification. However, there is a cumulative migration trend in the effect of cracks on GW features.

Establishment and Stability of MdUI-GMM
Considering that the traditional DI method has limitations in crack quantification under time-varying conditions. In this section, MdUI-GMM is used to model the probability distribution of GW features and suppress the influence of uncertainty.

The Construction Process of MdUI-GMM
For T 1 and T 2 , all the DIs extracted from channel 2-1, channel 2-3, and channel 4-1 are used as the coordinate values of the three dimensions to draw the spatial scatter diagram, which is shown in Figure 11. It can be seen that DIs of different crack lengths are overlapped seriously and the spatial distribution of the GW feature sample set has a cumulative migration trend under time-varying dynamic load. which is shown in Figure 11. It can be seen that DIs of different crack lengths are overlapped seriously and the spatial distribution of the GW feature sample set has a cumulative migration trend under time-varying dynamic load.
(a) (b) The first 30 multi-dimensional GW features extracted from 30 dynamic GW signals are used as a GW baseline sample set to construct a baseline MdUI-GMM. Once a GW signal is acquired during the online damage monitoring process, new GW features can be obtained, and the GW feature sample set is updated first. Then, the online monitoring MdUI-GMM can be constructed by the EM algorithm with a uniform initialization process. Finally, PMI can be calculated. The number of GCs is set to be c = 2 to cover the GW feature sample set. The total number of GW feature sample set used to construct MdUI-GMM is R = 30. The baseline MdUI-GMM of the typical specimen is shown in Figure 12a. Based on the migration method, some typical on-line monitoring MdUI-GMMs accompanying the crack propagation are shown in Figure 12b, c, d, respectively. The yellow regions indicate the PDF of GC 1 and the purple regions indicate the PDF of GC 2. These yellow and purple ellipsoids represent contour plots of 3-dimensional Gaussian distributions. The parameters of the ellipsoids are determined by the mean and covariance of the Gaussian distribution which account for 95% of the sample set (two standard deviations from the mean). The blue points indicate the GW features acquired from the multi-channel under dynamic loading. The x-axis coordinates of the blue points are DINCM of channel 2-1. The y-axis coordinates are DINCM of channel 2-3. The z-axis coordinates are DINCM of channel 4-1. DINCM can be calculated by Equation (1). The first 30 multi-dimensional GW features extracted from 30 dynamic GW signals are used as a GW baseline sample set to construct a baseline MdUI-GMM. Once a GW signal is acquired during the online damage monitoring process, new GW features can be obtained, and the GW feature sample set is updated first. Then, the online monitoring MdUI-GMM can be constructed by the EM algorithm with a uniform initialization process. Finally, PMI can be calculated. The number of GCs is set to be c = 2 to cover the GW feature sample set. The total number of GW feature sample set used to construct MdUI-GMM is R = 30. The baseline MdUI-GMM of the typical specimen is shown in Figure 12a. Based on the migration method, some typical on-line monitoring MdUI-GMMs accompanying the crack propagation are shown in Figure 12b, c, d, respectively. The yellow regions indicate the PDF of GC 1 and the purple regions indicate the PDF of GC 2. These yellow and purple ellipsoids represent contour plots of 3-dimensional Gaussian distributions. The parameters of the ellipsoids are determined by the mean and covariance of the Gaussian distribution which account for 95% of the sample set (two standard deviations from the mean). The blue points indicate the GW features acquired from the multi-channel under dynamic loading. The x-axis coordinates of the blue points are DI NCM of channel 2-1. The y-axis coordinates are DI NCM of channel 2-3. The z-axis coordinates are DI NCM of channel 4-1. DI NCM can be calculated by Equation (1).
When the crack propagates, the form of the distribution of the GCs in the GMM undergoes a cumulative migration change. That is, the PDF of monitoring GMM undergoes a cumulative migration change compared to the PDF of the baseline GMM. The PMI between monitoring GMM and baseline GMM can be calculated by Equation (13). In combination with the measured crack length, a calibration can be performed, which will be discussed in Section 3. When the crack propagates, the form of the distribution of the GCs in the GMM undergoes a cumulative migration change. That is, the PDF of monitoring GMM undergoes a cumulative migration change compared to the PDF of the baseline GMM. The PMI between monitoring GMM and baseline GMM can be calculated by Equation (13). In combination with the measured crack length, a calibration can be performed, which will be discussed in Section 3.4.

Stability Comparison of Different Initialization Methods
To analyze the stability of the proposed MdUI-GMM, a standard deviation (SD)based statistical analysis method is applied. First, MdUI-GMM is performed 100 times with the training dataset of the T1. In the case of multi-dimensional (3D) GMM, each GC corresponds to a weight, and the mean value of each GC is a vector of 3 × 1 and the covariance is a matrix of 3 × 3. The SDs of the MdUI-GMM parameter sets are then calculated according to Equations (15) to (17).

Stability Comparison of Different Initialization Methods
To analyze the stability of the proposed MdUI-GMM, a standard deviation (SD)based statistical analysis method is applied. First, MdUI-GMM is performed 100 times with the training dataset of the T 1 . In the case of multi-dimensional (3D) GMM, each GC corresponds to a weight, and the mean value of each GC is a vector of 3 × 1 and the covariance is a matrix of 3 × 3. The SDs of the MdUI-GMM parameter sets are then calculated according to Equations (15)- (17).
where c represents the cth GC, c = 1, C. where w c , µ c,1 , and Σ c,11 are the average of weight, mean vector, and covariance matrix, respectively, of the cth component among all of the 100 parameters sets of MdUI-GMM. i represents the ith modeling process, i = 1, 2, . . . , 100. Table 2 shows the calculated SD of a uniform initialization process. The number of GCs is set to be C = 2 to cover the GW feature sample set. It is observed that the values of SD are all less than 10 −15 , which can be considered to be zero. Hence, it can be concluded that the modeling results of the MdUI-GMM are always stable.  Table 3 shows the calculated SD of a k-means initialization process. The number of GCs is set to be C = 2 to cover the GW feature sample set. It is observed that the values of SD are higher than those of the MdUI-GMM, which means that the k-means initialization method is relatively less stable than the uniform initialization method.

Calibration between PMI and Crack Length Using Prior Specimens
The first three specimen T 1 -T 3 are used as prior samples for quantitative calibration, and T 4 is used for testing. For T 1 -T 3 , the GW features of channel 2-1, channel 2-3, and channel 4-1 are fused to construct MdUI-GMM. During the monitoring process, PMI between monitoring MdUI-GMM and baseline MdUI-GMM is calculated to realize quantitative calibration of front crack length. The variation of PMIs with the length of the front crack for T 1 -T 3 are shown in Figure 13.
For T 1 -T 3 , the relationship between PMI and crack length is smoothed by a secondorder polynomial and then averaged. After that, the quantitative calibration relationship between PMI and the front crack length is established.
As shown in Figure 14, the equation between PMI and the front crack length is obtained as Equation (18) where l represents the front crack length. For T1-T3, the relationship between PMI and crack length is sm order polynomial and then averaged. After that, the quantitative cal between PMI and the front crack length is established.
As shown in Figure 14, the equation between PMI and the fron tained as Equation (18) = × + × +   For T1-T3, the relationship between PMI and crack length is smo order polynomial and then averaged. After that, the quantitative calib between PMI and the front crack length is established.
As shown in Figure 14, the equation between PMI and the front tained as Equation (18) = × + × +

Online MdUI-GMM Crack Propagation Monitoring Result
During the loading process, the GW signal is obtained online. Th MdUI-GMM can be established. Next, PMI between monitoring MdUI-MdUI-GMM can be calculated. Through the calculated PMI, combine tion results above, the front crack length of T4 can be output in real-tim ation of PMI during fatigue loading cycles is shown in Figure 15a. Re GMM measurement of the crack length are in good agreement with the as illustrated in Figure 15b.

Online MdUI-GMM Crack Propagation Monitoring Result
During the loading process, the GW signal is obtained online. Then, the monitoring MdUI-GMM can be established. Next, PMI between monitoring MdUI-GMM and baseline MdUI-GMM can be calculated. Through the calculated PMI, combined with the calibration results above, the front crack length of T 4 can be output in real-time. For T 4 , the variation of PMI during fatigue loading cycles is shown in Figure 15a. Results of the MdUI-GMM measurement of the crack length are in good agreement with the actual crack length, as illustrated in Figure 15b. The maximum absolute error of MdUI-GMM online crack length measurement is 3.4 mm and the minimum absolute error is 0.1mm, are illustrated in Table 4. It can be noted that on-line crack quantification can be realized by measuring PMI between monitoring MdUI-GMM and baseline MdUI-GMM. MdUI-GMM based SHM method can be used for on-line crack monitoring without aircraft grounding or structural disassembly and inspection.

Comparison Between Single-Channel and Multi-Channel Fusion
The calibration relationship between PMI and the crack length using single-channel GMM is shown in Figure 16a. Results of the single-channel GMM based online crack quantification and actual crack length of T4 are illustrated in Figure 16b. The actual crack length, MdUI-GMM based online crack quantification results, and the absolute error of diagnosis are illustrated in Table 5. The maximum absolute error of MdUI-GMM online crack length measurement is 3.4 mm and the minimum absolute error is 0.1 mm, are illustrated in Table 4. It can be noted that on-line crack quantification can be realized by measuring PMI between monitoring MdUI-GMM and baseline MdUI-GMM. MdUI-GMM based SHM method can be used for on-line crack monitoring without aircraft grounding or structural disassembly and inspection.

Comparison between Single-Channel and Multi-Channel Fusion
The calibration relationship between PMI and the crack length using single-channel GMM is shown in Figure 16a. Results of the single-channel GMM based online crack quantification and actual crack length of T 4 are illustrated in Figure 16b. The actual crack length, MdUI-GMM based online crack quantification results, and the absolute error of diagnosis are illustrated in Table 5.  Compared with Figures 14 and 16a, the scatter of the relationship between PMI and crack length of single-channel is higher than multi-channel. In addition, as shown in Table  5, the square sum of the error of multi-channel GMM measurement is smaller than that of single-channel GMM. From this, it is clear that multi-channel fusion improves the quantification compared to single-channel measurement. This approach is quite promising for use in engineering applications.

Conclusions
In this paper, a MdUI-GMM is proposed to improve the accuracy and stability of crack quantification under uncertainty. This method integrates multi-channel GW features to establish multi-dimensional GMM considering that different GW channels have different sensitivity to crack. Moreover, a uniform initialization method is adopted to improve the constructing stability of multi-dimensional GMM. In addition, the relationship between the PMI and crack length is calibrated with fatigue tests on prior specimens. Finally, the proposed method is applied for online crack quantification on notched speci-  Compared with Figures 14 and 16a, the scatter of the relationship between PMI and crack length of single-channel is higher than multi-channel. In addition, as shown in Table 5, the square sum of the error of multi-channel GMM measurement is smaller than that of single-channel GMM. From this, it is clear that multi-channel fusion improves the quantification compared to single-channel measurement. This approach is quite promising for use in engineering applications.

Conclusions
In this paper, a MdUI-GMM is proposed to improve the accuracy and stability of crack quantification under uncertainty. This method integrates multi-channel GW features to establish multi-dimensional GMM considering that different GW channels have different sensitivity to crack. Moreover, a uniform initialization method is adopted to improve the constructing stability of multi-dimensional GMM. In addition, the relationship between the PMI and crack length is calibrated with fatigue tests on prior specimens. Finally, the proposed method is applied for online crack quantification on notched specimens of an aircraft spar with complex fan-shaped cracks under uncertainty. By integrating the multi-channel GW features, the square sum of the error is reduced from 38.57 to 23.30. The results show that multi-channel fusion improves the quantification compared to the single-channel measurement of the crack length. Moreover, the standard deviations of the MdUI-GMM parameter sets are all less than 10 −15 , which is much smaller than the traditional k-means initialization method when performing 100 times. The results show that the construction results of the uniform initialization GMM are always stable.