**Author Contributions:** Conceptualization and methodology, M.A.-M.S. and S.A.; software and validation, G.K. and M.A.; formal analysis, M.A. and M.A.-M.S.; writing—original draft preparation, M.A.-M.S.; writing—review and editing, G.K., M.A., and B.A.; supervision, M.A.-M.S.; project administration, S.A.; funding acquisition, S.A. and B.A. Reviewed by all authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received external funding for publication by the Deanship of Scientific Research at Majmaah University, which is funding this work under project RGP-2019-29.

**Acknowledgments:** The authors would like to thank Mamdouh Mahfouz, Professor of Radiology at Kasr-El-Aini School of Medicine, Cairo University, Egypt, for his help in revising the introduction in order to correct any misused medical terms.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Al-Shaikhli, S.D.S.; Yang, M.Y.; Rosenhahn, B. Automatic 3d liver segmentation using sparse representation of global and local image information via level set formulation. *arXiv* **2015**, arXiv:abs/1508.01521.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **HEMIGEN: Human Embryo Image Generator Based on Generative Adversarial Networks**

#### **Darius Dirvanauskas 1, Rytis Maskeliunas ¯ 1, Vidas Raudonis 2, Robertas Damaševiˇcius 3,4,\*and Rafal Scherer 5**

Received: 26 July 2019; Accepted: 14 August 2019; Published: 16 August 2019

**Abstract:** We propose a method for generating the synthetic images of human embryo cells that could later be used for classification, analysis, and training, thus resulting in the creation of new synthetic image datasets for research areas lacking real-world data. Our focus was not only to generate the generic image of a cell such, but to make sure that it has all necessary attributes of a real cell image to provide a fully realistic synthetic version. We use human embryo images obtained during cell development processes for training a deep neural network (DNN). The proposed algorithm used generative adversarial network (GAN) to generate one-, two-, and four-cell stage images. We achieved a misclassification rate of 12.3% for the generated images, while the expert evaluation showed the true recognition rate (TRR) of 80.00% (for four-cell images), 86.8% (for two-cell images), and 96.2% (for one-cell images). Texture-based comparison using the Haralick features showed that there is no statistically (using the Student's t-test) significant (*p* < 0.01) differences between the real and synthetic embryo images except for the sum of variance (for one-cell and four-cell images), and variance and sum of average (for two-cell images) features. The obtained synthetic images can be later adapted to facilitate the development, training, and evaluation of new algorithms for embryo image processing tasks.

**Keywords:** deep learning; neural network; generative adversarial network; synthetic images

#### **1. Introduction**

Deep neural networks (DNN) have become one of the most popular modern tools for image analysis and classification [1]. One of the first accurate implementations of DNN, AlexNet [2], was quickly bested by a deep convolutional activation feature (DeCAF) network, which extracted features from AlexNet and evaluated the efficacy of these features on generic vision tasks [3]. VGGNet demonstrated that increasing the depth of convolutional neural network (CNN) is beneficial for the classification accuracy [4]. However, the deeper neural network is more difficult to train. ResNet presented a residual learning framework to ease the training for very deep networks [5]. The ResNet architecture, which combined semantic information from a deep, coarse layer with appearance information from a shallow network, was adapted for image segmentation [6]. Region CNN (R-CNN) was proposed as a combination of high-capacity CNNs with bottom-up region proposals in order to localize and segment images [7].

The breakthrough of DNNs in quality led to numerous adaptations in solving medical image processing problems such as the analysis of cancer cells [8] and cancer type analysis [9]. Deep max-pooling CNNs were used to detect mitosis in breast histology images using a patch centered on the pixel as context [10]. The U-Net architecture with data augmentation and elastic deformations achieved very good performance on different biomedical segmentation applications [11]. A supervised max-pooling CNN was trained to detect cell pixels in regions that are preselected by a support vector machine (SVM) classifier [12]. After a pre-processing step to remove artefacts from the input images, fully CNNs were used to produce the embryo inner cell mass segmentation [13]. A set of Levenberg–Marquardt NNs trained using textural descriptors allowed for predicting the quality of embryos [14].

Image analysis methods were also applied to embryo image analysis. Techniques for extracting, classifying, and grouping properties were used to measure the quality of embryos. Real-time grading techniques for determining the number of embryonic cells from time-lapse microscope images help embryologists to monitor the dividing cells [15]. Conditional random field (CRF) models [16] were used for cell counting and assessing various aspects of the developing embryo and predicting the stage of the embryonic development. In addition to grading, embryonic positioning was achieved using linear chain Markov model [17]. Different solutions were developed for segmenting and calculating the number of cells: ImageJ, MIPAV, VisSeg [18]. Moreover, cell segmentation by marking its center and edges makes it possible to determine their shapes and quantities more quickly and accurately [19]. Two-stage classifier for embryo image classification was proposed in [20].

The rising progress of generative adversarial networks (GANs) applications on medical imaging state that most of the research is focused on synthetic imaging, the reconstruction, segmentation, and classification and proves the importance of the sector. Wasserstein-based GANs were applied for the synthesis of cells imaged by fluorescence microscopy capturing these relationships to be relevant for biological application [21]. The quality of generated artificial images also improved due to the improvements to deep neural networks [22] as well as on progressive growing of GANs for image data augmentation [23]. Variational autoencoders (VAE) and generative adversarial networks (GANs) are currently the most distinguished according to the quality of their results. The VAE models are more often used for image compression and recovery. Both methods have drawbacks, as recovering a part of the information loses some data and often introduces image fading or blurring. This effect could be reduced by matching the data as well as the loss distributions of the real and fake images by a pair of autoencoders used as the generator and the discriminator in the adversarial training [24]. Utilization of both generator and discriminator growing progressively from a low-resolution standpoint and adding new layers for an increase in fine details allowed for achieving a current state-of-the-art inception score of 8.80 on CIFAR10 dataset [25]. Network modifications can be applied to minimize fading of the resulting images [24]. Generative stochastic networks (GSN) were used to learn the transition operator of a Markov chain whose stationary distribution estimates the data distribution [26]. Deep recurrent attentive writer (DRAW) neural network architecture [27] was used to generate highly realistic natural images such as photographs of house numbers or other digits in the MNIST database. Dosovitskiy et al. [28] introduced an algorithm to find relevant information from the existing 3D chair models and to generate new chair images using that information. Denton et al. [29] introduced a generative parametric model capable of producing high-quality samples of natural images based on a cascade of convolutional networks within a Laplacian pyramid framework. Radford and Metz [30] introduced a simplification of training by a modification called deep convolutional GANs (DCGANs). GANs can also recover images from bridging text and image modelling, thus translating visual concepts from characters to pixels [31]. Zhang et al. [32] proposed stacked GANs (StackGAN) with conditioning augmentation for synthesizing photo-realistic images. GANs could also be applied for generating motion images, e.g., the motion and content decomposed GAN (MoCoGAN) framework for video generation by mapping a sequence of random vectors to a sequence of video frames [33]. The GAN for video (VGAN) model was based on a GAN with a spatial-temporal convolutional architecture that untangles the scene's foreground from the background and can be used at predicting plausible futures of static images [34]. Temporal GAN (TGAN) was used to learn a semantic representation of unlabeled videos by using different types of generators via Wasserstein GAN, and a method to train it stably in

an end-to-end manner [35]. MIT presented a 3D generative adversarial network (3D-GAN) model to recreate 3D objects from a probabilistic space by leveraging recent advances in volumetric CNNs and GANs [36]. Li et al. [37] used multiscale GAN (DR-Net) to remove rain streaks from a single image. Zhu [38] used GANs for generating synthetic saliency maps for given natural images. Ma et al. [39] used background augmentation GANs (BAGANs) for synthesizing background images for augmented reality (AR) applications. Han et al. presented a system based on a Wasserstein GANs for generating realistic synthetic brain MR images [40].

All these results are achieved by using large expert-annotated and ready-made databases and also exposing the problem of lacking good, core training datasets. In this work, we aimed to develop a method to generate realistic synthetic images that could later be used for classification, analysis, and training, thus resulting in the creation of novel synthetic datasets for research areas lacking data such as human embryo images. Here we used human embryo images obtained during cell development processes for training a DNN. We propose an algorithm for generating one-, two-, and four-cell images (the selection was based on the initial dataset provided by our medical partners) to increase the overall number of unique images available for training. For generating images, we have developed a generative adversarial network for image recovery, filling, and improvement. The significance of our approach would be that the method was applied to embryonic cell images as a type of image. It was very important for the GAN to accurately reproduce the outline of the cell (often poorly visible even in microscopy picture), since the whole image was almost the same (gray) and the cell itself is translucent. Our focus was not only to generate the generic image of a cell as such, but to make sure that it has all necessary attributes of a real cell image to provide a fully realistic synthetic version. We believe that as a large number of real embryo images required for training neural networks are difficult to obtain due to ethics requirements, the synthetic images generated by the GAN can be later adapted to facilitate the development, training, and evaluation of new algorithms for embryo image-processing tasks.

#### **2. Materials and Methods**

#### *2.1. Architecture of the Generative Adversarial Network (GAN)*

The generative adversarial network (GAN) consists of two parts: A generator *G* and a discriminator *D* (Figure 1).

**Figure 1.** A pipeline of a generative adversarial network (GAN).

The discriminator tries to distinguish true images from synthetic, generator-generated images. The generator tries to create images with which the discriminator could be deceived. The discriminator is implemented as a fully connected neural network with dense layers that represents an image as a probability vector and classifies it into two classes, a real or a fake image. The generator is a reverse model that restores the former image from a random noise vector. During the training, the discriminator is trained to maximize the probability of assigning the correct class to the training images and images generated by the generator. At the same time, the generator is trained to minimize classification error

log(1 − *D*(*G*(*z*))). This process can be expressed in terms of game theory as a two-player (discriminator and generator) minimax game with cost function *V*(*G*, *D*) [41]:

$$\max\_{G} \min\_{D} V(D, G) = E\_{x \sim P(x)}[\log(D(x))] + E\_{z \sim P(z)}[\log(1 - D(G(z)))].\tag{1}$$

The discriminator evaluates the images created by the training database and the generator. The architecture of the discriminator network is shown in Figure 2. The network consists of six layers. The monochrome 200 × 200 pixel image in the first layer is transformed and expanded into a single line vector. The dense layer is used in the second, fourth, and sixth layers. The LeakyReLU [42] function, where α = 0.2, is used in the third and fifth layers as it allows a small gradient when the unit is not active:

$$\text{LReLu}(x) = \max(x, 0) + a \min(x, 0). \tag{2}$$

**Figure 2.** Architecture of the discriminator network.

At the network output, we get a one if the network guesses that the image is real or zero if the network guesses that the image is fake.

The generator is composed of 11 layers (Figure 3). The input layer contains 1 × 100 vector of randomly generated data. These data are transformed using the dense, LeakyReLU, and batch normalization layers. The ReLU layer uses α = 0.2. The batch normalization layer was used to normalize the activations of the previous layer at each batch, by applying a transformation that maintains the mean activation close to 0 and the activation standard deviation close to 1 [43]:

$$\mathfrak{X}\_{i} \leftarrow \frac{\mathfrak{x}\_{i} - \mu\_{\beta}}{\sqrt{\sigma\_{\beta}^{2} + \varepsilon}}.\tag{3}$$

**Figure 3.** Architecture of the generator network.

The batch normalization momentum value was set to 0.8, while all other values were used as per default Keras model values. The output produces a 200 × 200 pixel sized black-and-white image.

#### *2.2. Training*

We trained one (single) network to generate all types of cell images. For the training, we have used the Adam optimization algorithm [44] with a fixed training speed of 0.00002, and β1 = 0.7 for training the generator. During training, the Adam algorithm calculates the gradient as follows:

$$
\xi\_l g\_t \leftarrow \Delta\_\theta f\_l(\theta\_{t-1}).\tag{4}
$$

The biased first moment estimate is updated as follows:

$$
\mu m\_t \gets \beta\_1 m\_{t-1} + (1 - \beta\_1)\mathbf{g}\_t. \tag{5}
$$

The biased second moment estimate is update as follows:

$$
\upsilon\_t \leftarrow \beta\_2 \upsilon\_{t-1} + (1 - \beta\_2) g\_t^2. \tag{6}
$$

Then we compute the bias corrected first moment estimate as follows:

$$
\hat{m}\_t \leftarrow \frac{m\_t}{1 - \beta\_1^t}.\tag{7}
$$

The bias corrected second moment estimate is computed as follows:

$$
\psi\_t \leftarrow \frac{\upsilon\_t}{1 - \beta\_2^t}.\tag{8}
$$

This gives the update parameters as follows:

$$
\partial\_t \leftarrow \frac{\partial\_{t-1} - a\mathfrak{n}\_t}{\sqrt{\mathfrak{d}\_{t}^{\prime\prime} + \varepsilon}}.\tag{9}
$$

The binary cross-entropy function is used to evaluate the error as follows:

$$-\frac{1}{N}\sum\_{i=1}^{N} \left[ y\_i \log(\mathfrak{H}\_i) + (1 - y\_i)\log(1 - \mathfrak{H}\_i) \right].\tag{10}$$

#### *2.3. Evaluation of Image Quality*

The problem with using the GAN model is that there are no clear ways to measure it qualitatively. There are few most commonly used criteria to evaluate GAN network: Average log-likelihood, classifier classification [41], and visual fidelity of samples [45]. These methods have both advantages and disadvantages. Log-likelihood and classifiers attribute the image to a certain class with a matching factor. However, they only indicate whether the generated image is like the average of one of the classes but does not value a qualitatively generated image. Such a qualitative assessment requires a human evaluator, but with a large amount of data, such estimates can change independently and be biased [46]. Likewise, a different expert person could evaluate the same data somewhat differently.

We use a combined method to evaluate the results obtained in our study, which includes: (1) Human expert evaluation, (2) histogram comparison, and (3) texture feature comparison. The generated images are evaluated by human experts (three medical research professionals) to determine if they are being reproduced qualitatively. If we can determine its class from the received image, the image is quality-restored. The expert based scoring was calculated using a visual Turing test [47] as this method has already been proven effective in the evaluation of synthetic images generated by GANs [48].

The histogram comparison method checks whether the histogram of the generated images corresponds to the histogram of the images in the training database. For comparison, we used the normalized average histogram of training data *H*1, and the normalized average histogram of generated images *H*2. Further, we apply four different methods for comparison as it is recommended in [49]. The correlation function determines the likelihood of two image histograms. Its value of two identical histograms equals one:

$$\mathbb{C}(H\_1, H\_2) = \frac{\sum\_{I} \left( H\_1(I) - \overline{H}\_1 \right) \left( H\_2(I) - \overline{H}\_2 \right)}{\sqrt{\sum\_{I} \left( H\_1(I) - \overline{H}\_1 \right)^2 \sum\_{I} \left( H\_2(I) - \overline{H}\_2 \right)^2}}. \tag{11}$$

The Chi-square function takes the squared difference between two histograms. The squared differences are divided by the number of samples, and the sum of these weighted squared differences is the likelihood value:

$$\chi^2(H\_1, H\_2) = \sum\_{I} \frac{\left(H\_1(I) - H\_2(I)\right)^2}{H\_1(I)}.\tag{12}$$

The histogram intersection calculates the similarity of two discretized histograms, with a possible value of the intersection lying between no overlap and identical distributions. It works well on categorical data and deals with null values by making them part of the distribution.

$$I(H\_1, H\_2) = \sum\_{I} \min(H\_1(I), H\_2(I))\tag{13}$$

The Bhattacharyya distance approximates the normalized distance between the histograms using the maximum likelihood of two object histograms as follows:

$$B(H\_1, H\_2) = \sqrt{1 - \frac{1}{\sqrt{H\_1 H\_2 N^2}}} \sum\_I \sqrt{H\_1(I) H\_2(I)}.\tag{14}$$

The texture-based comparison uses the texture analysis features based on grey level co-occurrence matrix (GLCM), which are related to second-order image statistics that were introduced by Haralick [50] as follows: (1) Angular second moment (energy), (2) contrast, (3) correlation, (4) variance, (5) inverse difference moment (homogeneity), (6) sum average, (7) sum variance, (8) sum entropy, (9) entropy, (10) difference variance, (11) difference entropy, (12) information measure of correlation, (13) information measure of correlation II, and (14) maximal correlation coefficient.

#### *2.4. Ethics Declaration*

The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of Faculty of Informatics, Kaunas University of Technology (No. IFEP201809-1, 2018-09-07).

#### **3. Experiment and Results**

#### *3.1. Dataset and Equipment*

We have used the embryo photos taken with Esco Global incubator series, called Miri TL (time lapse). The embryo image sets were registered in the German, Chinese, and Singapore clinics. No identity data were ever provided to the authors of this paper. The Esco's embryo database consists of three image classes that together have 5000 images with a resolution of 600 × 600 pixels. Number of different images in classes: One-cell—1764 images, two-cell—1938 images, four-cell—1298 images. Images are obtained from 22 different growing embryos in up to five stages of cell evolution. The database was then magnified by rotating each photo at 90, 180, and 270 degrees. Resulting number of different images in classes: (one-cell images—1764 × 4 = 7056; two-cell images—1938 × 4 = 7752; four-cell images—1298 × 4 = 5192). The images were taken in the following culturing conditions: Temperature of 37 ◦C, stable level of 5% CO2, and the controllable values of nitrogen and oxygen mixture. The embryos were photographed in a culture coin dish made from polypropylene with a neutral media pH = 7. The inverted microscopy principle was used, with 20x lenses, without zoom, with focusing and with field of view of 350 um. A camera sensor used was IDS UI-3260CP-M-GL (IDS Imaging Development Systems GmbH, Obersulm, Germany). An example of embryo images in the Esco dataset is presented in Figure 4. The evaluation and training of GANs was done on an Intel i5-4570 CPU with a GeForce 1060 GPU and 8 GB of RAM.

**Figure 4.** Sample images from Esco embryo image dataset.

#### *3.2. Results*

The training of GAN was carried out in 200,000 iterations, with a sample size of 256 images per batch. The duration of training using GeForce 1060 GPU was about 15 h 45 min.

From the final generated images, one can easily count the number of embryonic cells and determine the class (Figure 5). During subsequent iterations, the generator restored the embryo image and reduced the amount of noise in the generated images.

**Table 1.** Evolution of embryo cell images during GAN training.

**Figure 5.** Example of final generated embryo cell images using the proposed algorithm. From left to right: One, two, and four cells. Images were filtered from "salt and pepper" noise using a median filter. See Table 1 for raw outputs.

During the training, embryo cell images were generated at every 25,000 iterations. An example of these images is shown in Table 1.

After 25,000 iterations, we can see that he managed to restore the plate's flair and release the background, but the embryo itself was not restored. After 100,000 iterations, the generator has been able to clearly reproduce the image of the embryo cell.

The error function values for both generator and discriminator are shown in Figures 6 and 7. During an initial training up to 20,000 iterations, the error function value of the generator increases rapidly. Further increase of the error function results in a slower pace. Once the training number reaches 75,000 iterations, the value of the generator error function begins to increase with the training of two-cell and four-cell images. The generator becomes able to produce more complex images, which is also visible from the generated sample images (see Table 1). The one-cell image can be identified from 50,000 iterations, whereas two-cell and four-cell images can be recognized from 75,000–125,000 iterations. From the discriminator error feature, we can see that starting from 20,000 iterations it becomes difficult to separate real images from artificially generated images, and the value of the error function is less than 0.15. A value of discriminator log loss increases, as discriminator is unable to differentiate between a real image and a generated image.

**Figure 6.** Loss log graph for the generator network.

**Figure 7.** Loss log graph for the discriminator network.

For the expert-based evaluation (see Table 2), 1500 images (500 in each class) were generated. The number of artificial one-cell images, where one cell was clearly recognized, was 96.2%. The quality of two-cell and four-cell image generation of more complex images deteriorated. Of the two-cell images, the number of cells could be clearly determined from 86.8% of the images. When evaluating the four-cell images, 80% accuracy was obtained, i.e., one out of five images was generated inaccurately (as decided by an expert).

**Table 2.** Evaluation of the generated embryo image cells by human experts.

Additionally, the image similarity was evaluated by comparing their histograms. In Figure 8, one can see the comparison of image histograms, where the blue curve is the average histogram of the training image dataset and the red is the average histogram of the artificially generated images. In all three classes, we can see that the generated image is brighter than the real one. This could be explained due to the slight "salt-and-pepper"-type noise in the generated image. The highest histogram match was obtained in the single-cell-generated images. This shows that the generator was better able to reproduce images of a simpler structure. Note that Figure 8 shows an evaluation of a generated cell image with a maximum epoch number of 200,000. The particular histogram cannot be identical as a single network generates images of different cells. This comparison was done by evaluating the average histogram of all (training and generated) image histograms.

**Figure 8.** Comparison of histograms of real vs. generated images (real images—blue, generated images—red): (**a**) one-cell images, (**b**) two-cell images, (**c**) four-cell images.

Table 3 shows the results of a comparison of normalized histograms using correlation, Chi-square, intersection, and Bhattacharyya distance formulas. In all four cases, the largest histogram coincidence was shown by the one-cell generated images, while the two-cell and four-cell images performed worse.

**Table 3.** Evaluation of generated artificial embryo cell images using histogram comparison criteria.

We also compared the values of the textural (Haralick) features for the original and generated embryo images and evaluated the statistical difference between the values of features using the Student's t-test. The results were significant (*p* < 0.01) for all Haralick features except the sum variance feature (one-cell images), variance and sum average (two-cell images), and the sum variance feature (four-cell images).

To compare the Haralick features from both (original and generated) image datasets, we also used the Pearson correlation between the principal components of feature sets extracted using principal component analysis (PCA). For both feature sets, the first principal component (PC1) explains more than 98% of variance, so we use only PC1 for further comparison. The results are presented in Table 4 and show a very high correlation between the PC vectors.

**Table 4.** Results of principal component analysis (PCA) of Haralick features from original and synthetic embryo images.

#### **4. Discussion**

As the correct comparison with other algorithms not possible due to very different source datasets used, we tried to compare within the misclassification rate, by comparing percentage of the images from all generated groups were assigned to correct class or not at all. Our HEMIGEN method was rated at 12.3%.

Table 5 provides a comparison of the misclassification rate of synthetic images when compared to the results obtained by other authors. Gaussian mixture deep generative network (DGN) demonstrated 36.02% misclassification rate [51]. DGN with auxiliary variables and two stochastic layers and skip connections achieved the 16.61% misclassification rate [52]. Semi-supervised classification and image generation with four-layer generator demonstrated 8.11% misclassification rate on house number image generation [48]. The adversarial learned inference (ALI) model, which jointly learns a generation network and an inference network using an adversarial process, reached a misclassification rate of 7.42% using CIFAR10 test set (tiny images) [53]. The WGAN-based approaches ranged from 6.9% [21] to 50% [40] depending on the application. The methods indicated are only "loosely" comparable, taking into account the differences (importance) of features of synthetic image targeted and scopes of the works by other researchers.

We also checked for the problem of mode collapse where the generator would produce the same image over and over, thus fooling discriminator, while no new image would be generated in real-life. We have not found this case in our approach, possibly due to an adequate number of different embryos used. Some generated images had embryos, which overlap one another, but this was not considered a failure, but a realistic case of cell division.

#### **5. Conclusions**

We used generative adversarial networks trained on real human embryo cell images to generate a dataset of synthetic one-, two-, and four-cell embryo images. We have achieved the highest quality of generated images for single-cell embryo images, where 96.2% of the synthetic embryo images were recognized as accurate and usable by human experts. The worst accuracy was achieved for the synthetic four-cell images, of which only 80% could be identified correctly. These results were confirmed by the histogram comparison, which achieved the highest scores for synthetic single-cell images (an average correlation of 0.995 was achieved when comparing histograms of real and synthetic one-cell embryo images), as well as by comparison of image textures analyzed using the Haralick features.

As our algorithm allows us to manipulate the size, position, and number of the artificially generated embryo cell images, these images can then be used to train and validate other embryo image processing algorithms, when the real embryo images are not available, or the number of available real embryo images is too small for training neural networks.

**Author Contributions:** Conceptualization, R.M. and V.R.; data curation, D.D.; funding acquisition, R.S.; investigation, D.D, R.M., and V.R.; methodology, R.M.; software, D.D; supervision, R.M.; validation, D.D and R.M.; visualization, D.D; writing—original draft, D.D, R.M., and R.D.; writing—review and editing, R.M. and R.S.

**Funding:** The project partially financed under the program of the Polish Minister of Science and Higher Education under the name "Regional Initiative of Excellence" in the years 2019–2022 project number 020/RID/2018/19, the amount of financing 12,000,000.00 PLN.

**Acknowledgments:** The authors also would like to thank Esco Global for kindly provided embryo image dataset. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### **Image Thresholding Improves 3-Dimensional Convolutional Neural Network Diagnosis of Di**ff**erent Acute Brain Hemorrhages on Computed Tomography Scans**

#### **Justin Ker 1, Satya P. Singh 2, Yeqi Bai 2, Jai Rao 1, Tchoyoson Lim 3and Lipo Wang 2,\***

Received: 1 April 2019; Accepted: 4 May 2019; Published: 10 May 2019

**Abstract:** Intracranial hemorrhage is a medical emergency that requires urgent diagnosis and immediate treatment to improve patient outcome. Machine learning algorithms can be used to perform medical image classification and assist clinicians in diagnosing radiological scans.  The dataset used consists of 399 volumetric CT brain images representing approximately 12,000 images from the National Neuroscience Institute, Singapore. We used a 3D CNN to perform both 2-class (normal versus a specific abnormal class) and 4-class classification (between normal, SAH, IPH, ASDH). We apply image thresholding at the image pre-processing step, that improves 3D CNN classification accuracy and performance by accentuating the pixel intensities that contribute most to feature discrimination. For 2-class classification, the F1 scores for various pairs of medical diagnoses ranged from 0.706 to 0.902 without thresholding. With thresholding implemented, the F1 scores improved and ranged from 0.919 to 0.952. Our results are comparable to, and in some cases, exceed the results published in other work applying 3D CNN to CT or magnetic resonance imaging (MRI) brain scan classification. This work represents a direct application of a 3D CNN to a real hospital scenario involving a medically emergent CT brain diagnosis.

**Keywords:** 3D convolutional neural networks; machine learning; CT brain; brain hemorrhage

#### **1. Introduction**

Intracranial hemorrhage is a medical emergency that can have high morbidity and mortality if not diagnosed and treated immediately. This condition affects 40,000 to 67,000 patients in the United States annually and up to 52% of patients die within one month [1]. Three commonly-encountered sub-types of intracranial hemorrhage are subarachnoid hemorrhage (SAH), intraparenchymal hemorrhage (IPH), acute subdural hemorrhage (ASDH). In a severe brain trauma, various permutations of SAH, IPH, and ASDH can be seen, which we have termed brain polytrauma hemorrhage (BPH) in this work. The common causes of SAH are trauma and cerebral aneurysmal rupture, while IPH can be caused

by hypertension, amyloid angiopathy, brain tumor hemorrhage or trauma. ASDH and BPH appear because of head trauma.

When patients with intracranial hemorrhage present to the emergency department, a computed tomography (CT) scan of the brain is done to diagnosis intracranial hemorrhage, so that medical or surgical treatment can follow. CT scans work by exploiting the differential absorptive properties of body tissues placed between an x-ray emitter and detector. Brain tissue, blood, muscle, and bone give rise to different levels of x-ray attenuation, expressed in Hounsfield units. By moving the x-ray emitter and detector circumferentially around a subject, a three-dimensional image of the subject's internal tissues is obtained.

Limitations in the availability or experience of clinicians, especially in rural or resource-strapped health systems, to diagnose CT brains quickly can cause treatment delays. Automating the diagnosis of CT brain scans, or assisting the clinician in triaging critical from normal scans, would help patients by expediting their treatments and improve outcomes. Figure 1 shows the axial slices of five CT brain scans showing a normal brain, SAH, IPH, ASDH, and BPH. The outer rim of uniform white skull bone surrounds the dark grey brain tissue, and areas of acute hemorrhage, appear as patchy, light-grey areas of varying shapes.

**Figure 1.** Computed tomography (CT) brain scans. From left, normal (N), subarachnoid hemorrhage (SAH), intraparenchymal hemorrhage (IPH), acute subdural hemorrhage (ASDH), brain polytrauma hemorrhage (BPH). Each image represents an individual image slice. One patient's complete stack of CT images contained between 24 to 34 image slices in our dataset.

The artificial neuron was first described by McCulloch and Pitts in 1943 [2]. This has evolved through the symbolic, rule-based artificial intelligence (AI) paradigms, to manual feature-handcrafting algorithms, and to modern multi-layered or "deep" neural networks, which perform feature detection and classification automatically. Convolutional neural networks (CNN) owe their inception to Fukushima's Neocognitron model in 1982 [3], and their popularity to Lecun et al. [4] and Krizhevsky et al. [5] The latter employed a CNN to win the 2012 Imagenet Large Scale Visual Recognition Challenge, and since then CNNs have been used for many image classification tasks. The advantage of CNNs in image classification is the ability to perform feature-extraction and learn high-level image features automatically without feature-handcrafting, leading CNNs to become the dominant machine learning architecture in image recognition tasks. CNNs have been widely used in machine vision to perform a variety of tasks, such as image classification, object detection, and semantic segmentation. In the medical image analysis space, CNNs have been used for the classification and diagnosis of 2-dimensional medical images, such as chest x-rays, retinal photographs, skin dermoscopic images, and histology images, with performance comparable to or exceeding human clinicians [6–9].

In analyzing volumetric magnetic resonance imaging (MRI) or CT data, various machine learning approaches including 2D CNN have been attempted. These efforts have involved manual slice selection, extensive manual pre-processing, feature hand-crafting, and segmentation, before classification [10–12]. A number of these approaches classify individual images from a volumetric image stack one image at a time, turning a 3D classification problem into a 2D task. Other authors have used employed a combination of 2D CNN in the three planes (axial, coronal, sagittal) that define a 3D volume [13–15]. Roth et al. [14] detected abnormal lymph nodes on thoracic CT scans, by decomposing a 3D volume of interest into re-sampled 2D views, and then training their CNN on augmented variations of these 2D views. Their CNN achieved a satisfactory sensitivity of 90%, at six false positives per patient. The disadvantages of these previous strategies include the manual time and effort in hand-crafting, feature segmentation, and stripping, and the potential loss of spatial contextual information when a 3D volume is analyzed using 2D slices.

Three-dimensional CNNs are an emerging architecture, and have been used mainly in analyzing video or 3-dimensional volumetric medical images. Previously, the use of 3D CNN was limited as it was computationally expensive and lengthy to process 3-dimensional kernels and entire volumes of images. However, more papers on 3D CNN have appeared in the scientific literature with their adoption likely aided by decreasing computational hardware costs. In medical image analysis, 3D CNNs have been applied to detecting abnormalities (tumors, hemorrhage, ischemia) in brain, heart, lung, and liver organs on CT or MRI imaging [16–18]. As a measure of the interest in the clinical problem of automatically detecting intracerebral hemorrhage on MRI or CT, there is also a growing number of publications on this topic using 3D CNNs [19,20]. Dou et al. [20] analyzed cerebral microbleeds on MRI brain scans using a two-stage 3D CNN that screened and then detected microbleeds. Their method achieved 93% sensitivity, at a cost of 44% precision and 2.7 false positives per patient.

We propose a 3D CNN that classifies CT brain scans into normal (N), subarachnoid hemorrhage (SAH), intra-parenchymal hemorrhage (IPH), acute subdural hemorrhage (ASDH), and brain polytrauma hemorrhage (BPH).

In this work, we aim to create a 3D CNN that can automatically detect and diagnose SAH, ASDH and IPH on CT brain scans, and distinguish them from normal scans. This work would have direct clinical application in emergency departments and acute-care hospitals worldwide. We propose a simple and fast 3D CNN that is effective and accurate. It is hoped that the straightforward implementation of this 3D CNN will lead to its widespread adoption. Specifically, we modify well-known 2D CNN architectures into 3D CNN. The novel aspects of our work are as follow:

(1) To our knowledge, this is the first demonstration of a 3D CNN on volumetric CT brain data that classifies patient scans into different acute hemorrhagic variants, which impacts subsequent medical treatment. Previous work has been limited to detecting normal versus abnormal scans. We also demonstrate that our network performance gives state-of-the-art results in classification accuracy.

(2) We present a novel image thresholding method optimized for the detection of acute hemorrhage on CT brain scans, which mimics the visual analysis of radiological images by human radiologists. We demonstrate that the application of our method improves the classifier performance.

This paper is organized into the following sections. Section 2 describes our methods, with details on the dataset, network architectures, and training protocols. Section 3 reports our results and experimentation with various network parameters. Section 4 analyzes the impact, limitations and future work stemming from our results. Section 5 summarizes and concludes this paper.

#### **2. Methods**

#### *2.1. 3-Dimensional Convolutional Neural Networks*

The feature map of a convolution layer is formed by convolving the feature map of the previous layer with learnable kernels, adding a bias term, and then applying an activation function [21]. Specifically, we can define *h p k* as the *kth* feature map of the *pth* layer, and *h p*−1 *j* as the *j th* feature map of

the previous layer, *p* − 1. *Wp j*,*k* is a learnable kernel, and *bp* is the bias term. σ is the activation function, commonly a rectified linear activation unit (ReLu). This is written as:

$$h\_k^p = \sigma \left( \sum\_j h\_j^{p-1} \times \mathcal{W}\_{j,k}^p + b\_k^p \right) \tag{1}$$

A 2-dimensional convolution can be defined for convolving an input image *I* with a kernel *K*. Extending the equation for a 2-dimensional convolution [22] into a 3-dimensional convolution, we obtain:

$$S(\ell, m, n) = (\mathbb{K} \times I)(l, m, n) = \sum\_{a} \sum\_{b} \sum\_{c} I(\ell - a, m - b, n - c) \text{ K}(a, b, c) \tag{2}$$

Equation (2) may be expressed as:

$$S(\ell, m, n) = \sum\_{a, b, c} h\_j^{p-1} \left( \ell - a, m - b, n - c \right) \mathcal{W}\_{j,k}^{p} (a, b, c) \tag{3}$$

In Equation (3), *h p*−1 *j* is the *j th* 3-dimensional feature map of the previous layer *p* − 1, and *Wp j*,*k* is the 3-dimensional kernel. Substituting Equations (1) and (3), we obtain the equation for *h p k* , which is a 3-dimensional feature map:

$$\mathcal{H}\_k^p = \sigma \left( S\_{j,k}^p + b\_k^p \right) \tag{4}$$

These 3-dimensional convolution layers were stacked with max-pooling and fully-connected layers. Kernels were initialized with Gaussian distribution and parameters were tuned with standard stochastic gradient descent and cross-entropy loss minimization.

#### *2.2. Image Thresholding to Detect Acute Hemorrhage*

The CNN extracts features from the input images, and most of the features are associated with edges, shape, and curves present in the input images. It can also be seen by visualizing the different layers of CNN that the first layer mostly picks up edges present in the images, while the second layer picks up curves, and the third layer picks up shapes. The human eye is similar in that it also observes an image for its constituent edges, curves, and shapes, as suggested by the presence of ocular dominance and orientation columns in the primary visual cortex of the mammalian brain. In detecting anomalies on medical images, a region of interest (ROI) may be subtle or not apparent to visual inspection by either human or CNN. Just as thresholding aids a human radiologist to emphasize possibly abnormal areas, we propose a threshold operator to accentuate individual sharp edges, such that this improves the likelihood of a CNN detecting a feature, and, therefore, performing a correct classification.

To improve the discriminatory ability of our model, we propose a novel thresholding method to detect acute hemorrhage. We build on the work of Zhang et al. [23] who used spatial histograms to detect cars in images. There is a range of pixels intensities which is common in both normal and abnormal CT scans. Therefore, we can discard these intensities before feeding the image to the CNN without any loss of information. This is similar to the dimensionality reduction process. In our proposed method, we generate the average pixel intensity histograms of 2 classes of CT scans (such as Normal and SAH), to overlap and search for an optimal pixel intensity threshold that accentuates their difference. Pixel intensities below this threshold are then disregarded. This process creates sharp edges and shapes around a ROI which helps CNN in feature extraction. An added benefit is that a CNN will require less time for training.

The intuition underlying this approach is two-fold. First, we observed that normal CT brain scans, even across different patients are largely homogeneous, which makes a histogram representation of the class meaningful. Abnormal scans can be thought of as the addition of extraneous blood signal to normal scans. Second, we also observed that radiologists adjust image contrast levels when reading CT scans, to accentuate subtle amounts of blood, and to downplay the appearance of normal surrounding

tissue. Our method attempts to model this behavior in human visual cognition. Using our method improves the overall performance of the classifier across all classes of acute brain hemorrhage.

#### *2.3. Dataset and Pre-Processing*

The dataset consists of non-contrast CT brain images from the National Neuroscience Institute at Tan Tock Seng Hospital, Singapore. After institutional review board approval, a search of electronic discharge summaries with the diagnoses matching head injury or intracranial hemorrhage was performed. The resultant list of patient identifiers was used to query and retrieve relevant scans from the hospitals' Picture archiving and communication system (PACS) servers. Each scan was anonymized and manually checked by the authors to ensure ground truth. The final dataset consisted of 399 unique patient volumetric CT brain scans representing approximately 12,000 images, summarized by the different classes in Table 1. These scans had varying numbers of image slices and slice thickness, owing to variability in CT scanner model and scanning protocols. We prepare the data for 5 five-fold cross-validation. Training, validation, and test CT scan images were augmented eight-fold by flipping along the vertical axis, and rotation up to 45 degrees.

**Table 1.** Number of original unique patient computed tomography (CT) scans.

#### *2.4. Network Architecture*

Table 2 shows the model architecture used in our 3D CNN. We experimented with various model architectures including VGGNet and GoogLeNet, to optimize for the trade-off between classification accuracy and computational time. We aimed to have a model with a straightforward design for easy trouble-shooting and to facilitate real-world implementation.  However, while they opted for a two-part ensembled screening and discrimination stages, we opted for a single throughput architecture for simplicity and performance. Figure 2 is a pictorial diagram of our proposed 3D CNN. After the necessary preprocessing steps, input volumetric data of 3D CT scans are passed through a pre-defined threshold operator as discussed in Section 2.2 and becomes the input to the 3D CNN.

**Table 2.** Model architecture of the 3-dimensional convolutional neural networks (3D CNN) used in this work.

**Figure 2.** Proposed architecture for binary and multi-class classification of CT scans. The features are visualized using 3D deconvolution visualization methods at each pooling layer.

#### *2.5. Training and Implementation*

Training was performed on a computer with two Intel Xeon E5-2630 CPU processors, four NVIDIA GTX 1080 Ti Graphical processing units, and 128 GB of DDR4 Random access memory. The project was implemented using the Python programming language and the Google Tensorflow library. We used the rectified linear unit (ReLu) as the activation function, the Adam optimizer, and cross-entropy as the loss function. We used the grid search approach for optimizing the learning rate, dropout, and kernel size of the convolution and pooling layers. We varied learning rates from (0.1, 0.001, 0.0001), and dropout from (0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6) after each convolution layer and fully connected layer. The kernel sizes of convolution layers were varied from (1 × 1 × 1, 2 × 2 × 2, 3 × 3 × 3, 5 × 5 × 5), and pooling kernel sizes ranged from (1 × 1 × 1, 2 × 2 × 2). Eventually, we set the learning rate at 0.001 and set the dropout after the fully connected layer at 0.2. We employed a convolutional kernel size of 3 × 3 × 3 at each layer, with a 2 × 2 × 2 pooling kernel after each convolutional layer. We set β1 = 0.001, β2 = 0.999; ε = 10−8.

#### Metrics

To evaluate our network performance, we measured the Sensitivity (S), Precision (P) and F1 scores across each classification task. TP, FP, FN refer to true positive, false positive, and false negative, respectively. The F1 score is defined as the harmonic average of Sensitivity and Precision and is a measure of a test's accuracy. An F1 score of 1 indicates perfect Sensitivity and Precision, while a score of 0 indicates the opposite.

$$S = \frac{TP}{TP + FN}, \quad S = \frac{TP}{TP + FP} \tag{5}$$

$$F\_1 = \frac{2}{\frac{1}{5} + \frac{1}{7}}\tag{6}$$

The confusion matrix for this 4 Class classification problem is tabulated in Table 3. Table 4 summarizes the Sensitivity, Precision, and F1 scores, re-casting the multi-class problem into a 2-class problem (Normal versus a specific class) to calculate these metrics.

**Table 3.** Multi-Class Classification for Normal and Abnormal CT Scans.

Highest F1 scores are in bold.

#### **3. Results**

We performed experiments involving binary classification (Normal versus SAH, IPH, ASDH, BPH) and multi-class classification. In the latter, BPH was left out as BPH contains features of Normal, SAH, IPH, and ASDH. For the binary classification experiments, we implemented the experiments with and without the thresholding method described in Section 2.2.

For the multi-class classification experiments, the model discriminated between four classes with an overall F1 score of 0.684. Table 3 presents the confusion matrix for the multi-class classification. The actual number for each class represents the augmented test set after the original test set was augmented eight-fold. The respective F1 scores for each class are Normal: 0.819, SAH: 0.639, IPH: 0.427, ASDH: 0.829. Thresholding was not implemented for the multi-class classification as it is optimized for binary classification. The model performed well for Normal and ASDH scans, but only moderately well for SAH. Interestingly, although the training dataset for ASDH was the smallest, the model was able to discriminate this the best. This may be due to the fact that ASDH images are visually grossly asymmetrical compared to the other classes. This is due to brain compression from the significant subdural hemorrhage (see Figure 1), and may be a strongly activated discriminatory feature. Surprisingly, a significant number of IPH scans were misinterpreted as SAH, and we hypothesize that this may be due to subtle SAH traces that may appear on some of the IPH scans.

Table 4 summarizes the results from the 2-class classification experiments. Overall, our model was able to discriminate between normal and each of the classes with satisfactory results, matching or exceeding previously published results for similar work (Table 5). We demonstrate that for every class, the implementation of the thresholding technique improves all the evaluation metrics. The largest increase was seen in the ASDH class, where the F1 score increased from 0.706 to 0.952, despite a small training set. Various model architectures, input image, and filter sizes were modified to optimize for

accuracy and computational time. The original image size of the CT scans was 512 × 512 × 28 pixels (at 5 mm slice thickness), and the input size to our model was 50 × 50 × 28 pixels. Further decreasing the input size resulted in deteriorating model performance. We posit that the thresholding technique improves the signal-to-noise ratio for each input image by downplaying the extraneous image features, thereby accentuating the presence of acute hemorrhage. We also found that the thresholding technique decreased training time for the respective models significantly. For example, in the Normal versus SAH classification task, the model training time was decreased from over 10 h to 1 h 32 min. 3D CNNs are computationally intensive to train, and any decrease in computational cost and time is advantageous in clinical deployment.

**Table 5.** Comparison of results involving brain hemorrhage detection in volumetric brain scans.

Dou et al. expressed their evaluation metric as Sensitivity, Precision, and False positives per subject.

#### **4. Discussion**

Acute brain hemorrhage is a common neurosurgical emergency which can result in severe patient morbidity and mortality. It is the result of myriad causes, including trauma, hypertension, cerebral aneurysm rupture, and the treatment of each is different. Depending on the clinical condition, the patient may require close observation in a high dependency or intensive care setting, or immediate neurosurgical operation. It is imperative to minimize the time from diagnosis to treatment, to give the patient the best chance of recovery.

We propose an automated 3D CNN to classify volumetric CT brain data into various hemorrhagic variants, to assist doctors in expediting patient treatment.  These classes were chosen as the neurosurgical treatment for each class is different.

We also proposed and implemented a novel pixel thresholding method to detect acute hemorrhage on CT brain scans. This method improved classifier performance on our dataset and can be conceivably exported for use in other datasets for other anatomical regions, where acute hemorrhage detection is required. Potentially, aside from detecting acute blood, this method can also be generalized for the detection of other abnormalities such as tumors. Figure 3 demonstrates an image slice of a CT brain with acute subdural hemorrhage. The hemorrhage is the white crescent on the left of the image, which is putting pressure on the grey areas of brain, and pushing it to the right of the image. To visualize the activations in the 3D CNN better, we used the deconvolution technique described by Yosinski et al. [25] to visualize this single image. The top row of the image (boxes B and C) represent convolution layer 1 and pooling layer 1, respectively. The bottom row (boxes D and E) represent the same layers with thresholding applied. The difference with thresholding is that the edges of the target hemorrhagic lesion appear more distinct and sharper, which may account for why thresholding improves classification accuracy.

**Figure 3.** To demonstrate the effect of thresholding, a single slice of an ASDH CT Brain scan is shown, with the corresponding activations of convolution and pooling layers. **A**, original image. **B**, 1st convolution layer. **C**, 1st pooling layer. **D**, 1st convolution layer with thresholding applied. **E**, 1st pooling layer with thresholding applied. D and E appear sharper than B and C, demonstrating how thresholding can accentuate abnormal areas, and improve classifier performance.

In this paper, we demonstrate state-of-the-art 3D CNN classifier performance for different classes of acute brain hemorrhage. In addition, the application of our thresholding method for acute hemorrhage enables our model to achieve further additional improvement in classification. In the 4-class classification task, our model achieved an overall F1 score of 0.684. To our knowledge, there has been no published work involving multi-class classification of different classes of brain hemorrhage for comparison. In the 2-class classification task, the best performance was achieved in differentiating normal from ASDH scans, with a F1 score of 0.952 using our proposed thresholding method. The largest improvement with thresholding was also seen in the ASDH class, as the initial F1 score was only 0.706.

There has been a long history in the attempt to analyze hemorrhage in volumetric CT and MRI data. Before the widespread use of CNNs in image analysis, methods used included Hopfield networks [26], support vector machines [27], and segmentation masking with logistic classifiers [28]. Although these simple classifiers performed well, they were often applied to single image slices with obvious pathology, that were often manually chosen, which, therefore, represents an unrealistic problem scenario. Hybrid 2D CNN methods exemplified by Roth et al. [14] were a bridge to the 3D CNN training of networks. Fully 3D CNN model local and contextual spatial dependencies and extract features in all three dimensions of image voxels. Kamnitsas et al. [16] exploited the dense inference technique and small kernels, to segment lesional areas on brain MRI scans. Of note, they used a dual 11-layer 3D CNN pathway to process images at multiple scales, and Conditional Random Fields to decrease their false positive rate.

There are two points worth noting in our network architecture design. First, we deliberately kept our network architecture simple with relatively few layers, and did not leverage on other techniques, such as ensembling and transfer-learning. This was for both practical and theoretical considerations. The application of full 3D CNN to volumetric medical images is nascent, and by purposely keeping the model architecture straightforward, we are able to assess the effect size of various hyperparameter-tuning quickly, experiment with various architectures efficiently, and to troubleshoot network errors expediently. Running 3D CNN is computationally intensive, and a simpler network mitigates lengthy run-times.

From a more theoretical perspective, a straightforward architecture also allows us to grasp a sense of the baseline performance of a 3D CNN, and to describe our experimentations with a clear mathematical description. We believe that a firm theoretical framework will help in directing further areas for exploration, rather than blindly tuning hyperparameters, or implementing network boosting ensembling architectures without understanding.

The second point is that with a straightforward network architecture design, we were able to demonstrate the clear improvement that our thresholding method brought for detecting acute hemorrhage. We arrived at this thresholding method by observing radiologists as they scrolled through various patient CT scans. We noted that in situations where areas of acute brain hemorrhage were subtle, the radiologists would increase and decrease the image contrast to accentuate the appearance of the abnormality, which in this case was blood. We did the same thing while studying our dataset, and, therefore, explored how this optimization of human visual cognition or analysis could be implemented algorithmically. We took inspiration from the work of Zhang et al. [23] who used pixel intensity 'spatial histograms' for object detection within an image. Our underlying assumption was similar, that objects of a certain class would have similar patterns of pixel intensities. Where we differ is that while Zhang et al. were concerned with the spatial location of the object, we are concerned more with the actual presence of a particular pattern of pixel intensity, and the point of divergence from the pattern denoting a normal scan.

The proposed architecture in this work has important clinical significance. The different abnormal diagnoses studied in this work all exert a significant epidemiological and social-economic toll on healthcare systems and societies. SAH, IPH, ASDH, and BPH are all neurosurgical emergencies that require immediate but different treatments to maximize the likelihood of patient survival and to achieve a good long-term functional outcome. Our work has a role in helping clinicians minimize the time between diagnosis and treatment, especially in hospitals that may not have a radiologist after hours, or in remote rural settings where no clinicians are available. Even in large tertiary care hospitals with 24 h-radiologists, our proposed architecture framework can assist radiologists by triaging important abnormal scans from the large numbers of normal scans that are read sequentially as patients are scanned.

One limitation of this work is the relatively small and un-balanced dataset that we worked with, which is a common issue in medical image analysis, with a bias towards normal samples. At 399 images, our dataset size is comparable to those used in many other works [20,24] but smaller than the almost 40,000 CT scans used by Jnawali et al. [19]. Despite this, data augmentation techniques have resulted in performance comparable to Jnawali et al.'s much larger dataset. This may be due to the fact that medical data is relatively homogeneous in appearance, compared to natural image processing tasks. It would be interesting to study what is the optimal dataset size for processing specific volumetric image classification tasks, in CT, MRI, 3D ultrasound images for different lesions. We addressed concerns of overfitting with known mitigating techniques, such as dropout, which were applied to the early layers in our network.

There are many avenues for further exploration in volumetric medical image analysis. CNNs and 3D CNNs have been the dominant network architecture in image analysis, but unsupervised learning methods for medical image analysis are emerging for 3D object generation [29,30], and they have been largely unexplored in the context of 3D medical image analysis. Generative architectures, such as variational autoencoders and generative adversarial networks, have not been applied to volumetric medical data, and these techniques may potentially mitigate the need for large well-labelled datasets. Specific to our work, we intend to explore if adding a screening stage to a 3D CNN, or multi-scale receptive fields can improve performance, as some authors have demonstrated [16,20] for MRI brain scans. The addition of a memory or attention-based component to model long term dependencies in 3D medical image analysis [24,31,32] is also interesting for further investigation, as there is evidence for a strong biological correlate with the mammalian visual system [33].

#### **5. Conclusions**

This work presented the implementation of a 3D CNN to classify and diagnose volumetric CT brain data. Normal CT brains and a variety of abnormal scans constituting different types of brain hemorrhage were classified by our 3D CNN. We also implement a novel optimization method to detect acute hemorrhage on CT scans. The proposed 3D CNN can automatically detect important normal and abnormal features of cerebral anatomy without handcrafting, or significant data pre-processing. Computational costs were also modest, which will add to straightforward implementation. Results from classification experiments demonstrated that the 3D CNN outperforms previously published methods in detecting abnormal brain scans with hemorrhage, with higher sensitivity and recall. Our 3D CNN can be applied to other volumetric medical data and can be used to expedite and improve patient care.

**Author Contributions:** Conceptualization, L.W., J.R., T.L., J.K.; software S.S. and Y.B.; method J.K. and S.S.; validation, S.S., J.K. and Y.B.; writing J.K.; supervision L.W., J.R. and T.L.

**Funding:** This research was funded by the National Neuroscience Institute-Nanyang Technological University Neurotechnology Fellowship Grant.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### **Rapid Multi-Sensor Feature Fusion Based on Non-Stationary Kernel JADE for the Small-Amplitude Hunting Monitoring of High-Speed Trains**

#### **Jing Ning 1,2,\*, Mingkuan Fang 1, Wei Ran 1, Chunjun Chen 1 and Yanping Li 1**

Received: 29 March 2020; Accepted: 16 June 2020; Published: 18 June 2020

**Abstract:** Joint Approximate Diagonalization of Eigen-matrices (JADE) cannot deal with non-stationary data. Therefore, in this paper, a method called Non-stationary Kernel JADE (NKJADE) is proposed, which can extract non-stationary features and fuse multi-sensor features precisely and rapidly. In this method, the non-stationarity of the data is considered and the data from multi-sensor are used to fuse the features efficiently. The method is compared with EEMD-SVD-LTSA and EEMD-JADE using the bearing fault data of CWRU, and the validity of the method is verified. Considering that the vibration signals of high-speed trains are typically non-stationary, it is necessary to utilize a rapid feature fusion method to identify the evolutionary trends of hunting motions quickly before the phenomenon is fully manifested. In this paper, the proposed method is applied to identify the evolutionary trend of hunting motions quickly and accurately. Results verify that the accuracy of this method is much higher than that of the EEMD-JADE and EEMD-SVD-LTSA methods. This method can also be used to fuse multi-sensor features of non-stationary data rapidly.

**Keywords:** high-speed trains; hunting; non-stationary; feature fusion; multi-sensor fusion

#### **1. Introduction**

Hunting motion is a self-excited vibration that is a serious obstacle to the safety of high-speed trains [1]. Monitoring systems are designed to detect hunting only after it has developed to a specific degree. Besides, in most cases, the recognition result is only obtained using a single observation. Therefore, the accuracy and real-time performance of monitoring systems need to be further improved [2]. With the increasing performance requirements of high-speed trains, it is important to establish an accurate and rapid feature extraction method for hunting detection through multi-characterizations before the phenomenon has developed to any significant degree.

The structure of high-speed trains is very complicated and their working conditions are very poor, resulting in non-stationary vibration signals [3]. In the existing feature extraction research, a variety of data extraction algorithms are utilized. These algorithms can be roughly divided into two categories: those designed for stationary data and those for non-stationary data. Feature extraction methods for stationary data include Singular Value Decomposition (SVD), Linear Discriminant Analysis (LDA) [4,5], Principal Component Analysis (PCA) [6], Locality Preserving Projection (LPP) [7], and so on. However, it is hard to extract features of non-stationary data using these methods. For non-stationary data, manifold learning is a good feature extraction method [8–11]; however, it incurs a high computational cost and has extremely long calculation times, which means that the diagnostic information cannot be

fed back to the system in time. Therefore, a rapid yet precise method to extract the non-stationary signal in practical engineering applications is needed. Cardoso [12,13] proposed a method named Joint Approximate Diagonalization of Eigen-matrices (JADE) in the field of blind source separation, which is used to quickly separate multiple features. Because the method is simple and effective, it is also widely used in pattern recognition. Liu [14] improved the JADE method and achieved a good feature fusion performance by simplifying the calculations, and then developed a method based on kernel JADE to identify the rolling bearing fault [15]. The method based on JADE was adopted to diagnose the faulty parts of the rolling bearing in [16], and a similar method was applied to predict coaxial bearing performance degradation in [17]. However, in this application, the non-stationary and rapid character of the signal is not considered.

The structure of high-speed trains is complex. The vibration signals from high-speed trains are affected by many factors, such as the propagation path and the measuring point location, which lead the signal to couple with different characteristic information. Pattern recognition based on single sensor vibration signals has limitations, which cannot completely render the evolutionary characteristics of high-speed train vibrations [18]. Therefore, it is necessary to utilize a multi-sensor fusion method to extract the features of the train's operating state.

Furthermore, the evaluation parameters of lateral stability of railway passenger trains in different countries are different. Lateral force on the rail and on the wheel axis, lateral acceleration of the bogie frame, acceleration of the axle-box, and lateral acceleration of the vehicle body can be used as evaluation parameters, respectively [19–22]. However, only one parameter is used in each standard, and this parameter has different limits. For example, the influence of frequency is not considered when lateral acceleration of the bogie frame is used. In fact, the acceleration of the axle-box is a particularly important index in hunting monitoring. Therefore, in this paper, we try to fuse the signal from the bogie frame and axle box to extract the hunting motion features. To monitor the state of hunting of high-speed trains, four basic states are classified: normal, small convergence, small divergence, and hunting. Our aim is to identify the small divergence state before hunting occurs. Thus, in this paper, the monitoring of small amplitude hunting enables us to rapidly distinguish between these four basic states online [23].

Moreover, for real-time classification, it is extremely important to ensure that once the small amplitude appears to diverge, the entire calculation process can be completed immediately. Therefore, the time of the whole calculation process must be short enough for real-time classification.

Considering the above problems, a rapid multi-sensor feature fusion method based on the Non-stationary Kernel JADE method is proposed. The JADE method is a fast and accurate feature fusion algorithm, but it is generally used in stable environments. In order to use the algorithm in a non-stationary environment, the whole time series is divided into *M* time periods [24] and the kernel function is introduced. Then, *M* kernel matrices are obtained by the *M* time periods, which are jointly decomposed. After this, Jacobian rotation is used to obtain the unitary matrix by diagonalizing multiple kernel matrices simultaneously to extract the non-stationary fusion features. In addition, in order to visualize the data features, the extracted fusion features are expressed in three dimensions. The between-class indicator and within-class indicators [25] are also employed to describe the clustering performance of the features quantitatively.

In this paper, a multi-sensor data feature extraction framework is provided, in which a rapid feature fusion method using Non-stationary Kernel JADE (NKJADE) is proposed. This framework consists of the following series of steps. First, the Ensemble Empirical Mode Decomposition (EEMD) method is utilized to decompose the preprocessed signals to Intrinsic Mode Functions (IMFs). Then, the energy matrices are obtained using the IMFs and the fusion features are obtained through NKJADE, followed by inputting the extracted features to an LSSVM [26,27] for training and recognition.

Data from Case Western Reserve University (CWRU) are used to verify the performance of this method against that of the SVD and JADE methods. In this paper, the proposed method is also applied to identify the evolutionary trend of the hunting motion quickly and accurately. The results verify that the accuracy of this method is much higher than that of the JADE method.

The remainder of this paper is organized as follows: In Section 2, the theoretical backgrounds of Ensemble Empirical Mode Decomposition (EEMD) [28] and a separability metric are introduced. The proposed method of non-stationary Kernel JADE is also described. The framework of the proposed method is introduced in Section 3. In Section 4, data from Case Western Reserve University (CWRU) are used to test the proposed method, and the operational data of multiple states of high-speed trains are used to verify the accuracy and rapidity of the proposed method. Finally, the conclusion is presented in Section 5.

#### **2. Theoretical Background**

#### *2.1. EEMD*

The basic idea of the EEMD method is to use the sifting process to decompose the signal into several intrinsic mode functions, when Gaussian white noise is added. For the original signal, the specific steps to use EEMD to decompose it into IMFs are as follows:

(1) Obtain the overall signal by adding Gaussian white noise to the original signal:

$$a'(t) = a(t) + a(t) \tag{1}$$

(2) The overall signal is decomposed to obtain the IMF components of each order, where *i* represents the *i*-th component, *r* is the residual term, and *n* is the number of IMFs:

$$a'(t) = \sum\_{i}^{n} c\_i + r \tag{2}$$

(3) Repeat step (1) and (2), each time adding different white noise sequences of the same amplitude:

$$a\_j'(t) = \sum\_{i=1}^n c\_{ij} + r\_j \tag{3}$$

In Equation (3), *cij* is the *i*-th IMF component of the decomposition to which white noise at for the *j*-th time, while *rj* is the residual value of the decomposition.

(4) Using the zero-mean principle of the Gaussian white noise frequency, the effect of white noise can be eliminated, and the IMF component corresponding to the original signal can be expressed as:

$$c\_i(t) = \frac{1}{N} \sum\_{j=1}^{N} c\_{ij}(t) \tag{4}$$

where *n* represents the number of times the white noise is added and *ci* represents the *i*-th IMF component obtained by the EEMD decomposition of the original signal.

*2.2. IMF Energy Matrix*

For the IMF component *c i* = [*c i* (1), ...... , *c i* (*M*)], the energy feature can be expressed as:

$$c\_i = \sum\_{j=1}^{M} \left[ c'\_i(j) \right]^2 \tag{5}$$

All IMF components of a vibration signal form a feature energy vector *e* = [*e*1,*e*2, ......*en*], where *n* is the dimension of vector e and represents the number of IMF components.

#### *2.3. Proposed Non-Stationary Kernel JADE*

Since the original JADE algorithm is based on stationary signal analysis, considering the non-stationary nature of high-speed train signal, the kernel JADE method [15] is applied to a non-stationary environment.

The main idea of the kernel is to map the input matrix into the nonlinear space Φ Suppose *X* = {*x*1, *x*2, ... , *xm*}; then, the processing can be defined as follows:

$$\{\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_m\} \to \{\Phi(\mathbf{x}\_1), \Phi(\mathbf{x}\_2), \dots, \Phi(\mathbf{x}\_m)\} \tag{6}$$

During implementation, we need to calculate the inner product of two eigenvectors which have been mapped into a nonlinear space using a kernel function, and a kernel matrix will be calculated using Equation (7):

$$K\_{\dot{i}\dot{j}} = k(\mathbf{x}\_i, \mathbf{x}\_{\dot{j}}) = \prec \Phi(\mathbf{x}\_i), \Phi(\mathbf{x}\_{\dot{j}}) > \tag{7}$$

where *xi*, *xj* are vectors. The commonly used kernel function includes the following [29]:

$$k(\mathbf{x}\_i, \mathbf{x}\_j) = \exp(-\frac{\left\|\mathbf{x}\_i - \mathbf{x}\_j\right\|^2}{2\delta^2})\tag{8}$$

$$k(\mathbf{x}\_{i\prime}, \mathbf{x}\_{j}) = \left(\alpha \mathbf{x}\_{i}^{T} \mathbf{x}\_{j} + \mathbf{c}\right)^{d} \tag{9}$$

$$k(\mathbf{x}\_i, \mathbf{x}\_j) = \tanh(a\mathbf{x}\_i^T \mathbf{x}\_j + \mathbf{c}) \tag{10}$$

Similar to KPCA [30], the centered kernel matrix *K* can be calculated through Equation (11):

$$K\_{ij} = K\_{ij} - \frac{1}{M} \sum\_{r=1}^{M} K\_{ir} - \frac{1}{M} \sum\_{r=1}^{M} K\_{rj} + \frac{1}{M^2} \sum\_{r,s=1}^{M} K\_{rs} \tag{11}$$

The whole time series is divided into *M* segment intervals *T***1**, ...... , *TM*, and consequently *M* covariance matrices **S***T*1 , ...... , **S***TM* can be generated. These *M* covariance matrices are jointly diagonalized to find a unitary matrix *U* which can diagonalize *M* covariance matrices simultaneously; then, the energy features can be extracted.

**Step 1:** For *M* segment intervals *T*1, ...... , *TM*, the covariance matrices of signal *x*(*t*) can be expressed as:

$$\text{Cov}\_{T\_k} = \mathbf{S}\_{T\_k} = \frac{1}{|T\_m|} \sum\_{t \in T\_m} \left[ k(\mathbf{s}\_{t\prime} \mathbf{s}\_t) \right] \tag{12}$$

where *st* = *xt* − *E*(*xt*), *E*(*xt*) denotes the mean of *xt*.

**Step 2:** The most common method to diagonalize matrices **S***T*1 , ...... , **S***TM* is to diagonalize the first matrix, and then transform the remaining *M*−1 matrices into diagonalization. **W** is the diagonalized matrix of the covariance matrix **S***T*1 :

$$\mathbf{W} = \mathbf{S}\_{T\_1}^{-1/2} = \mathbf{V}^H \boldsymbol{\Lambda}^{-1/2} \mathbf{V} \tag{13}$$

In Equation (13), **V** is the eigen-matrix of **S***T*1 , while **Λ** is the eigenvector of **S***T*1 . For the remaining *M*−1 matrices, the diagonalization matrix can be defined respectively:

$$\mathbf{S}\_{T\_m}^\* = \mathbf{S}\_{T\_1}^{-1/2} \mathbf{S}\_{T\_m} (\mathbf{S}\_{T\_1}^{-1/2})^H, m = 2, \dots, M \tag{14}$$

**Step 3:** The approximate joint diagonalization problem is equivalent to finding an orthogonal matrix **U** that minimizes:

$$\sum\_{m=2}^{M} \left\| \text{off}(\mathbf{U} \mathbf{S}\_{T\_m}^\* \mathbf{U}^H) \right\|^2 = \sum\_{m=2}^{M} \sum\_{b \neq d} \left( \mathbf{U} \mathbf{S}\_{T\_m}^\* \mathbf{U}^H \right)\_{bd}^2 \tag{15}$$

where off(**US**∗ *Tm* **U***H*) has the same off-diagonal elements as **US**∗ *Tm* **U***H*, and the diagonal element is zero, while *b* and *d* represent the *b*-th row and *d*-th column of the matrix, respectively.

Since the sum of squares remains the same when multiplied by an orthogonal matrix, the problem is equivalent to maximizing the sum of squares of the diagonal elements:

$$\sum\_{m=2}^{M} \| \text{diag} (\mathbf{U} \mathbf{S}\_{T\_m}^\* \mathbf{U}^H) \|\ \ \ = \sum\_{m=2}^{M} \sum\_{b=1}^{p} (\mathbf{U} \mathbf{S}\_{T\_m}^\* \mathbf{U}^H)\_{bb}^2 \tag{16}$$

where *p* represents the dimension to which the feature is extracted.

**Step 4:** Givens rotation is used to transform the set of matrices to a more diagonal form, two rows and two columns at a time. The Givens rotation matrix is given by:

$$G(i,j,\theta) = \begin{pmatrix} 1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & & \vdots & & \vdots \\ 0 & \cdots & \cos(\theta) & \cdots & -\sin(\theta) & \cdots & 0 \\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ 0 & \cdots & \sin(\theta) & \cdots & \cos(\theta) & \cdots & 0 \\ \vdots & & \vdots & & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1 \end{pmatrix} \tag{17}$$

where

$$\theta = \frac{1}{2} \text{arccot} \left[ \frac{(\mathbf{S}\_{T\_m}^\*)\_{22} - (\mathbf{S}\_{T\_m}^\*)\_{11}}{2(\mathbf{S}\_{T\_m}^\*)\_{12}} \right] \tag{18}$$

The initial value for the orthogonal matrix *U* is the identity matrix *I*. Matrices **ST**2 **,** ... **,STm** are then updated using:

$$\begin{array}{c} \mathbf{S}\_{T\_m}^\* \leftarrow \mathbf{G}(1,2,\theta) \mathbf{S}\_{T\_m}^\* \mathbf{G}(1,2,\theta), m = 2, \dots, M\\ \mathbf{U} \leftarrow \mathbf{U} \mathbf{G}(1,2,\theta) \end{array} \tag{19}$$

When the values of all non-diagonal elements are less than a given threshold ε, the iteration is completed, and the joint approximation diagonalization is achieved. The unitary matrix **U**ˆ is obtained so that multiple matrices are diagonalized jointly.

**Step 5**: The transform matrix *A* can be calculated as **A**=**UW**ˆ **#** , where superscript # denotes the pseudo-inverse.

As such, the original signal *s*(*t*) can be expressed as:

$$\mathbf{s}(t) = \mathbf{A} \cdot \mathbf{K} \tag{20}$$

Since the NKJADE method can extract the nonlinear relationships hidden in the high-dimensional feature space, the fusion features can be estimated through the joint feature decomposition using multiple inputs. The fusion features can express the nonlinear and non-stationary relationships hidden in the inputs well, so they can represent the characteristic relationship of different states, and then distinguish different feature states quickly and accurately.

#### *2.4. Separability Evaluation*

To illustrate the merit of our proposed algorithm, the separability *J* is utilized to demonstrate the algorithm's ability to form distinct classes. The capability of the feature extraction in pattern classification can be described quantitatively using between-class scatter *Sb*, within-class scatter *Sw* [31], and separability [32]. Assuming that the data have a total of C classifications, the vector of the *i*-th classification is:

$$\mathbf{x}^{i} = \left(\mathbf{x}\_{1'}^{i}, \mathbf{x}\_{2'}^{i}, \dots, \mathbf{x}\_{n\_i}^{i}\right) \tag{21}$$

where *ni* is the number of *i*-th classification.

Between-class scatter *Sb* and within-class scatter *Sw* can be calculated as follows:

$$\mathcal{S}\_{\mathbf{b}} = \sum\_{i=1}^{\mathcal{C}} p\_i (m\_i - m) \left( m\_i - m \right)^T \tag{22}$$

$$\mathbf{S}\_{\mathbf{i}\sigma} = \sum\_{i=1}^{\mathbb{C}} \left[ \frac{p\_i}{n\_i} \sum\_{k=1}^{n\_i} \left( \mathbf{x}\_k^i - m\_i \right) \left( \mathbf{x}\_k^j - m\_i \right)^T \right] \tag{23}$$

where *pi* = *ni*/ *C j*=1 *nj*, *mi* = *mean*(*xi ni* ), *m* = *C i*=1 *pimi*.

The between-class scatter *Sb* describes how far different classes are separated, and the within-class scatter *Sw* indicates how compactly each class of samples is distributed. Based on between-class scatter and within-class scatter, separable evaluation *J* is introduced to describe the clustering ability of different methods. Separable evaluation *J* could be calculated as follows:

$$J = \text{trace}(\mathbb{S}\_b/\mathbb{S}\_w) \tag{24}$$

where function *trace* refers to the sum of diagonal elements.

#### **3. Methodology**

For the different characteristics of vibration signals, a multi-sensor data feature extraction framework is provided, in which a rapid feature fusion method using NKJADE is proposed. The framework of the feature extraction method is shown in Figure 1. The method can quickly and accurately extract different features of information contained in non-stationary vibration signals. The main steps of this framework are as follows:

**Figure 1.** Method framework.

#### **4. Experiment Results and Analysis**

In order to verify the valid of the method, we applied it to bearing fault identification (Case I) and small amplitude hunting monitoring of high-speed trains (Case II). Some traditional approaches, such as Ensemble Empirical Mode Decomposition Joint Approximate Diagonalization of Eigen-matrices (EEMD-JADE) and Ensemble Empirical Mode Decomposition Singular Valuable Decomposition Learning Technology Systems Architecture (EEMD-SVD-LTSA), were compared with the proposed method.

#### *4.1. Case I—Case Western Reserve University Data*

#### 4.1.1. Data Description

The bearing test data for normal and faulty bearings were from the Case Western Reserve University (CWRU). The test signal contained normal state, ball fault, and inner and outer race faults; for the latter three, the fault diameters were 0.07, 0.14 and 0.21 inches, respectively.

The data used in this paper were collected from the drive end. The sampling frequency was 12 KHZ and the speed of the shaft was 1725 r/min, corresponding to 400 points collected per revolution. In order to reduce the influence of equipment fluctuations, each sample contained 800 points.

As shown in Table 1, two datasets were selected. Dataset A had the same fault location (inner race fault), but the fault size was different. Dataset B had different fault locations, but the fault size was the same, and the outer race fault was at the location of 3 o'clock. Each of the datasets was divided into three categories.

#### 4.1.2. Signal Processing Results

The scatter plots obtained by applying the EEMD-SVD-LTSA, EEMD-JADE method and EEMD-NKJADE methods on dataset A are shown in Figure 2a–c, while the corresponding scatter plots for dataset B are shown in Figure 2d–f. The results of selecting parameters for kernel are shown in Figure 3. The results of the application of the three algorithms in dataset A are shown in Table 2, while the results for dataset B are shown in Table 3.

**Figure 2.** Scatter plots of faults. (**a**) EEMD-SVD-LTSA on dataset A; (**b**) EEMD-JADE on dataset A; (**c**) EEMD-NKJADE on dataset A; (**d**) EEMD-SVD-LTSA on dataset B; (**e**) EEMD-JADE on dataset B; (**f**) EEMD-NKJADE on dataset B.

**Figure 3.** The results of selecting parameters for the Gaussian kernel. (**a**) Parameter selection results of Gaussian kernel of EEMD-NKJADE method on dataset A; (**b**) Parameter selection results of Gaussian kernel of EEMD-NKJADE method on dataset B.

**Table 2.** Results from dataset A.

#### **Table 3.** Results from dataset B.

#### 4.1.3. Discussion

We found that the generalization ability of the Gaussian kernel is better than that of the other kernels. Therefore, we focused on the parameters of the Gaussian kernel.

According to the results of selecting parameters in [15], we tried to set the range of σ to (1, 10). The parameter was incremented step by step (parameter step is 0.5, shown in X-axis), and the separability *J* was calculated (shown in Y-axis). The larger *J* is, the better the classification result will be. Therefore, for the parameter selection of EEMD method, the optimal parameter of the Gaussian kernel on dataset A and dataset B is 2.

From Figure 2 and Table 2, we see that all three methods could extract the features with satisfactory results. The accuracy of the three methods was nearly 100% in all cases. The result may be attributed to the CWRU bearing fault data having great differences, which are quite easy to classify.

However, the separability *J* in Table 2 is different, which is consistent with (a)–(c) in Figure 2. In this respect, the classification ability of EEMD- NKJADE algorithm is obviously better than other algorithms. In addition, the results of the running time in Table 2 show that the time required for EEMD-JADE calculation was relatively short (the PC configuration was as follows: Intel Core i5-4460, 12GB of memory, NVIDIA GeForce GT720).

In order to further verify the robustness of the algorithm, we also tested the results of the algorithm on dataset B.

In Table 3, the accuracy of the three methods in dataset B was almost the same as that in dataset A. However, the separability of dataset B was the highest after being processed using the EEMD-SVD-LTSA algorithm, while in dataset A, the separability achieved using EEMD-SVD-LTSA was the lowest, which indicates that the selection of dataset had a greater impact on EEMD-SVD-LTSA. Compared with EEMD-SVD-LTSA, the results from the EEMD-NKJADE and EEMD-JADE algorithms were less affected by the different dataset, and therefore seem to be more robust. Therefore, the method of EEMD-NKJADE offers superior performance with respect to the classification effect and the robustness, but its calculation time is longer than EEMD-JADE.

#### *4.2. Case II—Small Amplitude Hunting Monitoring of High-Speed Trains*

#### 4.2.1. Problem Description

The stability of hunting has always been a key problem in the study of vehicle lateral motion stability [33]. Small amplitude hunting is a sign of hunting instability. In China, hunting phenomena are considered to occur when the amplitude of lateral acceleration signals from the bogie frame reaches or exceeds 8–10 m/s2 more than 5 times continuously after a 10 Hz low-pass filter [34]. In Figure 4, the lateral acceleration of the bogie frame signals will sometimes go through a normal operation/small amplitude hunting/normal operation periodic cycle, which is a gradually convergent

process. The hunting amplitude in this situation is small and convergent. Sometimes, the signals will go through a normal operation/small amplitude hunting/critical hunting process, which is a gradually divergent process. The hunting amplitude in this situation starts small and then diverges.

**Figure 4.** Lateral acceleration signal of the bogie frame when small-amplitude hunting motion occurs during an online test. (**a**) Bogie frame acceleration; (**b**) Small amplitude convergent hunting; (**c**) Small amplitude divergent hunting.

Therefore, it is necessary to extract the features of different states rapidly and accurately, especially the small amplitude divergent hunting states, to guarantee the system is alerted in time to ensure the safe operation of the train.

#### 4.2.2. Data Acquisition

The data used in this paper were lateral acceleration signals of the bogie frame and axle box from an online tracking experiment. The CRTS II ballastless track and seamless rail were used on the whole line. The speed of the train was 320–350 km/h. The sampling frequency was 2500 Hz. All the data were acquired in accordance with China's Railway Passenger Traffic Safety Monitoring Standard [35]. Although in China the amplitude of lateral acceleration signals from the bogie frame is used as the testing parameter to monitor hunting motion, research has proven that many other testing parameters are also important for hunting monitoring. As such, in this paper, acceleration signals of the bogie frame and the axle box were used.

The installation locations of the accelerometers are shown in Figure 5, where two accelerometers are installed in the diagonal direction on the H-shaped bogie frame. The lateral accelerometers from the bogie are respectively denoted as S1 and S2. Also, considering that the vibration of the axle box is important for hunting [36], a sensor located on the axle box was used. In Figure 6, the lateral accelerometer on the axle box is denoted as S3. Figure 5 shows a photograph of the installation site.

**Figure 5.** Installation locations of the accelerometers. S1, S2: Accelerometer on the bogie frame; S3: Accelerometer on the axle box.

**Figure 6.** Site installation drawing.

Considering the applications of other researchers, the sampling frequency of the original signal was set to 2500 Hz. However, for hunting, the characteristic frequency of the lateral acceleration of the frame is only 3–7 Hz. Therefore, in the preprocessing stage, a 250 Hz resampling method was applied [37], and a low-pass filter of 0–10 Hz was applied. Then, a zero-mean smoothing process was used for preprocessing to eliminate trend terms. In accordance with commonly followed practices, parameters such as the amplitude of the lateral acceleration signals from the bogie frame and axel box were used to obtain a synthetic assessment of the lateral stability of the high-speed train tested. In this paper, the filtered lateral acceleration signals were divided into four states: normal, small amplitude convergent hunting, small amplitude divergent hunting, and hunting. Ten groups of sample data were used for each of the four states and each sensor, yielding a total of 120 sample groups. The length of each sample was 500 points, corresponding to a sample time of 2 s.

Nonlinear factors [38,39] have been proven to affect the bifurcation evolution of small amplitude hunting and, according to [23], all the values of the Lyapunov exponent of the lateral acceleration are greater than 1, which means that the lateral acceleration signals from the bogie frame have non-stationary characteristics.

#### 4.2.3. Feature Fusion

First, EEMD was applied on each signal. Then, 7 IMFs and 1 residue were obtained using the EEMD method. Figure 7 shows an EEMD decomposition view of the three signals at different points during the same time period.

**Figure 7.** Ensemble Empirical Mode Decomposition (EEMD) diagram of signals at different positions. (**a**) S1 position (bogie frame); (**b**) S2 position (bogie frame); (**c**) S3 position (axle box).

The samples of the three sensors were processed, and the three corresponding energy matrices E1, E2, and E3, with a size of 40 × 8 were obtained. The NKJADE method was used to fuse the three high-dimensional energy matrices and the data dimensionality was reduced.

#### *4.3. Result and Discussion*

#### 4.3.1. Single Sensor Classification Using NKJADE

A scatter plot of the features extracted using a single sensor and NKJADE is shown in Figure 8. The separability *J* and accuracy of the NKJADE method are shown in Table 4.

**Figure 8.** Scatter plot of feature extraction from single sensor using Non-stationary Kernel JADE (NKJADE). (**a**) S1 (bogie frame); (**b**) S2 (bogie frame); (**c**) S3 (axle box).

**Table 4.** Separability *J* running time and accuracy of NKJADE method.

From Table 4, the accuracy achieved using S1 data only was 97.23%, which was greater than that attained at the other sensor locations. However, when the three sensors were used together, the accuracy rate became 100%, and *J* reached a much greater value.

#### 4.3.2. Multi-Sensor Fusion Using NKJADE

The EEMD-SVD-LTSA and EEMD-JADE methods were used to compare the identification accuracy and calculation speed of the feature fusion method in multi-sensor conditions. The results of the methods are shown in Figure 9. The separability *J* and accuracy obtained using the different sensor fusion methods with data from all three sensors are shown in Table 5.

The JADE and SVD-LTSA methods were compared with the proposed method in the multi-sensor feature fusion conditions. From Table 5, the accuracy rate of the SVD-LTSA method was 93.75%, which used the non-stationary method LTSA. The accuracy rate of the JADE method was only 30.12%, which used the stationary method. The accuracy rate of NKJADE (Gaussian kernel, δ = 0.6) was 100%, while the separability *J* of NKJADE was greater than those of the other methods. Considering that data from high-speed trains have more typical no-stationary characteristics than the bearing fault data from CWRU, this result shows that the proposed method may more be suitable for non-stationary data. The separability *J* of the JADE method was only 26.67, which is even lower than that achieved using single sensor S2 or S3 data. The reason for this is probably that in the JADE method, the non-stationary condition is not considered.

Besides, because real-time processing is a significant factor for the small amplitude hunting monitoring of high-speed trains, the calculation time is a very important factor. If the calculation time is too long, the diagnostic information cannot be fed back to the system in time. Compared with the SVD-LTSA method, fusion features can be extracted very quickly using the proposed NKJADE method. The run time of the NKJADE method was nearly the same as that of the JADE method, and the separability of *J* of NKJADE was greater than JADE. This shows that NKJADE is a rapid multi-sensor feature fusion method based on non-stationary condition, which outperforms the SVD-LTSA and JADE methods. It can be applied to the small amplitude hunting bifurcation evolution monitoring in high-speed trains.

**Figure 9.** Scatter plot of feature extraction from multi-sensor data. (**a**) EEMD-SVD-LTSA; (**b**) EEMD-JADE; (**c**) EEMD-NKJADE.

**Table 5.** The separability *J*, running time, and accuracy using different sensor fusion methods.

As shown in Figure 10, the parameter is incremented step by step (parameter of step is 0.01, shown in X-axis), and the separability index *J* is calculated (shown in Y-axis). The larger *J* is, the better the classification result will be.

**Figure 10.** Parameter selection results of Gaussian kernel of EEMD-NKJADE method.

Because the raw data collected are mainly distributed in the range of −1–1 g (gravity). In Equation (8), σ is the width parameter of the kernel function, so we tried to set the range of σ to (0, 1). As shown in Figure 10, with the increase of δ*, J* increases at first and then decreases. It indicates that the range of parameter selection is reasonable. When δ = 0.6, the classification result is the best.

#### **5. Conclusions**

In this paper, a rapid multi-sensor feature fusion method based on NKJADE is proposed, with which the features of multiple sensors can be extracted quickly and accurately. In order to use the algorithm in a non-stationary environment, the whole time series is divided into *M* time periods, and the kernel function is introduced. Then, Jacobian rotation is used to obtain the unitary matrix by diagonalizing multiple kernel matrices simultaneously to extract the non-stationary fusion features.

Our main findings are that:

**Author Contributions:** Data analysis and funding acquisition, J.N.; Writing—original draft preparation, M.F.; data collection: W.R.; experiment design, C.C.; detection equipment, Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundations of China, grant number 51975486 and 51975487.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A Cognitive Method for Automatically Retrieving Complex Information on a Large Scale**

#### **Yongyue Wang 1, Beitong Yao 1, Tianbo Wang 2, Chunhe Xia 1,3 and Xianghui Zhao 4,\***

Received: 19 April 2020; Accepted: 26 May 2020; Published: 28 May 2020

**Abstract:** Modern retrieval systems tend to deteriorate because of their large output of useless and even misleading information, especially for complex search requests on a large scale. Complex information retrieval (IR) tasks requiring multi-hop reasoning need to fuse multiple scattered text across two or more documents. However, there are two challenges for multi-hop retrieval. To be specific, the first challenge is that since some important supporting facts have little lexical or semantic relationship with the retrieval query, the retriever often omits them; the second challenge is that once a retriever chooses misinformation related to the query as the entities of cognitive graphs, the retriever will fail. In this study, in order to improve the performance of retrievers in complex tasks, an intelligent sensor technique was proposed based on a sub-scope with cognitive reasoning (2SCR-IR), a novel method of retrieving reasoning paths over the cognitive graph to provide users with verified multi-hop reasoning chains. Inspired by the users' process of step-by-step searching online, 2SCR-IR includes a dynamic fusion layer that starts from the entities mentioned in the given query, explores the cognitive graph dynamically built from the query and contexts, gradually finds relevant supporting entities mentioned in the given documents, and verifies the rationality of the retrieval facts. Our experimental results show that 2SCR-IR achieves competitive results on the HotpotQA full wiki and distractor settings, and outperforms the previous state-of-the-art methods by a more than two points absolute gain on the full wiki setting.

**Keywords:** information retriever sensor; multi-hop reasoning; evidence chains; complex search request

#### **1. Introduction**

Most previous work about retrievers focus on searching the relative lexicon from a single paragraph, and some easy search queries may be satisfied reluctantly using a single paragraph or document by the traditional retrieval methods of search engines, such as matching lexical or semantic similarity or relatedness. However, a very challenging retrieval task for complex search requests, called multi-hop reasoning retrieval, requires combining evidence from multiple sources, as shown in Figure 1. Single-hop retrieval methods are almost incompetent for intricate search tasks. The objective of a multi-hop retriever is to obtain and fuse scattered information from the internet step by step [1,2]. Figure 1 provides an example of a search request requiring multi-hop reasoning. In response to the search request, only when one first infers from the first context in the search request can the Internet information corresponding to the query be inferred with any confidence from the second context.

**Figure 1.** An example of a multi-hop retrieval task. The method in this paper mimics the human reasoning process while searching for information online. The retrieval sensor will finally provide users with reasoning chains, which will help obtain structured retrieval information.

The multi-hop IR is of great difficulty, because it requires reading and aggregating information over multiple pieces of textual evidence. Although the Reading Comprehension model is applied to candidates of the paragraph for prediction (e.g., QFE [3], MINIMAL [4]), the aggregation of evidence from sources for supporting fact prediction is beset with critical difficulty. One of the challenges is that since some pivotal supporting facts have little lexical or semantic relationship with the original retrieval query, the retriever often omits these crucial facts, as illustrated in Figure 2. The reason for this challenge is that multi-hop IR is open domain retrieval on a large scale, where search requests are given without any accompanying contexts. In this case, one is required to locate relevant contexts to queries from a large internet source before extracting the most relevant information. Retrieving a fixed list of documents independently does not capture the relationships between evidence documents through bridge entities required for multi-hop reasoning. They often fail to retrieve the required evidence for the multi-hop query. The recent open domain IR methods [5] do not capture lexical information in entities, but are encountered with challenges for entity-centric problems. The complex cognitive process in humans seems to provide more inspiration for this retrieval process. Recently, some researchers began to mimic the human reasoning process to search for complex information. Although Ding's Cognitive Graph [6] takes the entity relations, this method cannot find the reasoning path automatically. Multi-Step Reasoner introduces a multi-step reasoning method [7], but it only repeats the retrieval process for a fixed number, resulting in ineffectively integrated information. Qi et al. propose GoldEn IR [8], which cannot use the information of relations between documents during the iterative process. Most state-of-the-art approaches [2,9–11] for the open-domain retrieval leverage non-parameterized models to retrieve a fixed set of documents, where an information span is extracted by a neural model. However, they often fail to retrieve the required evidence needed by multi-hop reasoning. Some research on the reasoning process do not separate the information from the cognitive graph. If helpful information for users is not selected on the cognitive graph, the retriever may fail. In addition, the explosion of online textual content with unknown verification raises a vital

concern about misinformation, such as fake news, rumors, and online "water army" opinions [12–14]. While some misinformation may have semantic relations with the retrieval query, they cannot be efficaciously identified by some methods [15–17] with ideal performance. If the scattered information does not have stronger interactions with neighbors on the cognitive graph, a retriever is apt to select this misinformation in comparison with the whole nodes. While RE3 and Multi-Passage [11,18] introduce an extract-then-select framework to re-rank candidates to improve these interactions, they overlook the information between different evidence candidates.

**Figure 2.** A complex example of an open-domain retrieval task from HotpotQA. Paragraph 2 is difficult to be retrieved using traditional retrievers due to little lexical overlap with the given query.

In this paper, we study the problem of multi-hop retrieval methods, which requires multi-hop reasoning among evidence scattered around multiple raw documents on a large scale. Despite the given query utterance and a set of accompanying documents, not all of them are relevant. The information can be obtained by selecting two or more pieces of evidence from documents and reasoning among them. The models are expected to search the most useful information for search request in open domains. We propose a sub-scope cognitive reasoning information retrieval (2SCR-IR) sensor, a novel method to address the above-mentioned difficulties for multi-hop retrieving, as shown in Figure 3 (Nodes are fundamental entities, with the color denoting the identical entity in paragraphs. The blue edges are the relationship between the different entities in the same paragraph, while the red ones are the relationship between the same entities in different paragraphs. One search request and three paragraphs are given. Our 2SCR-IR sensor conducts multi-step reasoning over the supporting evidence by constructing a cognitive graph from the supporting information, automatically adjusting the sub-scope to select a subgraph, propagating information through the graph, verifying the reliability of the information, and finally producing supporting evidence chains. The circles, denoting the sub-scope, are chosen by 2SCR-IR in each step). For the first challenge, the 2SCR-IR sensor constructs a self-adjusting sub-scope cognitive graph based on entities referred to in the search request and web documents, which is iterated to accomplish multi-hop reasoning. In each step, the 2SCR-IR sensor is produced and reasons on a dynamically adjusted graph, with unrelated entities left out and reasoning information exclusively preserved in a scope prediction process. To solve the second difficult problem, we use an information fusion process in the 2SCR-IR sensor to eliminate misinformation. The process of iteratively expanding with clues in sub-scope can discover the weak paragraphs relative to the query in our framework, which can also play a pivotal role in filtering out misinformation. For the purpose of further verifying the authenticity of retrieved information, we introduce pageview verification to our sensor. Overall, our work incorporates a four-fold contribution:

#### **2. Related Works**

#### *2.1. Search Engine and Information Retrieval*

The search engine is a system that assists users in retrieving information that they wish to obtain after submitting searching queries. One search engine, in response to a retrieval query, typically compares keywords of the query with an index generated from a sea of web sources, such as text files, image files, video files, and other content items. Based on this comparison, it will provide users with the most relevant content items [19]. In our paper, only textual information retrieved from the Internet is taken into account. With the rapid development of deep learning, researchers have an increasing interest in bringing neural networks into IR tasks. Although IR tasks frequently use sentence-similarity, as discussed by Guo et al. [20], they pay more attention to relevance-matching, in which the match of specific content plays a significant part. Some researchers [21] have confirmed that IR is more about retrieving sentences with the identical semantic meaning, while traditional methods based on relevance matching show a serious deficiency in semantic understanding. Therefore, we used natural language reasoning techniques, instead of relevance-focused IR methods, to solve this issue. In many cases, the paragraph containing the information corresponding to the searching query has great lexical overlap with the query, adding to its difficulty for the retrieval of common search engines from a large open scope. For instance, the accuracy of a BM25 retriever for finding all supporting evidence for a query diminishes from 57.3% to 25.9% on the "easy" and "hard" subsets, respectively [2].

(**a**)

**Figure 3.** *Cont*.

**Figure 3.** An example of a multi-hop text-based information retrieval task. (**a**). Search request and corresponding information source. (**b**). Reasoning process.

#### *2.2. Single-Hop and Multi-Hop QA*

Relying on the complexity in underlying reasoning, information retrieval tasks can be categorized into single-hop and multi-hop ones. The former only needs one piece of evidence extracted from the underlying information, e.g., "which city is the capital of China". On the contrary, a multi-hop retrieval task requires recognizing multiple relevant evidence and reasoning about them, e.g., "what is the capital city of the largest province in China". Many IR techniques that are able to reply the single-hop searching query [22] are hardly introduced in multi-hop tasks, since single evidence can only partially satisfy its searching query.

#### *2.3. Open-Domain Retrieval Task*

Open-Domain retrieval task refers to the setting where the search scope of the supporting evidence is extremely large and complex. Recent open domain IR systems based on deep learning methods follow a two-step process, namely selecting potentially relevant documents from a large corpus and extracting useful sources from these selected documents. As the first to successfully apply this framework to Chen et al. [23], IR [24], introduced new benchmarks, [11,25] improved this framework by introducing feedback signals, and [26] proposed a method that can dynamically adjust the number of retrieved documents. However, these mainstream methods are weak in mining the information loosely related to the searching queries and fail to identify the misinformation.

#### *2.4. Similarity Matching Based Method*

Traditional frequency-based approaches, such as the "bag of words" model [27], have been employed to extract features in the search engine. However, these frequency-based techniques do not preserve the text sequence, bringing about the lack of understanding of the context's full meaning [28]. While other text-similarity-matching approaches, such as tree kernels [29] and high order n-grams [30], may have the perception of the word order and semantics, they cannot master the context meaning completely, thus heavily affecting the accuracy of recognition.

#### *2.5. Multi-Hop Reasoning for Retrieval Task*

Previous research on popular GNN frameworks has shown promising results in retrieval tasks requiring multi-hop reasoning. Coref-GRN aggregates information in disparate references from scattered paragraphs [31]. The different mentions of the same entity are combined with a graph recurrent neural network (GRN) to produce aggregated entity representations [32]. Based on Coref-GRN, MHQA-GRN [33] refines the graph construction procedure with more connections, which shows further improvements. Entity-GRN [34] proposes a method to distinguish different relations in the cognitive graphs through the convolutional neural network. Besides, Core-GRN [35] explores the cognitive construction problem. Nonetheless, it is rarely investigated by researchers about how to effectively reason on the constructed cognitive graphs, which is the primary problem studied in this paper. Other frameworks, such as Memory Networks [36], deal with multi-hop reasoning. These frameworks develop representations for queries and retrieval sources, between which interactions are then made in a multi-hop manner. The IR-Net [37] generates an entity state and a relation state at each step, computing the similarity degree between all entities and relations given by the dataset. However, these frameworks perform reasoning on simple datasets with a limited number of entities and relations, which is quite different from our work with the large-scale retrieval and intricate search queries.

#### *2.6. Pageview*

The pageview is a collection of Web objects or resources representing a specific user event, such as clicking on a link or viewing a product page [38]. Besides, pageview frequency has been shown to help improve the performance of evidence retrieval results [39].

#### *2.7. Peasoning Chains*

A reasoning chain is a sequence that logically connect users' search requests to supporting facts [40]. Reasoning chains should be intuitively related: they should exhibit a logical structure [41], or some other kind of textual relation that would allow human readers to quickly obtain what they really need.

#### **3. System Architecture**

#### *3.1. Overview*

In this section, we describe our new graph-based reasoning method that learns to find supporting evidence as reasoning paths for complex retrieval tasks. Our inspiration is motivated by the human cognitive process for searching an open question online, and our method starts from selecting the critical entity in the retrieval query. In order to make better use of the current entity's surrounding information, our method is to connect the start entity to some related entities either searched in the neighborhood or connected by the same surface mentioned. The abovementioned process is repeated to construct reasoning chains.

The main part to implement the 2SCR-IR framework is comprised of two components in our proposed system (Figure 4): a system to construct the cognitive graph, on which the other system reasons explicitly. Our method uses a strong interconnection between evidence to prevent important documents from being omitted in the reasoning paths.

**Figure 4.** Model architecture of the proposed 2SCR-IR sensor.

We used Wikipedia to retrieve relevant information for an open-domain search request, in which each article falls into paragraphs, resulting in billions of paragraphs in total. Any paragraph *p* can be perceived as our retrieval target. Given a query *q*, our framework aims at deriving its relevant textual sources to make up evidence chains by retrieving and reasoning paths, each of which is represented with a sequence of paragraphs.

#### *3.2. Gold Paragraph Selection*

To ease the computational burden of following steps, we first applied paragraph selection to the given context. Documents with titles matching the whole query were our priority. The selector will calculate the relevance score of the query and each paragraph between 0 and 1. We can set a threshold parameter τ, and paragraphs with predicted scores greater than τ are chosen, which are concatenated together as the context C. Since not every piece of text is relevant to the query, we further search for paragraphs that contain entities appearing in the searching query. If the two previous ways fail, a BERT-based paragraph ranker [42] will be used to select the paragraph with the scores higher than τ. Q and C are further processed by upper layers. After the first-hop paragraphs are identified, the next step is to search for evidence within paragraphs leading to other relevant paragraphs. Unlike the dependence on entity linking, we use hyperlinks in the first-hop paragraphs to find second-hop paragraphs to avert introducing noise. After all the links are chosen, we add edges between the evidence with these links and hyperlinks. When the two-hop paragraph is selected, we can obtain some candidate paragraphs. In order to further lower noise, we use the paragraph ranking model based on the BERT encoder to select the top-*N* paragraphs in each step. Besides, a retrieval query *q* and each of the 10 paragraphs *pi* are introduced into a BERT model [42]. Furthermore, a soft-max function is utilized to calculate the probability *P* (*q*, *pi*) of *pi*. A paragraph *pi* is chosen as a gold paragraph for the query *q* if *P* (*q*, *pi*) is larger than the setting τ. To collect most of the gold paragraphs for each retrieval query, the threshold value of 97.0% for recall and 69.0% for precision was set.

#### *3.3. Encoder and Attention Module*

We introduced a BERT model to acquire the representation of the query *Q* = [*q*1, *q*2 ... *qk*] ∈ R*k*×*d* and that of the context *C* = [*C*1, ... *Cl*] ∈ R*l*×*d*, where *k* and *l* are the lengths of the searching query and the context, and *d* is the size of the BERT hidden layer.

To gain more semantic information, following [43], we use the bidirectional attention module to update the representation for each word and a weighted self-aligned vector to represent the query.

The representation of the context *C* is

$$\mathcal{C} = [\mathcal{C}; \mathcal{C}Q; \mathcal{C} \odot \mathcal{C}Q; \mathcal{QC} \odot \mathcal{C}Q] \in \mathbb{R}^{l \times 4d\_b} \tag{1}$$

where is elementwise multiplication, and *dh* is the hidden dimension.

The representation of the searching query is

$$
\varphi = \operatorname{softmax}(w\_{\varphi}h) \tag{2}
$$

$$q = w\_a \\ \tanh\left(\mathcal{W}\_a \sum\_k q \rho\_k l\_k\right) \tag{3}$$

where *h* represents the query embedding; *hk* denotes the *k*-th word in the query; and *w*ϕ, *wa*, and *Wa* are linear projection parameters.

#### *3.4. Construction of Cognitive Graph*

To perform the multi-hop reasoning, we first needed to construct a graph of the information source to cover all the Wikipedia paragraphs. The Stanford CoreNLP Toolkit [44] was employed in an attempt to generate semantic graphs from the source text. The cognitive graph was constructed with entities as nodes and edges, while the number of extracted entities was denoted as *N*. We constructed the entity graph based on Entity-GCN [34] and DFGN [16], signifying that all mentions of the candidates searched in the documents are used as nodes in this cognitive graph. We used hyperlinks in Wikipedia to develop the direct edges, while the undirected edges were defined in line with the positional properties of every pair of nodes. The edges in this graph were split into two categories:

It needs to be emphasized that we do not apply the co-reference resolution to pronouns because it will introduce complex and unnecessary links.

#### *3.5. Multi-Step Reasoning*

After context encoding, the 2SCR-IR sensor performs reasoning over the cognitive graph. With the embedding process of the query Q and the context C, there is a huge challenge: to identify support entities and the text span of potential evidence paragraphs, as well as capturing relationships between evidence documents with little lexical overlap or semantic relationship to the original query. Based on the Cognitive Graph [6] and the DFGN method [16] for the QA system, an advanced knowledge-fusion module is adopted to mimic humans' step-by-step thinking and reasoning behavior—starting from *Q*0 as well as *C*0 and finding one-step supporting entities. The module achieves the following:

After obtaining the entity embeddings from the input context *Ct*, we applied a graph neural network to propagate the node information to its neighbors. Furthermore, a dynamic graph attention mechanism, which is a hierarchical and top-down process, was utilized to imitate humans' step-by-step exploring and reasoning behavior. It first attentively reads all knowledge graphs and then all triples in each graph for the final word generation. Yan et al. [45] show that the higher relevance given to the query can help the neighbor entity obtain more information from nearby. In order to calculate the degree of relation between the context and the query, 2SCR-IR calculates the attention between them. If one entity is more relative to the query, the entity is more pivotal. We followed [16] to multiply this relevance score:

$$\mathbf{v}\_{i}^{(t)} = \operatorname{sigmoid}\left(\frac{\tilde{\mathbf{q}}^{(t-1)}\mathbf{v}^{(t)}\mathbf{e}\_{i}^{(t-1)}}{\sqrt{d\_{2}}}\right) \tag{4}$$

$$
\tilde{\mathbf{E}}^{(t-1)} = \mathbf{y}^{(t)} \odot \mathbf{e}^{(t-1)} \tag{5}
$$

where V(*t*) is a linear projection matrix; *d*2 represents the dimension of each entity; and γ(*t*) *i* is the relevance score of the entity *ei* in step *t*. The entity graph is updated in a loop.

As a result, this step of the information propagation is confined to a dynamic scope on the cognitive graph. The next step is to disseminate information across the dynamic scope.

We introduce how to integrate the output of the entity into the context. To identify which sentences are more crucial to answer the question, we used a weighted self-aligned vector to represent the context. For each sentence *si*, we compute weight ε*i* by

$$\varepsilon\_{i} = \operatorname{softmax}(w\_{\psi}s\_{i})\tag{6}$$

The context vector is computed by

$$\mathbf{C}\_{V} = w\_{a} \tanh\left(\mathcal{W}\_{a} \sum\_{i} \varepsilon\_{i} \mathbf{s}\_{i}\right) \tag{7}$$

where *w*ϕ, *wa* are linear projection parameters, the context vector *CV* ∈ R*d*.

After the context vector is concatenated to each word in the context, the *i*-th word *ci* in the context C is updated as follows:

$$c\_i = [c\_i; C\_V] \tag{8}$$

After completing the contextual information, we used an LSTM layer to fuse the output of the entity graph and the context:

$$\mathbb{C} = LSTM([\mathbb{C}; E]) \tag{9}$$

The structure layer we used is the same as that used by Yang et al. [2]. This prediction process has three functions, including support sentences as well as the start and end position of the information needed, denoted as *Osup*, *Ostart*, *Oend*, respectively. We put three RNNs *Pi* layer by layer to solve the problem of output dependency. The information of the fusion block is sent to the first RNN *P*0, and the output of this module is as follows:

$$\mathbf{O}\_{\text{sup}} = P\_0 \left( \mathbf{C}^{(t)} \right) \tag{10}$$

$$\mathbf{O}\_{start} = P\_1(\left[\mathbf{C}^{(t)}, \mathbf{O}\_{sup}\right]) \tag{11}$$

$$\mathbf{O}\_{end} = P\_2(\left[\mathbf{C}^{(t)}, \mathbf{O}\_{sup}, \mathbf{O}\_{end}\right]) \tag{12}$$

Each *Pi* outputs a logit **O** ∈ R*M*×*d*2 .

In order to optimize the combined effect of these three sections, we compute these three cross entropy losses:

$$L = \mu\_{\text{start}} L\_{\text{start}} + \mu\_{\text{end}} L\_{\text{end}} + \mu\_{\text{support}} L\_{\text{support}} \tag{13}$$

Our method learns to extract reasoning paths within the constructed cognitive graph as evidence chains. Besides, evidence resources for a complex searching query do not necessarily have lexical overlaps. We adopted an automatically adjusting sub-scope to collect those entities with some sort of relation. Once one of the resources is retrieved, its entity mentions and the query often entail another resource.

#### *3.6. Producing Plausible Evidence Chains*

Finally, the 2SCR-IR sensor will provide users with reasoning chains, which first verifies every reasoning path, and then extracts a retrieval span from the most reasonable paths using a common method [42]. Our retriever model learns to predict plausible reasoning paths by capturing the paragraph interactions through the BERT (CLS) representations, after independently encoding the paragraphs along with the searching query. This paragraph interaction is crucial for multi-hop reasoning [46], especially when faced with open-domain information. However, as there is invariably some noise or misinformation, the best evidence chain is not always enough to fully cover all information to be retrieved. Hence, we re-computed the probability of the path with the increase in the pageview to lower the uncertainty and make our framework more robust. In the end, we select the top-*m* evidence chains for users.

#### *3.7. Example*

After *N* paragraphs and the query are put in, the selector module will filter out irrelevant paragraphs. Then, BERT is used as the encoder module to represent the filtered content in line with Formulas (1)–(3). After the process of knowledge representation, a classification layer is applied to the prediction of relevance scores between the paragraphs and the query. One paragraph will be signed as 1 if it contains supporting facts. Otherwise, it will be signed as 0. In the inference process, the paragraph will be chosen if the relevance score exceeds τ. In the process of reasoning, in order to make inferences on information, we calculated the degree of the relationship between the entities and the query in the light of Formulas (4) and (5) and built the entity graph. After we got the output of the entity graph, it was fused into the context. We used Formulas (6)–(9) to recognize the relevance of the search request to sentences. Next, Formulas (10)–(12) were utilized to predict the supporting sentences, as well as the start and end position corresponding to each search request. Moreover, Formula (13) was employed to optimize the combined effect. Finally, we will get supporting facts and construct reasoning chains in correspondence with the search request.

#### **4. Experiments**

#### *4.1. Datasets*

We evaluate our method on an open-domain Wikipedia-sourced dataset named HotpotQA [2], which is a human-annotated large-scale dataset that needs multi-hop reasoning and provides annotations to evaluate the prediction of supporting sentences. Almost 84% of queries require multi-reasoning. We used the full-wiki and the distractor setting of HotpotQA to conduct our experiments, with our primary target being the full wiki setting. Due to its open-domain scenario, we used the distractor setting to evaluate the performance of our method in a closed scenario where the evidence candidates are provided. For the sake of optimization, we used the Adam Optimizer [47] with an initial learning rate of 0.0005 and a mini-batch size of 32.

In addition, a good information retriever should be robust enough to prevent noise from disturbing themselves. IR can be armed with the well-structured and large-scale knowledge graph DBPedia [48] to offer some robust knowledge between concepts. However, the biggest disadvantage of this approach lies in that our cognitive model cannot be trained end-to-end and the errors may be cascaded. Inspired by [23], our retriever is trained to recognize relevant and irrelevant paragraphs at each step. We therefore mixed the noise information (negative examples) and generally right paragraphs together; to be more specific, we used two types of negative examples: Term Frequency–Inverse Document Frequency-based (TF–IDF-based) and hyperlink-based ones. When it comes to the single-hop retriever, merely the former type is used. Regarding the multi-hop retriever, we used both sorts, with the latter type carrying more weight to prevent our retriever from being distracted by reasoning paths without correct answer spans. Generally, the number of negative examples is set to 50.

#### *4.2. Implementation Details*