A New Framework for Synthetic Aperture Sonar Micronavigation
Salvatore Caporale Yvan Petillot
cs.SY cs.RO 93C85 I.2.9










Caporale et al.: A New Framework for Synthetic Aperture Sonar Micronavigation
Caporale et al.: A New Framework for Synthetic Aperture Sonar Micronavigation

























A New Framework for
Synthetic Aperture Sonar Micronavigation
    Salvatore Caporale and Yvan Petillot
S. Caporale  and Y. Petillot are with the Institute of Sensors, Signals and Systems, Heriot-Watt University, Edinburgh, Scotland, UK (e-mail: s.caporale@hw.ac.uk; y.r.petillot@hw.ac.uk).
    December 30, 2023
===================================================================================================================================================================================================================================




Synthetic aperture imaging systems achieve constant azimuth resolution by coherently summating the observations acquired along the aperture path. At this aim, their locations have to be known with subwavelength accuracy. In underwater SAS, the nature of propagation and navigation in water makes the retrieval of this information challenging. Inertial sensors have to be employed in combination with signal processing techniques, which are usually referred to as micronavigation. In this paper we propose a novel micronavigation approach based on the minimization of an error function between two contiguous pings having some mutual information. This error is obtained by  comparing the vector space intersections between the pings orthogonal projectors. The effectiveness and generality of the proposed approach is demonstrated by means of simulations and by means of an experiment performed in a controlled environment.




Synthetic Aperture Sonar, Micronavigation.






§ INTRODUCTION



SAS systems share with SAR many practical and theoretical aspects, as they were originally introduced by moving the synthetic aperture paradigm from radar to sonar. Therefore, most image formation algorithms which have been conceived for SAR have been also considered in SAS literature<cit.>. Despite the overlapping concepts such as range migration and range invariant resolution, underwater SAS systems are operated in a much more challenging environment than SAR as (i) the navigation in water is affected by non-negligible errors and cannot always rely on an external accurate positioning system, (ii) motion errors are comparable to the wavelength hence their effect on image formation algorithms is remarkably destructive (iii) due to relatively small sound velocity in water with respect to the range of interest, the desired along-track resolution cannot be achieved by means of an AUV provided with a single Tx and Rx and moving at a reasonable along-track speed <cit.>.

For the above mentioned reasons, SAS systems are usually equipped with an accurate INS to recover the real navigated trajectory. Moreover, each ping consists of a single Tx and an array of Rx, allowing for a higher along-track sampling rate, hence a higher along-track speed <cit.>. However, from ping to ping a certain degree of redundancy is imposed to perform data based motion estimation known as DPCA which might be coupled with the INS for a more accurate motion compensation <cit.>.

A typical SAS system and the way it operates are illustrated in Fig. <ref>. Other approaches based on autofocus techniques have been also presented ni the literature <cit.>.

In a stripmap imaging system where the AUV moves along a straight trajectory at constant speed, major issues are caused by cross-range displacements. In fact, the sensitivity of the SAS PSF is remarkably higher along the cross-track than along the along-track <cit.>. The cross-track motion has to be known with subwavelength accuracy, whereas the along-track motion has to be in the order of portion of the sampling step. Hence, micronavigation algorithms focus on finding cross-track errors whereas it is generally assumed that the along-track locations are either correct or can be estimated by the navigation system <cit.>. 
The theoretical foundation of several currently employed approaches relies on (i) PCA, (ii) the hypothesis that the AUV rotation angle is small such that it can be characterized as a linearly increasing delay along the array elements and (iii) contiguous pings have some superimposed phase centers, such that the mutual delays between the corresponding tracks can be estimated by finding their correlation peak. Many challenges in SAS micronavigation still have to be addressed <cit.>. Those include (i) dealing with non straight trajectories, such as circular ones <cit.>, (ii) reducing the hardware cost by relying on less accurate navigation and (iii) analysing theoretical limits in performing unsupervised motion compensation.


Following the work in <cit.>, we here exploit the superimposition of some phase centers as for DPCA. However, we assume that they can be arbitrarily interlaced rather than exactly superimposed.

Then, the vector space intersection between the space corresponding to the two pings is described as a functional of the hypothetical displacement. Each ping is employed to compute a pair of displacement-dependent outputs being the projections on the intersection subspace. A proper convex error function between those outputs features a minimum in correspondence of the ping-to-ping displacement, so that it can be identified by means of an optimization. This approach  is computationally demanding due to the evaluation of an error function involving the recomputation of the projecting operators. However, it requires less a priori knowledge and is remarkably less restrictive in comparison to DPCA.



The paper is organized as follows. Section <ref> reviews the acoustic model and SAS, while Section <ref> illustrates the proposed motion compensation procedure. Results from synthetic data and from a real experiment are shown in Section <ref> and <ref> respectively and some conclusions are drawn in Section <ref>.



§ ACOUSTIC AND SAS MODEL




In this Section we first review the observation model for a single Tx and Rx acoustic system when a generic signal is given as input, then we show how this can be employed to retrieve the scene reflectivity when multiple observations are combined together. A simplified 2-D model is considered by means of the assumption that (i) the seabed reflectivity can be represented as a function on a plane and (ii) the 2-D scene reflectivity can be obtained as the projection of the 2-D seabed reflectivity on a slanted plane. Moreover, all the sensors are assumed to be still. Finally, SAS is described as a collection of such single Tx and Rx systems.



 §.§ Acoustic Model


Given a Tx at position z_t=(x_t,y_t) of length D_t and orientation ϑ_t and a hydrophone at position z_r=(x_r,y_r) of length D_r and orientation ϑ_r, the system response to an input complex passband function s(t) with respect to a scene whose complex reflectivity is expressed by ρ(z), with z=(x,y) can be obtained by:


    r(t)=∫_zρ(z) α(z_t,z_r,z,ϑ_t,ϑ_r) s(t-τ(z_t,z_r,z)) dz


where τ is the delay corresponding to the propagation from z_t to z and from z to z_r, α is an attenuation factor taking into account (i) the attenuation due to the propagation distance  from z_t to z and from z to z_r and (ii) the transmission and reception radiation patterns. With regard to the attenuation, we consider the exploding sources model where the α^2∝δ(z_t,z)+δ(z_r,z)∝τ(z_t,z_r,z) rather than α^2∝δ(z_t,z)δ(z_r,z), where δ is the distance operator. The model (<ref>) is represented in Fig. <ref>.












When s(t) is a narrowband impulse modulated at frequency f_0, the above model can be used for observing ρ(z). A suitable sampling grid z_n=(x_n,y_n) matching the available bandwidth is introduced so that the reflectivity is represented as the input vector ρ(z_n). By setting G(z_n) as the Green's function of the system


    G(z_n)=α(z_t,z_r,z_n,ϑ_t,ϑ_r)e^-j2πf_0τ(z_t,z_r,z_n),




the following output time-sampled signal ϕ(t_m) can be obtained as


    ϕ(t_m)=A(t_m,z_n) G(z_n)ρ(z_n)


where the matrix A(t_n,z_n) performs the discrete-space integration of equation (<ref>), which is basically done by an interpolation kernel and a summation. From a numerical point of view, A(t_m,z_n) is usually defined by its transpose one, since defining an interpolation along the intrinsic one dimensional time axis t_m is easier than doing it over the fictitious space axis z_n. We consider a NUFFT based time interpolation <cit.>, having an approximately linear complexity and high accuracy.

The integral in (<ref>) and its space discrete equivalent in (<ref>) clearly perform a dimensionality reduction from the two-dimensional space z_n to the one dimensional space t_m.  In fact, all points lying on the ellipsis having z_t and z_r as foci are intrinsically ambiguous with respect to the Tx Rx pair. This is also represented in Fig. <ref> where we highlighted a propagation path (dotted line) having the same delay as the one relative to the generic point z_n. As a consequence, inverting (<ref>) recovers the reflectivity ρ(z) affected by defocusing along ellipsis. 
Due to the radiation patterns, this defocusing increases with range.



 §.§ SAS model


The working principle of a synthetic aperture consists in performing multiple observations at prescribed Tx and Rx locations along a (usually straight) path, called along-track or cross range. By doing so, a longer aperture is emulated and a narrower beam is obtained in the direction orthogonal to the track, called range. The emulation is performed by signal processing and allows for getting a constant resolution along range.

Given a set of Tx and Rx z_l,t and z_l,r, l∈Z, the observation model of SAS can be written as the column vector of ϕ_l(t_m)=A_l(t_m,z_n) G_l(z_n) ρ(z_n).























In case each observation is a monostatic system, i.e. z_l,t=z_l,r=z_l, and z_l are uniformly spaced along the cross range direction by D/4, with D=max(D_t,D_r), then ρ(z_n) can be approximately recovered with a cross range resolution equal to D/2 by the transpose observation operator, also know as backprojection


    ρ(z_n)≈∑_l∈ZG^*_l(z_n)A_l^†(z_n,t_m)ϕ_l(t_m)


which has been conveniently rewritten as a sum of backprojections of each single Tx and Rx system.

As for getting a cross range resolution equal to D/2 a sampling on the along-track by D/4 has been considered, the oversampling ratio is equal to 2. Many variations of this setup have been considered in the literature in order to decrease the oversampling ratio and employ different reconstruction techniques <cit.>.



Practical SAS systems are mounted on an AUV. Hence, the speed v of the vehicle has to be set according to the wanted along-track sampling and the desired maximum range. For instance, in case of an AUV equipped with a single Tx and Rx, v≤D/4max(τ(z_l,t,z_l,r,z_n))=Dc/8R, where R is the maximum range. For R=[150]m and D=[5]cm, we get v=[6.25]cm/sec. A feasible speed is obtained by using an array of N equispaced Rx at distance L=D/2 and a single Tx. By doing so, a single transmitted ping allows for N observations, thus gaining a N factor on the speed. Such a SAS system is a collection of bistatic single Tx and Rx systems. If the range of interest is sufficiently large with respect to the wavelength, it can be approximately modelled as a collection of monostatic systems such that z̅_l=(z_l,t+z_l,r)/2. This model is referred to as PCA.



§ MOTION ESTIMATION



The motions of an AUV in a 3-D space are represented by means of 3 linear motion parameters being heave, sway and surge, and 3 rotation motion parameters being pitch, roll and yaw. As the acoustic model has been described in 2-D, only the effective projections of those on the slanted plane are needed, i.e. surge and sway as effective linear motions and yaw and roll as effective rotations. Roll only affects the radiation pattern over range and can be neglected. Without loss of generality, we here refer to surge as the motion error with respect to the prescribed linear trajectory at constant speed. Surge, sway and yaw are illustrated in Fig. <ref>.


The SAS image formation requires an accuracy by approximately λ/10 on z̅_l, where λ is the wavelength. Underwater navigation is not accurate enough to allow this accuracy for the employed frequency band, usually [100-300]kHz. As a consequence, a motion estimation has to be performed to guarantee that the acquired observations are coherently combined. This is usually achieved by combining navigation information coming from physical inertial sensors place on the AUV, that is the INS, with information extracted from the acquired data by introducing redundancy from ping to ping. Instead of moving the AUV at the maximum allowed speed, a reduced speed is adopted in order to have a certain number of K<N superimposed phase centers, as represented in Fig. <ref>. This approach is usually referred to as DPCA. Under the assumptions that (i) PCA holds, (ii) surge is negligible and (iii) yaw is small, the tracks from two contiguous pings referred to the same phase center differ for a time delay and the delays of the K-N superimposed phase centers are an affine function of the sway and yaw. The delays can be estimated by correlating the corresponding tracks, thus also sway and yaw can be obtained. Interested readers are referred to <cit.> for more details. An approach to tackle surge estimation is presented in <cit.>.



We here propose a technique loosening the three assumptions above, more specifically (j) PCA is replaced with a less restrictive approximation, (jj) the along-track motion is such that an unknown non-integer number of superimposed phase centers are present and (jjj) yaw is only small enough to guarantee that the area illuminated by two contiguous pings has a non-null intersection.











 §.§ Ping to Ping Displacement Estimation



Given that the receiving array is composed by N equispaced Rx, the raw data obtained by collecting the reflections from the scene of a single impulse will be referred to as ping. For the sake of the stripmap SAS under consideration, the AUV is moving along a straight trajectory affected by motion errors. The scene is observed by means of the start-and-stop approximation. Given the trajectory, each ping can be represented as a collection of bistatic systems identified by the array midpoint and orientation and the Tx position and orientation:


    σ_r^(p)   =   (x_a^(p),y_a^(p),ϑ_r^(p))
    σ_t^(p)   =   (x_t^(p),y_t^(p),ϑ_t^(p)).


where x_a and y_a are the coordinates of the midpoint of the receiving array elements. We here assume that the relative position of σ_t^(p) with respect to σ_r^(p) is known. In a real scenario, the displacement between the transmitter and the array midpoint is time-dependent as the AUV moves between the transmission instant and the reception relative to the considered range interval (see <ref>). 
We also introduce the operator T^(p) representing the observation model relative to ping p, so that T^(p)ρ gives the raw data relative to ping p. We also set




    ρ̆^(p)=T̃^(p)-1 T^(p) ρ


where ρ̆^(p) represents the recovered reflectivity at ping p by means of the pseudoinverse operator T̃^(p)-1, where T̃ is the available computable approximation of the real observation operator T^(p). Mathematically, T̃^(p) is the column matrix of the real observation matrices A_l,γG_l,γ, l=1,…,N (see (<ref>)).

Then, we define the orthogonal projector on the subspace identified by ping p


    Q^(p)=(T̃^(p))^-1T̃^(p).


Given another generic ping q and its orthogonal projector Q^(q), the projection on the intersection of the subspaces identified by Q^(p) and Q^(q) can be obtained by


    ψ^(p,q)=ψ^(q,p)=lim_i→∞(Q^(q)Q^(p))^iρ=lim_i→∞(Q^(p)Q^(q))^iρ.


By taking advantage of (<ref>), an approximation of ψ^(p,q) can be obtained by starting either from ρ̆^(p) or ρ̆^(q). In fact


    ρ̆^(p)   ≃   T̃^(p)-1T̃^(p)ρ=Q^(p)ρ
    ρ̆^(q)   ≃   T̃^(q)-1T̃^(q)ρ=Q^(q)ρ


hence by highlighting the expressions (<ref>)-(<ref>) in (<ref>) we get


    ψ^(p,q)   =   lim_i→∞(Q^(p)Q^(q))^iQ^(p)ρ
       ≃   lim_i→∞(Q^(p)Q^(q))^iρ̆^(p)=ψ^(p)
    ψ^(p,q)   =   lim_i→∞(Q^(q)Q^(p))^iQ^(q)ρ
       ≃   lim_i→∞(Q^(q)Q^(p))^iρ̆^(q)=ψ^(q)


and finally we get the relationship


    ψ^(p)≃ψ^(q).


The meaning of (<ref>) is that if the intersection operator is known, it is possible to identify the same intersection image starting from the raw data obtained by two different pings. In case the operator T^(p) and T^(q) are orthogonal, both ψ^(p) and ψ^(q) are null, hence their computation is not relevant. Conversely, in case they are not orthogonal, the fact that ψ^(p)-ψ^(q)≃0, where · is a generic norm, can be exploited as a property to identify the intersection operator.

In order to apply the principle expressed in equation (<ref>), we introduce operator Q^(p)_0 being a shifted/rotated version of Q^(p) such that the center of mass of the PCA is the origin of the axes. In more detail, by setting

    x̅^(p)   =   (x_a^(p)+x_t^(p))/2
    y̅^(p)   =   (y_a^(p)+x_t^(p))/2


we get


    σ_r,0^(p)   =   (x_a^(p)-x̅^(p),y_a^(p)-y̅^(p),0)
    σ_t,0^(p)   =   (x_t^(p)-x̅^(p),y_t^(p)-y̅^(p),ϑ_t^(p)-ϑ_r^(p))


and we set


    σ̅^(p)=(x̅^(p),y̅^(p),ϑ_r^(p)).


By introducing the shift and rotation operator S_σ, where σ specifies a triplet (x,y,ϑ), we finally get


    Q^(p)=S_σ̅^(p)Q^(p)_0S_-σ̅^(p).


Ping q can be specified by Q_0^(q) and σ̅^(q), or, alternatively, by Q_0^(q) and σ̅^(p) and the differential displacement between σ̅^(p) and σ̅^(q), which will be referred to as σ̅^(q,p)


    Q^(q)   =   
    	S_σ̅^(q)Q^(q)_0S_-σ̅^(q)
       =   
    	S_σ̅^(p)S_σ̅^(q,p)Q^(q)_0S_-σ̅^(q,p)S_-σ̅^(p)


with


    σ̅^(q,p)   =   (x̅^(q)-x̅^(p),y̅^(q)-y̅^(p),ϑ_r^(q)-ϑ_r^(p))
    	
       =   (x^(q,p),y^(q,p),ϑ^(q,p))


hence the following is verified


    σ̅^(p,q)=-σ̅^(q,p).



It is worth noting that the differential displacement which has been highlighted is equal to the one which would be obtained by replacing the bistatic situation with its PCA, nevertheless, this model maintains the capability of representing a bistatic system where the transmitter and the array have different orientations.



















So, the vector σ̅^(q,p) fully identifies the displacement between ping p and q and represents the goal of the motion estimation procedure. In order to estimate it, we rewrite equations (<ref>) and (<ref>) by taking into account that ρ̆^(p) and ρ̆^(q) are not available, whereas their rotated ans shifted version by S_-σ̅^(p) can be computed


    S_-σ̅^(p)ψ^(p)   =   lim_i→∞(Q_0^(p)S_σ̅^(q,p)Q^(q)_0S_σ̅^(p,q))^iS_-σ̅^(p)ρ̆^(p)
    
    	S_-σ̅^(p)ψ^(q)   =   lim_i→∞(S_σ̅^(q,p)Q^(q)_0S_σ̅^(p,q)Q_0^(p))^iS_-σ̅^(p)ρ̆^(q).















Since the displacement between p and q is unknown, the above equations can be rewritten with respect to a hypothetical displacement σ=(x,y,z)


    S_-σ̅^(p)ψ^(p)_σ   =   lim_i→∞(Q_0^(p)S_σQ^(q)_0S_-σ)^iS_-σ̅^(p)ρ̆^(p)
    
    	S_-σ̅^(p)ψ^(q)_σ   =   lim_i→∞(S_σQ^(q)_0S_-σQ_0^(p))^iS_-σ̅^(p)ρ̆^(q)































In general, for non-null ρ̆^(p) and ρ̆^(q) and non-null intersection operator, we have


    ψ^(p)_σ=ψ^(q)_σ  ⇔  σ=σ̅^(q,p).



In practice, they will not be equal because of approximation errors. However,  ψ^(p)_σ and ψ^(q)_σ can be employed for building an error function. Since they represent wavefield reconstructions, they feature periodic-wise variations with period equal to half the wavelength. An error function based on those would also feature such oscillations, hence it is expected to be convex on a small hyperinterval around the displacement triplet (x^(q,p),y^(q,p),ϑ^(q,p)):


    η^(q,p)(x,y,ϑ)=
    	 ψ^(p)_σ-ψ^(q)_σ _2.


Conversely, their absolute values represents reflectivity amplitudes so an error function based on them is expected to be convex on a larger hyperinterval although it does not necessarily have its minimum in (x^(q,p),y^(q,p),ϑ^(q,p))


    ζ^(q,p)(x,y,ϑ)=
    	 |ψ^(p)_σ|-|ψ^(q)_σ| _2.


Functions ζ and η will be referred to as modulus and phase error function respectively. The interval of convexity of the above functions will be discussed later in the paper. Normalized version of the above functions are more conveniently employed, where the normalization can be done with respect to either ψ^(p)_σ, ψ^(q)_σ or the maximum between the two.

Hence, the estimation procedure can be performed by:

	
  * initialize σ: σ←(0,0,0)
	
  * minimize ζ^(p,q) starting from σ: σ←minζ^(p,q)
	
  * minimize η^(p,q) starting from σ: σ←minη^(p,q)


The characteristics of ζ and η cannot be analytically identified. A practical assessment of this proposed procedure will be provided in Section <ref>. The following considerations hold:

	
	
  * the superimposition between phase centers is not required as the along-track displacement is also part of the search;
	
  * the estimation relies on global similarities in the image domain rather than on track-by-track correlations, thus more robust;
	
  * rotation estimation relies only on non-null vector space intersection between pings, hence large rotation can be estimated.

However, there are some major computational burdens:

	
  * each optimization iteration requires the update of the observation model of one of the two pings;
	
  * both the orthogonal projector and the vector space intersection requires infinite iterations which must be properly truncated.

To decrease the computational load, the range of interest can be conveniently limited.


























 §.§ Approximate Bistatic Model




















In the previous subsection we assumed each ping to be obtained as the output of a collection of N bistatic systems under the start-and-stop approximation. As the estimation approach do not rely on the PCA, we here introduce a more general approximation allowing for representing (i) the orientation of the Tx with respect to the orientation of the array and (ii) the effective displacement along the cross-track between the Tx and the array which is happening because the AUV motion. The employment of bistatic systems has been shown to carry some benefits <cit.>, hence the capability of including them is beneficial.

With the assumption that the transmitted signal is ideally instantaneous (no pulse compression) and that the observed range interval tends to zero, σ_r^(p) and σ_t^(p) can be determined by sampling the effective AUV trajectory on the slanted plane. According to PCA are approximated by their average. Infinite different bistatic systems sharing the same PCA can be identified. Among these, we approximate σ_r^(p) and σ_t^(p) by σ̃_r^(p) and σ̃_t^(p), such that



    σ̃_r^(p)   =   σ̅^(p)+ς_r
    	
    σ̃_t^(p)   =   σ̅^(p)+ς_t+(0,0,ϑ^(p)_t-ϑ^(p)_r)


where ς_r=(x_a,y_a,0) and ς_t=(x_t,y_t,0) such that ς̅=(0,0,0). In practice, the real bistatic system is replaced by another bistatic system having the same PCA but conveniently chosen to be a more accurate representation according to the considered range. For instance, if the Tx is positioned on the array midpoint, y_a=y_t=0 and x_a=-x_t would be the distance covered by the AUV while the impulse reaches the target range.

In a real scenario, the whole illuminated range is observed. Hence, the range is practically split in subranges where the above approximation is adopted. The splitting in subranges must be anyway performed as the AUV ping-to-ping motions also change along range.


















































































































































§ SIMULATION RESULTS



In order to assess the performances of the proposed approach, we here present results from a simulated scenario. An experiment performed in a controlled environment is presented in Section <ref>.

The simulation focuses on highlighting the generality, the accuracy and the convergence of the proposed estimation procedure with no side information provided. No restrictive assumptions are made on surge and yaw.














 §.§ Interval of Convergence



A SAS system consisting of 16 elements is considered, both the Tx and the Rx elements have D=[5]cm aperture, hence the whole array is [80]cm wide, whereas the phase centers are distributed uniformly along a [40]cm interval. As the number of the expected ping-to-ping superimposed phase centers, we considered K=4, so that each ping is [20]cm apart from its adjacent ones. With this setup, an AUV moving at [2]m/sec would be capable of covering a [75]m range. The Tx is physically located in the array midpoint, whereas it is virtually located [10]cm apart from the array midpoint position as an effect of the vehicle motion. The range observed with this assumption would be around half the maximum available range, i.e. [37.5]m. Nevertheless, being a simulation with the start-and-stop approximation, we deliberately considered [10]m as midrange in order to have  a configuration where the PCA would be less accurate, than the actual model (<ref>)-(<ref>). The system is operating at [300]kHz, while the bandwidth is [30]kHz, thus resulting in [2.5]cm range resolution. The cross-range resolution is approximately equal to D/2, i.e. also [2.5]cm. 



An instance of the objective functions ζ and η has been plotted in Fig. <ref> with respect to surge, sway and yaw separately. In more detail, the surge has been varied in the interval [-0.2,0.2], i.e. the array has been moved from having its phase centers totally superimposed to the ones of the previous ping to not having any superimposed. On the other axis, a range of hypothetical motion errors has been considered. The surfaces on top qualitatively show how both the error function ζ and η have minima on the plane diagonal, meaning that the surge can be correctly estimated. Also, ζ evaluated along the x axis appears mostly convex, such that the minimum can be found by starting from x=0. The error function η is steepest around the minima, but features a more remarkable non-convex behavior which makes it suitable for optimization only if a good starting point is given. Similar considerations can be done for the sway, which has been varied also in the remarkably large interval [-0.1,0.1], and for the yaw. Their η functions are not showed as they become very steep around the minima and cannot be properly visualized. 
Those results are partial as surge, sway and yaw have been considered separately, i.e., it is not showed whether those functions are convex and having the right minima when all the three motions are considered at the same time.





 §.§ Ping-to-ping Accuracy


As a second experiment, we consider a single ping q=9 featuring motion errors with respect to its previous p=8, both relative to the scene depicted in Fig. <ref> (top). It features a noisy background emulating sand ripples plus a phantom circular object whose diameter is [40]cm. Pings 8 and 9 happen in the middle of the considered cross-range interval, thus facing the circular object. The motion errors are set as follows: surge x^(q,p) equal to [-3.45]cm, sway y^(q,p) equal to [2.17]cm and yaw ϑ^(p,q) equal to [4.10×10^-2]rad, i.e., approximately [2.35]deg. Those values are random, but their magnitude is such that other approaches such DPCA would fail. In fact, the along-track motion error is larger than along-track sampling D/2 and not a multiple of D/2. 

The rotation error is such that the effect on raw data track cannot be described by simple delays. In Fig. <ref>, raw data of tracks 9-12 of ping 8 are compared to tracks 1-4 of ping 9 . 

In Fig. <ref> we plotted the error functions ζ and η when the considered motion errors are applied one at a time, i.e., the same as in Fig. <ref>. The left column shows function ζ in solid line and η in dotted line, whereas the right column shows the opposite. Real displacements are identified by the dashed vertical lines. The figures also show the minima which are found by the optimizations. With regard to ζ it can be seen that by starting with no a priori knowledge, i.e. σ=(0,0,0), the optimization is convex. For surge ans sway, the size of the convexity interval is proportional to the along-track sampling step, whereas for yaw it is also inversely proportional to the number of superimposed phase centers K. From the figure referred to sway (middle-left) the oscillatory nature of function η can be observed. Its period is roughly proportional to half the wavelength making it essential that the optimization on ζ has an accuracy being higher than this. 
The minima found on ζ exhibit some deviations with respect to the real displacements, which are then corrected by the optimization on η function. The minima of η are not equal to 0 as the algorithms for computing the orthogonal projections and the intersection have been truncate to 5 iterations for the sake of implementability.





From a qualitative point of view, it is interesting to have an understanding about what the projections on the intersection space represent. In Fig. <ref> the evolution of the projections and their difference along the optimization process are shown. In the initial situation (top) obtained for (x^(q,p),y^(q,p),ϑ^(q,p))×1/2, the projections feature differences in shape, position and amplitude. The intermediate situation (middle) obtained for (x^(q,p),y^(q,p),ϑ^(q,p))×3/4, exhibits higher correlation between the projections but their displacements still causes the error to have the same degree of magnitude as the projections. Finally, the situation at convergence still shows some residual error and some shape differences between the projections. This example highlights the estimation robustness versus punctual differences in reflectivity observed from the two pings, as the the error functions takes into account global similarities. 

As a final step, we consider the whole optimization in the case the given displacements are happening at the same time. The optimization on function ζ gives surge equal to [-2.66]cm, sway  equal to [2.15]cm and yaw ϑ^(p,q) equal to [4.08×10^-2]rad, whereas the subsequent optimization on η gives [-3.49]cm, [2.17]cm and [4.12×10^-2]rad respectively. The optimization process is illustrated in Fig. <ref> where the surge, sway and yaw component of the vector along the steepest descent are shown together with the function values. 



From a qualitative point of view, the resulting estimation error is highly compatible with SAS operations, which usually require the error on sway lower than λ/8. 
From a quantitative point of view, the assessment of the performance is quite hard as it is difficult to separate the dependency on the specific scene reflectivity from the dependency on the system parameters. Hence, estimating the error variance and bias could be only approached statistically. Nevertheless, it could be argued that both theoretical and/or statistical results on simulated data might focus on aspects which could be remarkably overtaken from other issues when moving to real data. For this reason, we believe it is not worth pursuing a further statistical estimation error analysis before assessing the working principle on real data. 









 §.§ Trajectory Accuracy


As a final simulation, we provide the estimated trajectory along a [3]m path covered by 16 pings. The considered scene reflectivity is again the one depicted in the top of Fig. <ref>. The reconstructed scene with no navigation is plotted in the middle figure. 
To make the simulation more realistic, the generation of raw data takes into account the different location and rotation of the Tx with respect to the Rx, whereas the reconstruction is performed by means of the proposed bistatic model (<ref>)-(<ref>) assuming that the expected Tx to Rx distance along the cross-range direction is known and equal to [10]cm. 
The considered random trajectory is partially shown in Fig. <ref> (top-left), together with a graphical explanation of the employed bistatic model. In the top-right figure the position of the Tx and array midpoint are explicitly shown.

The result of the trajectory estimation are shown in Fig. <ref>. The differential and cumulative errors on surge and sway might suggest small biases, which can be neglected as the final sway error is still only about λ/25. Finally, the obtained trajectory is employed to reconstruct the scene in the top figure of Fig. <ref>. The obtained image is qualitatively comparable to the original image and the whole estimation procedure fulfils the initial goal of estimating motion errors with no priors.



§ EXPERIMENTAL RESULTS




Given the peculiarity of the proposed estimation technique and its complex algebraic structure, assessing its validity on real data is a major concern. An experiment in a controlled environment has been designed at this aim. Nevertheless, the a priori knowledge of the sensor locations could not be provided by means of the available equipment. For this reason, a specific manner to validate the experiment has been identified and designed.






 §.§ System Setup


For this purpose, a water tank being [3]m large, [4]m long and [2]m deep was available. In order to avoid reflections from the concrete wall and bottom, a target consisting of 4 metal reflectors lying on an octagonal plastic surface parallel to the surface and mounted on a wooden support at approximately [1]m from the bottom has been considered. The geometry of the object is detailed in Fig. <ref>. A SAS system consisting of 8 Rx of [5]cm aperture and spaced by [5]cm and a [5]cm aperture Tx positioned at the midpoint of the array has been provided by Hydrason Solutions. The SAS system is capable of operating between [20]kHz and [170]kHz. The tank is provided with an industrial plotter which has been employed in order to move the SAS system along a straight line and perform the acquisitions by the start-and-stop approach. 
According to the specifications, the plotter provide a submillimetre accuracy. However, as the plotter engine produces a non-negligible noise, it had to be switched on and off throughout the whole acquisition process. This operation has caused the displacement accuracy to be inaccurate. 
Other causes of inaccuracy might have come from other mechanical non-idealities, such as submillimetre deviations in the supporting guide.



The SAS system has been attached to the plotter by means of a pole and mounted on a rotating plate capable of providing roll and yaw. Roll has been adjusted in order to make the beam illuminate the object before touching the tank bottom, whereas yaw has been adjusted to be null. However, these setups have been manually arranged, so no apriori accurate information can be assumed for those parameters. Hence, the positioning of the array throughout all its displacement along the guide can feature both yaw and roll. Finally, the relatively closeness of the target object makes the PCA assumption weak and the reflectivity to vary along the trajectory, thus causing a lack of coherence from ping to ping.

The SAS raw data have been collected by moving the system by steps equal to the spacing between phase centers, i.e. [2.5]cm, being also the expected along-track resolution, i.e. D/2. 
An interval of [80]cm plus the length of the array PCA, i.e. [20]cm, has been covered by 32 pings, thus getting a synthetic aperture of [1]m. As far as the physical parameters are concerned, a [30]kHz chirp of duration equal to [3]msec, central frequency equal to [105]kHz and sampled at [1]MHz has been employed. After pulse compression, demodulation and filtering, a downsampling by 33 has been performed, thus getting a cross-track resolution equal to [2.45]cm.






 §.§ Validation Methodology




By means of M pings with K=N-1, MN different observations are collected. So, for this experiment 256 tracks are acquired, which can be used in various ways. For instance, given an array whose PCA covers uniformly [20]cm, the considered synthetic aperture of [1]m can be equivalently covered with 5 pings and K=0 or 10 pings and K=4. By means of the considered 32 pings, two main tasks can be accomplished: (i) the 32 pings can be subsampled in order to produce P=1,…,N-K different motion estimates for K=1,…,7; (ii) by selecting one of the 8 Rx at a time, 8 new synthetic aperture can be obtained as a collection of the corresponding SISO systems. Both these tasks are illustrated in Fig. <ref>. The simple case K=7 allowing for the single phase P=1 is rearranged in order to obtain the case K=4,5,6 with all the corresponding different phases. Also, the case K=7 is reorganized in virtual arrays of SISO systems.

The benefit of the SAS made up of SISO systems is the following. Given that the elements are moved on a straight line, no sway can be present. A small rotation of each element is not actually impacting on the imaging task, so surge is the only motion error to be identified. Hence, the 8 images which can be obtained the SISO approach should be similar apart from a slight warping in the track direction. Hence they can be employed as initial ground truth. Conversely, the multiple motion estimates at different sampling steps can be exploited as follows. For each value of K, each of the P=1,…,N-K estimates represents a subsampling of the full sampled trajectory with different initial phase. Hence, they can be integrated in order to obtain a single motion estimate for each value of K. Then the obtained K trajectories can be applied separately on each of the 8 SAS based on the SISO approach, or to each one of the 8 SAS based on the proper  array, which can be referred to as SIMO system. This would lead to K×8 images for the SISO and for SIMO case respectively. Ideally, if the motion estimation procedure works accurately, all those reconstructed complex reflectivities should be equal.





 §.§ Motion Errors


With regard to surge, it was known that in correspondence of ping 12 a major misalignment occurred. In Fig. <ref> all 8 tracks corresponding to a single location have been plotted. Tracks 8-11 and 12-15 are not correlated as a consequence of this surge event. Moreover, a finer comparison of the tracks would reveal that other differences are present even within each group. Those can be attributed to a lack of calibration among the array elements and also to the fact that the 8 SISO systems are supposed to share the same phase center but practically look at the scene from different angles. 

When considering the SIMO systems, surge is the same as for the SISO ones, whereas sway can appear as an effect of the initial orientation of the array. This is explained in Fig. <ref>, where it is shown that the midpoint of a rotated array moves diagonally with respect to the reference system aligned with the array. Hence, the angle of this diagonal line is a measurement for the array orientation with respect to the scene reference system.







 §.§ Estimated Trajectory


Before performing the actual trajectory estimation by means of the optimizations, it is worth analysing the error functions in order to assess that they feature the same characteristics as in Fig. <ref>. Error function ζ and η with respect to pings 11 and 12 are plotted in Fig. <ref>. All the macroscopic features highlighted for emulated data are replicated for real data. The functions are convex on comparable intervals and η with respect to sway oscillates with period equal to half the wavelength (the oscillation is much more visible here as the central frequency is [105]kHz rather than [300]kHz). For the considered ping, surge, sway and yaw are equal to [1.71]cm, [-3.1]mm and [1.35×10^-3]rad respectively. The yaw motion error is less than 1/10 of one degree, so it can be considered negligible.

Moving from one ping to all the tracks, the results are shown in Fig. <ref>. By properly combining the differential displacements obtained for all the different phases illustrated in Fig. <ref>, a trajectory estimation can be obtained for each considered value of K, thus resulting in 4 trajectories plotted in Fig. <ref> (top). The differential surge and sway are also shown (left-bottom). It is worth noting that all the phases different from 1 have an undetermined starting point which has been numerically tuned to maximize the match with the trajectory obtained for K=1. Also, the differential motion errors represented in the figure are interlaced. As the trajectory lays  along a diagonal path with respect to the track, the cumulative displacements have been rotated by an angle equal to the opposite of the diagonal orientation, i.e. [1.27×10^-2]rad or [0.73]deg, and the same angle has been added to the yaw. The resulting cumulative displacements (right-bottom) give a hint about how the accuracy might decrease when the number of overlapping phase centers decreases. The overall deviation is enclosed in less than [5]mm and [0.5]mm for surge and sway respectively. 







 §.§ Reconstructed Reflectivity


As mentioned before in this Section, the absolute value of the reconstructed complex reflectivities obtained from the SISO wise SAS have to be equal apart slight deformation due to surge. Absolute values have been used as an initial way to assess the correctness of the estimation approach. Nevertheless, it is more proper to look at the complex reflectivities and verify that the both the real and the imaginary parts, 
are equal for both the SISO and the SIMO SAS for all the 8 set of pings having K=0.

Following this, in Fig. <ref> absolute value of the real part of the reconstructed reflectivities for the SISO SAS have been plotted for 4 of the 8 possible set of pings with no motion compensation. Despite the fact the their absolute values are fairly similar, the phase is different as the SISO element phase centers are located on different parallel lines. In Fig. <ref> the same plots are proposed with respect to the SIMO SAS. The motion compensated version of these cases have been represented in Fig. <ref> and Fig. <ref> respectively, each by means of a different estimated trajectory in order to show that the error among them is negligible. The areas where the coherent metallic reflectors are located have been highlighted to emphasize the consistency across the various cases. As expected, the rest of the area features some differences for the intrinsic lack of coherence due to the close range.




§ CONCLUSION



A novel motion compensation technique for SAS capable of identifying surge, sway and yaw with no restrictions and no prior information from inertial sensors has been introduced. The proposed approach is based on the comparison between the projections on the intersection subspace of different observations. The effectiveness of the technique has been proved by extensive analysis and simulations together with an experiment on real data.



§ ACKNOWLEDGEMENTS


This work was supported by the Engineering and Physical Sciences Research Council (EPSRC)
Grant number EP/J015180/1 and the MOD University Defence Research Collaboration
in Signal Processing.







10
	
	url@samestyle
	
	
	
	
	
	
	CafforioPratiRocca1991
	C. Cafforio, C. Prati, and F. Rocca, “Sar data focusing using seismic
	migration techniques,” IEEE Transactions on Aerospace and Electronic
		Systems, vol. 27, no. 2, pp. 194–207, Mar 1991.
	
	BillonPinto1995
	D. Billon and M. A. Pinto, “Some general considerations for synthetic aperture
	sonar design,” in OCEANS '95. MTS/IEEE. Challenges of Our Changing
		Global Environment. Conference Proceedings., vol. 3, Oct 1995, pp.
	1665–1670 vol.3.
	
	HayesGough2009
	M. P. Hayes and P. T. Gough, “Synthetic aperture sonar: A review of current
	status,” IEEE Journal of Oceanic Engineering, vol. 34, no. 3, pp.
	207–224, July 2009.
	
	BellettiniPinto2009
	A. Bellettini and M. Pinto, “Design and experimental results of a 300-khz
	synthetic aperture sonar optimized for shallow-water operations,” IEEE
		Journal of Oceanic Engineering, vol. 34, no. 3, pp. 285–293, July 2009.
	
	BellettiniPinto2002
	A. Bellettini and M. A. Pinto, “Theoretical accuracy of synthetic aperture
	sonar micronavigation using a displaced phase-center antenna,” IEEE
		Journal of Oceanic Engineering, vol. 27, no. 4, pp. 780–789, Oct 2002.
	
	CallowHayesGough2009
	H. J. Callow, M. P. Hayes, and P. T. Gough, “Motion-compensation improvement
	for widebeam, multiple-receiver sas systems,” IEEE Journal of Oceanic
		Engineering, vol. 34, no. 3, pp. 262–268, July 2009.
	
	gough2004
	P. T. Gough and M. A. Miller, “Displaced ping imaging autofocus for a
	multi-hydrophone sas,” IEE Proceedings - Radar, Sonar and Navigation,
	vol. 151, no. 3, pp. 163–170, June 2004.
	
	marston2015
	T. M. Marston and D. S. Plotnick, “Semiparametric statistical stripmap
	synthetic aperture autofocusing,” IEEE Transactions on Geoscience and
		Remote Sensing, vol. 53, no. 4, pp. 2086–2095, April 2015.
	
	herter2016
	U. Herter, H. Schmaljohann, and T. Fickenscher, “Autofocus performance on
	multi channel sas images in the presence of overlapping phase centers,” in
	OCEANS 2016 MTS/IEEE Monterey, Sept 2016, pp. 1–6.
	
	marston2012
	T. Marston, “A correlation-based autofocus algorithm for coherent circular
	synthetic aperture sonar,” in EUSAR 2012; 9th European Conference on
		Synthetic Aperture Radar, April 2012, pp. 66–69.
	
	caporale2016
	S. Caporale and Y. Petillot, “A novel motion compensation approach for sas,”
	in 2016 Sensor Signal Processing for Defence (SSPD), Sept 2016, pp.
	1–5.
	
	FesslerSutton2003
	J. A. Fessler and B. P. Sutton, “Nonuniform fast fourier transforms using
	min-max interpolation,” IEEE Transactions on Signal Processing,
	vol. 51, no. 2, pp. 560–574, Feb 2003.
	
	caporale2007
	S. Caporale, L. D. Marchi, and N. Speciale, “A svd-based algorithm for dense
	nonuniform fast fourier transform,” in 2007 15th European Signal
		Processing Conference, Sept 2007, pp. 2120–2124.
	
	hunter2016
	A. J. Hunter, S. Dugelay, and W. L. J. Fox, “Repeat-pass synthetic aperture
	sonar micronavigation using redundant phase center arrays,” IEEE
		Journal of Oceanic Engineering, vol. 41, no. 4, pp. 820–830, Oct 2016.
	
	lepage2002
	K. D. LePage and H. Schmidt, “Bistatic synthetic aperture imaging of proud and
	buried targets from an auv,” IEEE Journal of Oceanic Engineering,
	vol. 27, no. 3, pp. 471–483, Jul 2002.
	












































