doi:10.1155/2007/71528 Research Article Underwater Noise Modeling and Direction-Finding Based on Heteroscedastic Time Series

We propose a new method for practical non-Gaussian and nonstationary underwater noise modeling. This model is very useful for passive sonar in shallow waters. In this application, measurement of additive noise in natural environment and exhibits shows that noise can sometimes be significantly non-Gaussian and a time-varying feature especially in the variance. Therefore, signal processing algorithms such as direction-finding that is optimized for Gaussian noise may degrade significantly in this environment. Generalized autoregressive conditional heteroscedasticity (GARCH) models are suitable for heavy tailed PDFs and time-varying variances of stochastic process. We use a more realistic GARCH-based noise model in the maximum-likelihood approach for the estimation of direction-of-arrivals (DOAs) of impinging sources onto a linear array, and demonstrate using measured noise that this approach is feasible for the additive noise and direction finding in an underwater environment.


INTRODUCTION
A passive sonar generally employs array processing techniques to resolve problems such as localization of targets [1,2]. As a matter of fact, all the DOA estimation methods make a crucial assumption for the noise model, that have a great impact on the performance of DOA estimation. In the underwater environment, the measurements of additive noise show that we have non-Gaussian process [3][4][5]. Natural and manmade sources such as reverberation and industrial noise that cause additive noise distribution exhibit performances far away from the Gaussian model. These factors are more in coastal and shallow waters. Thus, the algorithms that are optimized for Gaussian distribution will degrade in actual experiments. All this mentioned factors give a stochastic and time-varying nature to the background noise. Thus, a proper model presentation which could best and simply describe the different features of the realistic background noise affecting the desired signal is an important part of a sonar signal processing. In the last decade, after the seminal works by Engle [6] and Bullerslev [7] there has been a growing interest in time series modeling of changing variance or heteroscedasticity. These models have found a great number of applications in nonstationary time series such as financial records. Generalized autoregressive conditional heteroscedasticity; for example, GARCH [7], is a time series modeling technique that uses past variances and the past variance forecasts to forecast future variances. GARCH models account for two main characteristics: excess kurtosis; that is, heavy tailed probability distribution, and the volatility clustering; that is, large changes tend to follow large changes and small changes tend to follow small ones, compatible to a large extent to the additive noises in a natural environment. We suggested this more realistic dynamic model for additive noise modeling in array signal processing [8]. Now, we offer this model for the underwater noise in passive sonar due to the facts that the commonly used model for environmental additive noise exhibits heavier tail than the standard normal distribution [9], and the conditional heteroscedasticity suggests a time series model in which time-varying variances are presented, that is, a more logical modeling for the dynamic of the additive noise [7]. Hence, in this paper, we propose to assume a conditional heteroscedasticity-based time series for underwater noise modeling and that can be used in the direction-finding approach for passive sonar. This paper is organized as follows. In Section 2 we present the GARCH time series. The proposed noise modeling as the underwater noise and the DOA estimation based on the new noise model is provided in Section 3 and the simulation results of the proposed method come in Section 4. Some concluding remarks are provided at the end.

GARCH TIMES SERIES
The exploitation of time series properties has been extensively used in signal modeling and parameter estimation. For example, ARMA time series models have wide applications in signal processing such as sonar signal processing and noise modeling [10,11]. One of them that has been used in the past decade, conditional heteroscedasticity time series, was first introduced by Engle [6] in the context of modeling United Kingdom inflation as known autoregressive conditional heteroscedasticity (ARCH). Such models are characterized by being conditionally Gaussian, additionally represented by a nonconstant and state-dependent variance. However, in [6,7,12,13], it is shown that a time-varied variance over time is more useful than a constant one for modeling non-Gaussian and non-stationary phenomena such as economic series. Generalization of ARCH that is proposed in [7] is called GARCH. Generally speaking, in heteroscedasticity we consider time series with time-varying variance; the conditional implies a dependence on the observation of the immediate past, and autoregressive describes a feedback mechanism that incorporates past observations into the present. GARCH then is a mechanism that includes past variance in the explanation of the future variance. GARCH models account for heavy tailed PDF as excess kurtosis and volatility clustering a type of heteroscedasticity. Now, we let (k) denote a real-valued discrete-time stochastic process, the GARCH (p, q) process is then given by [7] (k) ∼ G 0, σ 2 (k) , where G denotes the Gaussian probability density function and α 2 0 , α 2 i , and β 2 i are GARCH model coefficients. For example, Figure 1 shows some realizations of the GARCH(1, 1) with different coefficients. The flexibility of GARCH process is displayed in this figure, so that some different time series such as impulsive data can be modeled. This ability can be obtained due to the complex coefficients structure of GARCH modeling. The estimation of orders p and q has an important role in the GARCH modeling of the time series. Otherwise, because of the importance of the orders in the computation of coefficients, we should have the proper consideration for it. In this way, some methods are proposed such as likelihood ratio tests [14], Akaike information theory criterion (AIC), and Bayesian information criteria (BIC) [15]. The likelihood ratio tests would be used to determine supporting the use of a specific GARCH model for a time series. Using the following order selection methods, AIC, BIC, and likelihood ratio test, the order that provides the best model for data fitting is selected.
In this way, we use AIC and BIC information criteria to compare alternative models such as GARCH(1, 1), GARCH(2, 1) and others. Since information criteria penalize models with additional parameters, the AIC and BIC model order selection criteria are based on parsimony [15].
The AIC and BIC statistics are defined as where L g is optimized log-likelihood objective function values associated with parameter estimates of the GARCH models to be tested so that the following is obtained: N p is the number of GARCH parameters and K is the number of observations. In the following section, we consider GARCH-based model for the underwater noise modeling and DOA estimation in a passive sonar. Figure 2 shows a general block-diagram of a passive system such as sonar so that it has as the input process (propagated source) the underwater channel, the additive noise, and the observed data in receivers. In this way, we consider the additive noise comprised of the additive noise and interferences as follows:

Underwater noise modeling
where n(k) is the received additive noise and interference in time, n P (k) is the interference part, n G (k) is the additive Gaussian noise part, and k stands for the snapshot index. Due to natural and manmade sources in the underwater environment, the measurements of noise shows that we have non-Gaussian process [3][4][5]. These factors such as reverberation and industrial noise are more in shallow waters. However, in practice the noise model is not known because of the time-varying characteristic of system and non-Gaussian behavior of noise source. These are two major factors that can limit the performance of general methods in the practical experiments. In different applications such as sonar, the timevarying characteristic is generally due to time-varying nature of the medium channel, environment, noise, and interferences [16][17][18][19]. For example, underwater acoustic channel is a time-varying and multipath channel specially in shallow water. It varies due to different season, area, and situation of sea face. The channel variations can be due to the spatial movement of the source and/or changes in the propagation conditions such as sound speed profile. All this mentioned factors give a stochastic and time-varying nature to the background noise. Hence, we accept a model in which some kind of changing variance in time is included. As a result of the above time-varying events, it can be assumed that the additive noise has time-varying variance in the receiver. Moreover, measurements of the additive noise in related application such as underwater environment shows that the noise can sometimes be significantly non-Gaussian and it can be shown for the widely accepted model of additive noise and interference excess kurtosis can be observed [5,9]. For this purpose, the narrowband Middleton class-A model [20] is used. This model is a general physical and statistical model for received additive noise and interference. With (A.5) in Appendix A, the excess kurtosis is determined as shown in Figure 3. In this figure the excess kurtosis is shown for the Middleton class-A model for general values of model parameters. Thus, the assumed noise model that covers the properties of additive noise such as time-varying variance and heavy-tail PDF is more attractive. Under the above assumptions and important features of the GARCH time series model, we use this model for the additive noise modeling in the underwater acoustics applications such as sonar: where k = 1, 2, . . . , K and K is the number of snapshots. At the start of the modeling technique, we need to the estimation of orders of proposed model, that is, p and q. In this way, we used both AIC and BIC and so the results of our simulation almost always reached GARCH(1, 1) using recorded noise. The results of this approach are given in the simulation section. Therefore, Generally, the unknown coefficients (α 2 0 , α 2 1 , and β 2 1 ) are estimated using maximum-likelihood method [7]. This model exploits time-varying variance and heavy-tail PDF for approaching the realistic properties of the additive noise in practical applications. It is well known that kurtosis is an important parameter for analysis of non-Gaussian random processes. For this purpose, the kurtosis is given by [7] It can be shown that if (α 2 1 + β 2 1 ) < 1 and 1 − (α 2 1 + β 2 1 ) 2 − 2α 4 1 > 0, then the kurtosis is greater than 3 and GARCH(1, 1) can conclude heavy-tail PDF. Figure 4 shows the ability of GARCH modeling for heavy-tail PDF with excess kurtosis.
As it is well known, the assumed noise models can have an important roles in the signal processing methods such as parameters estimation. In the following, we propose a new DOA estimation approach using GARCH noise modeling. This approach could encompass the DOA estimation not only in non-Gaussian environment but also it could handel heavy tailed and nonstationary processes.

Direction-finding approach
A proper model presentation which could best and simply describe the different features of the realistic background noise affecting the desired signal is an important part of a sonar signal processing, and so the algorithms that are optimized for Gaussian distribution degrade in actual experiments. The performance of the source localization and estimation of DOA in passive array applications such as sonar heavily rely upon the particular array signal processing algorithms used in practice. In these methods, the key assumption is the noise model; that is, additive noise covariance, that is used in estimation of unknown parameters. Generally, additive noise is assumed to follow a Gaussian distribution. However, measurements of acoustic noise and interference in underwater environment show that the noise can sometimes be significantly non-Gaussian [3][4][5]9]. Consequently, we note that in this model, noise is not uniform across L sensors, which is a realistic modeling resting on the assumption of non-uniformity [21,22] and non-stationarity; that is, time-varying variance. Thus, the assumed noise model that covers the properties of background noise is more attractive. Under the above assumption we use the GARCH(1, 1) process for the additive noise in direction finding in array signal processing. Let us assume that a linear array of L omnidirectional hydrophones receives D (D < L) plane wave from unknown directions of arrivals. The incident plane waves are assumed to be narrowband with a center frequency. Under these conditions, the kth snapshot vector of array observation can be expressed as where s(k) is the D × 1 vector of the source waveforms, n(k) is the L × 1 vector of sensor noise, A(θ) is the L × D steering matrix a(θ i ) is the direction vectors, θ {θ 1 , . . . , θ D } T is the D × 1 vector of the unknown signal DOA, K is the number of snapshots, and (·) T stands for the transpose operation. We make the following assumptions: the signal waveforms are stationary; both temporally and spatially, and the signals and noise are statistically independent of each other. According to the previous noise modeling section, we propose using the multivariate GARCH(1, 1) for noise modeling in array sensors applications such as sonar. Thus, using (6) and (7) the additive array noise can follow as multivariate GARCH(1, 1) with zero mean and covariance matrix Q(k), so n(k) ∼ MG 1,1 n; 0, Q(k) , where MG 1,1 stands for the multivariate GARCH(1, 1), and In this approach, the additive noise model at every sensor is distributed similar to (6) and (7). Same as (9), another application of GARCH model can be seen in [23] that is used in the adaptive portfolio management based on maximum likelihood in state space method. Consequently (it is well known) one of the efficient methods in the estimation of parameters in array signal processing is the ML approach [11,24,25]. For this method, the key assumption is the noise model; that is, additive noise covariance that is used in Hadi Amiri et al.

5
the estimation of unknown parameters. In the following, we exploit the deterministic maximum-likelihood (DML) approach model so that the signal waveforms are deterministic unknown sequences. Thus, the joint PDF of the observed array snapshots using GARCH(1, 1) model is expressed as where where g is vector of GARCH(1, 1) coefficients, is the vector of the unknown signal. Therefore, by using (12) and (13) it can be shown that the following holds for loglikelihood: L p (·) stands for the proposed log-likelihood function to be maximized over the vector of unknown parameters ψ through ML approach. Therefore, due to complicated nature of problems, this estimation cannot be found analytically, and −L p (·) can be minimized through numerical procedures [1] and then unknown parameters are found. In this way, we use the gradient-based minimization, the Newton approach. These methods are based on multidimensional searching and minimization of log-likelihood function onto parameter space and have heavy computational burden. Generally, we would like to decrease the number of unknown parameters as much as possible, and also select their feasible initial values of them at the beginning of the process. As a matter of fact, the initial values of the bounded parameters are one of the most important factors in the rate of the convergence and the computational volume of the minimization process.
In the proposed method we have different unknown parameters in the process. The most important is the DOA of sources so that they are estimated in our approach. Without loss of generality, we assume that the number of sources is given and then we use the popular method such as music to have initial values of DOAs at the beginning of approach. About noise model parameters, we considered some constraints on their values such as α 2 1 ≥ 0, β 2 1 ≥ 0, and (α 2 1 + β 2 1 ) < 1. Continuously, if we do not have real signal waveforms, then, the known least square error approach is used for initial estimation [1]. For the statistical analysis of proposed method, the CRB is derived in Appendix B.

SIMULATION AND RESULTS
In this section, we demonstrate the performance of the proposed approach for modeling of the additive noise in passive sonar with two major experiments. In the first experiment, we use the recorded noise with one hydrophone in shallow water. In this scenario, order selection and estimation of PDFs of the real and simulated data are considered. Using GARCH noise modeling in the DOA estimation of the underwater targets are examined in the latter experiment that utilize uniform linear array (ULA). Root-mean-squareerrors (RMSE) of estimated DOA are considered versus SNRs and the number of snapshots.

Single hydrophone
For the performance analysis of the ability of the proposed model, we utilize the underwater additive noise that is recorded in the shallow water. Before modeling process, we exploit the available approaches [14,15] for the estimation of GARCH orders p and q and so find that p = 1 and q = 1 are sufficient orders for this experiment. A typical results of AIC and BIC are shown in Table 1 based upon underwater measured noise. We use (2) on the recorded data and conclude that GARCH(1, 1) is a feasible model in this application. After this model order selection, we can simulate the data with GARCH(1, 1) using log-likelihood approach with (3), (6), and (7). Figure 5 shows one of the time series of the measured noise and simulated noise with GARCH model. For the statistical comparison of proposed model, PDF is estimated for the real, Gaussian, and GARCH simulated noises. The results of estimation of PDFs are shown in Figure 6 that can verify the flexibility of GARCH process for the additive noise modeling.

Hydrophone array
After the analysis of the proposed method for real data modeling, we assume that passive sonar has a uniform linear array (ULA) with half-wavelength inter-element spacing, and equally powered narrowband sources with DOA = [5 • , 10 • ] relative to the broadside. This array has six omnidirectional sensors. The experiments consist of Monte Carlo trails, a total of 50 trails are run. In all examples, the DOA estimation RMSE of the proposed method have been compared with derived CRBs. We conduct the experiments to show the performance of our proposed method with respect to SNR and the number of snapshots. In our experimental results music and DML results are also compared against the proposed method, GARCH-ML. In this scenario, we use the 6 EURASIP Journal on Advances in Signal Processing underwater noise for the performance analysis of the proposed method. This data is collected in the shallow waters of Persian Gulf and includes the background noise. The detail of the experiment is given in Table 2. RMSE and CRB for this scenario are shown in Figures 7 and 8. In this experiment, we see that GARCH(1, 1) is an appropriate choice for the modeling of the underwater noise, and observe that the proposed method has resolved the targets better than the other methods, and the RMSEs of the proposed method are less than the others and asymptotically approaching the CRB limit.

CONCLUSION
In this paper we propose a new method for the underwater noise modeling and DOA estimation in passive sonar signal processing. We utilized GARCH(1, 1) noise modeling in the ML approach to estimate DOAs of sources. This model accounts for heavy tailed PDFs with excess kurtosis and timevarying variance of a type of heteroscedasticity. For evaluation of the proposed method, two experiments, namely, univariate and multivariate measured underwater noise, are used. We also computed CRB for studying the statistical performance of the proposed method. The results of these simulations verify that the proposed method is suitable for the noise modeling in the realistic underwater acoustic environment and so for the direction-finding approach in a passive sonar.

A. KURTOSIS IN MIDDLETON CLASS-A
It can be shown for the widely accepted model of additive noise, and interference excess kurtosis can be observed. For this purpose, the narrowband Middleton class-A model [20] is used. This model is a general physical and statistical model  for received additive noise and interference where n(t) is the received additive noise and interference in time, n P (t) is the interference part, and n G (t) is the additive Gaussian noise part. Due to Middleton class-A model, the characteristic function for the data is as follows [20]: and the PDF after normalization: where A is the overlap index that is a measure of "non-Gaussiannity," σ 2 m (m/A + Γ )/(1 + Γ ), Γ is the Gaussian factor, c 2 m = σ 2 G (m/K + 1), and K AΓ . Using characteristic function, the moments of n(t) are computed, and then Hence, the kurtosis is acquired using the following equation:

B. CRAMÉR-RAO BOUND
In order to understand the performance of estimation process using GARCH modeling we develop CRB [1]. If we denote the covariance matrix of the estimation errors by C(ψ), then the multiple-parameter CRB states that for any unbiased estimate of ψ. The J matrix is commonly referred to as Fisher's information matrix with the following elements: For kth single snapshot problem, the J is obtained from where the superscripts "R" and "I" denote the real and imaginary parts. In our DOA estimation method, we have J as a partitioned matrix: For the estimation of CRB, all of the blocks of the above matrixes should be computed, so that in the following these are given: Next J s s is a block diagonal matrix (2KD × 2KD), and each block is 2D × 2D and so has an identical structure. The kth block corresponds to kth snapshot. It has the following structure: , i = j, , i = j, A * θ i d θ j s j (k) σ 2 (k) ,