A Leonard

F.-X Dupé

J.-L Starck

A Compressed Sensing Approach to 3D Weak Lensing

Keywords: Gravitational lensing: weak, Methods: statistical, Techniques: image processing, Cosmology: observations, Galaxies: clusters: general, Cosmology: large scale structure of Universe

Context. Weak gravitational lensing is an ideal probe of the dark universe. In recent years, several linear methods have been developed to reconstruct the density distribution in the Universe in three dimensions, making use of photometric redshift information to determine the radial distribution of lensed sources. Aims. In this paper, we aim to address three key issues seen in these methods; namely, the bias in the redshifts of detected objects, the line of sight smearing seen in reconstructions, and the damping of the amplitude of the reconstruction relative to the underlying density. We also aim to detect structures at higher redshifts than have previously been achieved, and to improve the line of sight resolution of our reconstructions. Methods. We consider the problem under the framework of compressed sensing (CS). Under the assumption that the data are sparse or compressible in an appropriate dictionary, we construct a robust estimator and employ state-of-theart convex optimisation methods to reconstruct the density contrast. For simplicity in implementation, and as a proof of concept of our method, we reduce the problem to one-dimension, considering the reconstruction along each line of sight independently. We also assume an idealised survey in which the redshifts of sources are known. Results. Despite the loss of information inherent in our one-dimensional implementation, we demonstrate that our method is able to accurately reproduce cluster haloes up to a redshift of z cl = 1.0, deeper than state-of-the-art linear methods. We directly compare our method with these linear methods, and demonstrate minimal radial smearing and redshift bias in our reconstructions, as well as a reduced damping of the reconstruction amplitude as compared to the linear methods. In addition, the CS framework allows us to consider an underdetermined inverse problem, thereby allowing us to reconstruct the density contrast at finer resolution than the input data. Conclusions. The CS approach allows us to recover the density distribution more accurately than current state-of-the-art linear methods. Specifically, it addresses three key problem areas inherent in linear methods. Moreover, we are able to achieve super-resolution and increased high-redshift sensitivity in our reconstructions.

Introduction

Weak gravitational lensing has become a powerful tool to study the dark universe, allowing us to place constraints on key cosmological parameters, and offering the possibility to place independent constraints on the dark energy equation of state parameter, w [START_REF] Levy | [END_REF][START_REF] Hoekstra | [END_REF]Munshi et al. 2008;Albrecht et al. 2006;[START_REF] Peacock | ESA-ESO Working Group on[END_REF][START_REF] Schneider | Weak Gravitational Lensing[END_REF]Van Waerbeke & Mellier 2003).

Until recently, weak lensing studies considered the shear signal, and recovered the mass distribution, in two dimensional projection (see [START_REF] Schneider | Weak Gravitational Lensing[END_REF], for a review of weak lensing). However, with improved data quality and wideband photometry, it is now possible to recover the mass distribution in three dimensions by using photometric redshift information to deproject the lensing signal along the line of sight (Simon et al. 2009;Massey et al. 2007a,b;[START_REF] Taylor | [END_REF].

Under the assumption of Gaussian noise, various linear methods have been developed to recover the 3D matter distribution, which rely on the construction of a pseudoinverse operator to act on the data, and which include a levels, with no scope for improvement by judicious choice of inversion or denoising method.

In this paper, we consider the deprojection of the lensing signal along the line of sight as an instance of compressed sensing, where the sensing operator models the line-of-sight integration of the matter density giving rise to the lensing signal. For simplicity in implementation, and as a proof of concept of our method, we consider only the line of sight transformation, reducing the 3D weak lensing problem to a one-dimensional inversion, with each line of sight in an image treated as independent. We note, however, that the algorithm presented here can be cheaply generalised to three dimensions.

We adopt a sparse prior on our reconstruction, and use a state-of-the art iterative reconstruction algorithm drawn from convex analysis, optimisation methods and harmonic analysis set within a compressed sensing framework. This enables us to find a robust estimator of the solution without requiring any direct prior knowledge of the statistical distribution of the signal. We note that the compressive sensing framework allows us to consider an underdetermined inverse problem, thereby allowing us to obtain higher resolution on our reconstructions than that provided by the input data.

This method produces reconstructions with minimal bias and smearing in redshift space, and with reconstruction amplitudes ∼ 75% of the true amplitude (or better, in some cases). This is a significant improvement over current linear methods, despite the adoption of a simplified, onedimensional algorithm. In addition, our method exhibits an apparent increased sensitivity to high redshift structures, as compared with linear methods. Our reconstructions do exhibit some noise, with false detections appearing along a number of lines of sight. However, these tend to be localised to one or two pixels, rather than coherent structures, and we expect improved noise control with a full three-dimensional implementation.

We note that we do not include photometric redshift errors in our simulations and, consequently, the simulations shown here should be considered to be idealised. Ma et al. (2006) have presented a method to account for such errors in lensing measurements. This method has been used by various authors in studies on real data (see, e.g. [START_REF] Simon | [END_REF], and is straightforward to implement in the algorithm presented here. A full treatment of photometric redshift errors, which is essential in order for the method to be useful on real data, will be presented in an upcoming work. This paper is structured as follows. In § 2, we outline the weak lensing formalism in three dimensions, and outline several linear inversion methods to solve the 3D weak lensing problem. We introduce our compressed sensing framework and describe our proposed algorithm in § 3. In § 4, we discuss practical considerations in implementing the method, and describe our simulated dataset. In § 5 we demonstrate the performance of our algorithm in reconstructing simulated cluster haloes at various redshifts. We conclude with a discussion of our results and future applications in § 6.

Throughout the text, we assume ΛCDM cosmology, with Ω Λ = 0.736, Ω M = 0.264, h = 0.71, σ 8 = 0.801, consistent with the WMAP-7 results (Larson et al. 2011).

3D Weak Lensing

The distortion of galaxy images due to the weak lensing effect is described, on a given source plane, by the Jacobian matrix of the coordinate mapping between source and image planes:

A = 1 -κ -γ 1 -γ 2 -γ 2 1 -κ + γ 1 , ( 1
)
where κ is the projected dimensionless surface density, and γ = γ 1 + ifl 2 is the complex shear. The shear is related to the convergence via a convolution in two dimensions:

γ(θ) = 1 π d 2 θ D(θ -θ )κ(θ ) , (2)
where

D(θ) = 1 (θ * ) 2 ,
(3)

θ = θ 1 + iθ 2
, and an asterisk * represents complex conjugation.

The convergence, in turn, can be related to the threedimensional density contrast δ(r) ≡ ρ(r)/ρ -1 by

κ(θ, w) = 3H 2 0 ΩM 2c 2 w 0 dw fK (w )fK (w -w ) fK (w) δ[fK (w )θ, w ] a(w ) , (4)
where H 0 is the hubble parameter, Ω M is the matter density parameter, c is the speed of light, a(w) is the scale parameter evaluated at comoving distance w, and

f K (w) =    K -1/2 sin(K 1/2 w), K > 0 w, K = 0 (-K) -1/2 sinh([-K] 1/2 w) K < 0 , (5)
gives the comoving angular diameter distance as a function of the comoving distance and the curvature, K, of the Universe.

If the shear (or convergence) data is divided into N sp redshift bins, and the density contrast reconstruction is divided into N lp redshift bins (where N sp is not necessarily equal to N lp ), we can write the convergence κ (i) on each source plane as

κ (i) (θ) N lp =1 Q i δ ( ) (θ) , (6)
where

Q i = 3H 2 0 Ω M 2c 2 w +1 w dw W (i) (w)f K (w) a(w) , (7)
and

W (i) (w) = w (i) 0 dw f K (w -w ) f K (w ) p(z) dz dw z=z(w ) . (8)
Thus, for each line of sight, equation ( 6) describes a matrix multiplication, encoding a convolution along the line of sight. It is the inversion of this transformation:

κ(z) = Qδ(z) , (9)
that is the focus of this paper. We note that the inversion of equation ( 2) can be straightforwardly performed on each source plane in Fourier space.

Linear Inversion Methods

We focus here on the methods presented in Simon et al. (2009) andVanderPlas et al. (2011). For a review of other linear methods, the reader is referred to Hu & Keeton (2002).

The three dimensional lensing problem is effectively one of observing the density contrast convolved with the linear operator R, and contaminated by noise, which is assumed to be Gaussian. Formally, we can write

d = Rs + ε, ε ∼ N (0, σ 2 ) , ( 10
)
where d is the observation, s the real density and ε the Gaussian noise.

The general idea behind linear inversion methods is to find a linear operator H which acts on the data vector to yield a solution which minimises some functional, such as the variance of the residual between the estimated signal and the true signal, subject to some regularisation or priorbased constraints.

The simplest instance of such a linear operation is an inverse variance filter (Aitken 1934), which weights the data only by the noise covariance, and places no priors on the signal itself:

ŝIV = [R † Σ -1 R] -1 R † Σ -1 d , (11)
where Σ ≡ nn † gives the covariance matrix of the noise. This method proves problematic when the matrices involved are non-invertible, such as when there are degeneracies inherent in the allowed solution. In order to make the problem invertible, some regularisation must be introduced. Simon et al. (2009) opt to use a Saskatoon filter (Tegmark 1997;Tegmark et al. 1997), which combines a Wiener filter and an inverse variance filter, with a tuning parameter α introduced that allows switching between the two. This gives rise to a minimum variance filter, expressed as:

ŝMV = [α1 + SR † Σ -1 R] -1 SR † Σ -1 d , (12)
where S ≡ ss † encodes prior information about the signal covariance, and 1 is the identity matrix. VanderPlas et al. (2011) have recently proposed a filter based on the singular value decomposition (SVD) of the inverse variance filter of Equation (11). Under this formalism, we can write

ŝIV = VΛ -1 U † Σ -1/2 d , (13)
where we have decomposed the matrix

R ≡ Σ -1/2 R = UΛV † , U † U = V † V = 1 and Λ is the square diagonal matrix of singular values λ i = Λ ii .
Their filtering consists in defining a cutoff value σ cut , and truncating the matrices to remove all singular values with σ i < σ cut . This effectively reduces the noise by removing the small singular values, which translate into large values in the inversion.

VanderPlas et al. note that the SVD decomposition is computationally intensive, and while they do describe a method to speed up the process, it may not be practical to use this method on large images.

Similar considerations must be made when using any of the linear methods described above, as these all involve matrix inversion, which is an O(N 3 ) process. While optimisations can be found, these methods become excessively time-and computer-intensive when large datasets are considered. This, in effect, limits the resolution attainable using these methods.

Moreover, as discussed extensively in Simon et al. (2009) andVanderPlas et al. (2011), these linear methods give rise to a significant bias in the location of detected peaks, damping of the peak signal and a substantial smearing of the density along the line of sight. The Compressed Sensing (CS) theory, described below, allows us to address the lensing inversion problem under a new perspective, and we will show that these three aspects are significantly improved using a non-linear CS approach.

Compressed Sensing Approach

Linear methods are easy to use, and the variance of each estimator is rather direct to compute. Furthermore, these methods generally rely on very common tools with efficient implementation. However, they are not the most powerful, and including non-Gaussian priors is difficult, especially when such priors imply non-linear terms. Obviously, using better-adapted priors is required for building a more robust estimator. In this paper, we adopt a compressed sensing approach in order to construct an estimator that exploits the sparsity of the signal that we aim to reconstruct. The estimator is modelled as a optimisation problem that is solved using recent developments from convex analysis and splitting methods.

Compressed sensing theory

We consider some data

Y i ( i ∈ [1, .., m]) acquired through the linear system Y = ΘX , ( 14
)
where Θ is an m×n matrix. Compressed Sensing [START_REF] Candès | [END_REF]Donoho 2006) is a sampling/compression theory based on the sparsity of the observed signal, which shows that, under certain conditions, one can exactly recover a k-sparse signal (a signal for which only k pixels have values different from zero, out of n total pixels, where k < n) from m < n measurements. Such a recovery is possible from undersampled data only if the sensing matrix Θ verifies the Restricted Isometry Property (RIP) (see Candès & Tao 2006, for more details). This property has the effect that each measurement Y i contains some information about all of the pixels of X; in other words, the sensing operator Θ acts to spread the information contained in X across many measurements Y i .

Under these two constraints -sparsity and a transformation meeting the RIP criterion -a signal can be recovered exactly even if the number of measurements m is much smaller than the number of unknown n. This means that, using CS methods, we will be able to far outperform the well-known Shannon sampling criterion.

The solution X of ( 14) is obtained by minimizing min

X X 1 s.t. Y = ΘX (15)
where the 1 norm is defined by X

1 = i | X i |.
The 1 norm is well-known to be a sparsity-promoting function; i.e. minimisation of the 1 norm yields the most sparse solution to the inverse problem. Many optimisation methods have been proposed in recent years to minimise this equation.

Further details about CS and 1 minimisation algorithms can be found in [START_REF] Starck | Sparse Image and Signal Processing[END_REF].

In real life, signals are generally not "strictly" sparse, but are compressible; i.e. we can represent the signal in a basis or frame (Fourier, Wavelets, Curvelets, etc.) in which the curve obtained by plotting the obtained coefficients, sorted by their decreasing absolute values, exhibits a polynomial decay. Note that most natural signals and images are compressible in an appropriate basis.

We can therefore reformulate the CS equation above (Equation ( 15)) to include the data transformation matrix Φ:

min α α 1 s.t. Y = ΘΦα (16)
where X = Φ * α, and α are the coefficients of the transformed solution X in Φ, which is generally referred to as the dictionary. Each column represents a vector (also called an atom), which ideally should be chosen to match the features contained in X. If Φ admits a fast implicit transform (e.g. Fourier transform, Wavelet transform), fast algorithms exist to minimise Equation ( 16).

One problem we face when considering CS in a given application is that very few matrices meet the RIP criterion. However, it has been shown that accurate recovery can be obtained as long the mutual coherence between Θ and Φ,

µ Θ,Φ = max i,k | Θ i , Φ k , |, is low (Candes & Plan 2010).
The mutual coherence measures the degree of similarity between the sparsifying basis and the sensing operator. Hence, in its relaxed definition, we consider a linear inverse problem Y = ΘΦX as being an instance of CS when:

1. the problem is underdetermined, 2. the signal is compressible in a given dictionary Φ, 3. The mutual coherence µ Θ,Φ is low. This will happen every time the matrix A = ΘΦ has the effect of spreading out the coefficients α j of the sparse signal on all measurements Y i .

Most CS applications described in the literature are based on such a soft CS definition. CS was introduced for the first time in astronomy for data compression (Bobin et al. 2008;Barbey et al. 2011), and a direct link between CS and radio-interferometric image reconstruction was recently established in Wiaux et al. (2009), leading to dramatic improvement thanks to the sparse 1 recovery (Li et al. 2011).

CS and Weak Lensing

The 3D weak lensing reconstruction problem can be seen to completely meet the soft-CS criteria above. Indeed, 1. the problem is undetermined as we seek a higher resolution than initially provided by the observations, which are noise-limited, 2. the matter density is seen to be sparsely distributed,

showing clusters connected by filaments surrounding large voids, 3. the lensing operator encoded by the matrix Q in Equation ( 9) (or equivalently, the combination P γκ Q, where P flˇe ncodes the convolution in Equation ( 2)) spreads out the information about the underlying density in a compressed sensing way.

To highlight point (3) above, Figure 1 shows the rows (top panel) and columns (bottom panel) of the transformation matrix, Q, encoding weak lensing along the line of sight. The top panel shows the lensing efficiency kernel, which reflects the sensitivity of a given source plane to the shearing effects of lenses at various redshifts. This is the broad convolution kernel of Equation ( 4). The bottom panel shows the effect of a discrete lens at a given redshift on source planes, and demonstrates that localised lenses give rise to effects that are non-local and affect all sources at redshifts greater than the lens redshift.

Sparsity prior

Sparse priors have been shown to be very useful in regularising ill-posed inverse problems (see [START_REF] Fadili | Proceedings of the International Conference on Image Processing[END_REF], and references therein). In addition, a sparse prior using a wavelet basis has been used in many areas of signal processing in astronomy, such as denoising, deconvolution and inpainting to recover missing data in images [START_REF] Starck | Sparse Image and Signal Processing[END_REF] . The idea underlying such priors is that there exists a dictionary in which a given dataset is sparsely represented. The dictionary used should therefore match as closely as possible the shapes of the structures that we aim to detect.

Many experiments, such as N-body simulations, have shown that the matter density in the universe is largely distributed in localised clusters, which are connected by thin filaments. Because structures in the Universe appear to be physically sparse, we may therefore assume that the matter density is sparse in a domain adapted for cluster-and curvelike structures. Such domains exist and can be constructed by gathering several well chosen transforms inside a dictio-nary (e.g. a combination of wavelets and curvelets). The dictionary should be chosen in order to match as closely as possible the type of structure we aim to recover (e.g. isotropic structures for wavelets, filaments for curvelets), and may be adapted to include structures that are not sparse in the direct domain, but which may be sparsely represented in an appropriate dictionary.

Regularisation using a sparsity constraint can be understood through a Bayesian framework. We assume that the distribution of the solution in the sparsifying dictionary has a Laplacian distribution. This is equivalent to constraining the 1 norm, which promotes sparsity. A Gaussian assumption, in contrast, constrains the 2 norm, and leads to the standard Wiener filtering, a Tikhonov solution or an SVD solution, depending on the way we constrain the solution.

For specific classes of inverse problem, it has even been shown that the sparse recovery leads to the exact solution of the problem (compressed sensing). Such a behaviour does not exist with any other prior than sparsity. Obviously, the sparse recovery will be optimal when the signal is sparse, in the same way that the Wiener filter is optimal when the signal (and noise) is Gaussian. Because the Laplacian assumption (in a appropriate space such as wavelets) is more applicable than the Gaussian distribution, restoration of astronomical data is generally much more efficient using sparsity.

This explains why wavelets have been so successful for astronomical image restoration/detection. For the reconstruction of clusters along the line of sight, we are in a perfect situation for sparse recovery since clusters are not resolved due to the bin size, and they can therefore be modelled as Dirac δ-functions. We therefore take Φ to be a δ-function dictionary. Clearly, in this case, the pixel domain is especially appropriate for sparse recovery. We believe that the sparse prior is a much better model for this kind of data compared to previous methods with implicit Gaussian assumptions. Our results in this paper appear to support this claim.

We note that the method presented here is somewhat similar to the point-source reconstruction method described in Hu & Keeton (2002), though they use 2 minimisation combined with a strict prior on the number of haloes present along a line of sight. Our method, in contrast, places no priors on the number of structures along the line of sight.

Problem statement

Under the CS framework, the reconstruction of the matter density amounts to finding the most sparse solution that is consistent with the data. There are many different ways to formulate such an optimisation problem, and we opt for the following:

min s∈R n Φ * s 1 s.t. 1 2 d -Rs 2 Σ -1 ≤ s ∈ C , (17)
where the term to minimise is a sparsity-penalty function over the dictionary coefficients. The second term in Equation ( 17) above is a data fidelity constraint, with Σ being the covariance matrix of the noise and the allowed distance between the estimation and the observation, while the final term forces the solution to have values inside a given interval, usually C = [-1, +∞[ n for matter overdensity.

Note that this latter constraint, encoding a hard minimum on the signal to be recovered, is not possible with linear methods, and is therefore an additional strength of our method. Enforcing such physical constraints on our solution helps to ensure the recovery of the most physically compelling solution given the data.

In the compressed sensing literature, we can find equivalent writing of the problem ( 17), but with a focus on the data fidelity. We prefer to directly seek for the sparsest solution for a given freedom (i.e. ε) on the distance to the observation, as such distance is usually easier to tune than an equilibrium between a regularization and the data fidelity terms.

Note this optimisation problem can be equivalently expressed in a Bayesian framework, assuming Gaussian noise and a Laplacian distribution of the dictionary coefficients.

In order to solve (17), we use the primal-dual splitting method of Chambolle & Pock (2011), and is described in full in Appendix A. The algorithm is iterative, and effectively splits the problem into two parts, applying the two constraints in equation ( 17) separately.

On each iteration, the estimate of the reconstruction is compared with the data, and the data fidelity constraint (the second constraint in equation ( 17)) is applied by projection of the residual onto an 2 ball of radius . Independently, the estimate of the solution is projected onto the sparsifying basis, and sparsity is imposed through soft-thresholding, which minimises the 1 norm of the basis coefficients. The threshold level is set by the parameter λ, which aids in controlling the noise. The noise modelling is discussed in section 4.3, whilst the algorithm, and all practical considerations related to its implementation, are discussed in detail in Appendix B. Figure 2 shows a simplified schematic of the algorithm used.

Implementation of the Algorithm

We first consider the simulated data to be used in the remainder of this paper, before discussing some practical considerations important when implementing the algorithm.

Cluster Simulations

In order to test our method, we need to simulate a realistic data set. To this end, we consider here a fiducial survey with a background galaxy number density n g = 100 arcmin -2 distributed in redshift according to

p(z) ∝ z 2 e -(1.4z/z0) 1.5 , (18)
with z 0 = 1.0 (Taylor et al. 2007;Kitching et al. 2011), and the distribution truncated at a maximum redshift of z max = 2. Figure 3 shows this probability distribution, normalised arbitrarily.

We take the intrinsic dispersion in shear measurements to be σ γ = 0.2, and consider a field of 1 • × 1 • divided into a grid of 60 × 60 pixels. These parameters are chosen to mimic the data quality expected from next-generation surveys such as Euclid [START_REF] Universirty Refregier | [END_REF] and LSST (Ivezic et al. 2008;LSST Science Collaboration et al. 2009).

In Apply soft threshold at level λ.

Apply the inverse basis transform

Apply transform Σ -1/2 Q to solution estimate.

Apply data fidelity constraint: mean square error < .

Apply transpose transform Q † Σ 1/2 . Add the two components; apply positivity constraint. lensing signal within each source redshift bin, and Gaussian noise is added, scaled appropriately by the number density of galaxies within that bin.

Reconstructions in 1D

As noted in § 2, the 3D weak lensing problem can be reduced to a one-dimensional problem, by taking as our data vector the (noisy) lensing convergence along each line of sight, which is related to the density contrast through Equation ( 9). Therefore, we take d = κ ij (z) and R = Q, and consider each line of sight in our images independently. Further, as discussed previously, we take Φ to be a δfunction dictionary.

In our simulations, clusters are placed into a region where the mean density in the absence of the cluster is equal to the mean density of the Universe at that redshift. In other words, δ is constrained to be greater than zero in all our simulations. Therefore, the projection onto the convex set C in algorithm 2 applies a positivity constraint at each iteration.

Clearly, a one-dimensional implementation throws away information, as we do not account at all for the correlation between neighbouring lines of sight that will arise in the presence of a large structure in the image; however, reducing the problem to a single dimension is fast and easy to implement, and allows us to test the efficacy of the algorithm using a particularly simple basis function through which we impose sparsity. A fully three-dimensional treatment of the problem, with more accurate noise modelling (see below) will be the subject of a future work.

However, the algorithm used is entirely general; therefore, with appropriate choice of a three-dimensional basis set and taking d = γ(θ, z) and R = P flˇQ , one can implement this algorithm as a fully three-dimensional treatment of the data with no modification to the algorithm itself.

Noise Modelling and Control

Noise Model for the Data

We assume that the redshifts of the sources are known exactly, so there is no correlation between the noise in each source bin. Therefore, the covariance matrix of the noise along the line of sight is diagonal, with

Σ ii = σ 2 g (z i ) = σ 2 γ n g (z i )A pix , (19)
where A pix is the pixel area, n g (z i ) is the number density of sources in the bin at redshift z i , and σ γ is the intrinsic dispersion in galaxy ellipticity, taken throughout to be 0.2. This covariance matrix is used in the evaluation of the data fidelity constraint in our algorithm above. Note that the covariance matrix is only diagonal if the galaxy redshifts are known exactly. In practice, photometric redshift errors mean that each redshift slice in the data is likely to be contaminated with a few galaxies whose redshift error bars overlap with neighbouring redshift bins. In this case, the covariance matrix will have additional, non-diagonal elements that are non-zero. This is straightforward to model, however, for the chosen method of photometric redshift estimations, and our algorithm is entirely general regarding the form of the covariance matrix. Therefore, the problem of photometric redshift errors is readily tractable in our method, and will be presented in a future work.

Noise Control in the Algorithm

The noise in the reconstruction is controlled and suppressed by two parameters in the algorithm described in Figure 2 and Appendix B. The first, and most important of these parameters is the data fidelity control parameter, .

This parameter controls how well the data are fit by the reconstruction, with = 0 implying a perfect fit to the data, which is not possible in the presence of noise. Figure 4 demonstrates the effect of varying in the reconstruction of two lines of sight from our simulations. Clearly, when is small, the algorithm attempts to fit each data point more closely which, in the presence of noise, can result in overfitting of the data (as seen in line of sight 1) and hence false detections along the line of sight. On the other hand, a large may result in a solution that is not a good fit to the data (as seen in line of sight 2).

The second parameter used to control the noise is the soft threshold parameter λ, which is used in the algorithm to impose the sparsity prior. A threshold set excessively high will result in a null solution, whilst a threshold set fairly low will allow for more false detections of noise peaks along a given line of sight. The appropriate value for this threshold should be related to the expected fluctuations in the density contrast resulting from noise variations, and should scale with the signal to noise in the image.

Note that while strongly affects the accuracy of the estimation in reproducing the underlying density contrast, λ simply affects the sparsity of the solution. In other words, changing λ will not greatly affect the reconstructions of true density peaks, but may affect the number of false detections and noise peaks seen. Also note that a thresholding λ does not imply that density peaks with δ < λ will not be detected, as soft thresholding is only applied to one part of the estimate of the solution.

Results

Comparison with Linear Methods

Firstly, to demonstrate the effectiveness of our method, we compare our method directly with the linear methods of Simon et al. (2009) andVanderPlas et al. (2011). As these linear methods are only defined for N lp ≤ N sp , we consider the case of N lp = N sp = 20, and generate a single cluster halo at a redshift of z cl = 0.25 following an NFW halo profile with M 200 = 10 15 M and c = 3.

In the SVD method of VanderPlas et al. ( 2011), we take

v cut = 1- n i=1 σ 2 i / nmax i
σ 2 i = 0.01, and in the transverse and radial Wiener filtering methods of Simon et al. (2009), we take the tuning parameter to be α = 0.05 in both cases. For our CS approach, we take the soft threshold parameter λ = 8, and the data fidelity control parameter = 3.

Note that while the linear methods take d = γ(θ, z), R = P flˇQ , our method takes d = κ(θ, z), R = Q as before. The noise levels in each case are identical.

The results are presented in figures 5 and 6. Figure 5 presents the 2D projections of the reconstructions, computed by integrating the reconstruction along each line of sight, and the 1D reconstructions along the four central lines of sight. In the 3D renderings of Figure 6, the reconstructions from our method, the SVD method and the radial Wiener filter method are thresholded at δ = 3 (i.e. the plot only shows δ rec ≥ 3), and each is smoothed with a Gaussian of width σ = 0.7 pix in all three directions. The reconstruction from the transverse Wiener filter method is heavily damped with respect to the amplitude of the density contrast; a threshold of δ = 5 × 10 -6 is chosen in this case, and no smoothing is applied as the reconstruction already shows a very smooth distribution.

The SVD method appears, in the 1D plots, to identify the correct redshift of the cluster, with a small amount of line of sight smearing, but the plots show a prominent highredshift peak along the line of sight. We note that it may be possible to remove this false detection by raising v cut , but at the cost of increased line of sight smearing (see VanderPlas et al. 2011). The two-dimensional projection consisting of the integrated signal along each line of sight is seen to be more noisy than our results. The three-dimensional rendering shows that the SVD method does well at identifying and localising the cluster, but the resulting reconstruction does suffer from widespread noise at a moderate level.

The radial and tangential Wiener filter methods show very little noise in the 2D projections, and the radial Wiener filter shows very little smearing of the reconstruction along the line of sight. However, neither Wiener method recovers correctly the redshift of the cluster. While the transverse Wiener filter shows very little noise, it exhibits a broad smearing in both the radial and transverse directions, and a heavy damping of the amplitude of the reconstruction. The results obtained using the radial Wiener filter are considerably better, with very little noise seen in the reconstruction, and with the cluster seen to be well-localised in the radial direction. However, the amplitude of the cluster reconstruction is a factor of ∼ 2.5 smaller than the input density.

Our results are seen to suffer from several prominent, pixel-scale false detections along noisy lines of sight not associated with the cluster. However, along the four central lines of sight an excellent correlation between the input density contrast and our reconstruction is seen. One Fig. 6. 3D rendering of the reconstruction of a z cl = 0.25 cluster using our method (top), the SVD method (2nd row, v cut = 0.01), radial Wiener filter method (3rd row, α = 0.05) and the transverse Wiener filter method (bottom, α = 0.05) as described in the text. line of sight exhibits a prominent high-redshift false detection; however this does not appear in the remaining three lines of sight and the overall amplitude of the reconstruction is around 75% of the true value. The three-dimensional rendering demonstrates that the noise in our reconstruction shows very little coherent structure (i.e. tends to be restricted to isolated pixels), and is largely low-level. Moreover, the cluster is incredibly well-localised in redshift space, with the smearing seen in the figure primarily arising from the applied smoothing.

Reconstruction with Noisier Data

The simulations described in § 4.1 represent a rather optimistic set of data. Realistically, one might expect a much smaller number density of galaxies, so it is important to consider how our method fares with noisier data. To this end, weak lensing data were simulated for the same cluster as described above, but with n g = 30 galaxies per square arcminute. This represents a factor of ∼ 1.8 reduction in the signal to noise of our data. In order to account for this reduction in signal to noise, we increased our soft threshold parameter λ by a similar factor, taking λ = 15, whist we again used = 3.

The projections and 3D reconstruction obtained using this noisier data is shown in Figure 7. The 2D projection shows a noisier reconstruction, with more false detections, and the central cluster appears less extended. This is to be expected, given the lower signal to noise in the data. The line of sight plots again show more noise, with the peaks of some lines of sight being shifted slightly from the true redshift. However, they are well localised around the true peak, and their amplitudes match quite well with the true underlying density.

The 3D rendering again shows a well-localised peak at the location of the cluster, albeit slightly extended along the line of sight. Note also that the false detections continue to appear as single pixel-scale detections, rather than as extended objects. As before, we believe that a fully threedimensional implementation of our algorithm will reduce the number of such false detections by searching for coherent structures of larger angular extent than a single pixel.

While improvements can be made to the noise model used, it is clear that, with appropriate choice of parameters, effective noise control can be obtained with noisy data, and a very clean reconstruction obtained using our method. Here again, our method is seen to improve on the bias, smearing and amplitude damping problems seen in the linear methods presented above, and the results at this noise level are of a comparable or better standard than the results from linear methods at higher signal to noise, as presented above.

Improving the Line of Sight Resolution

Given the success of our method at reconstructing lines of sight at the same resolution as the input data, it is interesting to consider whether we are able to improve on the output resolution of our reconstructions. It is also worthwhile to test our ability to detect clusters at higher redshifts than that considered above, given that compressed sensing is specifically designed to tackle underdetermined inverse problems. Indeed, noise-free simulations suggest that a resolution improvement of up to a factor of 4 in the redshift direction may be possible with this method.

Therefore, we generate clusters as before, with our data binned into N sp = 20 redshift bins, but aim to reconstruct onto N lp = 25 redshift bins. We further consider clusters at redshifts of z cl = 0.2, 0.6, and 1.0. Given the changed reconstruction parameters, we modify our noise control parameters slightly and, for all the results which follow, take λ = 7.5, = 2.8.

Figure 8 shows the two-dimensional projection of the reconstruction and the 1D reconstruction of the four central lines of sight as before. Figure 9 shows the reconstructions of these haloes using our method as a three-dimensional rendering. The three-dimensional rendering is, as before, thresholded at δ rec = 3 and smoothed with a Gaussian of width σ = 0.7 pix.

Several features are immediately apparent. Firstly, all three clusters are clearly identified by our method. This is particularly impressive in the case of the z cl = 1.0 cluster, as linear methods have, thus far, been unable to reconstruct clusters at such a high redshift (see, e.g. Simon et al. 2009;VanderPlas et al. 2011). We note that such a detection is dependent on the redshift distributions of sources, however, and the lack of detection in Simon et al. (2009) and VanderPlas et al. ( 2011) may be due, in part, to their choice of probability distribution. However, we also note that the background source density in our sample is highly diminished behind the z cl = 1.0 cluster (32.5 arcmin -2 ).

Again, the three-dimensional renderings indicate that there is very little smearing of the reconstruction along the line of sight, in contrast with linear methods. This is further evidenced by the line-of-sight plots, in which the unsmoothed reconstructions show very localised structure. Furthermore, the reconstructions exhibit minimal redshift bias, and some lines of sight are seen to recover the amplitude of the density contrast without any notable damping.

However, again we see several prominent "hot pixels" or false detections along noisy lines of sight. Such detections are more evident as the cluster moves to high redshift, and may be significantly larger than the expected density con- trast of the cluster. This is expected, as the number density of sources behind the lens diminishes.

We note that these false detections often manifest at high redshift, arising out of the overfitting effect seen in Figure 4. We also note that the false detections are very well localised in both angular and redshift space, and do not form coherent large structures, making them easily identifiable as false detections; they tend to be localised to isolated pixels. Because of this lack of coherence, we expect that a fully three-dimensional implementation would suppress many of these false detections by aiming to detect coherent structure in three dimensions, and thereby seeking structures of larger extent than a single pixel.

More Complex Line-of-Sight Structure

Given the ability of our method to localise structure in redshift space, it is interesting to consider whether the algo-rithm is able to disentangle the lensing signal from two clusters located along the same line of sight. We consider three different cluster pairings -z cl = [0.2, 0.6], z cl = [0.2, 1.0], and z cl = [0.6, 1.0] -and the results obtained by our method are shown in three-dimensional rendering in Figure 10. We find that the reconstructions of individual lines of sight in this case is significantly more noisy than before.

The figure shows the z cl = [0.2, 0.6] cluster pairing to be reconstructed as a single, coherent structure smeared out between the two redshifts. In the other two cases, two distinct clusters are observed, but their redshifts are slightly biased, and a moderate amount of smearing in the redshift direction is seen. This smearing arises from the fact that individual lines of sight detect the structures at different redshifts; the aggregate effect is an elongation along the line of sight.

While these results are by no means as clean as the results obtained for single clusters, it is promising to note that we are able to detect the presence of more complex line of sight structure, despite the reconstruction of individual lines of sight being fairly noisy. It seems clear that in this case, a fully three-dimensional treatment of the data which takes into account the correlations between neighbouring lines of sight and which seeks to reconstruct coherent structures would offer improvements in this area.

Summary and Conclusions

Current approaches to 3D weak lensing involve linear inversion, where a pseudo-inverse operator is constructed incorporating prior constraints on the statistical distribution of the measurement noise and the underlying density. These methods are straightforward to construct and implement, make use of common tools, and are usually fairly fast. This makes them a convenient choice when approaching the 3D weak lensing problem.

However, reconstructions obtained in this way suffer from line-of-sight smearing, bias in the detected redshift of structures, and a damping of the reconstruction amplitude relative to the input. It has further been noted (VanderPlas et al. 2011) that the reconstructions obtained using these techniques may be fundamentally limited regarding the resolution attainable along the line of sight, due to smearing effects resulting from these linear methods. In addition, such methods are unable to treat an underdetermined inversion, and therefore are limited in their output resolution by the resolution of the input data which, in turn, is limited by the measurement noise.

We have presented a new approach to 3D weak lensing reconstructions by considering the weak lensing problem to be an instance of compressed sensing, where the underlying structure we aim to reconstruct is sparsely represented in an appropriate dictionary. Under such a framework, we are able to consider underdetermined transformations, thereby relaxing the constraints on the resolution of the reconstruction obtained using our method.

We have reduced the problem to that of one-dimensional reconstructions along the line of sight. Whilst this is clearly not optimal, as it throws away a lot of information, it allows us to simplify the problem and to employ a particularly simple basis through which we impose sparsity. We employ techniques recently developed in the area of convex optimisation to construct a robust reconstruction algorithm, and demonstrate that our method closely repro-duces the position, radial extent and amplitude of simulated structures, with very little bias or smearing. This is a significant improvement over current linear methods. In addition, we have shown that our method produces clean results that demonstrate an improvement over linear methods, even when the signal to noise in our data is reduced by a factor of ∼ 2 in line with that expected from current lensing surveys.

Furthermore, we demonstrate an ability to reconstruct clusters at higher redshifts than has been attainable using linear methods. Although our reconstructions exhibit false detections resulting from the noise, these noisy peaks do not form coherent structures, and are therefore well-localised and easily identifiable as noise peaks.

We have also tested the ability of our method to reconstruct multiple clusters along the line of sight. Whilst the reconstructions are noisy, and exhibit a stronger redshift bias, damping and smearing than that seen in the single cluster case, our method is seen to be sensitive enough to detect the presence of more than one structure along the line of sight. It is hoped that by improving our method, we will be able to more accurately reproduce these structures.

It is clear from our results that this method holds a lot of promise; even in the one dimensional implementation presented here, our results clearly show improvements to the bias, normalisation and smearing problems seen with linear methods. There is much room for improvement, however, and a fully three-dimensional implementation of the algorithm described here (which, as written, is entirely general) will be presented in future work. Such a treatment is expected to reduce the incidences of false detection in our reconstructions, the key to this improvement being the choice of an appropriate three-dimensional dictionary with which to sparsify the solution.

In addition, we note that the simulations presented here are idealised as compared to real data. Most notably, we do not take any account of photometric redshift errors, which will be present in any real dataset. This is an important source of error in lensing measurements, and in order for the method presented above to have any real value in reconstructing the matter distribution in a lensing survey, it must include treatment of these errors. This will be presented in a future work, included in the three-dimensional implementation of the algorithm presented here.

The quality of data available for weak lensing measurements continues to improve, and the methods by which we measure the weak lensing shear are becoming ever more sophisticated. So, too, must the methods we use to analyse the data in order to reconstruct the underlying density. While linear methods appear to be limited in resolution, and offer biased estimators, it is clear that a nonlinear approach such as ours does not, and -in a fully three-dimensional implementation -may therefore allow us to map the cosmic web in far greater detail than has previously been achieved.

Acknowledgments

We would like to thank Jalal Fadili for useful discussions, and Jake VanderPlas for supplying the code for the linear reconstructions and for helpful comments on an earlier version of this paper. We would also like to thank the anonymous referee, for their useful comments and suggestions. This work is supported by the European Research Council grant SparseAstro (ERC-228261).

which are the accumulation variables at iteration N -1. Chambolle & Pock (2011) have shown that the sequence defined by G pd (x N , ξ N ) N ∈N is bounded, and decreases at a rate of O(1/t), where t is the iteration number.

In order to use (B.2) in our context, the two indicatrice functions inside (A.4) are not considered, as they play little role. Therefore, in our case, we may rewrite G pd (x, y) as: We determine the algorithm to have converged when

∆G N pd ≡ G pd (x N , ξ N ) -G pd (x N -1 , ξ N -1 ) G pd (x N -1 , ξ N -1 ) < α, (B.7)
where the appropriate α can be determined from simulations, and is dependent on a tradeoff between the desired level of accuracy of the reconstructed data and the time taken to complete the reconstruction. Choosing α to be large will result in an estimate of the solution that may be some distance away from the solution that would be obtained if absolute convergence were reached, but which is obtained in a small number of iterations. pd as a function of iteration number N for the simulations described in § 5. This function is clearly a smooth, largely steadily decreasing function of iteration number, and thus an appropriate choice for defining convergence. Note that the curve shows some oscillations with iteration number. Such oscillations arise on lines of sight where the noise results in the algorithm having difficulty fitting the data within the constraints. With appropriate choice of parameters, these oscillations are relatively small and the curve eventually becomes smooth.

We find from experimentation that α = 10 -6 yields a solution that is sufficiently accurate for our purposes, and is sufficient to largely remove the oscillations, thus we set this to be our threshold for all reconstructions presented above. In addition, we require that each line of sight undergoes at least 1500 iterations, to avoid misidentification of convergence due to early oscillations of ∆G N pd , specifically a strong dip seen in this curve at the start of iteration.

It is possible, along certain lines of sight, for the estimate of the solution to remain constant at zero due to the soft thresholding while ∆G N pd still varies. To account for this, if the current cumulative estimate x N does not vary for 200 iterations, we assume that convergence has been reached for that line of sight.

We now consider several practical issues involved in the implementation of this algorithm. where Θ ≡ Σ -1 2 Q 2 + Φ * 2 is the sum of the 2 norms of the operators used in the algorithm. This sum is dominated by the second term, as the elements of the Q lensing efficiency matrix are small. We have chosen a δ-function dictionary, therefore the application of the transformation Φ * represents a multiplication by the identity matrix.

Thus we have Θ ∼ Φ * 2 ∼ 1, which implies that ωτ ∼ 1 is appropriate. For all the results that follow, we choose ω = τ = 1. Smaller values of these parameters may be used with little effect on the resulting solution. However, the smaller the values chosen, the smaller the steps taken in each iteration of the algorithm. This results in a slower convergence than seen with larger values of ω and τ . Choosing larger values than implied by Equation (B.8) is not advised, as convergence of the algorithm is not guaranteed in this case; the algorithm may give rise to a strongly oscillating estimator, and the resulting solution may not be the one that best fits the data given the prior constraints.

Fig. 1 .

 1 Fig. 1. Top Panel : Lensing efficiency kernel for given source planes as a function of lens redshift. Bottom Panel : Lensing efficiency of a given lens plane as a function of source redshift.

  each simulated image, one or more clusters are generated following an NFW density profile with M 200 = 10 15 M , c = 3 binned into N sp redshift bins. The effective convergence and shear are computed by integrating the

Fig. 2 .

 2 Fig. 2. Simplified schematic of the reconstruction algorithm described in the text.

Fig. 3 .

 3 Fig. 3. Redshift probability distribution p(z) of the sources as a function of source redshift z for the simulations described in the text.

Fig. 4 .

 4 Fig. 4. Reconstructions for two lines of sight obtained with varying . The top row shows the input data and reconstructed convergence vector, while the bottom row shows the reconstructed density contrast.

Fig. 5 .

 5 Fig. 5. Comparison of our method with the linear methods as labelled. The left column shows the 2D projection of the reconstruction, while the right column shows the 1D reconstructions along the four central lines of sight (dashed lines). Note that, due to the amplitude damping effect in SVD and Wiener reconstructions, the y-axis scaling is different in each of the line of sight plots.

Fig. 7 .

 7 Fig. 7. Reconstructions of a cluster at reshift z = 0.25 using noisy data with n g = 30 galaxies per square arcminute. Shown above are the 2D projection of the reconstruction (top left), 1D line of sight plots for the four central lines of sight (top right) and smoothed 3D rendering of the reconstruction, as before.

Fig. 8 .

 8 Fig. 8. Reconstructions of single clusters located at a redshift of z cl = 0.2 (top row), z cl = 0.6 (middle row) and z cl = 1.0 (bottom row). As before, the left column shows the two dimensional integrated projection of the reconstruction, while the right panel shows the input density contrast along the line of sight (solid line) and the 1D reconstruction along each of the four central lines of sight (dashed lines).

G

  pd (x, y) ≈ max y y , Ux -G * (y ) -G * (y) , (B.4) ≈ G(Ux) -G * (y) , (B.5) ≈ λ Φ * x 1 .(B.6)

Fig. B. 1 .

 1 Fig. B.1. The function ∆G N pd , as defined in the text, as a function of iteration number.

Figure B. 1

 1 Figure B.1 shows a characteristic example of the function ∆G N pd as a function of iteration number N for the simulations described in § 5. This function is clearly a smooth, largely steadily decreasing function of iteration number, and thus an appropriate choice for defining convergence. Note that the curve shows some oscillations with iteration number. Such oscillations arise on lines of sight where the noise results in the algorithm having difficulty fitting the data within the constraints. With appropriate choice of parameters, these oscillations are relatively small and the curve eventually becomes smooth.We find from experimentation that α = 10 -6 yields a solution that is sufficiently accurate for our purposes, and is sufficient to largely remove the oscillations, thus we set this to be our threshold for all reconstructions presented above. In addition, we require that each line of sight undergoes at least 1500 iterations, to avoid misidentification of

B. 2 .

 2 Choice of Step Sizes ω, τ ω and τ control the step size in the evolution of the algorithm, and are required to be positive and to satisfy the inequality ωτ Θ 2 ≤ 1 (B.8)

Fig. 9 .

 9 Fig. 9. Reconstructions of single clusters located at a redshift of z cl = 0.2 (top), z cl = 0.6 (middle) and z cl = 1.0 (bottom). The reconstructions are thresholded at δ = 3 and smoothed with a Gaussian of width σ = 0.7 pix.

Fig. 10 .

 10 Fig. 10. Reconstructions of two clusters along the line of sight, located at a redshift of z cl = [0.2, 0.6] (top), z cl = [0.2, 1.0] (middle) and z cl = [0.6, 1.0] (bottom). The reconstructions are thresholded at δ = 3 and smoothed with a Gaussian of width σ = 0.7 pix.

Appendix A: Solving the optimization problem

In this section, we explain the development of the reconstruction algorithm. First, we will present an overview of some key concepts and results from convex optimisation, before introducing the primal-dual scheme chosen to solve Equation ( 17), and finally discussing the convergence criterion for the algorithm. For a complete introduction to convex analysis, the reader is referred to [START_REF] Rockafellar | Convex analysis[END_REF]; [START_REF] Lemaréchal | Convex Analysis and Minimization Algorithms I[END_REF]; [START_REF] Boyd | Convex Optimization[END_REF].

Notation and terminology

Let H a real Hilbert space; in our case, a finite dimensional vector subspace of R n . We denote by . the norm associated with the inner product in H, I is the identity operator on H, and . p (p ≥ 1) is the p norm. A real-valued function f is coercive if lim x →+∞ f (x) = +∞, and is proper if its domain is non-empty:

Lastly, Γ 0 (H) represents the class of all proper lower semicontinuous (lsc) convex functions from H to (-∞, +∞].

A.1. Convex optimisation

The theory of convex optimisation aims to solve problems of the form

where F : H → R is a convex function. If F represents a potential function, for example, then (A.1) will seek for the state of equilibrium, i.e. the state that minimises the potential.

Many algorithms exist for solving such problems, and an excellent introduction can be found in [START_REF] Boyd | Convex Optimization[END_REF]. Recently, methods have been developed which exploit the decomposition of f into a sum of smaller convex functions. This can be very helpful if these functions show properties that are not preserved when the functions are summed.

For example, if we assume that F can be decomposed into a sum of K functions all in Γ 0 (H), then solving (A.1) is equivalent to solving

Such formulation is interesting when the F k functions show properties that are lost while considering the sum F , directly. In our case, we seek a formulation such that the functions we seek to minimise have proximal operators that are simple, and have closed form. While the sum F may not have this property, the functions F k can be chosen such that they satisfy this condition.

A.2. Proximity operator

We may now begin to describe the proximal splitting algorithm used in this work. At the heart of the splitting framework is the notion of the proximity operator :

Definition 1 (Moreau (1962)). Let F ∈ Γ 0 (H). Then, for every x ∈ H, the function y → F (y) + xy 2 2 /2 achieves its infimum at a unique point denoted by prox F x. The operator prox F : H → H thus defined is the proximity operator of F . Therefore, the proximity operator of the indicator function of a convex set C:

is merely its orthogonal projector. Some key properties are presented in the following lemma:

Lemma 2 (Combettes & Wajs (2005)).

-

A.3. Primal-dual scheme

Consider the optimization problem:

where G and B are two convex, proper and lower semicontinuous functions and U is a linear bounded operator. Equation (A.3) can be see as a special case of Equation (A.2) where K = 2 and one of the functions contains a linear bounded operator. Then the Algorithm 1, proposed in Chambolle & Pock (2011), will converge to the solution of (A.3), assuming that the proximity operators of G and B are easy to compute or known in closed form.

Algorithm 1: Primal-dual scheme for solving .

Parameters: The number of iterations Niter, proximal steps σ > 0 and τ > 0.

Initialization:

-Primal step:

-Update the coefficients estimate: xt+1 = 2x t+1 -x t

End main iteration Output: Final solution x = x N iter .

Adapting the arguments of (Chambolle & Pock 2011), convergence of the sequence (x t ) t∈N generated by Algorithm 1 is ensured.

Proposition 3. Suppose that G and B are two convex, proper and lower semi-continuous functions. Let ζ = U 2 , choose τ > 0 and σ such that στ ζ < 1, and let (x t ) t∈R be that defined by Algorithm 1. Then, (x t ) t∈N converges to a (non-strict) global minimiser of Equation (A.3) .

A.4. Application to solution of Equation (17)

First, we must rewrite the problem in Equation ( 17) in an unconstrained form by replacing the constraints by the indicator functions of the corresponding constraint sets:

where ı 2(ε) is the indicatrice function of the 2 -ball of radius ε and ı C the indicatrice function of a closed convex set C. Notice that (A.4) can be expressed in the form of (A.3) with,

We now need only to apply Algorithm 1 in order to compute the solution. This requires computation of the three proximity operators, which are given by the following proposition:

, then its associated proximity operator, prox λF , is the component-wise soft-thresholding operator with threshold λ as defined by Equation (B.1);

then its associated proximity operator, prox F , is the Euclidian projector to the convex set C, Prj C ;

then its associated proximity operator, prox F , is the projector onto the 2 -ball with radius ε defined as:

x else .

The application of Proposition 4 and Lemma 2 to the optimization problem in Equation (A.4), using the the primal-dual scheme in Algorithm 1, yields the reconstruction algorithm expressed in Algorithm 2.

Appendix B: Practical Considerations

Application of the method described in Appendix A to the specific case of 3D lensing leads us to Algorithm 2, where one can note the projection over the two constraints described above (data fidelity and minimum value of the solution), and the operator St λ , which imposes a soft threshold at a level of λ/ω as:

Algorithm 2: Nonlinear iterative algorithm for solving Equation ( 17).

Parameters: Choose ω, τ > 0 such that ωτ Θ 2 ≤ 1 where

, where R 2 is the spectral norm of the operator R. Initialization:

Main iteration: For n = 0 to Niter -2, -Initialise auxiliary variables:

-Sparsity-promoting penalty:

-Data fidelity term:

-Projection on the convex set C:

-Update current estimate:

End main iteration

B.1. Convergence criteria

The main difficulty with the primal-dual scheme described above -indeed, with any iterative algorithm -is to define an appropriate convergence criterion. In this case, the difference between two successive iterates x t and x t+1 is not bounded. However, the partial primal-dual gap G pd , defined by: Let us define two variables from the sequences (x t ) t and (ξ t ) t produced by Algorithm 1:


Nissrine Boujia

Franziska Schmidt
email: franziska.schmidt@univ-eiffel.fr

Christophe Chevalier

Dominique Siegert

Damien Pham Van Bang
email: damien.pham_van_bang@ete.inrs.ca

Damien Pham Van Bang

Using rocking frequencies of bridge piers for scour monitoring

Keywords: Scour, soil-structure interaction, pier-deck interaction, modal shapes, rigid modes, monitoring Pier diameter [m]

come   L'archive ouverte pluridisciplinaire

1 Background

Scour is the erosive action of water carrying away sediment of the riverbed. Around bridge supports, scour holes are created, undermining the foundations and decreasing the bearing capacity of the soil. As a consequence, scour is the main cause of bridge failures in many countries [START_REF] Anderson | Ground-penetrating radar: A tool for monitoring bridge scour[END_REF][START_REF] Azhari | Laboratory validation of buried piezoelectric scour sensing rods[END_REF]. Thus, in order to address this risk, various monitoring techniques have been developed and tested on field. In early days, geophysical tools such as radar and sonar [START_REF] Bao | Vibration-based bridge scour detection: A review[END_REF][START_REF] Ting Bao | Critical insights for advanced bridge scour detection using the natural frequency[END_REF][START_REF] Boujia | Distributed Optical Fiber-Based Approach for Soil-Structure Interaction[END_REF][START_REF] Boujia | Effect of Scour on the Natural Frequency Responses of Bridge Piers: Development of a Scour Depth Sensor[END_REF] were most commonly used. Various monitoring techniques use magnetic field properties and includes sensors such as the magnetic sliding collar [START_REF] Chen | Maximum scour depth based on magnetic field change in smart rocks for foundation stability evaluation of bridges[END_REF] and smart rocks [START_REF] Deng | Bridge scour: Prediction, modeling, monitoring, and countermeasuresreview[END_REF]. The main limitations of these monitoring techniques are their sensitivity to suspended sediments and their vulnerability in harsh environment which is major hurdle to continuous scour monitoring.

Recently, the use of dynamic behaviour-based techniques to measure scour depth has been confirmed as a promising alternative to the previous techniques. Two approaches can be identified, namely an indirect and direct approach. The indirect approach consists in placing rodsensors around the bridge piles and monitor their frequency [START_REF] Florens | Scour monitoring on bridge piermethodology and implementations[END_REF][START_REF] Foti | Influence of Foundation Scour on the Dynamic Response of an Existing Bridge[END_REF]. As scour increases, the first frequency of the sensor decreases [START_REF] Gorin | Use of surface-geophysical methods to assess riverbed scour at bridge piers[END_REF].

The direct approach consists in monitoring the frequency of the bridge structure itself. [START_REF] Lu | Field measurements and simulation of bridge scour depth variations during floods[END_REF] performed impact tests on railway bridges to assess their foundation conditions. The results indicate that the more the piers are scoured, the lower their frequencies, with definition of a limit numerical value of the frequency making it possible to discriminate damaged piers from those in good health. [START_REF] Millard | Assessing bridge pier scour by radar[END_REF] studied the effect of scour on the modal parameters of a five spans bridge in Italy. This study showed that the frequencies and modal shapes of the spans are highly affected by scouras also reported by various other authors [START_REF] Prendergast | Determining the Presence of Scour around Bridge Foundations Using Vehicle-Induced Vibrations[END_REF][START_REF] Prendergast | An investigation of the changes in the natural frequency of a pile affected by scour[END_REF]. [START_REF] Shinoda | Nondestructive evaluation of railway bridge substructures by percussion test[END_REF] studied the effect of scour on a single pile (slenderness 𝛼 > 10) and developed a 2D-model based on the Winkler theory [START_REF] Taricska | An Analysis of Recent Bridge Failures[END_REF] for Soil-Structure Interaction to estimate the variation of the pile frequency with scour depth. Thus, the soil-pile interaction was modeled by a beam exposed to horizontal independent linear springs which were progressively removed to take into account scour. The proposed model was able to give a good prediction of the evolution of the single pile frequency with scour. However, it is difficult to extend this model to the case of a bridge pier. In fact, on the one hand bridge piers are rarely slender as this given pile [START_REF] Shinoda | Nondestructive evaluation of railway bridge substructures by percussion test[END_REF]. One the other hand, the pier stiffness with respect to soil stiffness is usually very high. A recent experimental study of [START_REF] Wardhana | Analysis of recent bridge failures in the United States[END_REF] studied the effect of scour on the vibration of a laboratory-scale concrete column (slenderness 𝛼 = 2). The results show that, indeed, the first frequency of the column decreases with scour and can be used for scour monitoring. However, open points in this study can be identified : First, the modal shape of the pier corresponding to the identified frequency was not determined. Then, even though the characteristics of the tested single pier were closer to those of a real bridge pier, the pier-deck interaction was not taken into account. Finally, no physical model was proposed in order to estimate the frequency of this kind of structures for various scour depths.

The present paper addresses some open issues in this topic: First, we determine the physical meaning of the measured frequency : When monitoring the vibration frequency of slender and flexible structures, such as piles, it is easy to visualize the modal shape corresponding to their bending mode. For rigid structures, such as piers, it is more complex to understand exactly the kinematics of the modal shapes. Then, we will study the effect of the presence of the deck on the pier frequency and its sensitivity to scour. Finally, we will try to establish whether or not it is adequate to model the soil-pier by horizontal springs only, considering the aspect ratio between the lateral faces and the base of a cylindrical pier.

The paper starts with a description of the experimental campaigns that have been undertaken in a water flume. The experimental results are then presented in Section 2. In Section 3, a theoretical formulation for the soil-pier interaction is proposed and a sensitivity study of the dynamic response of the piers to the slenderness ratio and soil stiffness is performed. In Section 4, a comparison between the numerical and experimental results is presented. Finally, the main achievements of this work are summarized.

Experimental campaigns

The dynamic behaviour of reduced-scale bridge piers is analyzed in a water flume, with and without water, before studying a reduced-scale, one-span bridge.

Single bridge pier tests

A laboratory scale model bridge pier is made in order to assess the effect of scour on its frequency 𝑓. The concrete pier is cylindrical and has a total length 𝐿 = 55 cm and a diameter 𝑏 = 11 cm, which are the reduced dimensions of the piers of a prestressed concrete bridge near Orléans [START_REF] Winkler | Theory of elasticity and strength[END_REF].

Dry sand tests on a single pile

The experimental set-up is presented in Figure 1 : A tank of dimensions 1𝑚 × 1𝑚 × 1𝑚 is progressively filled with dry sand of Seine until it reaches a height of 0.7 m. The pier is then placed in the middle of the sand volume to avoid edge effects, with an embedded length 𝐷. A tri-axial, high-sensitivity accelerometer is placed at the top of the pier and records its response with a sampling frequency of 25.6kHz. It is connected to a data-acquisition card from National Instruments, controlled with LabVIEW. A symmetric scour hole is then created by progressively extracting a 5 cm-thick soil layer, so that the embedded length of the pier varies from 25 to 5 cm. For each embedded length 𝐷, an impact is applied with a calibrated impact hammer and the response of the pier in the horizontal direction 𝑥 is recorded by the accelerometer.

Single pier, water flume tests

The next series of experiment are conducted in a water flume : The flume is progressively filled with Hostun sand (𝐷 50 = 0.6 mm) until it reaches a final height of 25 cm. The pier is then placed with an initial embedded length 𝐷 = 20 cm (Figure 2). The flow velocity is progressively increased to generate a scour hole. The acceleration of the pier in the flow direction 𝑥 is recorded with two accelerometers (Figure 2), which will be used to identify the modal shapes corresponding to the measured frequencies of the pier. Two different testing protocols are applied : The first protocol consists in applying a hammer impact, as in the laboratory tests, at regular intervals of scour depth. For the second experimental protocol, the flow induced vibrations of the pier are recorded all along the scouring process.

Reduced-scale bridge model tests

For the final series of tests, a bridge configuration (Figure 3) is tested to assess the effect of the pier-deck interaction on the first frequency of the pier. Two bridge piers distant of 140 cm are place in the flume with an initial embedded length 𝐷 = 20 cm. Rubber bearings are then placed on the top of each pier and the model bridge deck is placed over the two rubber bearings. The length of the deck is 160 cm and its total weight is 𝑀 = 26 Kg. A hammer impact is applied at the middle of the deck at regular intervals of scour depth and the piers accelerations, in the flow direction 𝑥, are recorded with two accelerometers.

Experimental results

Effect of scour on the first frequency of the pier

The accelerations of the pier, measured by the accelerometers in all tested configurations presented in Section 2, contain a large number of high frequencies. Therefore, the signal is first filtered using a 4 th order low-pass Butterworth filter with a cut-off frequency of 100 Hz, which makes it possible to capture the structural response of the pier . Then, the first frequency 𝑓 of the pier is calculated using a Fast Fourier Transform (FFT). The variation of the first frequency 𝑓 of the pier with the maximum scour depth in dry sand, in flume without the deck and in the bridge configuration are as expected (Figure 4) : as scour depth increases, meaning as the embedded length 𝐷 decreases, the first frequency of the pier decreases. However, the results also show that the frequency of the pier in the flume, with and without the deck, remains steady for the first scoured centimeters. One possible reason for this insensitivity to scour is that the first layer of the soil is very loose : when conducting the test, the flow is first kept steady until the water reaches the required level of ℎ 0 = 7.5cm. During this initial phase and since the water level is very low at the beginning, the pier is scoured, up to 5cm in our configuration. Once the required water level is reached, sand is moved back to fill the hole and the impact tests are performed. However, since no compaction of the added soil is possible, this layer of the soil remains loose and may not insure much constraint to the pier.

Regarding the effect of the deck on the pier's frequency, the experimental results (Figure 4) show a significant decrease in the frequencies : the frequencies of the pier in the bridge configuration are less by about 35% compared to the frequencies of the single pier.

The sensitivity of the pier's frequency to scour in the two tested configurations is compared : to achieve that, the variation ratio of the frequency |Δ𝑓|/𝑓 0 (Δ𝑓 = 𝑓 -𝑓 0 where 𝑓 is the measured frequency for a given scour depth and 𝑓 0 is the frequency of the pier without scour) is calculated for each scour depth. The results (Table 1) indicate that the sensitivity of the pier to scour does not change significantly when taking into account the pile-deck interaction.

Identification of the modal shape of the pier

The tested pier has a slenderness ratio 𝛼 = 5 and its stiffness with respect to the soil stiffness is 𝐸 𝑝 /𝐸 𝑠 > 10 3 , where 𝐸 𝑝 is the Young modulus of the pier and 𝐸 𝑠 the Young modulus of the soil. Under impact tests, the pier vibrates following a rigid-body motion given its slenderness and stiffness. Thus, placing two accelerometers on the pier is sufficient to determine the mode shapes corresponding to the measured first frequency.

First, the first frequency of the pier is calculated using a Fast Fourier Transform (FFT). The signal is then filtered using a low pass filter (Figure 5 (a), (b), (c)) to make sure that the measured accelerations correspond to the first frequency only and therefore, the first modal shape. The modal shapes are normalized to the maximum displacement. The normalized modal shapes for different embedded lengths 𝐷 are plotted in Figure 5 (d). The results show that the governing motion of the pier is a coupled swaying-rocking mode with a center of rotation embedded in the soil. However, the variation of the modal shape with scour is insignificant.

Proposed theoretical model 4.1 Analytical formulation

In the following, we consider a bridge pier embedded at a distance 𝐷 in the soil and under free vibration in the lateral direction 𝑥 (Figure 6). So the model takes into account only the rigid body motion of the pier.

The following assumptions are made in the theoretical formulation:  The pier is a rigid beam with a total mass 𝑚, a moment of inertia around the gravity center 𝐽 𝐺 and 𝐿 1 the distance between the base and the gravity center.  The soil is assumed to be homogeneous and its weight is neglected.  The lateral soil reaction is modelled as horizontal independent linear elastic springs according to the Winkler theory [START_REF] Prendergast | Determining the Presence of Scour around Bridge Foundations Using Vehicle-Induced Vibrations[END_REF] with a stiffness 𝑘 𝑠 .  The base rotational springs 𝐾 𝑟 take into account the moment produced by normal pressure on the base of the pier, induced by base rotation. In the plane, the pier has two degrees of freedom: the displacement 𝑥(𝑡) and the rotation 𝜃(𝑡) of its gravity center.

The equilibrium of horizontal forces can be written:

𝑚𝑥̈+ ∫ -𝐿 1 +𝐷 -𝐿 1 𝑘 𝑠 (𝑥 + 𝜃𝑦) d𝑦 = 0. (1)
The rotation equilibrium around the gravity center can be expressed as:

𝐽 𝐺 𝜃 ̈+ 𝐾 𝑟 𝜃 + ∫ -𝐿 1 +𝐷 -𝐿 1 𝑘 𝑠 (𝑥 + 𝜃𝑦)𝑦 d𝑦 = 0. (2)
Equations ( 1) and (2) lead to the following system:

{ 𝑚𝑥̈+ 𝑘 𝑠 𝐷𝑥 + 𝑘 𝑠 𝜃 2 𝐷(𝐷 -2𝐿 1 ) = 0, 𝐽 𝐺 𝜃 ̈+ 𝐾 𝑟 𝜃 + 𝑘 𝑠 𝑥 2 𝐷(𝐷 -2𝐿 1 ) + 𝑘 𝑠 𝜃 3 [(𝐷 -𝐿 1 ) 3 + 𝐿 1 3 ] = 0, (3)
which can be written as:

[ 𝑚 0 0 𝐽 𝐺 ] { 𝑥θ ̈} + [ 𝑘 𝑠 𝐷 𝑘 𝑠 2 𝐷(𝐷 -2𝐿 1 ) 𝑘 𝑠 2 𝐷(𝐷 -2𝐿 1 ) 𝐾 𝑟 + 𝑘 𝑠 3 [(𝐷 -𝐿 1 ) 3 + 𝐿 1 3 ] ] { 𝑥 𝜃 } = 00 (4)
The solution of Equation ( 4) can be expressed as:

{ 𝑥(𝑡) = 𝑋cos(𝜔𝑡 + 𝜑) 𝜃(𝑡) = Θcos(𝜔𝑡 + 𝜑) (5)
The following system is then obtained by substituting equation (5) in equation ( 4):

[ 𝑘 𝑠 𝐷 -𝑚𝜔
The eigenfrequencies 𝑓 1 and 𝑓 2 of the two rigid modes of the pier are the solutions of Equation [START_REF] Chen | Maximum scour depth based on magnetic field change in smart rocks for foundation stability evaluation of bridges[END_REF].

The modal shapes (𝑋 𝑖 , Θ 𝑖 ) 𝑖⊂{1,2} corresponding to the two eigenfrequencies of the beam are calculated using the following equation:

(𝑘 𝑠 𝐷 -𝑚𝜔 𝑖 2 )𝑋 𝑖 + 𝑘 𝑠 2 𝐷(𝐷 -2𝐿 1 )Θ 𝑖 = 0 (8) Hereafter, we choose to normalize the modal shapes to the maximal displacement.

Theoretical relation between 𝒌 𝒔 and 𝑲 𝒓

In order to perform a sensitivity study, a first attempt is made to establish a simple relation between the horizontal stiffness of the soil 𝑘 𝑠 and the rotational stiffness 𝐾 𝑟 , an energetic approach is used. We make the assumption that the base of the pier is resting on uniformly distributed closely spaced springs having the same stiffness 𝑘 𝑠 as the horizontal ones as shown in figure 7. The elastic energy of deformation at the base of the pier can be written:

1 2 ∫ + 𝑏 2 - 𝑏 2 𝑘 𝑠 (𝜃𝑦) 2 d𝑥 = 1 2 𝐾 𝑟 𝜃 2 . ( 9
)
Therefore, the initial estimated rotational stiffness (𝐾 𝑟 ) 𝑖𝑛𝑖𝑡 is proportional to 𝑘 𝑠 according to the following equation:

(𝐾 𝑟 ) 𝑖𝑛𝑖𝑡 = 1 12
𝑘 𝑠 𝑏 3 .

(10)

Numerical validation of the theoretical model

In order to validate the theoretical formulation of the model and perform a sensitivity analysis to the soil stiffness and the slenderness ratio of the pier, a 2D-finite element model is developed. The pier is modeled as an Euler Bernoulli beam and the soil-structure interaction is modelled with linear elastic springs of a stiffness 𝐾 that verifies: 𝐾 = 𝑘 𝑠 × Δ𝑦, (11) where 𝑘 𝑠 is the modulus of subgrade reaction of the soil and Δ𝑦 is the spacing between the springs. A rotational spring is modeled at the base of the beam. The finite element model is then computed for different soil stiffness. The results of the parametric study are summarized in table 4. This numerical model show good agreement with analytical model, for various value of the soil stiffness 𝑘 𝑠 (Table 2 The results also indicate that as soil stiffness increases, the gap between analytical and numerical results increases : This is because, as soil stiffness increases, the assumption that the pier behaves as a rigid body with respect to the soil stiffness is no longer valid. In this case, the pier behaves as an elastic beam.

A second sensitivity study was performed to evaluate the effect of the slenderness ratio 𝛼 on the validity of the proposed model (Table 3 However, for a slenderness ratio 𝛼 = 50, the vibration of the pier is governed only by elastic modes. Therefore, the analytical model is no longer valid for this kind of structures.

The effect of the rotational spring

A sensitivity study is performed to assess the contribution of the rotational spring to the soil stiffness. The horizontal stiffness is 𝑘 𝑠 = 50 MN/m 2 and the stiffness of the rotational spring (𝐾 𝑟 ) 𝑖𝑛𝑖𝑡 is calculated using Equation [START_REF] Foti | Influence of Foundation Scour on the Dynamic Response of an Existing Bridge[END_REF]. The frequencies of the pier are computed, with and without the rotational spring, for an embedded length 𝐷 from 5 to 25 cm. Only the first rigid frequency 𝑓 1 is significantly modified by the rotational spring (Figure 8). This effect is more predominant as the embedded length 𝐷 decreases. This is mainly linked to the fact that the less the pier is embedded, the less the lateral soil is mobilized in comparison with the soil base. When the pier is embedded enough, the vibration of the pier in governed mainly by the lateral soil interaction and the effect of soil base is therefore negligible.

Proposed model validation 5.1 Identification of the soil stiffness parameters

A 3D-finite element model is developed using Code-Aster software. Both the soil and the model bridge pier are modeled using 3D elements and are meshed using 10 nodes tetrahedron elements. The model is based on the following hypothesis: 1) the soil medium and the pier are elastic, 2) the soil's Young modulus 𝐸 𝑠 remains constant with depth and 3) the soil and the pier are perfectly bounded at the interface. The soil reaches 2 m depth and 2.5 m width, therefore preventing any border effect on the pier response. The displacement of the lateral faces of the soil are blocked in the normal direction and the base is blocked in all directions. Scour is simulated by creating a cylindrical hole around the pier as shown in Figure 9 The material and geometrical parameters of the pier used are the experimental ones. The mass density and the Poisson ratio of the soil are 𝜌 𝑠 = 1700 Kg/m 3 and 𝜈 𝑠 = 0.33 respectively. Since the Young modulus 𝐸 𝑠 of the soil is the only unknown parameter, its value is calibrated in order to have the best fitting numerical results for the experimental single pier results. A good agreement is obtained for 𝐸 𝑠 = 4.2 MPa. A modal analysis is performed for each 5 cm scour depth to evaluate the frequency of the pier (Figure 8 Once the 3D model is calibrated, the numerical values of 𝑘 𝑠 and 𝐾 𝑟 can be determined following two steps : First, for 𝑘 𝑠 , a second 3D-model, with a gap between the pier's base and the soil, is developed (Figure 10). This makes it possible to eliminate the contribution of the base stiffness to the pier response, therefore this model is equivalent to the analytical model with 𝐾 𝑟 = 0.

This numerical model is computed for high embedded lengths 𝐷 = 25cm and 𝐷 = 20cm to avoid any base effect.

Considering the analytical model, the numerical frequency (𝑓 1 ) 𝑛𝑢𝑚 is solution of Equation ( 7), and thus verifies:

where 𝜔 𝑛𝑢𝑚 = 2𝜋(𝑓 1 ) 𝑛𝑢𝑚 .

The value of 𝑘 𝑠 can then be derived by solving Equation [START_REF] Lu | Field measurements and simulation of bridge scour depth variations during floods[END_REF], see Table 3. This analysis determines a range of possible numerical values for 𝑘 𝑠 , and allowing use to choose 𝑘 𝑠 = 25~MPa in the following.

𝐷 (cm)

𝑓 In a second step, the value of 𝐾 𝑟 is determined. Since the effect of 𝐾 𝑟 is predominant for low embedded length, Equation ( 7) is solved for small values of D (𝐷 = 5cm and 𝐷 = 10 cm). The input parameters are 𝑘 𝑠 = 25 MPa (as determined previously) and 𝜔 = 2𝜋𝑓 𝑛𝑢𝑚 (𝐷 = 5𝑐𝑚; 𝐷 = 10𝑐𝑚) where 𝑓 𝑛𝑢𝑚 the numerical frequency computed with the first 3D-model (Figure 8). The average value of the rotational spring obtained is 𝐾 𝑟 ≃ 4100 N.m/rad. It is worthy to mention that this value correspond to 3 times the value of the initial stiffness (𝐾 𝑟 ) 𝑖𝑛𝑖𝑡 , as obtained with Equation [START_REF] Foti | Influence of Foundation Scour on the Dynamic Response of an Existing Bridge[END_REF]. The values of 𝑘 𝑠 and 𝐾 𝑟 are kept fixed in the following.

Comparison with the experimental results of the single pier

In order to assess the accuracy of the proposed model, the experimental data is compared to the theoretical results. The theoretical frequencies are computed for the previously calibrated value of 𝑘 𝑠 (𝑘 𝑠 = 25 MPa) and two values of the rotational spring: (𝐾 𝑟 ) 𝑖𝑛𝑖𝑡 computed using Equation [START_REF] Foti | Influence of Foundation Scour on the Dynamic Response of an Existing Bridge[END_REF] and the calibrated value of (𝐾 𝑟 ) 𝑐𝑎𝑙𝑖𝑏 so that the relation between the two rotational spring stiffness was find to verify:

(𝐾 𝑟 ) 𝑐𝑎𝑙𝑖𝑏 = 3 × (𝐾 𝑟 ) 𝑖𝑛𝑖𝑡 = 1 8 𝑘 𝑠 𝑏 3 . (13)
Figure 11(a) indicates that the model computed with the initial value of the rotational spring (𝐾 𝑟 ) 𝑖𝑛𝑖𝑡 does not simulate correctly the single pier behaviour. In fact, as the embedded length 𝐷 of the pier decreases the gap between the theoretical and experimental frequencies increases, which is predominately related to an underestimation of the base rotational spring, as demonstrated in Section 4.4. This assumption is confirmed when the model is computed with a higher value of rotational spring stiffness (𝐾 𝑟 ) 𝑐𝑎𝑙𝑖𝑏 . The results show a good agreement between experimental results and the calibrated theoretical model. For the single pier in flume, the model does not take into account the loose stiffness of the first 5 cm of the soil and therefore overestimates these frequencies. The theoretical modal shape of the pier is also calculated. The results (Figure 11(b)) confirm that the first analytical modal shape correspond to a coupled swaying-rocking motion and validates the model calibration.

Comparison with the experimental results of the pier in bridge configuration

For the validation of the modelling of the pier in bridge configuration, we assume that each pier supports half the deck weight. The deck is modeled as a nodal mass placed at the top of the pier and having a total weight 𝑀′ = 𝑀/2 with 𝑀 the total mass of the deck (Figure 12(a)). In this case, the values of the weight 𝑚 and the moment of inertia around the gravity center of the pier 𝐽 𝐺 are modified to take into account the bridge deck. In that respect, the modified values are expressed by:

{ 𝑚′ = 𝑚 + 𝑀′, 𝐽 𝐺 ′ = 𝐽 𝐺 + 𝑀 ′ × 𝑑, (14)
where 𝑚′ and 𝐽 𝐺′ are the modified weight and moment of inertia around the gravity center respectively and 𝑑 the distance between the gravity center of the pier and the gravity center of half the deck. A good agreement (Figure 12(b)) between the theoretical and experimental results is achieved. Therefore, the model is capable of predicting the pier frequency variation with scour when the bridge deck is correctly taken into account.

Confrontation to literature results

In this section, the analytical model is used to model the tests conducted by [START_REF] Bao | Vibration-based bridge scour detection: A review[END_REF]. The tested single bridge pier had a diameter of 𝑏 = 15.3 cm and a total length 𝐿 = 30.6 cm. The embedded length was varied from 𝐷 = 22cm to 𝐷 = 7cm with 3 cm scour increment. The pier was tested in a low density sand.

The proposed theoretical model is computed using a lateral spring stiffness 𝑘 𝑠 = 8.5 MPa and the rotationnel spring 𝐾 𝑟 is calculated using the equation deduced after spring calibration expressed (Equation ( 13)). The results of the proposed model after calibration are in a good agreement with the experimental results (Figure 13).

Final remarks and conclusions

This paper investigates the effects of scour on the frequency of a bridge pier and proposes a simplified model for this kind of structures.

To this end, experimental campaigns were conducted in a water flume using a single pier in two types of soil, dry sand and saturated sand. The pier was also tested in a bridge configuration to evaluate the deck effect on the pier frequency. These experimental results provided validation data for both the first frequency and the corresponding modal shape of the pier.

A theoretical model was proposed to model the pier behaviour. The soil-structure interaction was modelled by both horizontal springs along the pier and a rotational spring at the base to take into account the vertical effect of the soil that become predominant as scour increases. This model has been proven to be in good agreement with the experimental results, for both the frequencies and the modal shape of the single piers and the piers in bridge configuration.

It has therefore been shown that the first rigid frequency of the pier corresponds to a coupled swaying-rocking mode. The bridge configuration has been modelled correctly by applying nodal masses at the top of the piers, equal to half the mass of the deck .

The proposed theoretical model has also been compared to literature results, good agreement has been shown.

The main limitation of the model is the direct identification of the lateral stiffness of the soil instead of fitting the experimental results. Future research will try establishing empirical relations between the lateral spring stiffness and the Young modulus of the soil for rigid piers.

Figure 1 :

 1 Figure 1: Experimental setup of single pier, dray sand tests.

Figure 2 :

 2 Figure 2: Experimental setup of single piers, water flume tests : (a) the initial configuration without scour (h the distance between the two accelerometers) (b) front view and (c) side view of the scoured configuration.

Figure 3 :

 3 Figure 3: Experimental setup of the bridge configuration tests.

Figure 4 :

 4 Figure 4: Variation of the frequency with the embedded length D, in dry sand and flume tests with and without the deck.

Figure 5 :

 5 Figure 5: The filtered accelerations recorded with the accelerometer A and B (labels « ACC A » and « ACC B ») for the embedded lengths (a) D=20cm, (b) D=14cm and (c) D=8cm and (d) the corresponding modal shapes.

Figure 6 :

 6 Figure 6: Rigid body vibration of a pier in soil : (a) Definition sketch (b), Proposed equivalent analytical model.

Figure 7 :

 7 Figure 7: Rotational spring 𝐾 𝑟 equivalence to vertical spring 𝑘 𝑠 .

Figure 8 :

 8 Figure 8: Effect of the rotational spring with varying embedded length (a) on the first frequency 𝑓 1 and (b) the second frequency 𝑓 2 of the pier.

  (a).

  Comparison of experimental and numerical results

Figure 9 :

 9 Figure 9:Finite elements model if the soil-pier system (a) mesh details and (b) results of the single pier in dry sand.

  (b)).

Figure 10 :

 10 Figure 10: Figure of the second numerical model neglecting the base effect.

  (a) Comparison between experimental and theoretical first frequency (b) Theoretical first modal shape

Figure 11 :

 11 Figure 11: Single bridge pier (a) theoretical frequency and (b) corresponding modal shape.

Figure 12 :

 12 Figure 12: Pier in bridge configuration (a) Simplified model and (b) comparison between experimental and theoretical first frequency.

Figure 13 :

 13 Figure 13: Comparison between theoretical and experimental frequencies of [3].

Table 1 :

 1 Frequency variation with scour depth of the single pier and the pier in a bridge configuration.

	Scour depth		|Δ𝑓|/𝑓 0 (%)
	(cm)	Single pier	Pier with the deck
	0	0	0
	3	0	4
	6	13	21
	9	31	35
	12	40	44

Table 2 :

 2 Comparison of analytical and numerical frequencies for different values of soil stiffness.

	).
	Lateral stiffness 𝑘 𝑠 MN.m -2	First frequency (Hz) Analytical Numerical	Second frequency (Hz) Analytical Numerical
	1	7.21	7.19	29.95	30.00
	5	16.12	16.11	66.99	67.01
	10	22.80	22.78	94.96	94.87
	100	72.13	71.01	300.29	299.38
	500	161.29	144.56	671.47	662.02

  ) : the numerical and analytical models are in agreement for low slenderness ratios 𝛼.

	Slenderness ratio 𝛼	First frequency (Hz)	Second frequency (Hz)
	-	Analytical	Numerical	Analytical	Numerical
	2	38.57	38.88	101.16	100.3
	5	26.83	26.80	96.28	96.17
	10	25.85	25.22	96.01	95.71
	50	25.71	5.35 (bending mode)	95.97	32.29 (bending mode)

Table 3 :

 3 Comparison of analytical and numerical frequencies for different slenderness ratio α for an embedded length of D=L/2and 𝑘 𝑠 = 10 MN/m 2 ..

Table 4 :

 4 Numerical values for pier frequencies, with and without the base effect.

		without a gap (Hz)	𝑓 with	𝑘 𝑠 (MPa)
			a gap (Hz)
	25	33.0	30.00	20
	20	23.5	22.5	25