**Adaptive Intrawell Matched Stochastic Resonance with a Potential Constraint Aided Line Enhancer for Passive Sonars**

## **Haitao Dong 1,2,†, Ke He 1,2,†, Xiaohong Shen 1,2, Shilei Ma 1,2, Haiyan Wang 1,3,***<sup>∗</sup>* **and Changcheng Qiao 4,5**


Received: 1 May 2020; Accepted: 1 June 2020; Published: 8 June 2020

**Abstract:** Remote passive sonar detection and classification are challenging problems that require the user to extract signatures under low signal-to-noise (SNR) ratio conditions. Adaptive line enhancers (ALEs) have been widely utilized in passive sonars for enhancing narrowband discrete components, but the performance is limited. In this paper, we propose an adaptive intrawell matched stochastic resonance (AIMSR) method, aiming to break through the limitation of the conventional ALE by nonlinear filtering effects. To make it practically applicable, we addressed two problems: (1) the parameterized implementation of stochastic resonance (SR) under the low sampling rate condition and (2) the feasibility of realization in an embedded system with low computational complexity. For the first problem, the framework of intrawell matched stochastic resonance with potential constraint is implemented with three distinct merits: (a) it can ease the insufficient time-scale matching constraint so as to weaken the uncertain affect on potential parameter tuning; (b) the inaccurate noise intensity estimation can be eased; (c) it can release the limitation on system response which allows a higher input frequency in breaking through the large sampling rate limitation. For the second problem, we assumed a particular case to ease the potential parameter *aopt* = 1. As a result, the computation complexity is greatly reduced, and the extremely large parameter limitation is relaxed simultaneously. Simulation analyses are conducted with a discrete line signature and harmonic related line signature that reflect the superior filtering performance with limited sampling rate conditions; without loss of generality of detection, we considered two circumstances corresponding to *H*<sup>1</sup> (periodic signal with noise) and *H*<sup>0</sup> (pure noise) hypotheses, respectively, which indicates the detection performance fairly well. Application verification was experimentally conducted in a reservoir with an autonomous underwater vehicle (AUV) to validate the feasibility and efficiency of the proposed method. The results indicate that the proposed method surpasses the conventional ALE method in lower frequency contexts, where there is about 10 dB improvement for the fundamental frequency in the sense of power spectrum density (PSD).

**Keywords:** adaptive stochastic resonance (ASR); matched intrawell response; nonlinear filter; line enhancer; autonomous underwater vehicles (AUVs)

#### **1. Introduction**

Passive sonars have been proven to have practical efficiency in detecting and recognising selfemitting underwater targets such as ships, submarines, and autonomous underwater vehicle (AUVs), etc. [1–3]. Generally, narrowband discrete components of ship-radiated noise, known simply as lines or tonals, have largely prevailed so far [4–6]. From the view of rotatory machinery, ship-radiated line signatures that are relevant to the engine and shaft/propeller rotation attract lots of attention, of which the long-term challenging problem is to address the weak signatures from heavy noise background in the far field scenario [3,7,8]. In this way, increasing the signal-to-noise ratio (SNR) is expected as a tenet to detection performance for applications. To enhance the narrowband discrete components, passive sonars usually employ an adaptive line enhancer (ALE) as a pre-processing step [9–12]. However, performance of the conventional ALE is limited, and advanced signal processing techniques dedicated to better denoising performance are of highly practical importance.

In general, more noise in the system often leads to worse detection performance and degrades the estimation accuracy. However, despite its disruptive character in nature, noise does play an important constructive role under certain conditions [13]. This phenomenon so called stochastic resonance (SR), in which the output periodic signals can be enhanced by the cooperative effect between "signal", "noise" and certain "nonlinear systems" [14]. Since proposed by Benzi et al. in 1981 [15], it has become a hot research topic in the field of nonlinear science. Previous studies have addressed a lot of research progress, both theoretical and experimental, in achieving some noteworthy contributions in practical applications [8,16–21]. As a result, taking SR for weak signal detection is regarded as a potential novel technique for weak signal detection, especially under low SNR conditions.

Classical SR theories generally refer to a noise enhanced phenomenon by means of adding an appropriate amount of noise [14]. This is restricted to applications of how to remove proper noise, especially under low SNR circumstances. Rather than adjusting the input noise levels, Xu et al. [22] proposed a parameter-induced stochastic resonance (PSR) which highlighted the effect of the nonlinear systems and promoted the flexibility of SR utilization in designing systems to deal with the noise. In view of this, utilizing SR in weak signal processing could be considered as a special nonlinear filter, which can achieve superior performance compared with the traditional noise-suppression-based filters [8,23–25]. Among the publications on this point, the simplest first-order overdamped Langevin equation (LE) model is largely adopted, while the filtering performance is limited [8,26,27]. To get a cleaner filtered signal with higher signal-to-noise ratio (SNR), the second order Duffing equation with an underdamped system is used for a secondary filtering effect, which reflects a superior filtering performance compared with the traditional state-of-the-art methods [23,24,28]. For the purpose of better improvements, some authors reported cascaded SR systems [29,30], coupled SR system [31,32], etc. with bistable or multi-stable potentials. As a matter of fact, superior performance generally required a high sampling frequency condition to fully take advantage of the nonlinear property by concentrating most of the noise energy into the low-frequency region [33]. Several efforts have addressed the problem of large parameter stochastic resonance (LPSR) by tuning the signal structure and (or) the system parameters [24,27,34]. However, all the methods of LPSR are based on the classical SR part, and the classical SR generally requires a large sampling frequency that is more than 50 times the driving frequency. All these processing approaches are essential to deal with the periodic signal individually, and hence, flexible limited and lack of computation efficiency. An open problem for practical applications is "Can the SR be realized under limited sampling rate?"

To address this problem, limited studies on intrawell SR can give us inspirations. Actually, intrawell SR can exist for a periodically driven noisy signal in terms of underdamped bistable nonlinear dynamic systems [14]. Alfonsi et al. [35] gave a phemonological nonadiabatic description of intrawell SR by adding Lorentz' colored noise. On the basis of this, Li et al. [36] further demonstrate that the

intrawell SR phenomenon usually occurs in a system with optimal system parameters. Since the system response speed of intrawell motion is much larger than that of interwell, it can release the limitation on system response speed which allows a higher input frequency. In this way, the large sampling condition can be eased in properly controlling the escape time scale (characterizing interwell jumps) together with the relaxation time of the particle. Back to the view of potential parameter tuning, the adjustable potential parameters can be equivalent to tuning potential depth and well separation [28]. There is lack of flexibility in matching the escape time scale and the relaxation time only with the damping effect. In our previous work [37], the nonlinear filtering effects of intrawell matched stochastic resonance is analyzed, which have shown a superior filtering performance as well as a wider range of frequency response. This can give us a guidance to ease the large sampling frequency limitation.

On the basis of the aforementioned analyses, this paper proposes a novel parameterized filter implementation on adaptive intrawell matched stochastic resonance (AIMSR) with a potential constrained Duffing system. As can be seen in [24,37] matching a high frequency signal requires extremely large system parameters. This refers to a large range of parameter searching interval to limit the computational efficiency and parameterized realization, especially in an embedded system. In this regard, we further eased the potential parameter *a* = 1 which makes it practically realizable. The nonlinear filtering effects are analyzed and evaluated, reflecting a superior performance in enhancing the lines especially under a limited sampling rate. Application verification is further conducted with a set of low sampled data fragment of autonomous underwater vehicles (AUVs). The output performance is further compared with the conventional ALE method, which shows an excellent filtering performance in enhancing the lines, especially for the lower frequency fundamental signature.

The rest of the paper is arranged as follows. After introducing the signal model and measurement in Section 2, in Section 3, the theory and implementation of the AIMSR method are detailed and described with a potential constraint Duffing system. Section 4 verifies the validity of utilizing AIMSR in dealing with the single and muli-harmonic line signature signals. Section 5 verifies the practicability of the proposed method by analyzing a set of low sampled AUVs data fragments. Besides, intensive discussions are made to give an insight into the principal results and inspired future investigations. Finally, concluding remarks are drawn in Section 6.

#### **2. Signal Model and Measurement**

Ship radiated noise have been studied for years. Ross and Urick have given an excellent description of the mechanisms of sound generation by large surface ships and submarines [4,5]. They have shown that the radiated noise from a ship is a combination of broadband noise and sinusoidal tonal signals that are generated by many sources. Theoreitical and experimental results have shown signatures related to the rotation of engines, shaft-line dynamics, propeller cavitations are mostly efficient for the purpose of passive detection, tracking, and classification [1,6,38,39]. From the view of rotatory machinery of the shaft and propeller, the harmonic signature is typically connected to the running state and can be measured by vibrational and acoustic approaches. For an underwater surveillance system, to address the shaft and propeller signatures, algorithms such as the Detection of Envelope Modulation on Noise (DEMON) [40], cyclic modulation coherence (CMC) [7], etc. are employed to check for the presence of periodic components to the broadband cavitation noise. This generally requires acoustic measurement of broadband noise in a high-frequency range of thousands of Hertz. Additional studies in recent years have shown that the harmonic signature as a propeller rotates can be observed by low frequency acoustic/vibrational measurement as well, which may propagate for a very long distance [41]. Such measurement techniques have been ubiquitously seen in sonobuoys, ocean bottom seismometers (OBS), acoustic vector hydrophone and other small-scale equipment [41–43].

As a matter of fact, these sinusoidal tonals are commonly considered to be the "acoustic fingerprint" of a moving vessel. Table 1 shows the major contributions and relations to the sinusoidal tonal signals

from ship engine and propeller. In the situation of practical applications, these sinusoidal components can be modeled as follows,

(i) Discrete line signature

$$\begin{split} r(t) &= s(t) + n(t) \\ &= A\_0 \cos(2\pi f\_0 t + \varphi\_0) + n(t) \end{split} \tag{1}$$

where *f*<sup>0</sup> represent the character frequency, *A*<sup>0</sup> and *ϕ*<sup>0</sup> are the corresponding amplitude and phase, respectively. *n*(*t*) represent the combination of radiated broadband noise and ambient noise.

(ii) Harmonic related line signature

$$\begin{aligned} r(t) &= s(t) + n(t) \\ &= \sum\_{h=1}^{M} A\_h \cos(2\pi h f\_0 t + q\_i) + n(t) \end{aligned} \tag{2}$$

in which *h* is the harmonic number, *Ah* and *ϕ<sup>h</sup>* are the amplitude and phase of the *h*th harmonic component, and *f*<sup>0</sup> is the fundamental frequency of shaft rotation, *n*(*t*) represent the combination of radiated broadband noise and ambient noise.

**Table 1.** Character frequencies of ship engine and propeller.


For small targets such as AUVs, speedboats, etc. with rotation marine engine and propeller, the signatures are in a periodic pattern as well. Their acoustic characteristics have been comprehensively studied by measurement tests and modeling [2,39,44], which indicate that the model of radiated noise can be the same with large ships. Since an AUV's prop does not rotate sufficiently fast to cause cavitation, it is hard to be detected with broadband noise in the high-frequency range. As is known, the low frequency harmonics are quite important for remote passive detection, and the experimental study in this work is set up with a very low frequency measurement by ocean bottom seismometer (samapling frequency *fs* =200 Hz).

#### **3. Adaptive Intrawell Matched Stochastic Resonance**

#### *3.1. Generalized Matched Stochastic Resonance with Duffing Oscillator*

The bistable SR phenomenon can be described as a particle driven by periodic force and random force in a quartic double-well system, where the periodic motion can be heightened with the assistance of moderate noise. This can be governed by a two-dimensional Duffing oscillator as below [24],

$$\frac{\mathrm{d}^2 \mathrm{x}}{\mathrm{d}t^2} + \gamma \frac{\mathrm{d} \mathrm{x}}{\mathrm{d}t} = -\frac{\mathrm{d}V(\mathrm{x})}{\mathrm{d} \mathrm{x}} + s(t) + \mathfrak{J}(t) \tag{3}$$

where *γ* is the damping factor, *s*(*t*) represents the input periodic signal, and *n*(*t*) stands for the noise item with noise intensity *D*. In particular, Equation (3) describes the movement of a particle in a quartic double well potential *V*(*x*),

$$V(\mathbf{x}) = -\frac{a}{2}\mathbf{x}^2 + \frac{b}{4}\mathbf{x}^4, \qquad a, b > 0 \tag{4}$$

in which *a* and *b* are barrier parameters of bistable SR system, and the potential barrier can be calculated by Δ*V* = *a*2/(4*b*). As the analytic form of the SR output cannot be easily obtained, it is generally numerically calculated through the discrete fourth-order Runge-Kutta method [23]. Mathematically, let <sup>d</sup>*<sup>x</sup>* <sup>d</sup>*<sup>t</sup>* <sup>=</sup> *<sup>y</sup>*, then Equation (3) can be separated into two first-order differential equations as,

$$\begin{aligned} \frac{dx}{dt} &= y\\ \frac{dy}{dt} &= ax - bx^3 - \gamma y + s(t) + \xi(t) \end{aligned} \tag{5}$$

and then the output discrete time series *x*[*n*] corresponding to Equation (3) can be calculated by,

$$\begin{cases} y\_1 = y[n] \\ x\_1 = -V'(\mathbf{x}[n]) - \gamma y\_1 + \mathbf{s}[n] + \tilde{\xi}[n] \\ y\_2 = y[n] + \mathbf{x}\_1 h/2 \\ x\_2 = -V'(\mathbf{x}[n] + y\_1 h/2) - \gamma y\_2 + \mathbf{s}[n] + \tilde{\xi}[n] \\ y\_3 = y[n] + \mathbf{x}\_2 h/2 \\ x\_3 = -V'(\mathbf{x}[n] + y\_2 h/2) - \gamma y\_3 + \mathbf{s}[n+1] + \tilde{\xi}[n+1] \\ y\_4 = y[n] + \mathbf{x}\_3 h \\ \mathbf{x}\_4 = -V'(\mathbf{x}[n] + y\_3 h) - \gamma y\_4 + \mathbf{s}[n+1] + \tilde{\xi}[n+1] \\ \mathbf{x}[n+1] = \mathbf{x}[n] + (y\_1 + 2y\_2 + 2y\_3 + y\_4)h/6 \\ y[n+1] = y[n] + (\mathbf{x}\_1 + 2\mathbf{x}\_2 + 2\mathbf{x}\_3 + \mathbf{x}\_4)h/6 \end{cases} \tag{6}$$

where *h* is the calculation step of the Runge-Kutta method.

From the perspective of parameter tuning, it can be considered as a special nonlinear filter, and the main issue is to achieve an "matched filter" output [24]. In general, a bistable system response is optimized by maximizing the signal-to-noise ratio improvement (SNRI) on system potential parameters *a* and *b*, time scaling parameter *h*, and damping factor *γ*. In this way, the nonlinear matched SR filter can be generally modeled with a generalized time-scale matching constraint as,

$$\begin{array}{ll}\textbf{max} & \text{SNRI} \\ \text{s.t.} & r\_{\text{K}} = \frac{\Omega}{\pi} \\ & \gamma \in (0, 1], \ h \in (0, 1] \end{array} \tag{7}$$

in which *rK* <sup>=</sup> *<sup>a</sup>* <sup>√</sup>2*πγ* exp(<sup>−</sup> *<sup>a</sup>*<sup>2</sup> <sup>4</sup>*bD* ) is the famous Kramers rate, <sup>Ω</sup> <sup>=</sup> <sup>2</sup>*<sup>π</sup> <sup>f</sup>*<sup>0</sup> represents the angular frequency of periodic forcing. The constraint of *rK* = Ω/*π* is the time-scale matching condition that represent the statistical synchronization of interwell transition [14]. It is necessary to point out that this generalized condition is not sufficient as other optimized factors are required to judge the occurance of SR.

For sinusoidal signals with additive noise, the SNRI corresponding to Equation (3) could have a maximum at *DSR* = Δ*V*, and can be equivalent to *D*ˆ = *aopt*/4*bopt* in according to [24]. And as a consequence, the mathematic model of SR based nonlinear matched filter can be rewritten as,

$$\begin{aligned} \left(\gamma\_{opt}^{\*}, h\_{opt}^{\*}\right) &= \underset{\gamma, h}{\text{argmax}} \text{SNRI} \\ \text{s.t.} &\qquad a\_{opt} = 2\sqrt{2}\pi f\_{0}\gamma \text{e} \\ &\qquad b\_{opt} = a\_{opt}^{4}/4\dot{D} \\ &\qquad \gamma \in \left(0, 16\sqrt{2}\pi f\_{0}\text{e}\right), h \in (0, 1] \end{aligned} \tag{8}$$

where e is the base of natural logarithms, and the optimal potential parameters *aopt* and *bopt* can be obtained with matched relationship that related to the frequency of periodic signal *f*0, damping factor *γ* and the noise intensity estimator *D*ˆ .

To adaptively adjust these system parameters, the measurement index of input-output SNR improvement (SNRI) is extensively utilized and can be estimated by measuring the power spectral of a time series. This can be further limited with a bandwidth Δ*B* for the local SNR ( SNR*local*) calculation as it is generally connected to the performance of signal detection [45]. Besides, evaluation indexes such as power spectrum kurtosis (PSK), peak SNR (PSNR), spectral correlation coefficient (SC), structural similarity index (SSI), etc. can be well adopted and fused [20].

#### *3.2. Framework of Intrawell Matched Stochastic Resonance with Potential Constraint*

The "classical" description of the SR phenomenon in a bistable system is that a particle can jump across the potential barrier back and forth in the assistance of proper noise. This is called interwell jumping behavior as the random-switching frequency *rK* (Kramers rate) is made to agree closely with the periodic forcing angular frequency [14,35]. From the view of parameter tuning, the adjustable potential parameters can be equivalent to tuning the potential depth Δ*V* and well location (stable focus) ±*xm* [28], where Equation (4) can be transformed into the function of well location and potential depth parameter as,

$$V(\mathbf{x}) = -\Delta V \left[ 2 \left( \frac{\mathbf{x}}{\mathbf{x}\_m} \right)^2 - \left( \frac{\mathbf{x}}{\mathbf{x}\_m} \right)^4 \right] \tag{9}$$

In this way, changing the potential parameters *a* and *b* is essential to adjust the potential depth Δ*V* and the separation 2*xm*. A demonstration of particles in a bistable model with different potential parameters is given in Figure 1. It can be seen that a large potential depth can lead to longer system response time for transition, and even intrawell response [35,36].

As is known, the time-scale matching constraint is directly connected to the system response time. The nonlinear filtering effect of matched stochastic resonance can be essentially regarded as a proper control of the particle's response. From this view, it is easy to know that the more the degrees of freedom for parameter tuning, the better the output response. Such results can be found in plenty of SR-related publications with multi-parameter optimization. The classical Duffing equation model with bistable potential utilizes parameter *a*, *b*, and *γ* are used to characterize the system response. It is intuitive to add new degrees of freedom for better characterizing and controlling of the particle's response. Hence, a potential constraint is further adopted here to ease the insufficient time-scale matching constraint [37]. By this means, the mathematical framework of matched stochastic resonance can be further generalized by intuitively adding a constraint of barrier factor *K* as below,

$$\begin{aligned} \left(\gamma\_{opt}^\*, h\_{opt}^\*, K\_{opt}^\*\right) &= \underset{\gamma, h, K}{\text{argmax}} \text{ SNRI} \\ \text{s.t.} &\quad a\_{opt} = 2\sqrt{2}\pi f\_0 \gamma \text{e} \\ &\quad b\_{opt} = a\_{opt}^4 / 4K\hat{D} \\ &\quad \gamma \in \left(0, 16\sqrt{2}\pi f\_0 \text{e}\right), h \in \left(0, 1\right], K \in \left(1, K\_{max}\right] \end{aligned} \tag{10}$$

in which *K* should be a postive real number to adjust the particle motion, and can be optimized from a proper searching interval.

**Figure 1.** A demonstration of bistable stochastic resonance with interwell and intrawell motion.

The merits of this potential constraint model can be summarized as follows: (1) it can ease the insufficient time-scale matching constraint so as to weaken the uncertain effect on potential parameter tuning; (2) the inaccurate noise intensity estimation can be eased as well; (3) it can release the limitation on a system response which allows us a higher input frequency in breaking through the large sampling rate limitation. All in all, the potential constraint model is anticipated to be superior in nonlinear filtering performance, especially under a low sampled circumstance that is highly concerned for engineering applications.

#### *3.3. Adaptive Strategy for Optimized Implementation*

As mentioned in Equation (10), it is a multi-parameter optimization problem, where the constraint of the matched potential relationship is a nonlinearity that can properly be solved by intelligent optimization algorithms such as genetic algorithm (GA), particle swarm optimization (PSO), ant colony optimization (ACO), grey wolf optimization (GWO), et al. The main problem is to determine the proper searching interval of each parameter.

To address this problem, it can be known from [24,37] that in matching high frequency, signals require extremely large system parameters. As for practical engineering fields, the frequency of received signals generally varies from tens to thousands of Hertz, and hence, the order of magnitude for a matched potential parameter should be extremely large, especially for parameter *b* as there is a relationship with *b* ∝ *a*4. This problem is the same to the barrier factor *K* as there is an experimental guidance that *K*∗ *opt* approximate to *f* <sup>4</sup> <sup>0</sup> [37]. This refers to a large range of parameter searching interval to limit the computational efficiency and parameterized realization. Such a problem will limit the utilization in a real-time online embedded system in confronting data overflow problem (such as commonly utilized 32 bits digital signal processor (DSP)). From this view, the frequency shift to low frequency should be a kind of efficient solution for this problem, while lack of flexibility in dealing with unknown frequency signal or multi-frequency signal. In the consideration of the degree of freedom, the potential constraint of a matched relationship can be relaxed by the new freedom of the barrier factor *K*. As a consequence, the matched potential relationship can be eased. From the aforementioned analysis, for this problem, one can use the biquadrate relationship between the potential parameter *a* and *b*. In this regard, we can assume a particular case here by easing the potential

parameter *aopt* = 1. As a result, the computation complexity is greatly deduced, and the extremely large parameter limitation is relaxed simultaneously. The simplified model can be expressed as,

$$\begin{array}{rcl} \left(h\_{opt}^{\*}, K\_{opt}^{\*}\right) &=& \underset{\gamma, h, K}{\arg\max} \text{SNRI} \\ \text{s.t.} & a\_{opt} = 1 \\ & b\_{opt} = 1/(4K\bullet) \\ & \gamma\_{opt} = (2\sqrt{2}\pi f\_{0}\mathbf{e})^{-1} \\ & h \in (0, 1], \ K \in (0, K\_{max}] \end{array} \tag{11}$$

where the *Kmax* can be further limitted, and the optimal damping factor *γopt* is underdamped as well. Here, we need to point out that in this simplification process, the insufficient time-scale matching constraint result in an insufficient relation of *γopt* = (2 <sup>√</sup>2*<sup>π</sup> <sup>f</sup>*0e)−1. Therefore, this constraint should be eased by a generalized underdamped condition to keep the degrees of freedom. Hence, the optimization model can be further expressed as

$$\begin{aligned} \left(\gamma\_{opt}^\*, h\_{opt}^\*, K\_{opt}^\*\right) &= \underset{\gamma, h, K}{\text{argmax}} \text{SNI} \\ \text{s.t.} &\quad a\_{opt} = 1 \\ &\quad b\_{opt} = 1/(4K\mathcal{D}) \\ &\quad \gamma \in (0, 1], \ h \in (0, 1], \ K \in (0, K\_{\text{max}}] \end{aligned} \tag{12}$$

#### *3.4. Implementation of Adaptive Intrawell Matched Stochastic Resonance*

In the last subsection, we have the adaptive strategy as an optimization problem with nonlinear constraints. Here, a signal processing strategy is proposed by jointly optimizing the parameters with a classical genetic algorithm (GA) method for the sake of global optimality. Consequently, the detailed implementation steps are summarized as follows,

(1) Signal pretreament. Data normalization and prewhitening are executed on the actual received noisy signals.

(2) Searching range initialization. Initialize the optimization searching range of parameters *γ* ∈ [0, 1], *h* ∈ [0.001, 1] and *K* ∈ [0.001, 200], respectively.

(3) Parameter optimization. Obtain the optimal parameters *γopt*, *hopt*, and *Kopt* according to the following objective function,

$$\left(\gamma\_{opt'}^\* h\_{opt'}^\* \kappa\_{opt}^\*\right)\_{\text{opt}} = \underset{\gamma, h, \mathcal{K}}{\text{arg}\max} \text{SNRI} \tag{13}$$

where SNRI corresponds to the fitness criteria of GA method. Here, the input–output SNR improvement (SNRI) is constructed by the superposition of global SNRI (SNRI*global*) and local SNRI (SNRI*local*) to better characterize the filtering performance both in the view of global and local conditions. The SNR calculation for the input and output time series can be defined as,

$$\text{SNR} = 10 \log\_{10} \frac{P\_s - P\_n}{P\_n} \tag{14}$$

in which *Ps* <sup>=</sup> <sup>∑</sup>*f*+Δ*Bs*/2 *<sup>i</sup>*=*f*−Δ*Bs*/2 *Si* represent the total power around the characteristic frequency *<sup>f</sup>*<sup>0</sup> within a small bandwidth of Δ*Bs*, and *Pn* = <sup>1</sup> <sup>Δ</sup>*<sup>B</sup>* (∑*f*+Δ*B*/2 *<sup>i</sup>*=*f*−Δ*B*/2 *Si* <sup>−</sup> *Ps*) represent the average power of background noise within Δ*B* bandwidth around the characteristic frequency *f*0. For the global SNR and the local SNR calculation, Δ*B* is setted by the full bandwidth and *fs*/100, respectively.

For a multi-harmonic signal circumstance, a modified measurement index of average SNRI is utlized, which can be defined as,

$$\text{SNRI} = \frac{1}{M} \sum\_{h=1}^{M} \text{SNRI}(h) \tag{15}$$

where *M* is the harmonic number of the received signal. The parameter optimization algorithm is summarized in Algorithm 1.

(4) Signal post-treatment. Output optimal waveform and the coresponding optimal fittest value.

#### **Algorithm 1 Parameter Optimization Algorithm.**

#### **Parameter Initialization:**

*Cr* = [*γstart* = 0, *hstart* = 0.001, *Kstart* = 0.001; *γend* = 1, *hend* = 1, *Kend* = 200]: the searching intervals; *Nc* = 3: number of chromosome;

*Np* = 100: Number of individuals in the population;

*χ* = 0.95: The fraction to be replaced by crossover in each iteration;

*μ* = 0.01: The mutation rate;

*M* = 10: The maximal iteration times;

*λstop* = 0: The threshold of stop condition.

**Initialize generation** 0:

*k*:=0;

*Pk*:=a population of *Np* randomly-generated individuals;

**Evaluate** *Pk***:**

```
Compute fitness criteria SNRI for each i ∈ Pk;
```

```
{
```
1: Compute the corresponded MSR output by fourth order Runge–Kutta (RK4) method according to


**do**

```
{
```
1: **Copy**: Select (1 − *χ*) ∗ *n* members of *Pk* and insert into *Pk*+1;

```
2: Crossover: Select χ ∗ n members, pair them up to produce offspring and insert the offspring into
  Pk+1;
```
3: **Mutate**: Select *μ* ∗ *n* members of *Pk*+1, and invert a randomly selected bit;

```
4: Evaluate Pk+1;
5: if Pk+1 − Pk ≤ λstop then break;
6: else
7: Increment: k:=k+1;
8: end if
```

```
}
```

```
while k -
        M;
```
**return** the optimal fittest individual from *PM*;

#### **4. Filtering Performance Analysis and Evaluation**

SR can be regarded as a specific kind of nonlinear filter, while generally requires a high sampling rate condition. As is noticed in the Section 3.2, the proposed method can release the limitation on system response so as to allow us a higher input frequency in breaking through the large sampling rate limitation. Hence, here we intuitively exhibit the filtering performance referred to the discrete line signature and harmonic related line signature under a low sampling rate condition.

#### *4.1. Discrete Line Signature Signal Analysis*

Without loss of generality of detection, we consider two circumstances corresponding to *H*<sup>1</sup> and *H*<sup>0</sup> hypotheses, respectively. The sampling frequency is set to *fs* = 100 Hz, and the data length *N* = 2048.

For *H*<sup>1</sup> hypothesis circumstance, the input is composed of a narrowband component with the frequency *f*<sup>0</sup> = 10 Hz, and the ambient noise that simulated by the Gaussian noise. The input SNR is set as −15 dB, and the normalized waveform and power spectrum are shown in Figure 2a. First, we employed the lofargram to qualitatively evaluate the performance. A 2.56 s short-time Fourier transformation (STFT) with the overlapping length of 0.1 s was used to generate the lofargrams. The lofargrams of the original input is illustrated in Figure 2b, where we can see that the narrowband component is not easy to be identified. By utilizing the proposed AIMSR method, the optimal output waveform can be seen in Figure 2c, where it reflects an intrawell response with large amplitude. After a response delay, it is stable and periodic at 10 Hz. The lofargram of the AIMSR output is demonstrated in Figure 2d. It clearly identifies the narrowband component *f*0. It is worth noting that there is a second-order harmonic with periodic driving. This can be attribiuted to the resonance phenomenon and anticipate benefiting harmonic signature extraction. A comparison of normalized power spectrum density (PSD) of input and output (Note a highpass filter is adopted with passband frequency 1 Hz to deal with the DC bias, and this is the same in the rest of the paper) with the Welch method is given in Figure 2f, where we can see the superior filtering performance. The frequency response of the AIMSR background noise has the form of Lorentzian distribution by concentrating most of the noise energy into the low-frequency region. This refers to the special property of nonlinear SR phenomenon. The optimal SNRI obtained at each iteration is shown in Figure 2e. It can be seen after three iterations that an optimal result is obtained, which reflects a convergence of the proposed algorithm.

For *H*<sup>0</sup> hypothesis circumstance, the input is simulated by the pure white Gaussian noise. As is assumed by prior unknown, we conduct the optimization with *f*<sup>0</sup> = 10 Hz for parameters initialization as well. The corresponding results are demonstrated in Figure 3a–f. By utilizing the proposed AIMSR method, the optimal output waveform is shown in Figure 3c, where it reflects an intrawell response as well. After a response delay, the amplitude looks more fluctuant compared with *H*<sup>1</sup> hypothesis circumstance as illustrated in Figure 2c. The lofargram of the AIMSR output is demonstrated in Figure 3d, where there seems to be a narrowband component near 10 Hz. To further compare the normalized power spectrum density (PSD) of the input and output in Figure 3f, we can see a frequency bias with Δ*f* that does not match the initialized frequency *f*0. Such a circumstance does not have a high-order harmonic that is regarded as a non-SR occurrence. The optimal SNRI obtained at each iteration is shown in Figure 3e, where the optimal result is obtained after two iterations. The optimal fitness value is larger than *H*<sup>1</sup> hypothesis circumstance. This should be affected by the calculation result of input as the input SNR is extremely smaller.

In summary, it is clear to see the superior filtering performance of AIMSR with a limited sampling rate condition. The difference in response between the *H*<sup>1</sup> and *H*<sup>0</sup> hypotheses can be identified, which indicates the detection performance fairly well.

**Figure 2.** A simulation of discrete line signature signal on *H*<sup>1</sup> assumption: (**a**) the input waveform and its normalized power spectrum; (**b**) lofargram obtained with the input; (**c**) the adaptive intrawell matched stochastic resonance (AIMSR) output waveform (*γopt* = 0.0322, *Kopt* = 79.05, *hopt* = 0.4413); (**d**) lofargram obtained with the AIMSR output; (**e**) the optimal fitness at each iteration with genetic algorithm (GA) method; (**f**) filtering performance comparison with normalized power spectrum density (PSD, Welch method).

**Figure 3.** A simulation of discrete line signature signal on *H*<sup>0</sup> assumption: (**a**) the input waveform and its normalized power spectrum; (**b**) lofargram obtained with the input; (**c**) the AIMSR output waveform (*γopt* = 0.2875, *Kopt* = 109.4507, *hopt* = 0.4465); (**d**) lofargram obtained with the AIMSR output; (**e**) the optimal fitness at each iteration with GA; (**f**) filtering performance comparison with normalized PSD (Welch method).

#### *4.2. Harmonic Related Line Signature Signal Analysis*

Harmonic vibration is a typical vibration generated by rotatory machinery. As is mentioned in Section 2, the harmonic-related signature is commonly referred to in a ship engine and propeller. The sampling frequency is set to be *fs* = 100 Hz, and the data length *N* = 2048. We consider two circumstances corresponding to *H*<sup>1</sup> and *H*<sup>0</sup> hypotheses as well.

For *H*<sup>1</sup> hypothesis circumstance, the tested harmonic signal is a combination of a fundamental sinusoid and its two high-order harmonics, whose frequencies are all harmonically related to two times and four times the integer multiples of the fundamental frequency *f*0. These three sinusoidal signals have the same amplitude of 0.1, and the same Gaussian noise background. The input SNR is set as −15 dB. The lofargram of the original input is illustrated in Figure 4b, where we can see the three harmonic-related narrowband components. By utilizing the proposed AIMSR method, the optimal output waveform still responds well with a response delay as shown in Figure 4c. In comparing with the single periodic driving, the amplitude of stable response is smaller. The corresponding lofargram of the AIMSR output is demonstrated in Figure 4d, where the three harmonic-related narrowband components can be identified. To further evaluate the filtering performance with normalized power spectrum density (PSD) of input and output in Figure 4f, we see the great local denoising performance with the fundamental frequency *f*0, which is gradually lost with the 4th order harmonic (4 *f*0). The optimal SNRI obtained after two iterations as shown in Figure 4e, where we can find the optimal value is so small compared with the discrete line signature signal circumstance. There is a negative improvement of global SNRI for the high order harmonics. In this way, the local SNRI is recommended to be the evaluation index of high-order harmonics.

For the *H*<sup>0</sup> hypothesis, the input is simulated by the pure white Gaussian noise. As is assumed, we conducted the optimization with *f*<sup>0</sup> = 10 Hz for parameter initialization as well. The corresponding results are demonstrated in Figure 5a–f. By utilizing the proposed AIMSR method, the optimal output waveform is shown in Figure 3c. The output is an intrawell response with a larger response delay, and the amplitude of stable response reflects larger effects. The lofargram of the AIMSR output is demonstrated in Figure 5d, where we can see a broadband harmonic energy in the first 5 seconds, and then the energy is focused to a periodic mode with a frequency bias Δ*f* that does not match the initialized frequency *f*<sup>0</sup> as shown in Figure 5f. This means no SR occurs as the noisy input can not be matched to the nonlinear system. The optimal fitness value is extremely large which is consistent with the *H*<sup>0</sup> hypothesis of the discrete line signature signal case. This indicates the detectability as well.

In summary, the proposed AIMSR method can be utilized to enhance the harmonic-related signature, especially for the fundamental frequency. Since the fundamental frequency estimation is a topic that spans many disciplines including passive sonars [39], speech recognition [46], biomedical signal processing [47], musical pitch estimation [48], and etc., this method is anticipated to be a potential new technique for the future. The difference in responses between the *H*<sup>1</sup> and *H*<sup>0</sup> hypotheses can be identified with proper measurement index to establish a detector. The computation is efficient and generally converges to an optimum within five iterations, which can be realized in the embedded system. This work will be a topic for further study in the future.

**Figure 4.** *Cont*.

**Figure 4.** A simulation of harmonic-related line signature signals on *H*<sup>1</sup> assumption: (**a**) the input waveform and its normalized power spectrum; (**b**) lofargram obtained with the input; (**c**) the AIMSR output waveform (*γopt* = 0.1132, *Kopt* = 70.4507, *hopt* = 0.4429); (**d**) lofargram obtained with the AIMSR output; (**e**) the optimal fitness at each iteration with GA; (**f**) filtering performance comparison with normalized PSD (Welch method).

**Figure 5.** *Cont*.

**Figure 5.** A simulation of a harmonic-related line signature signal on *H*<sup>1</sup> assumption: (**a**) the input waveform and its normalized power spectrum; (**b**) lofargram obtained with the input; (**c**) the AIMSR output waveform (*γopt* = 0.0462, *Kopt* = 159.2903, *hopt* = 0.5763); (**d**) lofargram obtained with the AIMSR output; (**e**) the optimal fitness at each iteration with GA; (**f**) filtering performance comparison with normalized PSD (Welch method).

#### **5. Application Verification and Discussion**

#### *5.1. Verfication on AUV's Low Frequency Propeller Harmonic Tonals*

To validate the feasibility and efficiency of the proposed method, an experiment was conducted with an AUV in the Zhanghe reservoir. The AUV moved away in a trajectory of straight line, and the radiated noise was measured by a low frequency seismometer with the sampling frequency *fs* = 200 Hz. Figure 6a shows the lofargram of the received signal, where we can see three harmonic lines. The conventional ALE that employs the least mean square (LMS) algorithm is conducted to enhance the line signatures as shown in Figure 6b, where the harmonic lines are more clearly identified. By utilizing the proposed AIMSR method, the optimal output is the intrawell as well. The corresponding lofagrams of the AIMSR output and the direct current (DC) offset filtered output are illustrated in Figure 6c,d, respectively. It can be seen that the background noise of the AIMSR output is much lower than that of the conventional ALE. A comparison of normalized power spectrum denisity (PSD) is shown in Figure 6f, where there is about 10 dB improvement for the fundamental frequency compared to that of the conventional ALE. For high order harmonics, the conventional ALE should be better. This indicates that a combination of the two methods gives a better output.

**Figure 6.** Filtering performance verification: (**a**) lofagram of the original received signal of an autonomous underwater vehicle (AUV); (**b**) lofagram of the adaptive line enhancer (ALE)output (LMS); (**c**) lofagram of the AMSR output (*γopt* = 0.1227, *Kopt* = 85.6350, *hopt* = 0.4707); (**d**) lofagram of the AMSR output with direct current (DC) offset filter; (**e**) the optimal fitness at each iteration with GA; (**f**) normalized power spectrum denisty (PSD) comparison.

#### *5.2. Discussion*

(1) From the aforementioned numerical and practical analyses on AIMSR with potential constrainted bistable model, the SR is implemented under a low sampling rate condition. Essentially, it is to add a new degree of freedom to the bistable potential which can have an effect on the limitation of the system response. As we consider the system parameter tuning SR as a special nonlinear filter, it is a tenet that the more the tuning parameters (system complexity), the better output response. Similar results can also be found like multi-parameters tuning [49], improved potentials [30,50,51], coupled systems [31,32], etc. Hence, one of the contribution of this paper is to explain and lead the further explore of nonlinear filter. In addition, from the view of system nonlinearity, there seems to be a connection to deep neural networks for the system complexity. This might lead to an inspiration and possibility in training the deep SR network for applications, and gives us a guidance to better understand the innate character of deep learning from the view of nonlinearity as well.

(2) The evaluation index of SNRI in this paper generally requires a prior knowledge of signal frequency. However, for engineering applications, it is commonly unknown. As is known, the SR will always force the energy focus to a periodic mode. It is hard to distinguish whether the correctness of output response is simply achieved by frequency searching. Although the difference of response between the *H*<sup>1</sup> and *H*<sup>0</sup> hypotheses can be fairily well identified, the problem of proper selection of the measurement index in the sense of detection is still a problem for further studies.

(3) The superior filtering performance is validated in this paper, which means the advantage of AIMSR is a potential. However, its intrinsic Lorentzian property leads to a preference in enhancing the low-frequency signals, and as a result, to a loss of the performance for higher frequencies. As is mentioned in the last subsection, a conventional ALE should be better in dealing with the high-order harmonics. Proper combination on these two methods is needed for better outputs. This work will probably be studied in the future.

(4) The proposed AIMSR method is of low computational complexity which enables implementation in an embedded system. This is an interesting and highly important engineering problem. However, related work can rarely be found. To further improve the applicablity of the SR in the engineering fields, we think this work is necessary for future studies.

#### **6. Conclusions**

In this paper, we proposed an adaptive intrawell matched stochastic resonance (AIMSR) method to break through the limitation of the conventional ALE by a nonlinear filtering effect. The problem of parameterized implementation of SR under a low sampling rate condition is firstly addressed by implementing a framework of intrawell matched stochastic resonance with potential constraint. To further promote its practicality in an embedded system, the large parameter limitation of matched relationship is eased to deduce the computation complexity. Simulation analyses are conducted with a discrete line signature and harmonic-related line signature that reflect the superior filtering performance as well as the computational efficiency fairy well. Besides, two hypotheses corresponding to *H*<sup>1</sup> (periodic signal with noise) and *H*<sup>0</sup> (pure noise) circumstances are further considered to reveal the feasibility of detection. Application verification was experimentally conducted in a reservoir with an AUV to validate the proposed AIMSR method compared with the conventional ALE method. Additional intensive discussions have been made to give an insight into the principal results and inspire future investigations. In the end, we anticipate that the proposed method can be a potential new technique for passive sonar detection in the future.

**Author Contributions:** Conceptualization, H.D.; Formal analysis, H.D. and S.M.; Funding acquisition, X.S. and H.W.; Investigation, S.M. and C.Q.; Resources, C.Q.; Methodology, H.D. and K.H.; Validation, H.D., X.S. and S.M.; Writing—original draft, H.D.; Writing—review & editing, K.H., X.S. and H.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Natural Science Foundation of China (Grant No. 61571365, 61671386, 61771401, 61571367) and National Key Research and Development Program of China (Grant No. 2016YFC1400200

**Conflicts of Interest:** The authors declare no conflict of interest. The funding sponsors have the role in the design of the study and 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/).
