3.1

Gaussian mixture modeling

The input variable x and the output variable m are combined into a new vector of variables y = [xT mT], then the training set is Y = [y₁, y₂, …, y_K]. A Gaussian mixture model of the joint probability density of the inputs and outputs is built [18]:

where M is the number of unknown Gaussian components, j represents the jth Gaussian component, p(y) represents the probability density function of y, and p(j) represents the probability density function of j, which is calculated as follows:

where 𝛾_𝑗 is the probability weight of the jth Gaussian component, 𝛿_𝑗 = {𝜂_𝑗, Σ_𝑗},𝑓(𝑗|𝛿_𝑗) is calculated as follows:

where e is the dimension of y, where 𝜂_𝑗 is the mean vector of the jth Gaussian component, and Σ_𝑗 is the covariance matrix of the jth Gaussian component. p(y/j) represents the conditional probability density function of y with respect to j:

3.2

Acquisition of model parameters

It is first necessary to determine the range of M, i.e., to determine the maximum value, M_max, and the minimum value, M_min, of M. The absolute increasing log-likelihood criterion (AIL) is chosen to select the optimal Gaussian component [19]. The AIL criterion is defined as

where log F and log F denote the maximum log-likelihood function of a Gaussian mixture model with M and M-1 Gaussian components, respectively, and the maximum log-likelihood function of a Gaussian mixture model with M Gaussian components is computed as follows:

The maximum log-likelihood function for a Gaussian mixture model with M-1 Gaussian components is similar. where K is the number of samples and is the total parameters in the complete Gaussian mixture model with M components [20].

Let M be equal to M_max and use the K-means algorithm to initialize the parameters of the Gaussian mixture model with M Gaussian components, denoted as .

3.2.1

EM algorithm [21] for obtaining estimated model parameters

E step, calculate the posterior probability of the ith training sample for the jth Gaussian component at the kth iteration:

M step, the probability weight , mean and covariance of the jth Gaussian component at the (k+1)st iteration are computed:

Let , denote the (k+1)th total model parameter. The total parameters obtained from each iteration are combined to obtain the sequence {α}, where each element is α. Maximize the log-likelihood function to find the model parameters :

3.2.2

Calculate the corresponding absolute increasing log-likelihood criterion AIL

Remove the least probable Gaussian component and merge it to the closest Gaussian component. Select the least likely component by finding the component with the smallest mixing probability:

Then the component s closest to the rth component is chosen.

where the B_s function is the symmetric KL scatter between the rth Gaussian component and the other components, and the symmetric KL scatter is a common measure of the difference between probability density functions:

The rth and sth Gaussian components are combined into a single component whose mixed probability, mean, and covariance are updated as

Obtain a Gaussian mixture model with M-1 Gaussian components with initial parameters .

Let M = M-1 and repeat EM algorithm flow if M ≥ M_min, otherwise proceed to the next step.

Obtain the optimal M_opt by maximizing the AIL.

The final parameter estimates of the model are obtained.

3.3

Computational model output

Assume that the new input data is x_new and the true value of the corresponding output is m_new. divide the mean vector and covariance matrix of each Gaussian component into input and output parts as follows:

Estimate the posterior probability of each Gaussian component of the new input data:

Compute the conditional distribution of the output of each Gaussian component of the new input data with respect to the input:

where and are the mean and covariance parameters of the jth conditional Gaussian distribution, calculated as follows:

The conditional distributions of the test set output with respect to the inputs are combined into a single Gaussian distribution. Finally, each component is weighted and combined to obtain the final prediction of the output:

Set the threshold function to

Tunnel cable equipment fault state diagnosis, online monitoring equipment, and inspection robots collected cable equipment running current, grounding current, partial discharge, cable network temperature, and other monitoring quantities as the input to the cable fault state diagnosis model, to get the model output value ŷ. When ŷ=0, the cable is in the normal state, when ŷ=1, the cable is in a weak fault 1 state, when ŷ=2, the cable is in a weak fault 2 state, when ŷ=n₁, the cable is in a permanent fault 1 state, when ŷ=n₂, the cable is in a permanent fault 2 state, and so on.

The model is updated online when new training points are added, and the process can be regarded as a process of adding new training points one by one. When the number of samples in the training set is increased to K+1, let the added samples be y_K+1=(x+1, m_K+1), then the mean and covariance parameter of the jth conditional Gaussian distribution are updated as:

where is Gaussian white noise with zero variance. The overall flow of the algorithm is shown in Fig. 1.

Figure. 1

General flow of the methodology

4.1

Evaluation indicators

The common metrics of the recognition model are accuracy, precision, recall, F1 comprehensive evaluation metrics, etc., among which accuracy is the most commonly used evaluation metric in fault recognition, which is calculated as shown below:

where N_t is the number of samples correctly predicted by the recognition model and N_f is the number of samples incorrectly predicted by the recognition model.

The accuracy can only describe the recognition results of the recognition model for all test samples, but cannot reflect the recognition performance of the model in specific categories. In order to comprehensively assess the performance of the recognition model, the F1 comprehensive evaluation index is introduced for comprehensive evaluation.

The accuracy is calculated as follows:

where N_Tk is the number of samples correctly predicted by the discriminative model in the kth category sample, and N_Fk is the number of samples incorrectly predicted by the discriminative model in the kth category sample.

The F1 composite evaluation index is calculated as follows:

The larger the value of the F1 comprehensive evaluation index, the stronger the comprehensive performance of the identification model.

4.2

Experimental setup

There are 8 cable operation states set up, which are coded as shown in Table 1. The data obtained through the cable and fault modeling simulation consists of the original dataset, which has 12000 samples, of which 2000 samples obtained by changing the states other than the eight cable operation states in the table are used as test set 2, 8000 samples are selected from the remaining 1000 samples as the training set, and the remaining 2000 samples are used as test set 1.

Table 1.

Cable operation status code

SVM, DT, RF, and KNN diagnostic models are used as a comparison of the proposed method. The simulation is carried out on MATLAB and the processor of the simulation device is 5800X with 32G RAM.

4.3

Simulation results and analysis

Five different methods of tunnel cable equipment fault state diagnosis models were trained using the training set, and the diagnostic results using test set 1 as the test set are shown in Fig. 2 and Table 2.

Figure. 2

Confusion Matrix for Cable Condition Diagnostic Results

Table 2.

Failure evaluation indicator values for each method

In order to evaluate the generalization ability of the tunnel cable equipment fault state diagnosis models trained by each method, test set 2 is used as the test set. The diagnostic results are shown in Fig. 3 and Table 3.

Figure. 3

Confusion Matrix for Cable Condition Diagnostic Results

Table 3.

Failure evaluation indicator values for each method

The simulation results of each method show that the proposed sparse online hybrid Gaussian multi-classification algorithm-based tunnel cable state fast identification method has higher identification accuracy than the SVM, DT, RF, and KNN-based tunnel cable state identification methods, especially in the new operating conditions, the cable fault identification algorithms based on MLRM, SVM, and ANN produce larger measurement errors, but the proposed recognition algorithm has a smaller increase in the cable fault state recognition error, which indicates that it is more capable of adapting to new operating conditions, i.e., the generalization ability of this soft sensor is stronger.

The hybrid Gaussian process is the union of multiple Gaussian components, and the number of Gaussian components is the number of different operating conditions of the system, so different from the general global identification, the hybrid Gaussian will automatically identify the number of operating conditions, and judge the current operating conditions during the actual measurement, so as to select the corresponding Gaussian components. Due to the huge difference between different operating conditions, global regression methods often cannot maintain high prediction accuracy over the whole operating domain, and the hybrid Gaussian process divides the whole operating domain and maintains good prediction performance in each operating condition.

Since it is generally not possible for the training set to contain all cable states, the online updating capability of the model becomes more important. In this paper, we introduce an online updating mechanism for the model, which can introduce new training data and update the model faster, without the need to start the whole model training process from scratch. Therefore, the proposed method can still maintain a good recognition accuracy under new operating conditions.

The model training time of each method is shown in Table 4:

Table 4.

Model training time for each method (unit: s)

From the above table, it can be seen that the model training time of the proposed method is the shortest. Since the training time of the hybrid Gaussian multi-classification algorithm has a linear relationship with the third power of the size of the training samples, the introduction of the sparse technique in the proposed method in this paper can greatly reduce the training time of the model in the case of a small increase in the error, which provides the possibility of online identification of the cable fault status.