Localization of autonomous underwater vehicles (AUVs) is a very important and challenging task for the AUVs applications. In long baseline underwater acoustic localization networks, the accuracy of single-way range measurements is the key factor for the precision of localization of AUVs, whether it is based on the way of time of arrival (TOA), time difference of arrival (TDOA), or angle of arrival (AOA). The single-way range measurements do not depend on water quality and can be taken from long distances; however, there are some limitations which exist in these measurements, such as the disturbance of the unknown current velocity and the outliers caused by sensors and errors of algorithm. To deal with these problems, an AUV self-localization algorithm based on particle swarm optimization (PSO) of outliers elimination is proposed, which improves the performance of angle of arrival (AOA) localization algorithm by taking account of effects of the current on the positioning accuracy and eliminating possible outliers during the localization process. Some simulation experiments are carried out to illustrate the performance of the proposed method compared with another localization algorithm.

1. Introduction

In a three-dimensional underwater environment, the problem of localization of autonomous underwater vehicles (AUVs) has been widely investigated in recent years, as they make it accessible to those untouchable areas for human beings and assist with complex and arduous underwater tasks, which has important theory and application value in various robotic applications, such as the underwater target detection [1], the underwater target tracking [2, 3], and underwater search and rescue missions [4, 5]. Localization is one of the key components that enable the autonomy of AUVs. AUVs need accurate localization for the accuracy of the gathered data. For a mobile robot to be truly autonomous, it needs to be able to operate and navigate without human intervention and in a nonspecially engineered environment [6, 7]. And thanks to the uncertainty and complexity of the underwater environment and technical restrictions, it is still challenging and needs more complementary researches that localize a moving AUV. Indeed, the research of robotic navigation in two-dimensional space has achieved great results and extended to the study of multirobot systems [8], while during underwater navigation, AUVs cannot rely on Global Navigation Satellite Systems (GNSSs) due to the attenuation of electromagnetic radiation in the water domain, and in the absence of specific positioning systems, they can exclusively relay on dead-reckoning techniques. As the latter approaches integrate noisy and biased measurements from Inertial Measurement Units (IMUs) and velocity sensors, they suffer from numerical drift that makes them usable only for relatively short periods [9]. And the inertial navigation system (INS) is also mostly applied on AUVs, but INS usually suffers from error accumulation.

Commercially available underwater positioning solutions are mainly based on acoustic devices that, through the measurement of the time of flight of acoustic signals, allow us to measures the ranges from source to receiving nodes. Acoustic localization system for AUVs mainly includes long baseline (LBL), short baseline (SBL), and ultra-short baseline (USBL). More recently, research efforts focused on position systems based on the use of range measurements to a single node with the aim of developing the solutions which is simple, cheap, and easy to operate. Such approaches, known in the literature as single beacon localization [9], single range localization [10], or range-only localization [11], are based on the fusion of range measurements to the single source with information from AUV’s onboard sensors as Inertial Measurement Unit (IMU), Doppler Velocity Logger (DVL), and depth sensors.

However, time-synchronization is hard to achieve in underwater environment, which baffles the accuracy of ranging and leads to inaccuracy of localization. Sound speed uncertainty enlarges the inaccuracy of distance estimation. Long propagation delay of acoustic signals, the impact of AUV mobility, and the drift term caused by the water current are the urgent problems to be solved, which will affect the accuracy of underwater acoustic positioning. Moreover, the severe multipath property in underwater environment caused by the scattering from the seabed and water surface is also needed to be overcome. In LBL, AUV communicates with one or more transponders fixed on the seabed or the surface of the water to measure distance and achieve underwater acoustic localization; it can also be classified as time of arrival (TOA), time difference of arrival (TDOA), and angle of arrival (AOA) according to the specific positioning mode. TOA needs strict time-synchronization, which is energy-intensive and hard to achieve. Multiple-node cooperation in TDOA [12, 13] increases cost and is computationally complex. Moreover, the LBL multiple-array acoustic localization [14] based on multiple-sensor data fusion mostly adopts the above two or more methods, so as to acquire more accurate location estimation. However, the positioning algorithms based on data fusion are relatively complicated and need to take much time to calculate and fuse the results of different positioning methods; thus, practical application is not desirable in view of low-cost requirements.

As mentioned above, various improvements have been proposed to deal with the problems about underwater localization of AUVs. For example, Tan and Li [15] proposed a centralized algorithm to overcome the severe multipath property of the underwater environment. Cheng et al. [16] presented a TDOA-based localization scheme for stationary underwater acoustic sensor networks, which does not require time-synchronization among network nodes. Gao et al. [17] proposed the round-trip ranging (RTR) technology frustrated by the long propagation delay, the drift term of the water current, multipath property of the underwater environment, and AUVs’ mobility. And Li et al. [18] presented a self-localization algorithm with accurate sound travel time solution (SL-STTS) which solved the problem of time-synchronization and sound speed uncertainty; however, the drift term of the water current was ignored and the measurements of multiple sensors are not fully utilized.

In this paper, an improved single-way ranging localization approach based on PSO of outliers elimination is proposed. In the proposed approach, the problem of vehicle autonomous localization under complex underwater environment is fully considered and some solutions are presented in the following paragraphs, including the drift term caused by the current and the outliers from multiple-array acoustic localization caused by sensors. Finally, various experiments are conducted under different noise conditions for the vehicle autonomous localization. The results show the efficiency and affectivity of the proposed approach. The main contributions of this paper are summarized as follows. (1) The effect of the water current on the localization for AUV is fully taken into account and the relative analysis is carried out. (2) An outlier detection approach of multiple-array acoustic localization is proposed to eliminate the abnormal results of multiple localization combination.

This paper is organized as follows. Section 2 presents the improved SL-STTS-based method for the vehicle autonomous localization method. Section 3 presents a method for detecting and removing outliers from sensor in multiple-array acoustic localization. Performance evaluation of the proposed algorithms is presented in Section 4 through several simulations. Section 5 gives out some discussions on the robustness and fault tolerance of the proposed approach. Finally, the conclusions are given in Section 6.

2. Single-Way Ranging Localization

Considering the 3D AOA localization problem of the underwater vehicle in Figure 1, where the vehicle at an undetermined location is to be located by communicating with acoustic transponders at the fixed inertial location (where is the number of the acoustic transponders), using the azimuth and elevation angles. Here, at least three transponders are required for 3D AOA vehicle localization. The orientation of the vehicle includes yaw and pitch angle (ignore the posture of the vehicle, so skip the roll angle). With no loss of generality, the vehicle and transponders are assumed to be at the Cartesian coordinate system. The vehicle velocity is given by a superposition of a drift term and a controlled input term . For underwater vehicle, the velocity term can be generated through a guidance controller exploiting an onboard navigation sensor as a Doppler velocity logger (DVL). The drift term models a constant unknown water current. The vehicle is also equipped with axis gyro and inclinometer for measuring its yaw angle and pitch angle and AOA antennas for measuring AOA of received acoustic signals. Similarly, the transponders are equipped with AOA antennas for measuring AOA of received acoustic signals. Conductivity-temperature-depth (CTD) instruments are preloaded onto the vehicle and transponders to help estimate sound speed by measuring salinity, temperature, and depth [19].

2.1. Problem Formulation

As the vehicle transmits a sound signal at time instance , meanwhile, a timer starts. Assume that there are transponders detecting the sound signal from the vehicle. For transponder , after receiving the sound signal at its local time , it puts its number, the location, the sound speed measured by CTD [20], and the AOA of received acoustic signal into a localization packet . At its local time , it adds the leaving time into as soon as it responds to the AUV with . The vehicle receives at its local time . When the timer is out, the vehicle stops receiving replies from transponders and processes the received localization packets; then the vehicle at time instance is located according to the location packet. In order to overcome the severe multipath property in underwater environment, the AUV will no longer receive acoustic signals from the same transponder before it completes one localization.

It can be observed that the vehicle may have moved to the next location when it receives at due to the vehicle’s motion and the drift term of the current as shown in Figure 2. Thus, the one-way travel time (OWTT) of an acoustic signal traveling from the AUV to a transponder may not equal the opposite case. The problem of the time is not synchronized between transponders and AUV is considered in [18], which measures the two-way travel time (TWTT) based on transponder with

Assume the vehicle has moved from to when the vehicle receives as shown in Figure 2. For simplicity we assume to be a constant in a short time, but is unknown but bounded; thus, the traveling distance of the vehicle is calculated by

Define the OWTT of a signal traveling from the vehicle (at ) to transponder as , from transponder to the vehicle as . Because and may not be equal; thus, the TWTT model is redefined as follows:where denotes the difference between and .

In Figure 2, let and , respectively, denote the angle between vectors and and the angle between vectors and . Define the elevation AOA and azimuth AOA of the received acoustic signal from vehicle as and and the elevation AOA and azimuth AOA of the received acoustic signal from transponder as and , which are measured from the positive -axis and -axis by AOA antennas, respectively. Since in shallow water, the ray trace of sound can be treated as straight line [21], the values of direction cosine for along -axis, -axis, and -axis are , , and , respectively. And the values of direction cosine for along -axis, -axis, and -axis are , , and , respectively [22].

The pitch angle and yaw angle of the vehicle can be measured by inclinometer and -axis gyro. The only difference between the vehicle’s local coordinates and the global Cartesian coordinates is that -axis is in the opposite direction of -axis. Thus, in the global Cartesian coordinates, the pitch angle of vehicle is and the yaw angle is still . The values of direction cosine for along -axis, -axis, and -axis are , , and , respectively [18].

Define the unit vectors of , , and as , , and , respectively. Then it can easily get their expression according to the angle measurements, resulting in the fact that the cosine value of angles and can be solved by using their unit vectors as follows.

It can be observed that and in Figure 2; thus, cosine value of and can also be denoted as follows:where is the sound speed. Thus, the following system of equations can be derived by substituting (7), (8) into (9), (10), respectively,

It easily leads to and by substituting (1), (2) into (11); then the OWTT of a signal traveling from the vehicle (at ) to transponder and the OWTT from transponder to the vehicle (at ) can be further solved by combining (3).

Note that the vehicle orientation may have changed when the vehicle receives the reply from transponder , as shown in Figure 2. Thus, in this case, we should reconsider the combined displacement according to displacement in each vector velocity. The specific strategies are as follows: the vehicle changes direction every time; its direction can be expressed as , and separately. Define the unit vectors of them as , and , respectively. Then it can be known thatwhere , and can be calculated with the same form as described in (5).

2.2. Multiple-Array Underwater Acoustic Localization

In underwater environment, the sound speed is not constant because of the different temperature, pressure (depth), and salinity. But with varying salinity, temperature, and depth, thus the sound speed is computed by exploiting the empirical formula [12]: namely,where is the temperature (in Celsius), is the depth (in meters), and is the salinity (in parts per thousand). Define the sound speed for the vehicle and transponder as and , respectively, which can be calculated after measuring the corresponding parameters.

The single-way distance measurements are obtained by the distance estimation method [18]; it is called accurate single-way distance measurement (ASDM) owing to the fact that the current is taken into account in this paper: namely,where denotes mathematical expectation of sound speed acquired by calculating the sound speed of different depths.

In a three-dimensional underwater environment, at least three transponders are required for 3D AOA target localization. Thus, transponders will have kinds of positioning combinations, which can solve the multiple-group target position according to single-way distance measurements as follows:where is the undetermined position coordinate of the vehicle, denote the determined position coordinates of the corresponding transponders, and and , , , and denote the distance measurements from the vehicle to the corresponding transponder, respectively.

Although multiple measurements of the undetermined vehicle location can be obtained according to various localization combinations, diverse measurements contain different credible information of the actual vehicle location due to the error of the algorithm and the complexity and uncertainty in the underwater environment. As a result of this difference, partial measurements acquired by the above method can belong to outliers. Therefore, a particle swarm optimization (PSO) algorithm based on the outliers elimination is proposed in this paper and the details will be given in the following paragraphs.

3. PSO Based on Outliers Elimination

Based on multiple-array positive sound source localization system, more accurate target position can be obtained by fusing multisource information; however, on the basis of the problem proposed in Section 2, it is necessary to eliminate the outliers that contain poor information of target position. The occurrence of outliers is mainly caused by the following factors, the measured values obtained by the sensor at a certain moment contain too much noise information thanks to the uncertain underwater environment, the packets can also be broken or damaged when the transponders respond to the vehicle with , and so on. In this paper, we adopt the anomaly detection algorithm called the Isolation Forest [23] to search for the outliers from multiple measurements of the undetermined vehicle location. Then the algorithm of particle swarm optimization based on simulated annealing is adopted to optimize measurements of the target location after the outliers have been eliminated.

3.1. Elimination of Outliers

The key to eliminate the outliers is to detect outliers, which is called the outlier detection and has wide application in the field of machine learning and research. An outlier is an observation which deviates so much from other observations as to arouse suspicions that it was generated by a different mechanism [24]. Isolation Forest (iForest) Algorithm [23] is an anomaly detector that measures without distance or density. It performs an operation to isolate each instance from the rest of instances in a given data set. Because outliers have characteristics of being ’few and different,’ they are more susceptible to isolation in a tree structure than normal instances. Therefore, outliers have shorter average path lengths than those of normal instances over a collection of isolation trees (iTrees). Therefore, based on this consideration, we adopt the approach of iForest to eliminate the possible anomalies in the multiple-array underwater acoustic positioning.

The whole algorithm of the iForest mainly contains two stages. The first is training stage; iTrees are constructed by recursively partitioning the given training set until instances are isolated or a specific tree height is reached that will result in a partial model [24]. When building an iTree, first we randomly extract a batch of samples from the given training set and then randomly select an attribute and a value between the maximum value and the minimum value of the attribute. The data in the sample that are smaller than the value of the data are divided into the left tree, and the data that are greater than or equal to the data are divided into the right tree. Finally, in the left and right branches, repeat the above steps until the following conditions are met: the data cannot be divided, that is only one data, or all the same data; the binary tree reaches the limit of the maximum depth. Details of the training stage are described in [23].

The second stage is evaluating stage; an anomaly score is derived from the expected path length for each test instance. are derived by passing instances through each iTree in an iForest; a single path length is derived by counting the number of edges from the root node to a terminating node as instance traverses through an iTree. When is terminated at an external node, where , the return value is plus an adjustment . Given a data set of instances, thus the path length of the instance in each iTree can be calculated as follows:where denotes the number of edges from the root node to a terminating node and is an adjustment and accounts for an unbuilt subtree beyond the tree height limit. is calculated bywhere is the harmonic number and it can be estimated by (Euler’s constant). As is the average of given , we use it to normalize . The anomaly score of an instance is defined aswhere is the average of from a collection of isolation trees. is the sample size of training sample in single iTree. is the average path length of the binary tree constructed by using instances.

According to formula (18) of the anomaly score, when , , then the entire sample does not really have any distinct anomaly; when , , then they are definitely anomalies; when , , then they are quite safe to be regarded as normal instances.

The result of multiple-array localization will be taken as the evaluation sample and input to the trained network of the binary tree. Then it can be searched for outliers from the multiple-array localization and remove them to get relatively accurate measurements for the undetermined vehicle location. Though it is still unknown which one is closer to the real vehicle location, PSO is a population-based stochastic optimization technique which can effectively solve complex optimization problems [25]. Thus, an improved PSO algorithm will be adopted to optimize the relatively accurate measurements in multiple-array positive sound source localization.

3.2. Improved PSO Localization

The basic idea of PSO approach based on outliers elimination is that the vehicle optimizes the result of multiple-array localization after eliminating outliers, and then it easily gets the suboptimal localization result for the real vehicle position by this optimization. In the interest of ensuring the convergence of PSO and eliminating the boundary, we use PSO to increase the contraction factor. The update equation of the speed and location is as follows:where and (nonnegative constants) are the learning factor. is the contraction factor, . and are the independent pseudo-random numbers obeyed uniform distribution of the . , , and represent the number and dimension of the particle, respectively. represents the best position that each particle has been experienced, and represents the best position that has been experienced globally.

Note that the particles are likely to move towards a local minimum by using the above formula to update. Then take the basic idea of the simulated annealing (SA) into account and determine a location from numerous as the global optimum to replace the actual global optimum with the strategy of roulette. Thus, update the speed as follows:

The probability for each individual can be selected as the global optimum is calculated bywhere is the size of species. is the constant. represents the fitness of the th particle and represents the fitness of the global optimal particle undergone.

The single-way distance measurements are assumed to be disturbed by noise with Gaussian distribution due to the measurement error of the time and anglewhere is the vector of real distances between the vehicle and transponders , , ; is the measurement of ; represents the measurement noise vector which cause ranging errors; is the covariance matrix of , where is the variance of the distance measurement noise based on transponder .

In order to get the accurate vehicle position , we take the sum square error function of distance measurements as the fitness function of PSO methodwhere represents the 2-norm.

Remark 1. Based on the proposed network model, the weight for distance measurement noise is comprised of time and angle measurement noise. In simulation experiments, their impact on localization accuracy will be further analyzed in the following paragraphs, respectively.

4. Simulation Experiments

In this paper, to test the performance of the proposed approach, some simulation experiments are conducted which were coded in MATLAB. In these experiments, a vehicle with some sensors moves in a predetermined trajectory, some single-way range localization algorithms are used to keep track of the robot position, and the vehicle travels freely with velocity with and constrained by the current with in a 3D underwater environment of 400 m × 400 m × 300 m. The space coordinate varies from to . Eight transponders are deployed in the space coordinate varying from to (see Figure 3). The noise model is unknown, for simplification without loss of generality, the noise is given out artificially in the simulation experiments. The sound speed is set to be 1500 m/s. The specific type of sensors and the noise function for a given sensor are ignored in this paper. To show the advantages of PSO localization approach based outliers elimination, it is compared with the RTR technology [17] and SL-STTS algorithm [18]. To remove the random effects on the single-way range localization and validate the performance of our algorithm, every experiment was conducted 20 times. The parameters in all of the experiments are the same, which are listed in Table 1.

ASDM noises are mainly caused by the noises of time and angle of arrival measurement. Thus, the distance measurement performance of each algorithm under the same distance measurement noise condition will be firstly analyzed in Section 4.1; the impact of distance measurement noise on the localization accuracy will be compared in Section 4.2; then this paper further investigates the impact of the time measurement noise and the angle measurement noise on the localization accuracy, respectively.

4.1. Performance Comparison of Range Measurement

In this experiment, the noise of the distance measurement is under Gaussian distribution condition; we set , (in Figure 4), and (in Figure 5). The results of the distance measurement performance based on the proposed approach in this experiment are shown in Figures 4 and 5. The absolute difference between the real distance and its measurement is regarded as error standard of distance measurement. In Figure 4, for the simplicity, we compare the distance measurement performance under the different vehicle speed due to the vehicle speed that can affect the time difference and increasing monotonically with . In Figure 5, and represent the unit vector of and , respectively, where the distance measurement performance of the proposed approach under the different current speed is compared to other methods, and several typical current directions relative to the vehicle speed are considered.

Remark 2. The average distance measurement error (ADME) is calculated bywhere and are the real value and the predicted value of single range from the vehicle to the transponder at the th step, respectively; is set as 20. represents the 2-norm.

The results of ADME in Figure 4 show that the distance measurements calculated by the ASDM are closest to the actual distance from the vehicle to the transponders, the time difference will simultaneously grow with increasing, and the error of the proposed approach is less than both the RTR and SL-STTS methods, while the ASDM method is less affected, which shows that the ASDM approach can deal with the single-way ranging measurements stably. Because the size and direction of the current are unknown but bounded, thus we adopt several typical directions for the current relative to , and the results of single-way distance measurement performance relative to the current of the different size and direction in Figure 5 show that the error of the proposed approach is less than both RTR and SL-STTS approaches. It is easy to observe that the proposed method is less affected by the current . The results in this experiment show that the proposed approach has a good performance to deal with the localization problem with increasing.

4.2. Performance Comparison of Localization

To further validate the performance of the proposed approach, this experiment is conducted. We adopt the root mean square error (RMSE), the maximum absolute error (MAXE), and the minimum absolute error (MINE) as performance evaluation indexes. In this experiment, we firstly analyze how the range measurement noise affects localization accuracy. It should be noted that the accuracy of time and angles dominates the precision of OWTT estimation and further affects localization performance. So how the measurement noise of time and angles impacts localization will be further analyzed, respectively. The results of the localization task based on the proposed approach are compared with the SL-RTR and SL-STTS algorithm in Figure 6, where the vehicle moves on the specified trajectory shown in Figure 3 and performs its self-localization, and the sampling time interval . Figure 7 is the root mean square error comparison of the vehicle under different noise levels, and the performance indexes of the localization task under different noise levels are listed in Tables 24. For simplicity, the variance of distance measurement noise based on transponder is set equal and independent of each other; the time measurement based on transponder is assumed to be disturbed by Gaussian distributed noise with ; the input angles are disturbed by Gaussian distributed noise with , and we set and .

Remark 3. The root mean square error (RMSE) in this study is calculated bywhere is the real position of the vehicle at time instant and is the measurement of the vehicle position. represents the 2-norm.

The results in Figure 6 show that the trajectory calculated by PSO-OE is closest to the actual trajectory; it means that the absolute error based on the proposed approach is less than both the SL-RTR and SL-STTS methods. And the fluctuation of errors in the proposed approach is very little, which shows that the proposed approach can deal with the localization problem stably. From the results shown in Figures 7(a) and 7(b), we can see that the values of RMSE for SL-RTR, SL-STTS, and PSO-OE have a growth, but the RMSE values of PSO-OE rise much slower (see Table 2). In Figures 7(c) and 7(d), the values of RMSE for all these three methods go up sharply. The reason is that the bias of time measurements in milliseconds may cause error of distance estimates in meters. The localization performance of PSO-OE is better compared with SL-RTR and SL-STTS under the same noise conditions (see Table 3). We can also observe that the RMSE of PSO-OE goes up slowly as increases quickly (see Figures 7(e) and 7(f), and Table 4). In other words, the localization accuracy is sensitive to the time and angle measurement noise, which verifies the significance of precise ranging, and PSO-OE works better than the other two methods.

5. Discussions

The results of the simulation experiments in Section 4 show that the proposed approach can satisfy the localization task under various situations. In this section, the robustness and fault tolerance of the proposed approach are discussed, where the problem of noise or faults acting on transponders will be taken into account [26].

At first, the noise acting on transponders of the proposed approach is discussed. Transponder time measurement noise has been taken into consideration in Section 4, where the experiments were based on the fact that transponders location is at the fixed inertial location. However, transponders location will inevitably drift with the accumulation of time caused by the current in the actual situation. Therefore, an experiment is conducted in the case of a drift in the transponder location, where the proposed approach is also compared with the SL-RTR and SL-STTS methods. For simplicity, is set equal to each other and distance independent in this experiment, and we set . Location noise on each dimension of transponder is assumed to be disturbed by Gaussian distributed noise with . The setting of other parameters is the same as the second experiment in Section 4. The experimental results are shown in Figure 8 and Table 5. The results show that the localization performance is sensitive to location noise of transponders, which verifies the significance of transponders location; however, the proposed approach contributes to accurate location estimates of vehicle and its sensitivity to transponders location noise is less than the other two methods.

To further discuss the robustness of the proposed approach under the condition of faults acting on transponders, a simulation experiment was conducted, where the ranging noise distribution and the parameters of the proposed approach are the same as the experiment in Section 4, except that the number of the transponders working well is different. In this experiment, it is assumed that two transponders will fail. The experimental results are shown in Figure 9 and Table 6. The results show that the accuracy of localization will decrease with the increase of ranging noises when two transponders fail. However, the localization accuracy of the proposed approach is higher than another two methods. Moreover, comparing it with the results in Table 2 (where all transponders work well), it is not difficult to find that the accuracy of the proposed approach decreased more slowly compared to the other two methods as two transponders malfunctioning: namely, the influence of transponder faults on the proposed approach is less than the other two methods.

6. Conclusions

Localization of autonomous underwater vehicles (AUVs) is investigated in this paper, and a single-way range localization approach based on PSO of outliers elimination is proposed. In the proposed approach, various conditions are considered in the localization algorithm based on SL-STTS, including time-synchronization free localization and long propagation delay. Furthermore, an accurate single-way range localization algorithm based PSO of outliers elimination is proposed to improve the localization performance of the SL-STTS, which can reduce the error caused by the current and make full use of the vehicle position information from different measurements. The proposed approach can improve the accuracy of the single-way range localization, without any a priori knowledge of the noise model. Also, the proposed approach can deal with the searching and exploring problem in unknown environments, which has broad applications, such as the perceiving and detecting for the underwater dam crack.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


The authors would like to thank the National Natural Science Foundation of China (61203365 and 61573128), the Jiangsu Province Natural Science Foundation (BK2012149), and the Fundamental Research Funds for the Central Universities (2015B20114) for their support of this paper.