An array virtual extension method for power equipment acoustic emission direction of arrival estimation based on high‐order cumulants

National Basic Research Program of China (973 Program), Grant/Award Number: 2017YFB0902705 Abstract The performance of array‐based power equipment acoustic emission (AE) source directional of arrival (DOA) estimation methods is decided by not only algorithm performance, but also the array hardware parameters, that is, array aperture, sensor number, and so forth. Improving the estimation performance by hardware approaches always leads to the increase in cost. In this study, an array virtual extension method based on high‐order cumulants is proposed which can improve the sensor array performance without the need for hardware change or extra investment. Compared to low‐order cumulants, the high‐order cumulants established from Taylor expansion contain higher order terms which are omitted in low order cumulants. Therefore, they have more detailed directional information about the signal source to improve the DOA estimation performance. Moreover, because high‐order cumulants contain information only for the non‐Gaussian components, the array after extension shows better anti‐Gaussian noise performance. Depending on the establishment method of the high‐order cumulant, two types of extension including aperture extension and sensor number extension are proposed for different application scenarios. Simulation and field tests results show that the array directional resolution is improved around 40% by aperture extension and the spectrum sidelobe level is reduced about 14 dB by sensor extension.


| INTRODUCTION
Power equipment is an important and fundamental part of the power system. The safe and reliable operation of power equipment is of great significance to the stability and safety of the power system [1,2]. Power equipment on-line condition monitoring based on acoustic emission (AE) directional of arrival (DOA) estimation can not only provide early warning for the power equipment defects, but also prevent the spreading of further failure [3]. Common power equipment AE sources include ceramic capacitors [4] with physical defects generating acoustic waves yielding to the mechanical stress [5,6], air around electrical discharges emitting AE due to the vibration caused by the unbalanced electromagnetic force [7,8] and other localized mechanical fracture may occur in power equipment [9].
In recent years, several numbers of sensor-array based AE source DOA estimation methods are proposed to improve the estimation performance by software approaches. The signal sparse decomposition method is utilized in Xie et al [10]. to improve the subspace decomposition accuracy which eventually improves the DOA estimation accuracy. A new searching algorithm is proposed in Wang et al [11]. for the estimation process of power transformer partial discharge detection. The searching method of steered response power (SRP) algorithm is studied and improved in Nunes et al [12] while Diaz-Guerra and Beltran [13] utilized neural network to reduce the algorithm computation complexity.
Without denying the merits of these algorithms, it should be noted that the performance of these software approaches is heavily related to the hardware parameters of the utilized acoustic sensor array, that is, array aperture, sensor number, sampling frequency, and so forth [14]. For example, the array directional resolution which is also the wideband signal beamforming ability, is mainly decided by the array aperture, that is, the distance between sensors [15]. The array anti-interference and anti-signal attenuation ability can be enhanced by increasing the number of sensors in the array [16]. The array DOA estimation accuracy and anti-noise ability are also highly related to the parameters of the utilized sensor array [14]. Therefore, many researchers focus on developing advanced and more accurate sensor array to improve the estimation performance by hardware approaches [17]. However, improving the estimation performance by hardware approaches would always lead to extra hardware and investment cost, which would probably limit its large-scale actual applications. Thus, it would be meaningful if we can upgrade the array hardware parameters with software approaches because it can not only improve the DOA estimation performance but also spare the investment cost.
In this study, an array virtual extension method based on high-order cumulants, which can improve the sensor array performance without the need of hardware upgrade is proposed. The main idea of design is to replace the original signal low-order cumulants with the high-order cumulants established from Taylor expansion. Compared to low-order cumulants, the established high-order cumulants contain higher order terms which are omitted in low-order cumulants [18]. Therefore, they have more detailed directional information about the signal source to improve the DOA estimation performance. Moreover, because high-order cumulants contain information only for the non-Gaussian components [19], the array after extension shows better anti-Gaussian noise performance. Depending on the establishment method of the high-order cumulant, two types of extension including aperture extension and sensor number extension are proposed for different application scenarios. The aperture virtual extension method can increase the main lobe width of the DOA estimation spectrum, which can increase the array directional resolution. The sensor virtual extension method can reduce sidelobe amplitude of the spectrum, thus reducing the signal attenuation. Both methods have also shown certain level of anti-noise ability.
Simulation tests results show that the array directional resolution is improved around 40% after the aperture virtual extension. The spectrum sidelobe level which indicates the signal attenuation is reduced around 14 dB after the sensor number extension. Both methods have also shown certain level of antinoise ability which is improved around 40% under high noise level environment. Field test results have also verified the effectiveness of the proposed methods. Our research may promote the application of sensor array-based power equipment online monitoring.
Notations: In this study, bold upper-and lower-case letters stand for matrices and vectors respectively. I N ∈ ℂ N�N the identity matrix. The symbol (⋅) T denotes transpose, (⋅) * denotes conjugate transpose, E[⋅] denotes expectation. The notation ≜ denotes definition.

| Acoustic signal and array model
We focus our discussion on far-field array which means the distance from the source to the array is much larger than the array's aperture size. Suppose the acoustic signal received by the reference microphone in the linear array is modelled with y (t) [10,19]: where t is the sampling point, s(t) is the signal emitted by the acoustic source, h(d s , t) is the impulse response related to the position and direction from the source to the array (denoted as d s ) and v(t) stands for the uncollated noise. Suppose there are M microphones in the array. Due to the wave path difference, there is a time difference τ m (m ¼ 0, 1, …, M-1) between different sensors receiving the same signal. Thus, the signal received by the m th sensor is denoted as y m (t) ¼ y(t-τ m ) with its Fourier transformation as [13]: where ω is the angular frequency, Y(ω) is the Fourier transform of y(t). Therefore, the signal received by the array can be denoted in the frequency domain as: where a → is denoted as the array response vector. The response vector is dependent on the time difference τ, which is decided by the array parameter and array type. Take the M element uniform linear array (ULA) in Figure 1 as an example.
For far-field array, we assume the incident angle θ for different sensors are the same. Thus, the time difference is path difference, as labelled in the figure divided by signal speed: τ m ¼ d ⋅ cosðθÞ ⋅ m=c where c is the wave propagation velocity, m denotes the mth sensor in the array. d denotes the array aperture which is the element spacing. The array response vector of ULA is: As for the M � M element uniform rectangular array (URA) with d as aperture, we set the bottom left sensor as a reference point. θ denotes the incident angle and ϕ denotes the pitch angle which is the angle between the signal source and the horizontal plane (x-y plane) as shown in Figure 1b. According to the folding angle formula in solid geometry, the relationship between θ, ϕ and β is cosðβÞ ¼ cosðθÞcosðϕÞ. Thus, the time delay of the mth sensor in the row of URA can be calculated as τ m ¼ d cosðθÞcosðϕÞm=c and that of the mth sensor in the column is τ m ¼ d sinðθÞcosðϕÞm=c.Eventually, the corresponding array response vector can be written as: where a → URA;r is the response vector of URA for rows and a → URA;c for columns. The URA response vector in matrix is:

| Wideband signal beamforming
The signal model introduced in Section 2.1 focus on narrowband signals, which means the signal bandwidth is much smaller than its central frequency. The acoustic signals of AE sources are typical wideband signals, which should be preprocessed before applying the DOA estimation methods [14]. This process is called the wideband signal beamforming, as shown in Figure 2.
The main idea of wideband beamforming is that we first divide the frequency domain wideband signals into small segments. These segments are called frequency bin and, in each bin, the signal can be viewed as a narrowband signal, which means the methods for narrowband signals can be utilized. The width of each bin is the sampling frequency divided by the discrete Fourier transform (DFT) length L.
To do so, we first transfer the time-domain signal into the frequency domain with DFT after receiving the sampled signal y m (t) at each of the sensors. In Figure 2, the l th bin of each signal is first multiplied with a corresponding time delay element e -jωτ to make sure they are time-aligned. This time delay is called steering delay which is equal to the values of the propagation delays as are explained in Equations (4) and (5). Then the signals are multiplied with a filter H(ω l ). These filters are often adapted during the beamforming process to enhance the desired source signal while attenuating others.
Eventually, they are added together to produce a single, focused signal Z(ω l ), where l ¼ 1, 2, …, L: The summation in Equation (7) is called the steered response. It focuses the multi-channel signals into a single spectrum, which can then be used to scan, or steer, over a predefined spatial region by adjusting its steering delays. Applying the inverse DFT (IDFT) to Z(ω) gives us the enhanced version of time-domain signal z(t).

| SRP-PHAT DOA estimation method
In this section, we introduce the SRP with phase transform (SRP-PHAT) DOA estimation method. Recall in Section 2.1 and 2.2, the direction of the signal source θ is assumed to be a known value in order to derive the steered response. In many cases, these values are unknown and should be estimated. To do so, we express the steered response as a function of both DOA and angular frequency as Z(ω, θ). The power of the steered response Z(ω, θ) is called the SRP as [12,13]: where (⋅) * denotes the complex conjugate. Substitute Equation (7) into Equation (8) where (ω l )) is the cross power spectrum where E(⋅) stands for the mathematical expectation operator.
By analysing Equation (9), we can see that the value of SRP P(θ) is related to the weighting function Ψ mn , received signals Y m , Y n , and the time delay element e -jωτ . Y m , Y n , and Ψ mn are all determined values which can be calculated with the received signals. Thus, the value of P(θ) changes only with the changes of time delay τ, which is essentially decided by the incident angle θ as we explained in Section 2.1. Thus, the DOA of the acoustic source can be estimated by: The selection of filter Ψ mn depends on the application aspect. Here we introduce the PHAT function as: This weighting function enhances the influence of signal phase which has shown well performance in low signal to noise ratio (SNR) environment.

| Fourth-order cumulants
We first introduce the definition of the moments and cumulants. Suppose the probability density function of variable x is f (x). We denote the Fourier transform of f(x) as its first eigenfunction Φ(ω): Expand it by applying the Taylor expansion theorem, we have: where O(ω n ) stands for the high order terms which are omitted. The relation between m k and the k th derivative of Φ (ω) is: where m k is called the k th -order moment of variable x. Define the second eigenfunction of x as Ψ(ω) ¼ In Φ(ω). Following the similar process, the k th -order cumulant c k of variable x can be obtained as: The difference between moments and cumulants is that the kth-order moment is only the kth coefficient of the Taylor expansion of Φ(ω). But the kth-order cumulant is composed of all the first k coefficients of this Taylor expansion. Thus, a higher order cumulants should contain more information about the variable than which are omitted in a lower order one.
We suppose that Y m , Y n , Y o , and Y p and are the complex random series with zero mean. We directly give their cumulant from second order to fourth order [18,19]: Because Y m and Y n are normalized to zero-mean, the second term in second order cumulant c 2 (m, n) is 0, which means c 2 (m, n) is equivalent to the cross-power spectrum of the two complex series.
Based on the cumulants introduced above, Figure 3 shows the road map of three algorithms: the original SRP-PHAT method, the array aperture virtual extension method and sensor virtual extension method based on SRP-PHAT. Depending on the establishment method of the high-order cumulant, the extension can be divided into aperture extension and sensor number extension.
The difference between the three algorithms lies in the calculation of power spectrum W mn in Equation (9). In SRP-PHAT, W mn is the power spectrum E(Y m ⋅Y n * ), which is also the second order cumulants between a pair of signals (Y m , Y n ). For aperture extension method, W mn is replaced with a variant of fourth order cumulant c 4 (Y m Y n * Y m Y n * ) as will be introduced in Equation (17). For virtual sensor extension, it is replaced with another variant of fourth order cumulant c 4 (Y r Y x * Y r Y y * ) as shown in Equation (21). After the time delay and adding weighting function process, DOA can be obtained by peak search of their respective power spectrum. The detailed explanation of the extension process is discussed next.

| Array virtual aperture extension
To discuss the array virtual extension, we should first look back to the SRP algorithm. The calculation of the cross power spectrum W mn (ω) in Equation (9) can be viewed as the calculation of second order cumulant between signal Y m and Y n once the signals are normalized to zero-mean. In most cases, second order cumulants cannot completely preserve the signal information because the higher-order non-zero terms are omitted in the Taylor expansion. However, they are preserved in the high-order cumulants. Thus, by substituting the second order cumulants in Equation (9) with higher order ones, we should improve the algorithm performance. A fourth order cumulant with expression below is used in the proposed array aperture virtual extension method: To explain Equation (17), we first recall the definition of cumulants and the array response vector introduced in Section 2. The original second order cumulant used in Equation (9) can be written in the form of: where σ m is the signal power received at sensor m. Similarly, the fourth order cumulant Equation (17) can be written as: where γ 4,m represents the fourth order cumulant of Y m .
Comparing the cumulants in Equations (18) and (19), an obvious difference is that the time delay value is extended from t to 2τ. Recall the explanation that the time delay is decided by the incident angle, sensor distance, and signal speed in Equations (4) and (5). Assuming that the source direction θ and sound speed are fixed values, the extension of time delay from τ to 2τ can be viewed that the distance between different sensors, that is, the array aperture is doubled. An illustration is given in Figure 4a, b where the solid line represents the original aperture and the dotted line represents the virtual extended aperture. Because the array aperture is one of the critical parameters determining the array performance, utilizing the fourth order cumulant can be viewed as the virtual array extension which can increase the array performance.
Note there are different variants of cumulants based on the different ways to configure the complex conjugate on the parameter to c 4 in Equation (17). By adjusting the different parameters, array virtual sensor extension can be achieved which will be introduced in the next subsection. Figure 4c illustrates the array sensor virtual extension where the solid circles represent the real sensors and the white circles represent the virtual sensors. r 1 (t), x(t), y(t), and v 1 (t) are respectively the signals received by reference sensor, two real sensors, and virtual sensor. Their complex form after Fourier transform are Y r1 (ω), Y x (ω), Y y (ω), and Y v1 (ω), respectively. The second cumulant between r(t) and v(t) writes as:

| Array virtual sensor extension
where k → is the propagation vector of the source signal. d dependent on the time difference, which is dependent on array arrangement. In Equation (20), because v 1 (t) stands for the signal received by virtual sensor, we cannot directly calculate the second order cumulant of c 2 . Instead, we use the fourth order cumulant of r(t), x(t) and y(t) to replace it: where γ 4,r1 represents the fourth order cumulant of r 1 (t). Note the positional relation between the real sensors and the virtual sensor is d Thus, the relation between c 4 and c 2 can be obtained as: This means the second order cumulant of the virtual sensor can be obtained from multiplying the fourth order cumulant of actual sensors with a constant term σ r1 2 /γ 4,r1 . Most importantly, multiplying a constant term would not interfere the peak search of DOA estimation in Equation (10). Thus, we can substitute the second order cumulant with fourth order cumulants in the SRP calculation. This can also be viewed as a virtual array extension by increasing the virtual sensors in the array, which should improve the array performance.
To obtain the second order cumulant from other virtual sensors, we only need to change the elements in the fourth order cumulant to satisfy the positional relation between these sensors. For example, to calculate the second order cumulant of virtual sensor v 2 (t), we can calculate c 4 ) as long as the positional relation is satisfied.

| Simulation steps and conditions
To verify the effectiveness of the proposed method, simulation tests are performed with the steps as follows: i. Determine the simulation parameters: Suppose the acoustic source emits broadband far-field acoustic chirp signals. A chirp is a signal in which the frequency increases or decreases with time, which is commonly used in acoustic signal simulations. Suppose the signals are received by a 16sensor URA with 7 cm array aperture. The sampling rate is 44.1 kHz/s. Direction angles are defined to be 0 degrees for a source located at the front of the sensor array and approaches �90 degrees as the source approaches the side. Gaussian white noise is added ii. Broadband signal processing: An oversampled uniform DFT filter bank is used to transform the wideband signal into sub band signals. The DFT length and filter length are both set to 512 with the increment size set to 256. Here the increment is the DFT increase from one bin to the next bin. The bin width is 86 Hz iii. Virtual array aperture extension: For a pair of sub band signals (take the mth and nth bin for example), calculate fourth order cumulant c 4 (Y m Y n * Y m Y n * ) with Equation (17). Substitute the corresponding power spectrum W mn in Equation (9) with the calculated fourth order cumulant. In this study, the aperture is extended twice as the actual aperture iv. Virtual array sensor extension: Select one actual sensor r 1 as the reference position and denote the virtual sensor as v 1 . Select two other actual sensors x and y which can satisfy with Equation (21). Substitute the power spectrum W r1v1 in SRP calculation Equation (9) with the calculated fourth order cumulant. The 4�4 array is extended to 6�6 array v. DOA estimation: The DOA is estimated with the estimator below: For aperture extension, cum 4 is in the form of Equation (17) and for virtual sensor extension is in the form of Equation (21). The direction corresponding to the maximum value of P (θ) is considered to be the DOA of the acoustic signal source.

| SRP spectrum analysis
To have an impression of the performance of the proposed methods, simulation tests are performed following the steps above. A 16 sensor URA with 7 cm array aperture is utilized as sensor array as shown in Figure 5. A simulated acoustic source is placed at the incident angle of 30 degree and the SNR is set to 30 dB. Other simulation parameters are as described in Section 4.1. Figure 6 shows the SRP spectrum obtained using the SRP-PHAT algorithm, the array aperture virtual extension (abbreviated as APE-EXT) and the array sensor virtual extension (abbreviated as SEN-EXT) method.
The red dash line with the circle mark is the spectrum of the original SRP-PHAT method. It has a main lobe and two obvious sidelobes. The dashed line with square mark and dashdot line are respectively the SEN-EXT and APE-EXT method. As discussed in our previous research, the directional resolution of a DOA estimation method is in inverse proportion to the width of its main lobe [15,20]. By applying the virtual extension method, the main lobe width in both of the two SRP spectrum is narrowed, which means the directional resolution of the SRP algorithm is improved. Here the directional resolution which is also the wideband signal beamforming ability indicates the ability of an algorithm to distinguish multiple sources and resist the influence of acoustic reverberation. In actual applications, improving the directional resolution is often realized by increasing the array aperture. The spectrum of APE-EXT method shows that the proposed virtual array aperture extension method can achieve the effect of array extension to a certain extent.
Moreover, the sidelobes of the SRP spectrum are reduced after applying the virtual array extension method. Sidelobes diffuse sound energy and increase attenuation and the sidelobes of a strong signal may interfere with the main lobe of other signals [16]. Results in Table 1 are the main lobe width and sidelobe level of the three methods averaged from 50 monte carlo tests. It can be seen that the array directional resolution is improved around 40% after the aperture virtual extension and the spectrum sidelobe level is reduced around 14 dB after the sensor extension. Meanwhile, the APE-EXT method is better at narrowing main lobe and the SEN-EXT method does better in reducing sidelobes.

| Directional error and array arrangement
Change the incident angle of the acoustic source from 0 degree to 90 degree and perform DOA estimation with the three methods respectively. The results are obtained with both the 16-sensor ULA and URA. SNR is set as 5 dB. Figure 7 shows the relation between the incident angle of the acoustic source and the DOA estimation root mean square error (RMSE). The SRP-PHAT method, aperture virtual extension and sensor virtual extension methods are tested with both the ULA and URA.
For methods tested on ULA, the more the incident angle approaches 90 degree, the higher the estimation error is. This is because the time delay between different sensors of ULA is t ¼ d•sin(θ)/c which is related to θ. By substitution, we have θ ¼ arcsin(cτ/d). Perform derivation of the relationship and we have: � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À ðcτ=dÞ which means the DOA estimation error is inversely proportional to cosθ. This is in accordance with our results in Figure 7.
To solve this problem, we rearranged the 16�1 linear array to a 4�4 rectangular array. Suppose 0 degree is located at the front of the sensor array and 90 degree is located at the side as in Figure 5. estimation RMSE is increased (mostly due to the decrease in sensor number), the error distribution is stable in all directions.

| DOA estimation under noise influence
To test the DOA estimation performance under noise influence, simulation tests are performed under different SNR conditions with the 4�4 URA. The SNR is calculated with SNR ¼ 20⋅log(A S /A N ) where A S is the signal amplitude and A N is the noise amplitude. The RMSE is calculated with: RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where θ i is the measured DOA and θ s is the actual DOA. N is the number of samples. Figure 8 shows the relationship between estimation RMSE and SNR obtained with different methods. It can be seen that the estimation error of all methods decreases as the SNR increases. The performance curve of the SRP method in Figure 8 is moved down after applying the PHAT weighting function, the aperture virtual extension and sensor virtual extension methods. The detailed data are presented in Table 2.
In Figure 8, we have also plotted a dotted line with the RMSE equals 1.5°which is a common standard for calculating DOA estimation success rate. It can be seen that the 1.5°R MSE is reached at SNR ¼ 4.1 dB for SRP method. The same level of accuracy is achieved with SNR down to 3.6, 2.8, and 2.4 dB for SRP-PHAT, sensor virtual extension and aperture virtual extension method respectively. This can also be viewed as applying the proposed array virtual extension have improved the anti-noise ability for the DOA estimation methods by 12%, 31%, and 41%, respectively. Moreover, the simulations results show that the estimation RMSE of all methods can be lower than 1°if the SNR is higher than 5 dB, which should meet the requirement of practical applications.

| Experiment steps and conditions
To further verify the effectiveness of the proposed array virtual extension method, field tests are conducted in an actual 110 kV substation site.
Different from the simulation tests, the AEs received in substation site contains not only white noise, but also the interference caused by signal refraction or reflection. Moreover, the received signals contain frequency components in a wide frequency range, which means the wideband signal beamforming is indispensable. A simplified illustration of the established AE source DOA estimation system, which consists of acoustic sensor array, pre-processing unit, 16-channel synchronous sampling unit and the computer is shown in Figure 9a. The acoustic sensors utilized in this study is the ICS-41352 MEMS acoustic sensor provided by Tokyo Denkikagaku Kogyo Electronics [21]. Its typical frequency response is 10 Hz to 15 kHz. In standard mode, it has sensitivity of À 26 dB, signal-to-noise ratio of 64 dB(A), total harmonic distortion (THD) of 0.2%, acoustic overload point of 120 dB sound pressure levels, and supply current of 430 μA. The sampling frequency for all the 16-channel sampling unit is set as 44.1 kHz.
The experiment steps are similar to the simulation steps in Section 4.1. The main difference is the signal acquisition method. A brief description of the experiment steps is described as below:  (23). Each test is conducted for 50 times with and without applying the array virtual extension method. The results are as shown in Figure 10 and Table 3

| Result analysis
We first calculate the SNR in the field environment. Take the signals in Figure 9b as example, the signals in [0s, 0.01s] are viewed as environment noise and those in [0.015s, 0.040s] are viewed as AE signals. The SNR calculated with the equation given in Section 4.4 is around 15 dB. Figure 9 gives the SRP spectrum calculated with the signals in Figure 8b before and after the virtual extension. Comparing the results in Figure 10 to Figure 6, it can be seen that the main lobe width in the spectrum is narrowed after applying the virtual extension method, which means the directional resolution of the SRP algorithm is improved. To be more specific, the main lobe width is calculated similar to the -3 dB bandwidth, which is the direction width from the peak value point of the power spectrum to the 1/2 peak value point. The main lobe widths of all the three methods have been plotted in Figure 10. By applying the sensor extension, the main lobe width is narrowed from 47.5 degree to 15.5 degree, and it can be further narrowed to 7.5 degree applying the aperture extension method. The directional resolution is improved 67% and 84%, respectively.
Note the main lode width of the aperture virtual extension method is smaller than that of the sensor virtual extension method. This is also in accordance with our results in the  simulation tests. However, for the sensor virtual extension method, it has shown better performance in reducing sidelobes. There is no local peak point other than the highest point. The relative value from the peak point to the bottom is increased by around 60% from 0.58 to 0.95. The DOA estimation error distribution and maximum error of the different methods conducted in field tests are given in Table 3 and Figure 11. It can be seen from Table 3 that not only the average estimation error is reduced from 2.67°to around 0.6°, but the maximum estimation error is also reduced from 3.75°to 1.5°as well. Thus, the performance of the DOA estimation algorithm is improved about 80% after applying the array virtual extension.
The cumulative distribution function (CDF) [22] of the DOA estimation error is plotted in Figure 10. Its definition is CDF(a) ¼ P(x<¼a) where P(⋅) is the probability function. The CDF value represents the proportion of the results smaller than the value in x-axis. Thus, the curves in Figure 11 represents the error distribution of the different methods.
We can see that the CDF curves of the aperture extension method and sensor extension method are both above that of the original SRP method. This rises of CDF curve mean the estimation errors, both the average value and the variance are reduced by applying the array virtual extension method. Moreover, the slopes of the three cures at x-axis ¼ 0.5 are also calculated and given out in the figure. It can be seen that the slope of SRP method is increased from 0.243 to 0.625 after the array extension. The slope of the CDF curve represents the scale of the error distribution, that is, the error variance. Thus, the results have verified that the array extension methods can indeed improve the DOA estimation performance of the SRP method.

| CONCLUSION
In conclusion, this study presents an array virtual extension method based on high-order cumulants for power equipment AE source DOA estimation. It can improve the sensor array performance without the need of actual hardware upgrade. The principle and process of the virtual extension method have been explained. Simulation tests and field tests are conducted to verify the effectiveness of the proposed method.
i. By substituting the second order cumulants in SRP algorithm with the fourth order cumulants, the sensor array can be virtually extended. The selection of different form of fourth order cumulants leads to the difference of array aperture extension and array sensor extension ii. The virtual aperture extension method can improve the directional resolution of the DOA estimation array while the virtual sensor extension method does better in sidelobe suppression. Both methods can improve the anti-noise ability of the DOA estimation algorithm which is verified by simulation and experiment results iii. Without loss of generality, the proposed array virtual extension method could also be applied to other DOA estimation methods to improve the sensor array performance ACKNOWLEDGEMENT This work was supported in part by the National Key Research and Development Program (2017YFB0902705).