Gianni Franchi
email: gianni.franchi@uni-siegen.de

Angela Yao
email: yao@cs.uni-bonn.de

Andreas Kolb
email: andreas.kolb@uni-siegen.de

Supervised Deep Kriging for Single-Image Super-Resolution

We propose a novel single-image super-resolution approach based on the geostatistical method of kriging. Kriging is a zero-bias minimum-variance estimator that performs spatial interpolation based on a weighted average of known observations. Rather than solving for the kriging weights via the traditional method of inverting covariance matrices, we propose a supervised form in which we learn a deep network to generate said weights. We combine the kriging weight generation and kriging process into a joint network that can be learned end-to-end. Our network achieves competitive super-resolution results as other state-ofthe-art methods. In addition, since the super-resolution process follows a known statistical framework, we are able to estimate bias and variance, something which is rarely possible for other deep networks.

Introduction

Super-resolution aims to transform low resolution (LR) images into images with high resolution (HR). In computer vision, super-resolution is relevant for applications where high-frequency information and detailing is desirable yet not always captured, e.g. medical imaging, satellite imaging, surveillance, etc. Our interest is in single image super-resolution (SISR), a special case where only one LR image is available. SISR is an ill-posed inverse problem with multiple solutions, i.e. multiple HR images could lead to the same LR image after applying a low-resolution filter.

SISR been solved with many different approaches. Early approaches relied on spatial relationships present in the image and include nearest neighbours [START_REF] Chang | Super-resolution through neighbor embedding[END_REF], bi-cubic interpolation [START_REF] Keys | Cubic convolution interpolation for digital image processing[END_REF]. Other techniques build upon the assumption that the underlying signal is sparse, e.g. compressive sensing [START_REF] Donoho | Compressed sensing[END_REF] and dictionary learning [START_REF] Yang | Coupled dictionary training for image super-resolution[END_REF]. These early approaches tend to yield HR images which are blurry and or noisy, due mainly to the fact that these methods are crafted heuristically. As a result, they cannot sufficiently recover enough information to yield sharp photo-realistic images.

Deep learning, and more specifically, convolutional neural networks (CNNs) [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], has afforded significant performance gains in SISR. The newest frameworks [START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF][START_REF] Kim | Deeply-recursive convolutional network for image super-resolution[END_REF][START_REF] Lim | Enhanced deep residual networks for single image super-resolution[END_REF][START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF][START_REF] Sajjadi | EnhanceNet: Single image superresolution through automated texture synthesis[END_REF] learn feature representations for super-resolution in a supervised end-to-end manner and are achieving ever-improving performance both in terms of speed and quantitative evaluation. However, despite their impressive results, most deep networks are being used as black boxes to make point estimates. The regressed outputs represent only a mean value and give no indication of model uncertainty. In the case of SISR, this means that we make (possibly very accurate) guesses for the unknown pixels, but we have no idea how good these guesses are. This can be highly problematic for applications such as microscopy or medical imaging, where model confidence is just as important as the actual output.

We would like to leverage the strong learning capabilities of deep networks and embed them within a known statistical framework for which we can derive model uncertainty. As such, we propose deep kriging to solve SISR. Kriging [START_REF] Matheron | Random sets and integral geometry[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF] is a geostatistics method used for spatial interpolation of attributes such as topography and natural resources. It has close relations to Gaussian process (GP) regression, although some sources oversimplify and describe the two as being equivalent. The main difference is that kriging makes stationary and ergodic assumptions on the random field representing the data, while GP regression considers it as a Gaussian process. 3Kriging interpolates unknown values by taking a weighted average of known values and results in an unbiased estimator with minimum variance. One can solve for the weights by inverting a covariance matrix computed on the known instances. Kriging makes for a natural extension to the application of SISR, where HR images are interpolated from pixels observed in the LR version of the image. The computational cost of the kriging is cubic with respect to the number of known instances; hence it is preferable to work locally. However, interpolated results are often overly smooth and depend on the choice of the covariance function.

Our contributions can be summarized as follows:

-We propose a novel deep learning-based statistical estimator with SISR results competitive with state-of-the-art. -We propose a deep kriging approach that solves for kriging weights and performs the kriging in a single network; this allows us to perform a supervised form of spatial interpolation which can be applied not only to SISR, but other geostatistical computations such as field reconstructions.

Related work

Prior to the use of deep learning, SISR approaches applied variants of dictionary learning [START_REF] Freeman | Example-based super-resolution[END_REF][START_REF] Donoho | Compressed sensing[END_REF][START_REF] Yang | Coupled dictionary training for image super-resolution[END_REF][START_REF] Peleg | A statistical prediction model based on sparse representations for single image super-resolution[END_REF][START_REF] Dong | Image deblurring and super-resolution by adaptive sparse domain selection and adaptive regularization[END_REF]. Patches were extracted from the low-resolution images and mapped to their corresponding high-resolution version which are then stitched together to increase the image resolution. Other learning-based approaches to increase image resolution include [START_REF] Zhao | Wavelet-domain hmt-based image super-resolution[END_REF][START_REF] Zhang | Single image super-resolution with non-local means and steering kernel regression[END_REF][START_REF] Anbarjafari | Image super resolution based on interpolation of wavelet domain high frequency subbands and the spatial domain input image[END_REF]. State-of-the-art SISR methods are deep-learning-based [START_REF] Lim | Enhanced deep residual networks for single image super-resolution[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF][START_REF] Kim | Deeply-recursive convolutional network for image super-resolution[END_REF][START_REF] Sajjadi | EnhanceNet: Single image superresolution through automated texture synthesis[END_REF][START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF]. The VDSR [START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF] and DRCN [START_REF] Kim | Deeply-recursive convolutional network for image super-resolution[END_REF] approaches showed the benefits of working with image residuals for super-resolution. The DRRN approach [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF] generalizes VDSR and concludes that the deeper the network, the better the super-resolved image. We also use a residual network in our approach, but unlike all the other deep SISR methods, we are solving for a set of filter weights to perform the super-resolution with our network rather than directly estimating the HR image itself.

Several unsupervised techniques have also been developed for SISR [START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF][START_REF] Marquina | Image super-resolution by tv-regularization and bregman iteration[END_REF], though their performance is usually poorer than supervised approaches. One work of note [START_REF] He | Single image super-resolution using Gaussian process regression[END_REF], uses GP regression to perform super-resolution and resembles ours in spirit in that we both model pixel intensities as a random process regressed from neighbouring pixels. However, [START_REF] He | Single image super-resolution using Gaussian process regression[END_REF] is unsupervised while our approach learn the weights from supervise examples.

Our proposed approach can be thought of as a form of local filtering [START_REF] Brabandere | Dynamic filter networks[END_REF], and is most similar in spirit to works which combine deep learning and GP regression [START_REF] Wilson | Deep kernel learning[END_REF][START_REF] Damianou | Deep gaussian processes[END_REF]. However, we differ from these techniques since we do not apply GP regression but a modified version of local kriging. These techniques are similar to us in that they also consider the data to follow a random process. However, we do not explicitly learn the covariance matrix, which offers us more flexibility in the relationships we wish to express. Furthermore, we treat each image as following a random process; the set of training images is then a set of random processes, while the aforementioned works consider all the training data to follow the same process.

Deep Kriging

We begin with a short overview on classical kriging in Section 3.1 before we introduce our proposed method of supervised kriging in Section 3.2 and elaborate on its statistical properties in Section 3.3. Finally, we show how the proposed form of deep kriging can be implemented in a deep network in Section 3.4.

Classical kriging

Consider an image f as a realization of a random field F . This random field is a collection of random variables, i.e. F = {F i } where each F i is a random variable and i is a position index. Unknown values on f at position x * , f * 4 , can be interpolated linearly from n known realizations f i , with normalized weights w i , i.e.

f * = n i=1 ω i (x * )f i , with n i=1 ω i (x * ) = 1. (1)
In classical kriging, all known realizations are used in the interpolation and n is the number of pixels in the image, though local variants have been proposed and studied in [START_REF] Pronzato | Bayesian local kriging[END_REF][START_REF] Meier | Local gaussian regression[END_REF].

The weights {ω i } are found by minimizing the true risk:

R(ω) = E (F * -F * ) 2 = var F * -F * . (2)
We can equate the E[•] with the var(•) term in Eq. 2 because the constraint that the weights must sum up to 1 implies that E[f * -f * ]=0. However, since we do not have access to different realizations of the random field, this variance cannot be solved directly.

To infer the variance from a single event, the theory of geostatistics replaces the classical statistics assumption of having independent and identically distributed (iid) random variables with assumptions of stationarity and ergodicity on the random field. The random field is assumed to be first and second order stationary. This means that E F i (f i (x)) does not depend on x and that

E F i (f i (x) × f i (x + τ ))
depends only on τ . In addition, the field is assumed to be second order ergodic, so a covariance or mean estimate in the probability domain is equivalent inthe spatial domain. This second assumption implies that an empirical covariance can be estimated and depends only on a distance τ between two realizations, i.e. τ = x i -x j . However, it is difficult to work directly with the empirical covariance, since it can be noisy and also is not guaranteed to form a positive semi-definite matrix. As such, in kriging, we use the true covariance by fitting a parametric model. While different models can be used, the most common one is Gaussian, with covariance C ij between points x i and x j defined as

C ij = C 0 exp -1 σ 2 x i -x j 2
, where C 0 and σ are parameters of the model. A Lagrange multiplier (with constant λ) can be used to minimize the risk function with the constraints on the weights w i , leading to the following cost function:

L(w) = var F * -F * + 2λ n i=1 w i (x * ) -1 , (3)
and the associated solution expressed in matrix form as:

     C11 . . . C1n 1 . . . . . . . . . . . . C1n . . . Cnn 1 1 . . . 1 0           w1(x * ) . . . wn(x * ) -λ      =      C1 * . . . Cn * 1      . ( 4
)
The process of kriging then is reduced to solving for the weights w i . This involves inverting the covariance matrix on the left. Since the matrix is of dimension (n + 1), where n is the number of observed pixels in the image (patch), it can be very computationally expensive. Even though one can limit n with local kriging, n needs to be sufficiently large in order to accurately capture and represent the spatial relationships. We aim to bypass the matrix inversion and still maintain sufficient generalization power by using a deep network to directly learn the weights instead.

Supervised kriging

Consider E, a subset of the discrete space Z 2 , as a representation of the support space of a 2D image. We denote an image as f ∈ R and f i as its value at pixel x i ∈ E. This assumes we work with greyscale images, which is a standard assumption made in most SISR algorithms [START_REF] Kim | Deeply-recursive convolutional network for image super-resolution[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF][START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF].

We assume that we are given a set of

n 1 training image pairs {( f i , f i )}, i ∈ [1, n 1 ]
, where f i is an up-sampled low resolution version of f i with the same size. Our objective is to super-resolve the low-resolution images of the test set

N test = { f i } n2
i=1 . We further assume that f i is a realization of a first and second moment stationary random process F i . For convenience we denote F i the distribution of the random process. In addition f is a realization of a random field F . We denote the training set as

{ f i , f i } n1 i=1 ∼ M,
where M is a joint meta-distribution where each training pair ( f i , f i ) ∼ F i × F i follows its own random process.

In classical kriging, the estimated f (x) is expressed as a linear combination of all the observed samples i.e. pixels of f . In the case of super-resolution, the observations come from the low-resolution image f (x). Furthermore, we estimate f (x) as a linear combination of only local observations of some window radius K, leading to the estimator:

f i * = k∈{ x k -x * 1≤K} ω k (x * ) f i k with k∈{ x k -x * 1 ≤K} ω k (x * ) = 1. (5)
When no confusion is possible we denote F i * = n k1=1 ω k (x * ) F i k , with n = (2K + 1) 2 . We found that a window of 7 × 7 provided the best results. To find ω, we want to minimize the true risk as given in Eq. 2, leading to:

R(ω) = n1 i=1 E F i × F i F i * -F i 2 . ( 6
)
This risk is a compromise between the geostatistical model and the deep learning model, where the covariance of each field is learned through supervised learning.

Furthermore, if we replace the true expectation with the empirical one and the random field with its realization, we arrive at

R(ω) = n1 i=1 Ni j=1   f i * - n k=1 ω k (x * ) f i k 2   , (7)
where N i is the number of pixels in image f i . Now, rather than solving for the weights ω by inverting the covariance matrix in Eq. 4, we propose learning this set of weights directly with a CNN. This allow us to capture and model more complex spatial relationships between the data than the classical model with covariances. As such, the estimator becomes f (x, f , g, θ), where w = g( f , θ), with g and θ representing the CNN network function and its parameters respectively.

Statistical properties

We can derive statistical properties for our proposed estimator F . Given {ω 1 , . . . , ω n }, the bias can be expressed as

bias(x) = E F |ω 1 ,...,ωn F (x) -F (x) = E F |ω 1 ,...,ωn F (x) -E F |ω 1 ,...,ωn F (x) = E F |ω 1 ,...,ωn F (x) - n k 1 =1 ω k 1 (xj)E F |ω 1 ,...,ωn F (x) .
Since n k1=1 ω k1 (x j ) = 1 and the random field is first-order stationary, we are left with zero bias according to the weights of the neural network, i.e. bias(x) = 0. In other words, our network does not add bias to the field F . As such, from F to F , if the two fields have the same first moment, then no bias added to F 5 . This point is critical since in classical statistics a good estimator should have zero bias and minimal variance.

By definition, the variance of our estimator F is

V F (x) = E F F (x) 2 -E F [ F (x)] 2 . (8)
Again, we estimate the covariance given {ω 1 , . . . , ω n }. Since we assume that F is second-order stationary, the covariance of F depends only on the distance between two points, i.e.

E F F (x k ) F (x k ) = C( x k -x k 2 ) -µ 2 . ( 9
)
By setting µ = E F [ F (x)], we arrive at:

V F |ω 1 ,...,ωn (x) = (k,k )∈[1,(2×K+1) 2 ] 2 ω k (x)ω k (x) C( x k -x k 2) ( 10
)
where the covariance C(•) is estimated from the low resolution image f . The variance depends on the position x. Note, however, since we directly estimate the weights ω(x * ) without making assumptions on the covariance function, we cannot compute the variance and an approximation is needed to estimate the covariance empirically from f . Eq. 10 however is particularly interesting since it shows that the estimator variance at x is directly proportional to the variance of f and the different values of ω near x. The bigger these values are, the more uncertain the estimator is.

Network Implementation and Training

5 Of course, we cannot account for the low resolution process that produced F So far, the theory which we have presented on supervised kriging is generic and does not indicate how one should minimize the risk in Eq. 7 to solve for the kriging weights. We solve for the weights with a CNN network which we show in Fig. 1. The weight estimation branch is composed of a residual network of 9 residual units; the residual unit itself has 2 convolutional layers, each preceded by a batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] and a ReLU.

This branch follows a similar architecture as [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF], with the difference that they apply this architecture to directly estimate the HR image, while we use it to estimate our kriging weights, i.e. as a local dynamic filter [START_REF] Brabandere | Dynamic filter networks[END_REF]. In addition, contrary to [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF] our network is not recurrent.

The network learns a total of 20 convolution layers to determine the kriging weights. The first 19 are of depth 128, the last prediction layer has depth (2 × K + 1) 2 and outputs the kriging weight vector ω(x) = [ω 1 (x), . . . , ω (2×K+1) 2 (x)]. The weights are applied to the repeated input image via a point-wise multiplication and summation along the depth dimension to yield the HR output. Overall, the network with K = 3 is very lightweight, inference on a 320 × 480 sized image takes only 0.10 seconds with a titan X.

Note our network is actually learning how to estimate the kriging weights ω, i.e. as a local dynamic filter. To derive this dynamic filter end-to-end, we adapt the following formulation of the Eq. 5 in convolution terms:

fi (x * ) = k∈{ x k -x * 1≤K} ω k (x * ) fi * h k (x k ) with k∈{ x k -x * 1 ≤K} ω k (x * ) = 1 (11) with h k (x) = δ(x -x k )
or the Dirac function. The h k filter is denoted by "repeat input" in Fig. 1, since applying the Dirac functions in each position x * is a way to extract the pixel neighbourhood at position x * . After having determined the weights, we normalize them so that they sum up to 1, as given by the constraint in Eq. 11. We point-wise multiply the normalized weights to the output of the repeat filter; to arrive at the super-resolved image, we simply sum along the depth. We implement our network in Tensorflow, minimizing the loss in Eq. 7 with the Adam optimizer, a batch size of 8 and a learning rate of 10 -4 . We apply gradient clipping and dropout between the last ReLU and the convolution layer that predicts ω, with a dropout rate of 0.3. The effective depth of our network is d = 2 + 2 × U , with U the number of times the residual units are applied. In our case, we use U = 9.

Experiments

Datasets, Pre-processing and Evaluation

For training, we follow [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF] and use the 291 images combined from [START_REF] Yang | Image super-resolution via sparse representation[END_REF] and the Berkeley Segmentation datasets [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] . We further augment the data by rotating (90 • , 180 • and 270 • ) and scaling the data (×2, ×3, ×4). We limit to these fixed grades of rotation and scaling to follow the same protocol as [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF]. Like [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF], we work with patches of size 31 × 31, sampled from the augmented dataset with a stride of 21. Testing is then performed on the four commonly used benchmarks: Set 5 [2], Set 14 [START_REF] Kim | Deeply-recursive convolutional network for image super-resolution[END_REF], B100 [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] and Urban 100 [START_REF] Huang | Single image super-resolution from transformed self-exemplars[END_REF]. To be consistent with [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF], we convert the RGB images to YCbCr and only resolved the Y component. This is then combined with an upsampled Cb and Cr component and converted back into RGB. Up-sampling is done via bi-cubic interpolation.

We evaluate the resulting super-resolved images quantitatively with the peak signal-to-noise ratio (PSNR): 2 , where f ref and f are the ground truth and the super-resolved images respectively and N refers to the number of pixels in the (high-resolution) image. A higher PSNR corresponds to better results. We also evaluate using Structural Similarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], which measures the similarity between two images. The closer the SSIM value is to 1, the more similar the super-resolved image to the ground truth high resolution image.

PSNR( f , f ref ) = log 10 255 2 N i=1 ( f (xi)-f ref (xi))

Comparison to state-of-the-art

We compare our PSNR and SSIM measures for the different scales against stateof-the-art in Table 1. The first two methods, bi-cubic and local kriging are unsupervised, while all others are supervised approaches. We use the Matlab implementation of bicubic interpolation.

For the local kriging, we use a neighbourhood of 90 × 90 and a stride of 81. At this sparse setting, unsupervised kriging does very well in comparison to bicubic interpolation. However, it is already extremely slow, since for each patch, we need to compute the empirical covariance, true covariance, and then invert the covariance matrix. In total, it takes approximately 1 minute for an image of size 320 × 480. In comparison, our proposed deep kriging, on the same image, applied densely is one magnitude faster, and takes only 0.25 seconds with a titan x.

Looking at our approach (reported in the second last column in Table 1) with respect to the supervised methods, our performance is comparable to the DRRN [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF] B1U9 setting, DRCN [START_REF] Kim | Deeply-recursive convolutional network for image super-resolution[END_REF] and VSDR [START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF]; all three have networks with depth similar to ours. The best results are reported by the DRRN B1U25 [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF], which has 50 convolution layers and is more than twice as deep as our network.

Given the trend of the community to work on deep learning-based approaches for super-resolution, as well as the fact that no labelled data is required for training, one can work with (infintely) large datasets [START_REF] Tong | Image super-resolution using dense skip connections[END_REF][START_REF] Lim | Enhanced deep residual networks for single image super-resolution[END_REF]. For fair comparison with state-of-the-art, however, we omit from this table the techniques which do not use the fixed 291 image training set as per the training protocol set by [START_REF] Tai | Image super-resolution via deep recursive residual network[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF].

Model uncertainty

One of the key strengths of our technique is the fact that we can estimate model uncertainty. We show in Figure 2 the estimated variance for each pixel based on Equation 10. To evaluate this equation, we need to apply covariance models for f and f . We do this by first estimating the empirical covariance c(τ ) = i f (x i )f (x i + τ ) -µ and then solving for C 0 and σ of the Gaussian model in section 3.1 that is the closest to the empirical covariance.

The main advantage of our model of uncertainty is that we can have information on the black box CNN. It gives us feedback about the reliability of results provided by an unknown image. One can see in Figure 2 our uncertainty estimate. The estimated variance has a higher value than the real PSNR this is due to the fact that we have a noisy estimation of the PSNR from a low resolution image. Holistically, however, the are similar, in that areas with high variance corresponds to high PSNR. Quantitatively, we find for Set 5 that 91.9% of the super-resolved pixel values fall within 3 standard deviations based on the estimated variance. In addition, we evaluate the similarity between the images resulting from the PSNR and the one resulting from the variance in set 5. We use the correlation to measure the similarity and find that these images have high correlation up to 0.8.

Conclusions

In this paper, we have proposed a joint deep learning and statistical framework for single image super-resolution based on a form of supervised kriging. We solve for the kriging weights in the manner of a local dynamic filter and apply it directly to the low resolution image, all within a single network that can be learned end-to-end. Since we work within the known statistical framework of kriging, we can estimate model uncertainty, something typically not possible for deep networks. More specifically, we show through derivations that the statistical estimator generated by our network is unbiased and we calculate its variance.

Fig. 1 .

 1 Fig. 1. Overview of the network structure. The input image goes through two branches.The first one calculates the weights ω; the second copies the input image for an efficient application of the weights via point-wise multiplication.

Input

  Image Bicubic local kriging DRRN ours variance RMSE real B1U9 estimate (a) Super resolution results of set 5 data set with scale ×3. Input Image Bicubic local kriging DRRN ours variance RMSE real B1U9 estimate (b) Super resolution results of set 5 data set with scale ×3.

Fig. 2 .

 2 Fig. 2. Super resolution results

-

  Our proposed deep kriging is a hybrid deep learning and statistical framework for which we can derive statistical bias and variance; such measures of model performance and uncertainty are not commonly available for other deep networks, making us the first to model and compute pixel uncertainty in deep SISR approaches.

Set 5

 5 ×2 33.66/0.930 35.68/0.923 36.66/0.953 37.53/0.959 37.63/0.959 37.66/0.959 37.74/0.959 37.65/0.960 ×3 30.39/0.868 31.72/0.864 32.75/0.909 33.66/0.921 33.82/0.923 33.93/0.923 34.04/0.924 33.94/0.923 ×4 28.42/0.810 29.83/0.814 30.48/0.862 31.35/0.883 31.53/0.885 31.58/0.886 31.68/0.889 31.56/0.886 Set 14 ×2 30.24/0.869 31.26/0.874 32.45/0.907 33.03/0.912 33.04/0.912 33.19/0.913 33.23/0.914 33.10/0.912 ×3 27.55/0.774 27.94/0.781 29.30/0.822 29.77/0.831 29.76/0.831 29.94/0.831 29.96/0.834 29.81/0.831 ×4 26.00/0.703 26.43/0.715 27.50/0.751 28.01/0.767 28.02/0.767 28.02/0.767 28.18/0.770 28.08/0.771 B 100 ×2 29.56/0.843 31.06/0.848 31.36/0.888 31.90/0.8960 31.85/0.894 32.01/0.897 32.05/0.897 32.01/0.896 ×3 27.21/0.8431 27.86/0.7357 31.36/0.8879 28.82/0.7976 28.80/0.7963 28.91/0.7992 28.95/0.800 28.82/0.800 ×4 25.96/0.668 26.55/0.665 26.90/0.711 27.29/0.725 27.23/0.723 27.35/0.726 27.38/0.728 27.36/0.727 Urban 100 ×2 26.88/0.840 28.32/0.843 29.50/0.895 30.76/0.914 30.75/0.913 31.02/0.916 31.23/0.919 30.95/0.921 ×3 24.46/0.735 25.06/0.715 29.50/0.799 27.14/0.715 27.15/0.828 27.38/0.833 27.53/0.838 ×4 23.14/0.658 23.70/0.627 24.52/0.722 25.19/0.752 25.14/0.751 25.35/0.757 25.44/0.764 25.10/0.749

Table 1 .

 1 Results of various super resolution techniques. The column 3th and 4th are unsupervised. The other columns are trained on the 291 data set[START_REF] Yang | Image super-resolution via sparse representation[END_REF][START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. The Urban 100 set having no images at scale 3 we did not evaluate it.

We refer the reader to our supplementary materials for a more detailed comparison of the two.

For concise notation, we use f * to denote f (x * ), fi to denote f (xi) and w * i = wi(x * ).

Acknowledgement This research was funded by the German Research Foundation (DFG) as part of the research training group GRK 1564 Imaging New Modalities.


Liang Yu

Weikang Jiang

Jerome Antoni

Quentin Leclere

Haijun Wu

BeBeC-2018-S10 ACOUSTIC BEAMFORMING BY NON-SYNCHRONOUS MICROPHONE ARRAY MEASUREMENTS

Acoustic beamforming consists in measuring the sound field by the microphone array and infers the source strength in conjunction with the propagation model. A fundamental limitation of acoustic beamforming is determined by the aperture size of the array and the microphone density. In order to extend the working frequency range of conventional phased array beamforming, the sound sources are scanned by moving sequentially a small prototype array, resulting in a large array and high microphone density measurements, which is referred to as non-synchronous measurements beamforming. The main issue of nonsynchronous measurements beamforming is that the phase relationships between consecutive snapshots of non-synchronous measurements are missing and result in missing entries in the spectral matrix, however, complete spectral matrix information is crucial for high quality image of beamforming. Thus, non-synchronous measurements beamforming boils down to a spectral matrix completion problem which is modeled by a structured low rank model and solved by designed fast iteration algorithm. In particular, the structured low rank model is constructed by two ingredients: (1) the low rank property of the spectral matrix and (2) the continuity of the acoustic field. The simulation results show the nonsynchronous measurements beamforming extended the working frequency range of given small phased array from 2 -2.8 kHz to 1.6 -4.5 kHz with only 9 times non-synchronous measurements.

INTRODUCTION

Acoustic imaging aims to obtain a visualized sound map from the radiating objects, which has a wide range of applications in source identification, vibration analysis and machine diagnosis [START_REF] Lu | A gearbox fault diagnosis scheme based on near-field acoustic holography and spatial distribution features of sound field[END_REF][START_REF] Pézerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF]. It is generally accepted that the beamforming [START_REF] Dougherty | What is Beamforming?[END_REF][START_REF] Michel | History of acoustic beamforming[END_REF] and Near-field Acoustical Holography (NAH) [START_REF] Williams | Fourier acoustics: sound radiation and nearfield acoustical holography[END_REF] are two most widely exploited methods dedicated to the aforementioned issues. A fundamental limitation of acoustic source imaging is imposed by geometry of the array. Specifically, the minimum frequency of reconstruction depends on the dimension of the array and the maximum frequency depends on the distance between neighboring microphones [START_REF] Kim | Optimal regularisation for acoustic source reconstruction by inverse methods[END_REF]. A solution to achieve a large array and/or high microphone density is to scan the object of interest by moving sequentially a prototype array (an arbitrary array of achieving the sequential measurements), which is referred to as non-synchronous measurements. Non-synchronous measurements have been addressed in [START_REF] Bi | Broadband acoustic holography from intensity measurements with a three-dimensional pressure-velocity probe[END_REF][START_REF] Loyau | Broadband acoustic holography reconstruction from acoustic intensity measurements. i: Principle of the method[END_REF][START_REF] Nejade | Reference-less acoustic holography techniques based on sound intensity[END_REF] in the context of NAH, which breaks through the restriction of using references by reconstructing the phase from the quadratic pressure and tangential components of the sound intensity, yet it requires the use of an array of intensity probe.

The non-synchronous measurements in the context of NAH was further investigated in Ref [START_REF] Antoni | Reconstruction of sound quadratic properties from non-synchronous measurements with insufficient or without references: Proof of concept[END_REF] by taking the incomplete measurements explicitly into account to formulate the acoustic inverse problem, which is recognised an intrusive method since it directly solves the inverse problem with the incomplete measurements and requires manipulation of the structure of current reconstruction methods. Non-intrusive method of non-synchronous measurements has been next investigated in Ref [START_REF] Liang | Acoustical source reconstruction from non-synchronous sequential measurements[END_REF], which is reformulated as a matrix completion problem, and the Cyclic Projection (CP) [START_REF] Yu | Spectral matrix completion by cyclic projection and application to sound source reconstruction from non-synchronous measurements[END_REF] and FISTA [START_REF] Yu | Acoustical source reconstruction from nonsynchronous sequential measurements by fast iterative shrinkage thresholding algorithm[END_REF] algorithm were proposed correspondingly to find a full spectral matrix subject to given constraint of hermitian symmetry, measurement fitting, reduced rank and spatial continuity of the sound field. The non-intrusive method aims to complete a spectral matrix, then feeds it to the spectral matrix based acoustic image method.

However, there remains a need for an efficient beamforming method of non-synchronous measurements. The spectral matrix completion in the context of beamforming is still dissimilar with the one in NAH. First, The beamforming demands a fast and robust implementation of spectral matrix completion, whereas the NAH requests a high accuracy implementation of spectral matrix completion. The different challenges of computational efficiency and accuracy will result in various choices of nuclear norm minimization algorithms. Second, the beamforming is usually implemented in the far field and the NAH in the near field, and the measurements distances are highly connected with the spatial basis choices in spectral matrix completion. Thus, the differences between the NAH and beamforming result in the different methods of spectral matrix completion of data missing from the aspects of interpolation basis choices and the nuclear norm minimization algorithms.

In this paper, the non-synchronous measurements in the context of beamforming is formally formulated as the problem of spectral matrix completion, which is solved by the Fast Iterative Shrinkage Thresholding Algorithm. The paper is organized as follows. The Forward model of sound propagation and acoustic beamforming is described in section 2, and the spectral matrix completion with nuclear norm minimization is given in section 3. Section 4 and 5 are separately the simulation and conclusion.

FORWARD MODEL OF SOUND PROPAGATION AND ACOUSTIC BEAMFORMING

Conventional beamforming

Assume the source plane is equally discretized into N grids at the known positions r n , n = 1, ..., N, and S incoherent broadband sound sources with coordinates r s , s = 1, ..., S are distributed on these grids (the position set {r s } belongs to the set {r n }). The radiated sound is measured by M microphones (N > M > S) with coordinates r m , m = 1, ..., M. The measured sound pressures in the time domain are modeled as stochastic processes, then a number of snapshots are produced by the Fourier transform of a short-time segment tapered with a smooth data window. The M dimensional vector of measured pressures p t = (p 1t , ..., p Mt ) T at a frequency f corresponding to the t-th snapshot (t = 1, ..., T ) is obtained by:

p t = G(r)s * t + n t . (1)
Here, the M × S matrix G = (g 1 , ..., g s , ..., g S ), in which a general term g ms describes the freefield sound propagation model from the s-th source to the m-th microphone:

g ms = G(r m |r s ) = e jk r m -r s 4π r m -r s , (2)
where k = 2π f c is the wavenumber and c is the sound speed. The vector s * t = (s * 1t , ..., s * St ) T contains the uncorrelated random amplitudes of the sound sources: it is assumed to follow an S dimensional complex-valued Gaussian distribution with 0-mean and diagonal covariance matrix

S * cc = E{s * t s * t H } = diag(α α α 2 ), (3)
where α α α 2 = (α 2 1 , ..., α 2 S ) T is vector of powers (variances of amplitudes) of the sound sources and diag(x) denotes the diagonal matrix with the elements of x on the main diagonal. The measurement noise vector n t = (n 1t , ..., n Mt ) T is also assumed to be complex Gaussian distributed with 0-mean and covariance matrix S nn = σ 2 I M . It is reasonable to further assume that s t 1 and n t 2 are independent for any two snapshots t 1 ,t 2 ∈ {1, ..., T }. Therefore, the sound pressure p t measured by a microphone array follows a complex-valued multivariate Gaussian distribution with 0-mean and spectral matrix of propagation is

S pp = E{p t p H t } = GS * cc G H + S nn = Gdiag(α α α 2 )G H + σ 2 I M . (4)
Conventional beamforming consists in estimating the S sources by the sound pressures measured by the M microphones, and the parameters to estimate are the locations r s and the powers α α α 2 s (s = 1, ..., S) of the sound sources. The conventional beamforming is based on the assumption that the signals are emitted by a single source of the discretized source plane r n . Thus the forward model of Eq. ( 1) can be further written as a discrete forward model of signal propagation:

p t = G(r)s t + n t . (5)
Here G(r) is M × N matrix with G = (g 1 , ..., g n , ..., g N ), in which a general term g mn describes the free-field sound propagation model from the n-th source to the m-th microphone:

g mn = G(r m |r n ) = e jk r m -r n 4π r m -r n . (6)
The idea behind each beamformer is to steer the microphone array to search for the source position, and the beamformer power output B(r n = r s ) at position r n of measurement p t is defined as

B(r n = r s ) = w H n Ŝpp w n , where Ŝpp = 1 T ∑ T t=1 p t p H t w n = g n g n 2 2 , n = 1, ..., N (7)
The weights

w n = g n g n 2 2
are designed to be independent of the array data or data statistics, which compensate the time delays and amplitude attenuation of the forward propagation.

Forward model of non-synchronous microphone array measurements

The non-synchronous microphone array measurements are achieved by scanning the sound sources (moving sequentially the prototype array), then the sound image is obtained by the conventional beamforming. Thus, the i-th measurement of non-synchronous microphone array measurements is modeled as the same propagation of Eq. ( 1),

p {i} t = G {i} (r)s * t + n t . (8)
where i = 1, ..., P and P is the number of non-synchronous array measurements. The spectral matrix model of propagation in Eq. ( 4) can then be correspondingly reformulated as the nonsynchronous measurements form

S (i) pp = E{p {i} t p {i} t H } = G {i} S * cc G {i} H + β i 2 I M , i = 1, ..., P, (9)
where the noise power β i 2 , i = 1, ..., P is estimated for each configuration i but equal for all sensors. The estimated spectral matrix is calculated by

Ŝ(i) pp = 1 T ∑ T t=1 p {i} t p {i} t H
. It is noted that the S * cc is invariant for measurements mode since the fundamental premise of the acoustic source distribution is considered as a stationary stochastic field in this paper. The spectral matrices Ŝ(i) pp , i = 1, ..., P of each with size M × M can be methodically rearranged in the block diagonal positions of an incomplete spectral matrix Ŝm pp of MP × MP (see left part of Fig. 1). Contrary to non-synchronous measurements, spectral matrix Ŝpp of simultaneous measurement (assumes the number of microphone in Eq. ( 1)) is MP) is a full matrix with size MP × MP since all the coherence information of array measurements are known. With the purpose of using non-synchronous measurements for acoustic beaforming, the problem is formulated as to seek an full spectral matrix from the spectral matrix Ŝm pp of data missing.

SPECTRAL MATRIX COMPLETION WITH NUCLEAR NORM MINIMIZATION

Low rank modeling of acoustic field

In acoustics, the rank of spectral matrix measures the stochastic complexity of the acoustic field that is produced by the uncorrelated sources (latent variables in stochastic terminology), and it indicates the number of uncorrelated sources. That is to say, although the measured sound field (usually with measurements noise) might look complicated, it is assumed that the observed sound field is only radiated from limited number of sound sources, and many piratical acoustic scenarios have justified this assumption. This limited number of uncorrelated sources produce a correlated sound field, and the rank of spectral matrix can be used to measure the correlation of the acoustical field. For instance, a rank one spectral matrix corresponds to a fully coherent field (the measurement in one position can be easily predicted from other positions in the field), whereas a spectral matrix with full rank corresponds to a diffuse field in general. Thus, the spectral matrix completion problem boils down to seeking a spectral matrix with low rank (minimized rank) subject to a given constraint of hermitian symmetry and corresponding

Liang Yu, Weikang Jiang etc.

measurements of missing, i.e. minimize S

Rank(S)

subject to A (S) -Ŝm pp F ≤ ε S * = S 0, [START_REF] Michel | History of acoustic beamforming[END_REF] where the unknown spectral matrix S to be searched is linked to the partial measurements Ŝm pp by the sampling operator A : C MP×MP → C MP×MP that extracts the elements from diagonal blocks, and the difference between A (S) and Ŝm pp in Frobenius norm is less than a given tolerance ε. S * = S 0, in which denotes positive semi-definite. The hermitian property is an inherent characteristic of spectral matrices, as well as the fact that its eigenvalues are non-negative. Unfortunately, it has no unique solution for block diagonal matrices of the type Ŝm

pp . An example is given here to illustrate this fact, and solid mathematical proof can be found in [START_REF] Candès | Exact matrix completion via convex optimization[END_REF]. Assume that the measured matrix Ŝm pp is a 2-by-2 complex hermitian matrix

1 a -1 a 1
, where a = a e jϕ is an unknown complex element and the diagonal elements are known measurements without error. By definition the objective matrix S is a rank-1 matrix and it is straightforward to find that a e jϕ = a -1 e -jϕ as a result of hermitian symmetry. Therefore, the solution of Eq. ( 10) is non-unique since it can be written {a : a = 1, ∀ϕ}. Thus, another constraint is needed as explained in the next section.

Enforcing spatial continuity of the acoustic field

The spectral matrix reflects the correlation between the Fourier coefficients of the measurements, however, it does not contain information on the position of the microphones. For example, two adjacent columns in the spectral matrix must become equal if the corresponding microphones are infinitely close to each other (nearly overlapping) due to the spatial continuity of the acoustic field; accordingly, this information must be incorporated in the current problem.

In this work, the spatial position of the microphones is encoded in the matrix completion problem by introducing an additional constraint. It is assumed that the measured pressure p with size MP × 1 can be decomposed onto a dimension reduced spatial basis Φ Φ Φ ∈ C MP×K p (K p < MP) with coefficients ϑ ∈ C K p ×1 as follows:

p = K p ∑ i=1 φ i (r)ϑ i = Φ Φ Φϑ , (11)
where φ i (r) is the i-th column of matrix Φ Φ Φ and ϑ i is the i-th entry of vector ϑ . In this thesis, the Fourier basis is chosen, meaning that the acoustical field is decomposed as a sum of "acoustical modes" [START_REF] Photiadis | The relationship of singular value decomposition to wave-vector filtering in sound radiation problems[END_REF]. In other words, a spatial structure is encoded. The coefficients ϑ are related to the measurements as:

ϑ = Φ Φ Φ † p, (12)
where † denotes the pseudoinverse of a matrix as a result Φ Φ Φ is generally not invertible. Then the smoothed pressure p can be represented as [START_REF] Pascal | Patch near-field acoustic holography: Regularized extension and statistically optimized methods[END_REF]

p = Φ Φ ΦΦ Φ Φ † p = Ψ Ψ Ψp, (13)
where Ψ Ψ Ψ = Φ Φ ΦΦ Φ Φ † is defined as a projection operator. Now, let S . = E{pp * } and S . = E{ p p * }. Therefore,

S = Ψ Ψ ΨS S SΨ Ψ Ψ * . (14)
Equation ( 14) may be rewritten as Ψ Ψ ΨS S SΨ Ψ Ψ * -S F ≤ ε where ε is the representation error between Ψ Ψ ΨS S SΨ Ψ Ψ * and S owing to the Φ Φ Φ truncation. It is important that Ψ Ψ Ψ imposes a specific structure that encodes the information on microphone positions into the matrix S. This can be seen as another constraint to be added to Eq. [START_REF] Michel | History of acoustic beamforming[END_REF]. It is noted that the Ψ Ψ Ψ is an orthogonal projection operator, which can be decomposed as U U U I K p 0 0 0 U * , where UU * = I, I K p is the identity matrix with n elements, and U is the modal matrix in the eigenvalue decomposition (EVD) of Φ Φ Φ. Then Eq. ( 13) can be written as

p = Ψ Ψ Ψp = K p ∑ i=1 (p * u i ) * u i . ( 15
)
This means that only the p * u i i = 1...K p part corresponding to the eigenvalues equal to one are kept by this projection operator, and this is exactly the property of orthogonal projection operator. It is important to emphasize that the role of the orthogonal projection operator Ψ Ψ Ψ is to smooth the spectral matrix and to ensure the spatial continuity of the acoustical field.With the intension of defining the spatial basis Φ Φ Φ in this projection basis, it is firstly denoted that K is the dimension of basis to represent the sound sources and MP is the total number of microphones by non-synchronous measurements. Then the Singular value decomposition (SVD) of the propagation function H ∈ C MP×K between the array of non-synchronous measurements and sound sources is H = UΣ Σ ΣV * ( * denotes conjugate transpose), where Σ Σ Σ, U ∈ C MP×MP , and V ∈ C K×K are separately the singular values (a diagonal MP × K matrix with non-negative real numbers on the diagonal), left-singular vectors and right-singular vectors of H. The physical meaning of the left-singular vectors are the "acoustical modes" to represent the acoustic field, and the right-singular vectors are the basis to represent the sources [START_REF] Leclère | Acoustic imaging using under-determined inverse approaches: Frequency limitations and optimal regularization[END_REF][START_REF] Photiadis | The relationship of singular value decomposition to wave-vector filtering in sound radiation problems[END_REF]. Eventually, spatial basis Φ Φ Φ is defined as Φ Φ Φ = Ext(U, K p ), where the operator Ext extracts the first K p columns of U. K p is chosen as K p = M √ P, which is a "rule of the thumb" as explained in Ref [START_REF] Antoni | Reconstruction of sound quadratic properties from non-synchronous measurements with insufficient or without references: Proof of concept[END_REF].

Low rank model of acoustic field with continuity constraint

The nuclear norm defined by a sum of the singular values of spectral matrix is an alternative relation of rank minimization. Thus, the data missing spectral matrix completion problem can

Liang Yu, Weikang Jiang etc.

be then formulated as the following constrained optimization problem:

minimize S S * subject to A (S) -Ŝm pp F ≤ ε Ψ Ψ ΨSΨ Ψ Ψ * -S F ≤ ε S * = S 0, (16)
This optimization problem can be solved by the following FISTA algorithm [START_REF] Yu | Acoustical source reconstruction from nonsynchronous sequential measurements by fast iterative shrinkage thresholding algorithm[END_REF]. The operator shrink is calculated as Udiag(max{σ 2 iλ k µ, 0})U * , where max{σ 2 iλ k µ, 0} takes the maximum value between 0 and σ 2 iλ k µ, U and diag(σ 2 i ) are the elements of the Eigenvalue Decomposition (EVD) of matrix G k . A continuation technique is used here to tune the regularization parameter λ . It can be seen as a method to solve a sequence of problems ( 16) defined by an increasing sequence {λ k }: (i) λ 0 as an initial regularization parameter and λ d as final regularization parameter are chosen firstly; (ii) a high starting value of λ implies stronger regularization and the optimal manifold has a lower dimension and is easier to identify; (iii) a coarse to fine search is achieved in the solution space by decreasing the regularization parameter λ step by step (with a ratio η). The calculation procedure of FISTA is illustrated in Algorithm 1. The algorithm is stopped by the following criterion:

[ Ŝm pp ] i j -[G k+1 ] i j F [ Ŝm pp ] i j F ≤ SC, i j ∈ ϒ ϒ ϒ, (17)
where SC is a constant value that is chosen by the user, [] i j denotes the (i, j) element of a matrix and ϒ ϒ ϒ the positions of fixed block diagonal entries.

Algorithm 1 Nuclear Norm minimization by FISTA 1: Starts with G 0 = S 0 = 0 ∈ C MP×MP , t 1 = 1; µ is step size, λ 0 is an initial regularization parameter and λ d is the final regularization parameter, and η is the radio to decrease λ k-1 in each step. 2: While λ k ≥ λ d do 3: For 1 : N m (N m is the maximum iteration steps for each regularization) In order to design a far field array, the array resolution and spatial aliasing must be taken into account. Array resolution specifies how well an array is able to resolve the direction of propagation. It is typically specified in terms of a 3 dB down point from the peak at which the spatial filter formed by an array processing algorithm. The width of the spatial filter at the 3 dB down point is the beamwidth, which determines the array resolution of conventional beamforming. Spatial aliasing is an undesirable side effect from spatially under-sampling the aperture of an array. The primary challenge with array designs is that they must work well for a broad range of frequencies. The lowest frequencies of interest drive the array aperture spiral array at an extreme range of frequency are shown in Fig. 3. Figure 3 (a) shows the sidelobe is 10 dB smaller than the peak at 2.8 kHz, therefore it can be ignored. Figure 3 (b) shows that the sidelobe appeared in dynamic range of 10 dB at 2.9 kHz, which causes the spatial aliasing in the imaging results of beamforming. To avoid spatial aliasing, the highest working frequency of current array is 2.8kHz. With the purpose of determine the working frequency of lowest, the beamwidth of the Dougherty spiral array is shown in Fig. 3. The relative positions between the sources and array is shown in Fig. 2(b), and with the sources at a spacing of 0.76m, the lowest frequency of working of the current array is 2kHz. This is verified by the beamwidth of 0.36 along kx(x) and 0.3 along ky(y) at 2kHz, which are shown in Fig. 4. The center positions of the non-synchronous measurements are shown in Fig. 2 (c), 9 times non-synchronous measurements are applied in the current setup. Figure 6 shows the results of acoustic imaging in non-synchronous measurements, simultaneous and single measurements respectively. It is observed that the non-synchronous measurements broaden the frequency range of working from 2 -2.8 kHz to 1.6 -4.5 kHz. The imaging results of single measurement (with only an array of 29 elements) are non-distinguishable when it is applied outside the frequency range of 2 -2.8 kHz. However, the imaging results of nonsynchronous measurements (with only an array of 29 elements and moving 9 times sequentially) keep distinguishable in the range of 1.6 -4.5 kHz, which are comparable with the simultaneous measurements of array of 261 elements.

4: G k = G k-1 -µ(A (G k-1 ) -Ŝm pp ) 5: Sk = shrink(G k , λ k µ) 6: S S S k = Ψ Ψ Ψ S S S k Ψ Ψ Ψ * 7: t k+1 = 1 2 (1 + 1 + 4t 2 k ) 8: G k+1 = S k + t k -1 t k+1 (S k -S k-
(c) (d) (e) (f) (g) (h) (i) (j) (k) (l)

CONCLUSION

In this paper, non-synchronous measurements beamforming is formulated as a spectral matrix completion problem which is solved by minimizing the nuclear norm of a spectral matrix subject to measurements fitting, Hermitian symmetry, and spatial continuity of the sound field. The proposed non-synchronous measurements beamforming extended greatly the working frequency range of a single array measurement. The simulation results show the non-synchronous measurements beamforming extended the working frequency range of given small phased array from 2 -2.8 kHz to 1.6 -4.5 kHz with only 9 times non-synchronous measurements. It is verified that the non-synchronous measurements improve the frequency range of working with limited number of microphones, which improve the capability of array and reduce the budget of measurements greatly in the aspects of number of microphones and channels of acquisition system.

Figure 1 :

 1 Figure 1: measurement mode and spectral matrix: non-synchronous measurements (left) versus simultaneous measurement (right).

Figure 2 :

 2 Figure 2: The positions of microphone array: (a) a 29-element Dougherty spiral array, (b) the relative position of sources and spiral array, (c) the center positions of 9 times nonsynchronous measurements.

Figure 3 :

 3 Figure 3: Far-Field array pattern: (a) 2.8kHz, (b) 2.9kHz.

Figure 4 :

 4 Figure 4: Far-Field array pattern at 2kHz: (a) profile representation along the kx(x) direction, (b) profile representation along the ky(y) direction.

Figure 5 :

 5 Figure 5: Sound source imaging result at different frequencies of single measurement of (a) 1.8kHz, (b) 2kHz, (c) 2.8kHz, (d) 3kHz.

Figure 6 :

 6 Figure 6: Sound source imaging result of (a) non-synchronous measurements of 9 times at 1.6kHz, (b) simultaneous measurements at 1.6kHz, (c) single measurements at 1.6kHz, (d) non-synchronous measurements of 9 times at 2kHz,(e) simultaneous measurements at 2kHz, (f) single measurements at 2kHz, (g) non-synchronous measurements of 9 times at 2.8kHz, (h) simultaneous measurements at 2.8kHz, (i) single measurements at 2.8kHz, (j) non-synchronous measurements of 9 times at 4.5kHz,(k) simultaneous measurements at 4.5kHz,(l) single measurements at 4.5kHz.

  1 ) 9: End 10: If Stopping criteria ≤ SC, 11: break 12: End for if 13: λ k = max(ηλ k-1 , λ d ) 14: Go to step 2.

th Berlin Beamforming Conference 2018 Liang Yu, Weikang Jiang etc.

ACKNOWLEDGEMENTS

The authors would like to thank the continuous discussion with Ms.Wenhuang Wang and Jianzheng Gao (Shanghai Jiao Tong University). This work was supported by the National Natural Science Foundation of China (Grant No. 11704248).


Paul Verlaine Gakne
email: pvgakne@ucalgary.ca

Kyle O'keefe

Tackling the Scale Factor Issue in a Monocular Visual Odometry Using a 3D City Model

His research focuses on vision

In open sky environment with good satellite visibility, the Global Navigation Satellite Systems (GNSS) provide seamless and accurate solutions. However, their performance significantly degrades in harsh environments such as urban canyons. To overcome this limitation, various integration strategies that combine GNSS with other sensors (e.g., inertial navigation systems -INS, odometers, radar, barometer, LiDAR, altimeter, cameras etc.) to ensure satisfactory position solutions in terms of accuracy, availability and continuity in almost all environments have been thoroughly studied with more or less success. In general, the use of camera systems (thanks to the concept of the Visual Odometry -VO) is able to provide an accurate position estimation that can in turn be integrated with solutions obtained from other sensors.

In fact, the VO aims at recovering the motion (position and orientation) of a platform (e.g., vehicle, robot as International Technical Symposium on Navigation and Timing (ITSNT) 2018 13-16 Nov 2018 ENAC, Toulouse, France shown in Figure 1) by exploiting the images captured by a camera rigidly attached to the platform. Both monocular (i.e., one camera) and stereo (i.e., two cameras) systems can be used for VO implementation. Vision-based positioning for pedestrian [START_REF] Marouane | Visual Odometry for Pedestrians Based on Orientation Attributes of SURF[END_REF] or vehicular applications [START_REF] Aumayer | Ultra-tightly Coupled Vision/GNSS for Automotive Applications[END_REF] has been previously studied by exploiting various concepts such as the rigid-body motion [START_REF] Ma | An Invitation to 3-D Vision: From Images to Geometric Models[END_REF], and the use of feature points [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF]. When more than one camera is used, the VO outputs both the (change in) rotation and the true translation magnitude directly by triangulating the feature points and obtaining the depth information. However, with a single camera system, the scale cannot be directly measured. This means that the translation of the platform can be estimated accurately but up to a scale. This is one of the primary sources of error in VO-based (change of) position estimation and the main bottleneck that prevents monocular-based VO from reaching accuracies comparable to stereo-based VO. The scale factor problem in monocular VO is illustrated in Figure 2 (results for a controlled indoor area). It is clear that the VO provides an accurate solution but up to a scale only (causing the mismatch observed at the "End"). To resolve the scale ambiguity and drift, the general trend is to fuse the image measurements with other sensors. Specifically, Inertial Navigation Systems (INS) are widely used [START_REF] Spaenlehauer | A Loosely-Coupled Approach for Metric Scale Estimation in Monocular Vision-Inertial Systems[END_REF]; [START_REF] Ji | Visual-Inertial Combined Odometry System for Aerial Vehicles[END_REF]; [START_REF] Gabriel | Fusion of IMU and Vision for Absolute Scale Estimation in Monocular SLAM[END_REF]) as well as GNSS sensors [START_REF] Soloviev | Integration of GPS and vision measurements for navigation in GPS challenged environments[END_REF]; [START_REF] Lim | Integration of GPS and monocular vision for land vehicle navigation in urban area[END_REF]; [START_REF] Gakne | Tightly-Coupled GNSS/Vision Using a Sky-Pointing Camera for Vehicle Navigation in Urban Areas[END_REF]). [START_REF] Spaenlehauer | A Loosely-Coupled Approach for Metric Scale Estimation in Monocular Vision-Inertial Systems[END_REF] combined inertial measurements with monocular odometry in a loosely coupled approach. Regarding their visual odometry, the tracking stage of the ORB-SLAM algorithm (Mur-Artal, Montiel and Tardos, 2015) was employed while the Euler forward integration was used for the inertial measurement processing. The main advantage of using inertial sensors is their high measurement rate that can be used to rapidly compute the pose estimate. The scale factor is computed as the ratio of L 2 -norm [START_REF] Horn | Norms for Vectors and Matrices Ch. 5 in Matrix Analysis[END_REF] of translation vectors of the camera position given by the integration of INS measurements in the world coordinate frame and the L 2 -norm of the translation vectors of the camera position given by the VO. [START_REF] Gabriel | Fusion of IMU and Vision for Absolute Scale Estimation in Monocular SLAM[END_REF] fused the visual and inertial data to determine the scale factor unobserved in the visual framework. Two approaches were used to estimate the scale factor in a monocular visual Simultaneous Localization and Mapping (SLAM). The first method consisted of making use of the spline fitting task [START_REF] Jung | Camera trajectory estimation using inertial sensor measurements and structure from motion results[END_REF] and the second made use of a multi-rate Extended Kalman Filter (EKF) that was embedded in a Parallel Tracking and Mapping (PTAM - Klein and Murray, 2007). Both approaches provided accurate results at reasonable processing times. [START_REF] Ji | Visual-Inertial Combined Odometry System for Aerial Vehicles[END_REF] proposed a method aiming at operating aircraft in GPS-denied environments using visual odometry (with a camera pointed downward) and a high-accuracy INS. It was observed that the noise contained in the INS measurements can significantly affect the platform motion estimation resulting drift in the aircraft position estimate. To reduce the position estimation error, their method parametrizes feature points with their depth direction perpendicular to the ground. The depth of feature points is obtained in two different ways (by using an altimeter and a 2D laser). This leads to a slower drift because the position error resulting from the INS orientation noise is partially cancel by this process.

Although methods employing inertial sensors are effective and have advantages of lightweight, high measurement rate, immunity to RF interference etc. they are subject to the inertial sensors errors. INS drift rapidly when they do not get updates from other sensors such as GNSS. Other research focused on directly combining GNSS with vision to solve for the scale factor. [START_REF] Soloviev | Integration of GPS and vision measurements for navigation in GPS challenged environments[END_REF] resolved the scale factor ambiguity by integrating the Global Positioning System (GPS) carrier phase measurements with a vision system. Their method consisted of integrating GPS/vision in order to estimate the position as well as the orientation changes of the camera. The integration is realized by combining GPS carrier phase measurements with feature points that are extracted from images. The state relating the position vector with the feature point ranges is first defined, then the changes in the carrier phase measurements between two consecutive images are defined as function of the position change vector. The combination of these two terms allowed to unambiguously resolve the position changes, as well as the range estimates. A similar approach was presented in Gakne and O'Keefe (2018) but the pseudorange measurements are used in lieu of the carrier phase. [START_REF] Lim | Integration of GPS and monocular vision for land vehicle navigation in urban area[END_REF] proposed a method for augmenting GPS with a monocular camera for accurate and reliable positioning by combining the GPS measurements and the relative attitude obtained from the vision system (by the use of the vanishing point). They show that their method performs well in various environments and has results in better solutions than loosely coupled GPS/INS integration approach.

International Technical Symposium on Navigation and Timing (ITSNT) 2018 13-16 Nov 2018 ENAC, Toulouse, France However, in urban environments, GNSS measurements are subject to multiple source of error (mainly multipath) that can lead to position errors up to hundreds of meters.

In the Vision/INS case, the scale factor correction can be heavily degraded. For this reason, some researchers have integrated GNSS/Vision/INS or more sensors. [START_REF] Chu | Monocular Camera/IMU/GNSS Integration for Ground Vehicle Navigation in Challenging GNSS Environments[END_REF] integrated a monocular camera measurement with IMU and GNSS for land vehicle navigation in harsh environments based on an EKF. In this case, the translation magnitude is determined from GNSS measurements using differential GNSS technique.

The acquired image data is synchronized with GNSS and the baseline of the position change between two GNSS epochs coincides with the camera translation magnitude.

Ben Afia, Escher and Macabiau (2015) integrated even more sensors namely GNSS, IMU, vision and wheel speed sensor (WSS) for vehicular applications. In their implementation the INS estimation errors are corrected by the mean of GNSS, vision and WSS as well as motion constraints for land vehicles (using an error state EKF based on a closed loop configuration). The inertial measurement scale factor was modelled as a Gauss-Markov process.

Finally, another path of research has dealt with the monocular scale factor ambiguity without adding any other sensors by defining the road as a plane and using the vehicle height to obtain the scale information [START_REF] Choi | Resolving Scale Ambiguity for Monocular Visual Odometry[END_REF]; [START_REF] Kitt | Monocular Visual Odometry using a Planar Road Model to Solve Scale Ambiguity[END_REF]). However, these approaches still contain errors because the assumption that the road is a plane is only true up to an extent.

The objective of this paper is twofold: (i). to combine the information obtained from a 3D building model by matching synthesized 3D building model images and camera images and also with visual odometry; (ii). To calculate the scale factor from the 3D building model to increase the translation magnitude accuracy for monocular VO.

The remainder of this paper is organized as follows: Section 2 briefly presents the existing works related to this paper; Section 3 details the methodology employed in this paper. Section 4 presents the experiments, results and analysis. Finally, Section 5 concludes the paper.

RELATED WORK

A 3D city model can be defined as a digital representation of areas that depicts buildings, terrain surfaces and other features of a given city. Nowadays, digital representation of most of the major cities around the world is available whether commercially (e.g., from 3dcadbrowser.com) or freely (e.g., google.com/earth). Several researchers have made prior attempts to combine the VO with a 3D building model and/or with other sensors for navigation purposes. This works can be classified into several categories.

The first category uses the vision sensor alone to compute the scale factor by making assumptions regarding the camera height relative to the ground, defined as a plane. [START_REF] Kitt | Monocular Visual Odometry using a Planar Road Model to Solve Scale Ambiguity[END_REF] compensated the scale factor drift in a monocular VO by making use of constraints and using the knowledge of the camera height (assumed constant during the data collection) as well as by making assumptions about the environment and the camera orientation. The road is assumed to be planar and the camera pitch and roll equal to zero. The tracked feature points are assumed lying on the road plane and then the translation's magnitude is obtained by considering the height of the camera above the ground. A similar approach has been used later by [START_REF] Zhang | Robust Monocular Visual Odometry for a Ground Vehicle in Undulating Terrain[END_REF]; [START_REF] Choi | Resolving Scale Ambiguity for Monocular Visual Odometry[END_REF]; [START_REF] Song | High Accuracy Monocular SFM and Scale Correction for Autonomous Driving[END_REF]. However, it is obvious that in real-world scenarios, the road is not necessarily planar and the pitch and roll are non-zero in high dynamic maneuver scenarios for example (which can be a typical case for vehicular applications). Moreover, in the high dynamics case, even the height of the camera can significantly vary. However, it should be acknowledged that errors caused by such assumptions can be smaller for robots driven in indoor environment scenarios for example. [START_REF] Aqel | Estimation of image scale variations in monocular visual odometry systems[END_REF] proposed a method dealing with the camera height variation in order to improve the accuracy of the scale factor estimation. A downward pointing camera installed under the vehicle was used. Their method consists of setting reference points on the images by the mean of two laser pointers. The scale variation is then estimated by monitoring the changes of distance between the two reference points after obtaining the height variation of the vehicle. [START_REF] Manolis | Accurate Scale Factor Estimation in 3D Reconstruction[END_REF] also note that measuring and/or estimating the height of the vehicle in order to estimate the scale factor is not sufficiently accurate and proposed three alternative techniques for estimating the translation magnitude pertaining to a three-dimensional reconstruction (two for stereo systems and one for monocular systems). In their approach, the camera pose is estimated from a single image by matching sets of 2D and 3D feature points and making the pose estimation refinement iteratively. To do this, they embedded a P3P solver [START_REF] Xiao-Shan | Complete Solution Classification for the Perspective-Three-Point Problem[END_REF] into the random sample consensus (RANSAC -Fischer and Bolles (1981)) framework.

Gräter, Schwarze and Lauer (2015) associated the vanishing point (VP) concept with the road definition as a plane (estimated using structure from motion techniques) to correct the scale drift for advanced driver assistance systems. The proposed method fits the reconstructed feature points and VPs and uses the least-squares optimization to refine the plane and obtain the scale.

A second group of research deals with the scale estimation in monocular systems by combining the vision system with GNSS. Gakne (2018); Gakne and O'Keefe (2018) modelled the scale in a tightly-coupled GNSS/Vision system and show that the scale factor drift can be reasonably modelled with a Gaussian random walk. The scale was observed from the previous and actual position estimates from VO. [START_REF] Li | Integration of Low-Cost GNSS and Monocular Cameras for Simultaneous Localization and Mapping[END_REF] later proposed an integration of GNSS and a monocular system for SLAM. Their approach fuses the measurements of each system using an optimization-based scheme. It

International Technical Symposium on Navigation and Timing (ITSNT) 2018 13-16 Nov 2018 ENAC, Toulouse, France outputs the absolute position and attitude of the driven car as well as the translation's scale factor.

The fusion of inertial systems and vision can be classified as the third research path for scale estimation in monocular applications. This is presented in [START_REF] Gabriel | Fusion of IMU and Vision for Absolute Scale Estimation in Monocular SLAM[END_REF] as described in the previous section. [START_REF] Weiss | Real-time metric state estimation for modular vision-inertial systems[END_REF] combined a monocular vision system with an inertial sensor equipped with a three-axis gyroscope and accelerometer. Their algorithm is built to operate in real time and is independent of the vision algorithm that is used to estimate the camera poses instead.

The fourth category includes works where more than two sensors are added. Ben Afia, Escher and Macabiau (2015) integrated the GNSS, IMU, monocular vision and the wheel speed sensors for navigation in urban canyons. The WSS velocity was particularly affected by an unknown scale that was then modelled as a constant (scale factor error). [START_REF] Won | Selective Integration of GNSS, Vision Sensor, and INS Using Weighted DOP Under GNSS-Challenged Environments[END_REF] evaluated the navigation performance of GNSS/INS/vision integration simulated data, where the number of observed satellites was decreased from three to one. Their method directly outputs absolute (scaled) solutions.

It is also worth mentioning here that a stereo system combined with GNSS sensors that allow direct observation of the scale factor has been extensively studied as presented in [START_REF] Fei | Tightly Coupled Stereo Vision Aided Inertial International Technical Symposium on Navigation and Timing (ITSNT) 2018 13-16 Nov 2018 ENAC[END_REF] and [START_REF] Aumayer | Ultra-tightly Coupled Vision/GNSS for Automotive Applications[END_REF]. Also, RGB-D cameras have been used for direct scale determination such as presented in Mur-Artal, Montiel and Tardos, (2015).

Finally, any other sensor or digital map that can provide useful information can be combined with a camera in order to address the scale factor issue. [START_REF] Zhang | Scale Estimation and Correction of the Monocular Simultaneous Localization and Mapping (SLAM) Based on Fusion of 1D Laser Range Finder and Vision Data[END_REF] proposed an integration strategy that fuses a monocular visual SLAM with a 1D-laser range finder to obtain the scale estimation and the drift correction. An analytical feasibility for estimating the scale factor by fusing the vision system outputs (derived based on the local dense reconstruction of image sequences) and the laser information is described in detail. [START_REF] Kai | Enhanced Monocular Visual Odometry Integrated with Laser Distance Meter for Astronaut Navigation[END_REF] proposed an integrated monocular VO/laser system that corrects the scale factor drift for astronaut navigation. The laser distance/range finder provides a distance to a single point that can be used for 3D scene reconstruction or scale estimation. The work presented in this paper is close to this category but instead uses an existing 3D city model to obtain the distance from a driven vehicle to a point.

One aspect presented in this paper is the skyline-based positioning. It consists of matching images obtained from a camera with images synthesized from a 3D city model. This technique has been presented in our previous works [START_REF] Petovello | Skyline Positioning in Urban Areas Using a Low-cost Infrared Camera[END_REF][START_REF] Gakne | Monocular-based pose estimation using vanishing points for indoor image correction[END_REF]. In this paper particular attention is given to a combination of a VO (realized from a special setup with a sky-pointing camera) and an existing 3D building model. The next section provides details on the method used for the scale factor drift correction proposed in this paper.

METHODOLOGY

This section elaborates the proposed method for scale factor estimation and drift correction. This is principally based on the use of a 3D city model. The proposed work is mainly subdivided in two subsections: the visual odometry and the skyline-based positioning.

Monocular Visual Odometry

VO can be defined as the computation of the camera motion from monocular or stereo image sequences.

Herein, a feature points-based VO is developed from image sequences from an upward-facing camera and follows the steps bellow:

• Image acquisition and correction: this step consists of acquiring images from a monocular camera and rectifying the image distortion via the calibration process. The calibration process also allows the determination of the camera intrinsic parameters such as the focal length, the principal point coordinates and the skew coefficient between the image axes; • Feature detection, description and matching: feature points are detected and described by using the Oriented FAST and Rotated BRIEF (ORB - [START_REF] Rublee | ORB: An Efficient Alternative to SIFT or SURF[END_REF]. The description of this algorithm as used in this paper is presented in Gakne and O'Keefe (2018); • Motion estimation: from the tracked feature points and their displacement on the images, the relative motion of the camera is computed.

An example of feature detection, matching and outlier rejection is given in Figure 3 where the red circles represent feature points that are detected but not used (outliers); Green circles are those that are detected properly matched between consecutive image frames. The established matches are illustrated by the green lines. The step-by-step VO algorithm formulation is given in Gakne (2018) and summarized in Figure 4. Having explained the visual odometry, the next section presents the 3D building model based positioning.

3D Building Model-based Positioning

A skyline is literally defined as an outline of land and buildings defined against the sky. If accurately generated, such information can be used as a fingerprint that can uniquely describe a given city/area. The skyline-based positioning used in this paper is presented in detail in [START_REF] Gakne | Monocular-based pose estimation using vanishing points for indoor image correction[END_REF] and summarized in Figure 5.

The localization problem in this case is solved by using the skylines obtained from a 3D city model and from an upward-facing camera. At each vehicle location, the skyline is computed from the 3D building model (see Figure 6) over a range of azimuth. The skyline is defined as:

( )   ,, , , 0.
.

1 p p j p j jN = = - ς ε h (2)
where j and 720 N = are the azimuth index and number respectively (i.e., an azimuth resolution of 0.5 is used). ε is the highest elevation angle of the obstructing surface and h the corresponding height.

One of the main challenges of the skyline-based positioning is to accurately segment the camera image.

The difficulty arises because of the lighting variations, weather conditions and the facade of the buildings (e.g., buildings with glass). In this work, the Flood-fill algorithm described in Gakne (2018) is used for the camera image segmentation. Readers can refer to the provided reference for details the image synthesization from the 3D city model.

International Technical Symposium on Navigation and Timing (ITSNT) 2018 13-16 Nov 2018 ENAC, Toulouse, France The comparison/matching between the camera images and the database images (obtained from the 3D building model) is done by the mean of a similarity metric. Herein, the cross-correlation coefficient ( CC ) is used. It is defined as: The similarity metric obtained as in Figure 7 gives the location of the vehicle on the travelled path. In order to improve the matching accuracy between the observed (camera) and the 3D building model images the vanishing points have been used to determine the camera pitch and roll and to rectify images accordingly as in [START_REF] Gallagher | Using vanishing points to correct camera rotation in images[END_REF] and detailed in [START_REF] Gakne | Monocular-based pose estimation using vanishing points for indoor image correction[END_REF]. In order to be succinct, this process will not be repeated in this paper. Having the skyline defined as in Equation ( 2) the next step is to compute the slant distance that will be used to compute the scale factor.

    ( ) () () 3 , 1 ( , ) ( , ) ( , )

Slant Distance and Scale Factor Computation

In this paper the highest point in view from the vehicle location is chosen and the slant distance to this point from the vehicle is computed. From Equation (2), the slant distance is computed as: slant sin( )

= h d ε (4)
where the parameters h and ε are defined as in Equation

(2). The computation of these parameters follows that presented in [START_REF] Petovello | Assessment of Skyline Variability for Positioning in Urban Canyons[END_REF] as implemented in Gakne (2018).

The slant distance computed in this way is similar to the distance obtained from LiDAR data for example. As such, if carefully used, this information can be used to compute the scale factor ambiguity observed from a monocular system. The global scale factor is thus computed as: frames. The relative scale factor is estimated based in the actual and previous estimated platform position [START_REF] Gakne | Tightly-Coupled GNSS/Vision Using a Sky-Pointing Camera for Vehicle Navigation in Urban Areas[END_REF]. This is given as:

( ) ( ) ( ) 1 1 1 2 2 2 kk k k k kk s x x y y z z - - - = - + - + - (6)
where ( )

,, x y z represent the estimated platform's position and k is the image frame number.

The steps for correcting the scale factor drift are given as follow:

• Compute the scale factor as given in Equation ( 6) from frame-to-frame; • Check if the 3D building model solution is available. If yes, compute the slant distance as defined in Equation ( 4) then use this information to compute the global scale factor as in Equation ( 5)

These steps are summarized in the Algorithm 1.

Algorithm 1 Scale factor computation and correction

Input: ς Output: s Initialization: 1 r s  , 1 s  for ( 0 : 1 k -NumberFrame ) ( ) ( ) ( ) 1 1 1 2 2 2 kk kk kk r s x x y y z z - - - = - + - + - ; if (3DBM_Solution_available)  slant sin( ) = h d ε ; slant 3DBM VO g = - d s PP ; r g s  s ;   r ss  ; return s

EXPERIMENT AND RESULTS

The performance of the proposed scale drift correction is evaluated in this section.

Experiment Setup

An experiment consisting on a monocular system rigidly mounted on the top of the car driven in downtown Calgary, AB, Canada was conducted to test this other method (and is the same as used in [START_REF] Gakne | Tightly-Coupled GNSS/Vision Using a Sky-Pointing Camera for Vehicle Navigation in Urban Areas[END_REF]). A sky-pointing camera and a reference system (a SPAN LCI from NovAtel) are both mounted on the vehicle roof as depicted in Figure 8.

Results and Analysis

With this experimental data, the performance of VO when aided by a 3D building model can now be compared to the case where the building model is not available. For the comparison, the scale factor is first computed as in Equation ( 6) after estimating the current position for every system (reference, VO and VO aided 3DBM).

Then, the scale factor error is computed by subtracting each from the reference scale factor. The time series results are depicted in Figure 9.

Figure 9: Scale errors

The scale factor error obtained from the VO is overall noisier and larger than the VO aided by the 3DBM. However, it is still clear that there are peaks observed in the aided VO scale. These arise when the 3DBM solution is not available or significantly degraded. The final position output by each system (VO aided by the 3D building model and the pure VO) is depicted in Figure 10.

International Technical Symposium on Navigation and Timing (ITSNT) 2018 13-16 Nov 2018 ENAC, Toulouse, France This justifies the degraded solution at the upper turn (eastwest) for the VO aided 3D building model. This issue can be solved by employing a recently created/updated 3DBM. Also, due to the scale inaccuracy obtained from the VO-only case, it can be seen that at the turns, the VO solution drifts (went over the reference at the lower turn) while the translation magnitude is smaller than the reference at the upper corner.

The cumulative distribution function (CDF) of the horizontal position error is depicted in Figure 11.

Figure 11: CDF of the horizontal errors

This figure shows that position error for the VO-only scenario is bounded between -30 m to 50 m while for the VO-aided 3DBM, the position error is bounded between -10 m and 15 m. This clearly suggests that the position error introduced by the scale drift degrades the VO position solutions compared to the VO aided by the 3DBM. Results are similar about 10% of the time, suggesting that the proposed method improved the final solution up to 90% of the time.

The assessment of the influence of the scale factor drift on other metrics such as the velocity are left for future work.

It is also worth mentioning that the skyline-based positioning performs better in environments with tall buildings. It is thus important to consider a camera fieldof-view (FOV) such as it is possible to capture as much as possible surrounding buildings.

CONCLUSIONS

This paper presents a method of resolving the scale factor drift observed in monocular visual odometry systems by using a 3D city model to obtain the slant distance, similar to the information provided by LiDAR data. The proposed method was assessed by using data collected in downtown Calgary and compared with monocular visual odometry alone. The results have shown that the proposed approach improved the final solution 90% of the time compared to the VO-only.

Figure 1 :

 1 Figure 1: Navigating robot (credit: jpl.nasa.gov)

Figure 2 :

 2 Figure 2: Illustration of the scale factor issue in monocular visual odometry (adapted from Gakne and O'Keefe (2017))

Figure 3 :

 3 Figure 3: Feature point detection, description (ORB algorithm), matching and outlier removal (RANSAC). Left: image frame at time t; Right: image frame at time t + ∆t. Given two sets of M feature points represented byM

Figure 4 :

 4 Figure 4: Feature points-based visual odometry

Figure 5 :

 5 Figure 5: Skyline-based positioning flowchart

Figure 6 :

 6 Figure 6: Location of the vehicle indicated on the 3DBM as well as the ideal skyline synthesized forward and upward skylines at the same location

  the total number of pixel per image; p is the database image's position; is the Hadamard product of two matrices; image synthesized from the 3D city model.

Figure 7 :

 7 Figure 7: Similarity metric indicating the location of the vehicle on the road

  obtained from the 3D building model.Whenever the 3D building model solution is available, g s is used to correct the scale drift that is introduced in the relative scale factor over a certain number of image International Technical Symposium on Navigation and Timing (ITSNT) 2018 13-16 Nov 2018 ENAC, Toulouse, France

Figure 8 :

 8 Figure 8: Experiment Setup. Left: the vehicle; Right: top view of the vehicle

Figure 10 :

 10 Figure 10: Trajectories comparison (3DBM stands for 3D building model)

ACKNOWLEDGMENTS

The authors would like to thank Dr. Mark Petovello for his suggestions during the early development of this work.