Yves Jannot

Didier Lasseux

A new quasi-steady method to measure gas permeability of weakly permeable porous media

A new quasi-steady method for the determination of the apparent gas permeability of porous materials is presented in this paper along with the corresponding interpretative physical model derived from the unsteady flow equations. This method is mainly dedicated to the measurement of very low permeability of thin porous media, although thicker but more permeable samples may also be analyzed. The method relies on quasi-steady flow resulting from a (quasi) constant pressure maintained at the inlet face of the sample. Gas flow-rate, as low as 3 × 10 -10 m 3 /s, is determined from the record of pressure increase in a reservoir connected to the outlet face of the sample. An estimate of the characteristic time, t c , to reach quasi-steady flow after imposing a constant pressure at the inlet is derived. It is validated by direct numerical simulations of the complete unsteady flow, clearly defining the required experimental duration for the method to apply. Experimental results obtained on rather permeable and thick rock samples are reported showing an excellent agreement of the measured permeability with that determined independently on the same sample whereas the experimental value of t c is also in very good agreement with the predicted one. The method is further employed on a composite material sheet allowing the identification of an apparent gas permeability of about 10 -23 m 2 .

I. INTRODUCTION

Accurate gas-permeability determination of porous materials having permeabilities lower than 1 mD (∼10 -15 m 2 ) remains a tricky task. Above this typical value, a classical measurement technique on a homogeneous 1D sample consists in performing a stationary Darcy experiment. 1,2 This can be carried out using an experimental device as the one schematically represented in Figure 1.

From the measurement of the inlet volumetric flow-rate Q 0 , pressure drop P, and inlet pressure P 0 under steady flow conditions, gas permeability k can be deduced from the following relationship: 3 μeQ 0 P 0

S P = k P m . ( 1
)
To arrive at this result, the gas is supposed to be ideal while the flow is assumed to be isothermal and slow enough for inertial effects to be negligible. Moreover, compressibility effects are assumed to be significant at the sample scale only, not at the pore scale. An alternative method, that makes use of a variable charge, was proposed to measure concrete permeability, 4 neglecting, however, Klinkenberg effects (see below for this physical mechanism). Results obtained with the method were successfully compared with a technique developed earlier. This method however neglects air compressibility and makes use of a water column in a capillary as a pressure gauge. The meniscus, above which a low pressure is imposed to perform the measurement, might be subjected to evaporation that was not taken into account in the description. 5 a) Author to whom correspondence should be addressed. Electronic mail: didier.lasseux@ensam.eu. Tel.: +33 556 845 403. Fax: +33 556 845 436.

While dealing with materials having low permeability, Klinkenberg effects may be significant. [START_REF] Klinkenberg | API Drilling and Production Practice[END_REF] They are the macroscopic signature of gas molecules collisions at pore walls which effects become non negligible compared to gas molecules/molecules collisions. [START_REF] Carman | Flow of Gases through Porous Media[END_REF] Their impact can be estimated with the Knudsen number K n = λ/l, λ being the mean free-path of gas molecules at the pressure and temperature under consideration and l, a characteristic size at the pore scale. Typically, Klinkenberg effects are significant when Kn is in the range 0.01 to 1. In this range, they can be taken into account using a classical slip flow description at the pore-scale while at the macroscopic scale, they are characterized by a Klinkenberg coefficient b having the unit of pressure. When Kn is extremely small, i.e., when the Klinkenberg coefficient b of the gas and porous material under concern is much smaller than P m , k in Eq. ( 1) can be identified as the intrinsic permeability k l . However, when Klinkenberg effects are significant, k must be understood as the apparent gas permeability at P m . The intrinsic permeability k l and Klinkenberg coefficient b may then be determined from several experiments performed under stationary conditions but at different values of P m using the analogue of Eq. (1), [START_REF] Wu | [END_REF][START_REF] Rushing | Klinkenberg-corrected permeability measurements in tight gas sands: Steady-state versus unsteady-state techniques[END_REF][START_REF] Blanchard | [END_REF][11] namely,

μeQ 0 P 0 S P = k l (P m + b). ( 2
)
From a practical point of view, when the above described steady-state method is used on weakly permeable materials (typically for k < 10 -18 m 2 for example), difficulties appear to measure the gas flow rate with enough accuracy. In fact, since this flow rate can be extremely small, less than few cm 3 per hour, conventional gas flowmeters are inappropriate and more sophisticated methods, that are more difficult to implement, FIG. 1. Experimental set-up used to perform steady-state gas permeability measurement. B.P.R. is a back pressure regulator allowing to impose an outlet pressure. must be envisaged. Time necessary to reach stationary flow, that can become extremely long due to its dependence on the inverse of the permeability as detailed below, has not been determined so far. An accurate estimate of this characteristic time is provided in this work.

The objective of the present work is to develop an alternative quasi-steady state method that is simple to implement, making use of an adequate flow-rate measurement technique. The corresponding interpretative model to determine gas permeability, derived from the complete unsteady flow model, is provided. An equivalent mass diffusivity is used yielding an accurate estimate of the characteristic time necessary to achieve quasi-steady flow in the proposed experimental procedure. This estimate is validated numerically from the solution of the complete unsteady flow model. As will be detailed below, the method is particularly well adapted, but not restricted, to thin and weakly permeable porous materials. Finally the method is illustrated with experimental results obtained on a thin composite material on the one hand and more permeable but much thicker rock samples on the other hand and on which both the intrinsic permeability and characteristic time of the quasi-steady regime are validated. The main advantages of the method are (i) the simplicity of the experimental device that can be carefully pre-designed by making use of the time limits and volume/pressure constraints that are detailed in this work; (ii) the simplicity of the interpretation of the measurement to derive the permeability (and eventually the Klinkenberg coefficient when the experiment is repeated at different pressure levels, just as with a pure steady-state experiment); (iii) the possible real-time control of the experiment by making use of the time limits and volume/pressure constraints required for the method to apply.

II. PHYSICAL MODEL AND THE NEW QUASI-STEADY METHOD

The basic idea of the technique developed in this work is inspired from non-stationary gas permeability measurement methods often referred to as "Pulse-Decay" or "Draw-Down" experiments and their variants. [START_REF] Bruce | Petroleum Transactions[END_REF][START_REF] Aronofsky | Method of Measuring Characteristics of Porous Materials[END_REF][START_REF] Jones | [END_REF][15][16][17][18][START_REF] Jannot | Symposium of the Society of Core Analysts[END_REF][START_REF] Jannot | Symposium of the Society of Core Analysts[END_REF][START_REF] Lasseux | [END_REF] In these methods, the inlet face of the sample is connected to a reservoir of volume V 0 while the outlet face is connected to a reservoir of volume V 1 which is finite for the Pulse-Decay or infinite (i.e., corresponds to the atmosphere) for the Draw-Down. A mass of gas is confined in the upstream reservoir at a larger pressure than the initial one in both the sample and downstream reservoir. At t = 0, the upstream reservoir is opened at the inlet edge of the sample letting the pressure pulse relax through the sample. From the record of the time evolution of the pressure in the upstream reservoir -or the pressure difference between the upstream and downstream reservoirs -the permeability and Klinkenberg coefficient can be identified by making use of the physical model detailed below.

A. Physical model

If we assume a homogeneous non-deformable sample and a 1D isothermal gas flow along the sample axis x at a very low Reynolds number -this is the case in practise -so that no significant inertial (or Forchheimer) effects are present, [START_REF] Blanchard | [END_REF]22 the mass and momentum conservation equations are

ε ∂ρ(x, t) ∂t + ∂(ρ(x, t)u(x, t)) ∂ x = 0, (3)
u(x, t) = - k μ ∂ p(x, t) ∂ x = - k l (1 + b/ p(x, t)) μ ∂ p(x, t) ∂ x . ( 4
)
In addition, if gas is considered as ideal ( p(x, t) = ρ(x, t) RT M ), which is a valid approximation for gases like N 2 or He at experimental operating pressures of the order of ten bars, these two equations can be combined using the relationship between the apparent and intrinsic permeabilities k

= k l (1 + b/p(x, t)) yielding ∂ 2 ( p(x, t) + b) 2 ∂ x 2 = ε μ k l ( p(x, t) + b) ∂( p(x, t) + b) 2 ∂t , (5)
or equivalently

∂ 2 ( p(x, t) + b) 2 ∂ x 2 = ε μ k p(x, t) ∂( p(x, t) + b) 2 ∂t . ( 6
)
When a completely unsteady experiment is considered, the upstream boundary condition (at x = 0) is such that the mass flow rate through the inlet face of the sample is exactly that out of the upstream tank. The mass flow rate through the inlet face is

-S p 0 (t)M R T k l (1+b/ p 0 (t)) μ ∂ p(x,t) ∂ x | x=0
, with p 0 (t) = p(x = 0, t), p 0 (t) being also the pressure in the upstream tank out of which the mass flow rate is

-V 0 M R T dp 0 (t) dt . Matching the two leads to dp 0 (t) dt = S μV 0 k l ( p 0 (t) + b) ∂ p(x, t) ∂ x at x = 0. ( 7
)
Following the same lines, the downstream boundary condition is given by

dp 1 (t) dt = - S μV 1 k l ( p 1 (t) + b) ∂ p(x, t) ∂ x at x = e, ( 8
)
where p 1 (t) is the pressure at x = e at t, i.e., in the downstream tank. The associated initial conditions are

p(x, 0) = P 1 for 0 < x ≤ e, (9)
p(0, 0) = P 0 , (10) in which P 1 = p 1 (t = 0) and P 0 = p 0 (t = 0). This model is a general one for 1D gas flow. Indeed, for a classical steady-state experiment, the stationary version of Eq. ( 5) can be integrated twice making use of the two bound-ary conditions p(x = 0) = P 0 and p(x = e) = P 1 yielding

p(x) = (P 1 + b) 2 -(P 0 + b) 2 e x + (P 0 + b) 2 1/2 -b. ( 11
)
When this solution is introduced back in Darcy's law (Eq. ( 4)), the volume flow rate Q 0 at the entrance x = 0 is readily obtained as

Q 0 = u (0) S = - k l μ (P 1 + b) 2 -(P 0 + b) 2 2 e P 0 S, (12)
which, when re-arranged, gives the result in Eq. ( 2).

In the case of a completely unsteady experiment, the model given by Eqs. ( 6)-( 10) can be used in an inverse procedure in order to estimate k l and b. This is achieved by minimizing the residual between the measured signal p 0 (t) (or p 0 (t)p 1 (t)) and the corresponding one computed from the model as illustrated in Lasseux et al. [START_REF] Lasseux | paper presented at the Pressure Vessels and Piping Conference[END_REF] The method avoids major difficulties of the steady-state technique pointed out in the introduction, although the underlying physical model involves the porosity, ε, of the material. Because the simultaneous identification of the porosity with k l and b in the inverse procedure is extremely difficult due to the poor sensitivity of the pressure decay signal to ε, this parameter must be known a priori and provided as an input for inversion. Moreover, it was shown in an earlier analysis that the input value of ε must be very accurate for a reliable estimate of k l and b. [START_REF] Jannot | Symposium of the Society of Core Analysts[END_REF] While dealing with weakly permeable materials having also frequently small porosity, the accurate determination of this parameter, which requires a separate dedicated experiment, is a difficult task. It must also be noted that the inverse process can be cumbersome so that the overall procedure makes a simpler method desirable. This is achieved with the method described below which takes advantage of both the steady and unsteady state techniques.

B. The new quasi-steady method

This method is inspired from the unsteady technique detailed above but, instead of letting the gas pressure pulse in the upstream reservoir relax through the sample, a constant pressure P 0 is maintained in V 0 . This can be achieved either by continuously adjusting the pressure or more practically by simply choosing V 0 large enough so that, as can be easily inferred from Eq. ( 7), the variation δp 0 of p 0 (t) over the whole period of measurement remains extremely small compared to P 0 . A criterion to satisfy this constraint is given in Sec. II C 2. Moreover, in this new method, we require the measurement of p 1 (t) to be performed at time t > t c (see Sec. II B 2 below for the estimation of t c ) so that the capacitive term in the right-hand side of Eq. ( 5) or (6) becomes negligible. This actually corresponds to the time necessary to reach steady flow when a constant pressure is applied at the inlet. Under these circumstances, the experiment is quasi-steady, the unsteady character remaining only through the time evolution of p 1 (t), as indicated by the boundary condition (8). The increase of p 1 (t) is actually recorded as an indirect measurement of the flow rate through the sample. A schematic repre-FIG. 2. Experimental configuration for the quasi-steady method.

sentation of the experiment is reported in Figure 2. To keep the method simple, the focus is on the determination of the permeability that can be achieved in a straightforward manner from the downstream pressure record restricted to a period t c < t < t f , corresponding to a pressure increase δp 1 that remains small compared to P 1 . The development of the method is proposed in the Sec. II B 1 and the associated constraint t < t f is made clear in Sec. II C 1.

Determination of the permeability

For t > t c , (see Sec. II B 2 below for the estimation of t c ), the flow can be described by a quasi-steady version of the initial boundary value problem ( 6)- (10), in which the accumulation term in Eq. ( 6) is removed and the upstream boundary condition ( 7) is replaced by a constant pressure P 0 at x = 0, hence yielding

∂ 2 ( p(x, t) + b) 2 ∂ x 2 = 0, (13)
p (0, t) = P 0 , ( 14
)
dp 1 (t) dt = - 1 2 S μV 1 k l ∂( p(x, t) + b) 2 ∂ x at x = e. ( 15
)
The analytical solution of this problem can be simply written as

( p(x, t) + b) 2 = - 2μV 1 S k l dp 1 (t) dt x + (P 0 + b) 2 0 ≤ x ≤ e, (16)
which provides an expression for the rate of pressure increase in the downstream reservoir given by

dp 1 (t) dt = k l S 2μV 1 e ((P 0 + b) 2 -( p 1 (t) + b) 2 ). ( 17
)
This form could be used in an inverse procedure to identify k l and b with the major drawback that the time derivative of the experimental record of p 1 (t) is required, making the technique inefficient. A much simpler and easy-to-use method can be derived, however, by noticing that, when p 1 (t) is recorded over a time interval [t c , t f ] such that the downstream pressure increment δp 1 is extremely small compared to P 1 , dp 1 (t)/dt remains constant. Under these circumstances, the determination of the permeability can be made from a simple linear regression on p 1 (t) on the time interval [t 0 , t 1 ], t 0 ≥ t c , t 1 ≤ t f , according to Equivalently, since

k l = 2 μV 1 S e (P 0 + b) 2 -(P 1 + b) 2 a = μV 1 S e (P m + b) (P 0 -P 1 ) a. ( 18
)
k = k l (1 + b/P m ) , ( 19
)
we have

k = 2 μV 1 S e P 2 0 -P 2 1 a = μV 1 S e P m (P 0 -P 1 ) a, ( 20
)
where a is the slope of p 1 (t) estimated on the interval [t 0 , t 1 ]. More precise forms can be used upon averaging Eq. ( 17) over the time interval on which a is estimated leading to

k l = 2 μV 1 S e (P 0 +b) 2 -( p1 +b) 2 a or k = 2 μV 1 S e P 2 0 -p2 1 a where p1 = 1 t 1 -t 0 t 1 t 0 p 1 (t) dt.
For simplicity, the forms given by Eqs. ( 18) and ( 20) are kept in the sequel of the paper. Apart from the simplicity of the signal interpretation, it shall be noted that an evident advantage of the present method is that the mass flow rate,

V 1 M RT dp 1 (t)
dt , which can be extremely small, is measured by recording p 1 (t) that can be performed accurately using rather standard pressure sensors.

As it is obvious from Eqs. ( 18) and ( 20), the method is such that a single experiment allows to estimate k only but not k l and b separately. To estimate both parameters, several experiments, carried out at different pressures, are required as in the case of a steady experiment. For this reason, the method is particularly well suited to thin samples for which experiments are shorter due to the dependence of t c on e 2 as shown in Sec. II B 2.

Characteristic time t c and validation

In this section, we shall first determine the characteristic time for the quasi-steady flow hypothesis to hold and for the solutions (18) or (20) to apply. This is achieved from scaling arguments on the dimensionless form of the quasi diffusion equation (6). To do so, it is convenient to start from equation (6) of the complete unsteady model and let φ(

x, t) = (p(x, t) + b) 2 yielding ∂φ * (x * , t * ) ∂t * = D m t c e 2 ∂ 2 φ * (x * , t * ) ∂ x * 2 , ( 21
)
where D m is a pseudo mass diffusivity that is pressure dependent. In Eq. ( 21), all quantities with the superscript * are dimensionless and obtained by scaling all variables by the corresponding characteristic reference value (noted below with the subscript c). In order to correctly identify the characteristic time t c , associated to the propagation of the pressure signal through the sample, proper scalings on φ(x, t), p(x, t) and x that make them O(1) must be chosen. These scalings are given by φ c = (P 0 + b) 2 , p c = P 0 and x c = e yielding D m = k P 0 p * ε μ . Since all the dimensionless variables in Eq. ( 21) are now O(1), the pre-factor of the diffusive term must also be O(1) yielding

D m ∼ k P 0 ε μ , ( 22
)
and hence

t c ∼ = ε μ e 2 k P 0 ∼ = ε μ e 2 k l (1 + b/P m ) P 0 . ( 23
)
This last relationship represents an estimate of the criterion for the quasi-steady hypothesis to hold. Indeed, when the flow is observed at a time larger than t c , the unsteady term in the left hand side of Eq. ( 21) becomes insignificant.

This prediction is now checked using direct numerical simulations performed on the complete unsteady model given by Eqs. ( 6)- (10). Details of the discrete procedure used to solve the problem are provided in Appendix A.

Numerical simulations were performed with values of the parameters reported in Table I using nitrogen (μ ∼ = 1.8 × 10 -5 Pa s) with P 1 = 10 5 Pa, a very large value of V 0 (V 0 = 10 3 m 3 ) ensuring a constant P 0 and S = 1.96 × 10 -3 m 2 . Values of the characteristic time t c obtained from Eq. ( 23) are also indicated in this table. The Klinkenberg coefficient b was correlated to k l according to an empirical relationship proposed in the literature, obtained for granular porous media in the range of permeability 10 -15 -10 -17 m 2 and given by 14

b = 0.189 k -0,36 l . ( 24
)
Case 5 is such that all parameters are kept identical to Case 4 except k l and b that were calculated to keep the term k l (1 + b/P m ) unchanged with respect to this later case. Under these conditions, the predicted value of t c from Eq. ( 23) remains the same in both cases.

The evolution of p 1 (t) obtained from the numerical simulations is reported in Figure 3, for each of the four configurations under consideration while the corresponding evolutions of dp 1 (t)/dt are represented in Figure 4 in which the prediction of Eq. ( 17) in the quasi-steady regime was also reported as open symbols.

Results clearly show the three distinct regimes characterized by (i) the "pore-filling" period for which the downstream pressure remains unchanged at the very early stage of the experiment; (ii) a rapid evolution of p 1 (t) during a fully unsteady regime at intermediate times, up to t = t c ; (iii) the quasi-steady regime of particular interest in this work for times larger than t c .

First, the evolution of dp 1 (t)/dt obtained from the numerical solution of the complete unsteady model is perfectly well predicted by the solution developed above (Eq. ( 17)) in the quasi-steady regime (t > t c ). Second, in this last regime ) 4, validating the use of Eqs. ( 18) and ( 20) to determine k l and k. Finally, the characteristic time of the transient behaviour is perfectly well predicted by Eq. ( 23) as shown by the vertical dashed lines materialising the estimated values of t c (see Figure 3). In particular, the excellent prediction of t c , that remains unchanged for Cases 4 and 5, must be noted. As expected, t c only depends on the characteristics of the sample and on the pressure levels P 0 and P 1 . In particular it does not depend on V 1 (see Table I where several values of V 1 were considered). This validation confirms the physical relevance of the estimate for t c which clearly highlights that the quasi-steady method developed in this work is better suited when ε e 2 /k l is small in order to keep the experimental duration small.

t c =276 s Case #3 t f =6330 Eq. ( 15

C. Time t f and upstream reservoir volume constraints 1. Time constraint t f

The simple method of determination of k l and k (Eqs. ( 18) and ( 20)) developed in this work requires that the downstream pressure increase δp 1 remains small compared to P 1 over the experimental observation of p 1 (t). This implies that the downstream pressure record must be restricted to the time interval [t c , t f ], corresponding to a quasi-linear evolution of p 1 (t). The time constraint t f is now made explicit by noticing that dp 1 (t)

dt ∼ = k l S
μV 1 e (P m + b) (P 0 -P 1 ) for t c < t < t f whereas dp 1 (t) dt is increasing with time for t < t c . As a consequence, δp 1 can be overestimated by writing

δp 1 < k l S μV 1 e (P m + b) (P 0 -P 1 ) t f . ( 25
)
The constraint δp 1 P 1 can hence be safely expressed as

t f μV 1 e k l S P 1 (P m +b) (P 0 -P 1 ) = μV 1 e k S P 1 P m (P 0 -P 1 ) = P 1 a . ( 26
)
As shown in Figure 4, the limit t f = 0.2 μV 1 e k S P 1 P m (P 0 -P 1 ) materialized by dashed lines is a satisfactory criterion for all the situations under investigation.

It must be noted that the estimate of t c and the constraint t f can be used as a real-time feedback control of the experiment. Indeed, k can be determined from Eq. ( 20) on an interval [t 0 , t 1 ] where p 1 (t) is apparently linear while carrying out the experiment. An estimate of t c (Eq. ( 23)) and the constraint (26) can be then post-compared to t 0 and t 1 to check the relevance of the permeability estimation and redefine the time interval [t 0 , t 1 ] if necessary while giving an indication on whether the experiment can be carried on. Note also that the expression of t c along with the constraint (26) can be employed to pre-design the experimental setup when the permeability and porosity ranges to be investigated are defined.

Upstream volume constraint

When the experiment is performed with a fixed upstream volume V 0 to ensure a constant pressure P 0 , this volume must be large enough for the simple quasi-steady method developed in this work to be effective which means that the pressure decrease δp 0 must remain small compared to P 0 for t up to t f . The constraint can be made clear by estimating the corresponding mass of gas, δm, flowing out of the upstream reservoir over t f . Over the time period t f , a fraction δm 1 of this mass of gas contributes to the pressure rise in the downstream reservoir while the remaining part, δm p , is accumulated in the pores to increment the pressure in the sample. An overestimate of δm 1 can be readily obtained from (25) and the assumption of an ideal gas leading to δm 1 < k l S μ e (P m + b) (P 0 -P 1 ) t f M RT . An overestimate of δm p may be obtained by considering that gas accumulated in the pores corresponds to a pressure increase from P 1 to P m yielding δm p < (P m -P 1 ) ε e S M RT so that an overestimate of the upstream pressure decrease δp 0 = δm 1 +δm p V 0

RT

M can be expressed as

δp 0 < 1 V 0 k l S μ e (P m + b) (P 0 -P 1 ) t f + (P m -P 1 ) ε e S . ( 27
)
In Appendix B, a more accurate estimate -but of less practical usefulness -of δm p (and hence of δp 0 ) is derived.

The upstream constraint δp 0 /P 0 1 ensuring quasisteady conditions can finally be derived from Eq. ( 27) by requiring the following sufficient condition:

P 0 -P 1 P 0 V 0 k l S μ e (P m + b) t f + ε e S 2 = P 0 -P 1 P 0 V 0 k S μ e P m t f + ε e S 2 1. (28)
Again, this last constraint can be used either a posteriori to check the validity of a measurement or as a design guideline for an experimental setup to be used on given permeability-porosity intervals.

III. EXPERIMENTS

To illustrate the method described above, experimental results obtained on two different types of materials (rocks and composite material) are reported in this section. We focus on the determination of the apparent permeability, k, at a given pressure P m , keeping in mind that k l and b could be determined from repeated identical experiments at different values of P m as documented elsewhere. [START_REF] Wu | [END_REF][START_REF] Rushing | Klinkenberg-corrected permeability measurements in tight gas sands: Steady-state versus unsteady-state techniques[END_REF][START_REF] Blanchard | [END_REF][11]

A. Materials and method

The device used for the experiments is represented in Figure 5 and is composed of the following parts:

• A gas cylinder (N 2 ) is connected to the upstream reservoir through a pressure regulating valve V D and a three-port valve V 3 . This valve enables to connect the upstream reservoir either to the gas supply at P 0 or to the atmosphere. • A digital manometer M1 (Keller LEO1) measures the pressure P 0 in the upstream reservoir. • A differential pressure sensor M2 (a Hottinger Baldwin Messtechnik 0-100 mbar for the experiment on the composite material and a Rosemount 3051S 0-620 mbar for experiments on rock samples) is employed to measure the pressure variation in the downstream reservoir. The high pressure input of the sensor is connected to the downstream reservoir and to a valve V 2 enabling the connection to the atmosphere. The low pressure input of the sensor is connected to a small chamber drilled in the sample holder, close to the downstream reservoir. The small chamber, that can be connected to the atmosphere through the valve V 1 , is used as a reference pressure cavity that remains at the same temperature as the downstream chamber so that temperature perturbation are (partially) filtered out from the pressure signal recorded in the downstream reservoir. Moreover, potential atmospheric pressure fluctuations are also filtered out since the cavity remains closed during measurement. • A digital data logger (Almemo 2290-5) records simultaneously the differential pressure, the temperature difference with respect to the enclosure, T 1 = T 1 -T 0 , of the tubing connecting the volume V 1 to the high pressure input of the pressure sensor and the temperature difference with respect to the enclosure, T 2 = T 2 -T 0 , of the sample holder, close to the reference pressure cavity. Temperatures are measured with type K thermocouples. Time sampling of the record can be varied in the range 0.1 s -1h.

The measuring cell, materialized by the dashed area in Figure 5, was adapted to the type of sample under consideration. For the measurement on rock samples reported in Sec. III B 1, a Hassler sleeve core holder, represented in Figure 6(a), was used as usually employed for this type of sample. [START_REF] Al-Hinai | [END_REF] In the case of the composite material sheet (see Sec. III B 2), the measuring cell schematized in Figure 6(b) was employed. In this last configuration, a plate of sintered porous material, which permeability (about 10 -12 m 2 ) was much larger than that of the composite material, was placed below the composite sheet in the downstream reservoir. This plate avoids significant deformation of the composite material which would result from the applied pressure in the upstream reservoir and would bias the pressure rise record p 1 (t). The operating process is the following:

• the sample is positioned into the device and a confining or clamping pressure is applied; • the device is placed in a temperature controlled enclosure (typically at 30 • C) while valves V 1 , V 2 and V 3 (as well as V 4 , (see Figure 6(a) for experiments on rock samples) are open to the atmosphere; • valves V 1 and V 2 are closed when thermal equilibrium is reached and the atmospheric pressure is measured separately in the enclosure. Although equilibrium with the enclosure is reached, some low frequency temperature fluctuations (for instance on a day/night cycle for a long experiment) may still persist as will be further shown below on experiments carried out on the composite material; • the pressure P 0 is set to the chosen value by adjusting the pressure regulating valve V D ; • the valve V 3 is switched to connect the gas supply to the upstream reservoir. For experiments performed on rock samples, the valve V 4 (see Figure 6(a)) is closed.

Recording of the downstream differential pressure and temperatures is started.

B. Experimental results

Experimental results were obtained on two different thick and relatively permeable rock samples. An additional result on a thin and very poorly permeable composite material sample is also reported to highlight the capability of the method to measure extremely small permeabilities.

Rock samples

The quasi-steady method was tested on two rock samples R1 and R2 which characteristics are reported in Table II. The samples were cylindrical cores of circular cross section and their characterization was made in the axial direction. For comparison purposes, the intrinsic permeability, k l and Klinkenberg coefficient, b, were pre-determined independently from the present method using a fully unsteady (Draw-Down) method described above (see Sec. II A) and further detailed elsewhere. [START_REF] Jannot | Symposium of the Society of Core Analysts[END_REF][START_REF] Lasseux | [END_REF] Since the sensitivity of the upstream pressure decay signal, p 0 (t), to the porosity, ε, was poor, this parameter was measured with a classical gas pycnometry (N 2 ) experiment. The two parameters k l and b were determined using an inverse procedure based on the physical model ( 6)- (10). Accuracy of the procedure is excellent since residuals on p 0 (t) (i.e., the difference between the measured signal and the simulated one at inversion convergence) are less than 100 Pa while the Draw-Down experiments were carried out with N 2 at P 0 up to 11 × 10 5 Pa and 16 × 10 5 Pa. Accuracy is corroborated by the relative standard deviation, σ k l and σ b , on k l and b reported in Table II.

Quasi-steady experiments were carried out with the parameters V 0 , V 1 , P 0 , and P 1 reported in Table III where we have also indicated the times t f at which experiments were ended. The expected apparent permeability, k, at these values of the pressures along with the corresponding expected characteristic time t c reported in this table were estimated from data of Table II and Eqs. ( 19) and ( 23) respectively.

Experimental results of the evolution of p 1 (t) are represented in Figure 7 where the transient periods are clearly appearing followed by the quasi-steady periods.

In order to check the validity of the characteristic time estimates t c reported in Table III, the time derivative dp 1 (t) dt of the downstream pressure was determined. This was done from a polynomial least square fit of degree 6 (sample R1) and 4 (sample R2) performed on the signal p 1 (t). Evolutions of dp 1 (t) dt are represented in Figure 8 in which we have also reported the predicted values of t c as vertical dashed lines. These predictions are in excellent agreement with the emergence of the linear evolution of p 1 (t) characteristic of the quasi-steady regime observed experimentally, confirming the validity of the estimate. Moreover, the time derivative remains constant up to t f for both experiments, confirming that permeability may be estimated following the quasi-steady model for t > t c .

Next, apparent permeabilities were estimated from linear fits performed on p 1 (t) for t c ≤ t ≤ t f . The fits are represented in Figure 9 showing the relevance of the linear approximations over these time intervals with excellent determination coefficients. The slopes a were found to be 4.0 Pa s -1 (sample R1) and 9.4 Pa s -1 (sample R2), which upon using Eq. ( 20), yields k ∼ = 5.15 × 10 -17 m 2 (sample R1) and k ∼ = 12.12 × 10 -17 m 2 (sample R2). Again these two results are in excellent agreement with predicted values reported in Table III and validate the approach proposed in this work. TABLE III. Experimental parameters for tests performed up to t = t f on rock samples R1 and R2. k and t c are the expected apparent permeability at (P 0 + P 1 )/2 and characteristic time of the transient period, both being evaluated from data of Table II.

Sample

V 0 (×10 -3 m 3 ) V 1 (×10 -3 m 3 ) P 0 (×10 5 Pa) P 1 (×10 Finally, the time and volume constraints expressed by the relationships ( 26) and (28) respectively may be checked a posteriori. According to our experimental parameters, the right hand side of the relationship ( 26) is roughly 2.53 × 10 4 s (sample R1) and 1.03×10 4 s (sample R2), values that are much larger than the values of t f used to estimate k for both samples, confirming that the time constraint is fully satisfied. This is in agreement with observations made on Figure 8. Similarly, the left hand side in relationship ( 28) is evaluated to be approximately 1.96 × 10 -2 (sample R1) and 9.01 × 10 -2 (sample R2), both values being smaller than 1 which indicates that the experimental parameters were conveniently chosen for the quasi-steady experiment. In Figure 10, we have reported the evolutions of the upstream pressure for samples R1 and R2.

As can be inferred from these figures, the relative upstream pressure decrease, δp 0 /P 0 , is roughly 2.0 × 10 -2 (sample R1) and 3.7 × 10 -2 (sample R2). It is to be compared to the two values calculated from the left hand side of (28) and mentioned above which should be upper bounds. This is the case for sample R2, while, although experimental noise on p 0 (t) should be considered, it is not for sample R1. However, one must keep in mind that this upper bound is a restrictive one and when better estimated according to the expression provided in Appendix B, one finds 2.06 × 10 -2 instead of 1.96 × 10 -2 for sample R1 which better corresponds to the expected overestimate of δp 0 /P 0 .

Composite material sheet

An additional experiment was carried out to determine the transverse permeability of a thin and very poorly permeable medium, namely a sheet of composite glass-fiber/epoxy material. The experimental parameters were such that e = 2 mm, S ∼ = 17.67 × 10 -3 m 2 , V 1 = 10 -4 m 3 , P 0 ∼ = 6.08 × 10 5 Pa, P 1 ∼ = 1.013 × 10 5 Pa. Since the upstream reservoir was directly connected to the gas cylinder through the valves V d and V 3 , V 0 can be considered as infinite in this experiment so that the constraint expressed by the relationship (28) is safely verified.

The time record of p 1 (t) obtained during this experiment is represented in Figure 11(a) while the evolution of the temperatures T 1 (t) and T 2 (t) are reported in Figure 11(b). As expected, the downstream pressure rise is very small (less than 5 mbars over roughly 18 hours of record) suggesting that the permeability of the material is indeed extremely small. As can be seen in Figure 11(b), the temperatures remained relatively constant for 5.5 ≤ t ≤ 15 h. As for the data on the rock samples, a least square polynomial fit of degree 6 was performed on p 1 (t) in the interval of time 5 ≤ t ≤ 10 h where the transition to the quasi-steady regime obviously occurs. The time derivative of this fit is depicted in Figure 12(a) indicating that t c can be reasonably identified as t c ≈ 8 h. From this value, the slope a of p 1 (t) was identified in the time interval 8 ≤ t ≤ t f = 15 h (see Figure 12(b)) where the temperature is free of significant fluctuations yielding a ∼ = 10 -2 Pa s -1 with a very satisfactory determination coefficient. When reported in the expression (20) for the apparent permeability determined from this quasi-steady regime, one finds k ∼ = 1.1 ×10 -23 m 2 which is far below the classical values of permeability determined with conventional measurement techniques. Note that the time constraint expressed by the relationship ( 26) is verified for this measurement. Indeed, the right hand side in ( 26) is evaluated to be ∼10 7 s (i.e., ∼2900 h), which is much larger than the value of t f used to estimate a, validating our experimental protocol for this material.

In the absence of data for the porosity of the material, no prediction of t c can be performed. However, if ε is postevaluated from the graphically identified value of t c (t c ≈ 8 h), this yields ε ≈ 0.3%, a value that remains in the range of porosity reported in the literature for this type of material. 25,[START_REF] Lin | the 17th World Conference on Non-Destructive Testing[END_REF]

IV. CONCLUSION

A method to determine the permeability of porous media, that is particularly well adapted to materials for which the ratio εe 2 /k is small (ε is the porosity, e the sample thickness and k the apparent permeability of the material), i.e., to very poorly permeable porous media and/or thin samples, was proposed in this work. The method is based on the quasisteady regime occurring during 1D flow when the sample is subjected to a (quasi) constant upstream pressure. The flow rate is indirectly measured through the pressure rise in the downstream reservoir, which is a net advantage of the technique. Indeed, since standard sensitive pressure sensors can be used, extremely small flow rates can be detected allowing the determination of very small permeabilities. The physical model to interpret the experiment was carefully derived from the complete unsteady model of flow showing that part of this quasi-steady regime is characterized by a linear evolution of the downstream pressure signal. A simple interpretation of this part of the signal can be made to determine the permeability and the Klinkenberg coefficient from repeated experiments at different mean pressure levels. Along with the fact that the experimental setup for this method is simple to implement, the simple interpretation makes the method of particular interest.

From scaling arguments, the characteristic time at which the quasi-steady regime is observed was explicitly derived. It was validated by numerical simulations of the complete unsteady model in a variety of operating conditions. This characteristic time defines the lower time-limit of the experimental downstream pressure signal to be considered with this method. The constraint on the upper time-limit was also determined, clearly defining the time interval of interest for the method to apply. Moreover, when the experiment is to be carried out with an upstream reservoir of finite volume, which is of practical common use, a constraint was provided for the method to remain valid. All these constraints can be advantageously used either as a set of guidelines to design an experiment and/or as a real time control tool to assert that the experiment is carried out under appropriate conditions.

Experiments carried out on two different rock samples were presented, demonstrating the validity of the method, both on the characteristic time of the quasi-steady regime and on the permeability estimates. The capability of the method to determine very low permeabilities (down to 10 -23 m 2 ) was illustrated by a measurement in the transverse direction of a composite material sheet. This fast and simple method, taking advantage of both fully steady and unsteady techniques, is a promising tool for permeability characterization of very poorly permeable materials like samples of tight-gas or gasshale reservoirs, cap rocks of gas or nuclear waste storage areas, sealing or composite materials, etc., all materials on which permeability measurement remained a challenging task so far and for which reliable and accurate methods are still scarce.

An alternative and more accurate expression of the constraint given by the relationship (28) is developed in this Appendix.

We start from a mass balance of gas between t = 0 and t = t f , given by P 0 V 0 + P 1 (V 1 + ε e S) = (P 0 -δp 0 ) V 0 + (P 1, the integral term in Eq. (B1) can be explicitly calculated as e 0 ε S p(x, t f )dx ∼ = 2 3 ε e S (P 0 +b) 3 -(P 1 +b) 3 (P 0 +b) 2 -(P 1 +b) 2 . When inserted back into (B1) while taking into account the upper bound of δp 1 in (B2), one obtains an overestimate of δp 0 as Although this form is more accurate than (28) and can be easily shown to be less restrictive, it is not of real practical usefulness. Indeed, since it involves the value of b which is not known while performing the measurement of k at a given P m , a real-time feedback control of the validity of the experiment is impossible with this form.

  1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 t

FIG. 5 .

 5 FIG. 5. Schematic representation of the experimental device. The measuring cell is one of those represented in Figure 6 depending on the sample type.

FIG. 6 .

 6 FIG.6. Two measuring cells used for (a) a rock sample (Hassler sleeve coreholder); (b) a composite material sheet.

FIG. 7 .FIG. 8 .

 78 FIG. 7. Downstream pressure evolutions (a) sample R1; (b) sample R2. For convenience of the representation, the signal was truncated at t = 500s.

FIG. 10 .

 10 FIG. 10. Upstream pressure (a) sample R1; (b) sample R2.

FIG. 11 .

 11 FIG.11. Experimental results on the composite material sheet. Evolution of (a) the downstream pressure; (b) temperature differences T 1 and T 2 .

FIG. 12 .

 12 FIG. 12. (Color online) (a) Time derivative of the downstream pressure for 5 ≤ t ≤ 10. (b) Least square linear fit on p 1 (t) for 8 ≤ t ≤ t f = 15: a ∼ = 36.9 Pa h -1 .

+ ε e S 2 3 (+ ε e S 2 3 (

 33 P 0 + b) 3 -(P 1 + b) 3 (P 0 + b) 2 -(P 1 + b) 2 -P 1 , (B4)and hence an alternative form to (28) developed in the paper given by 1 P 0 V 0 kS μ e P m (P 0 -P 1 )t f P 0 + b) 3 -(P 1 + b) 3 (P 0 + b) 2 -(P 1 + b) 2 -P 1 1. (B5)

Nomenclature a slope of p 1 (

 1 t) in the linear quasi-steady regime (Pa/s) b Klinkenberg coefficient (Pa) D m pseudo mass diffusivity (D m = k P 0 p*/ε μ) (m 2 /s) e sample thickness (m) k apparent gas permeability (m 2 ) k l intrinsic permeability (m 2 ) l characteristic length at the pore scale (m) M molar mass of gas (kg/mol) p(x, t) pressure at x and t (Pa) p 0 (t) absolute pressure in the upstream reservoir at t (Pa) p 1 (t) absolute pressure in the downstream reservoir at t (Pa) P 0 initial absolute pressure in the upstream reservoir (=p 0 (t = 0)) (Pa) P 1 initial absolute pressure in the downstream reservoir (=p 1 (t = 0)) (Pa) P m arithmetic mean pressure (P m = (P 0 + P 1 )/2) (Pa) Q 0 volumetric flow-rate (m 3 s -1 ) R ideal gas constant (R = 8.314 J/mol/K) S sample cross sectional area (m 2 ) t time (s) t c characteristic time to reach the quasi-steady regime (s) t f final time of the experiment (s) T temperature (K) u(x,t) filtration velocity (m/s) V 0 volume of the upstream reservoir (m 3 ) V 1 volume of the downstream reservoir (m 3 ) x coordinate along the sample axis (m) Greek δp 0 pressure variation in the upstream reservoir over t f (Pa) δp 1 pressure variation in the downstream reservoir over t f (Pa) P pressure drop ( = P 0 -P 1 ) (Pa) T 1 temperature difference between the temperaturecontrolled enclosure and the tubing connecting V 1 to the pressure sensor ( T 1 = T 1 -T 0 ) (K) T 2 temperature difference between the temperaturecontrolled enclosure and the sample holder, close to the reference pressure cavity ( T 2 = T 2 -T 0 ) (K) ε porosity φ pseudo potential (φ = (p + b) 2 ) (Pa 2 ) λ mean free-path of gas molecules (m) μ gas dynamic viscosity (Pa s) ρ gas density (kg/m 3 )

TABLE I .

 I Parameters employed in the numerical simulations.

	Case nos.	k l (m 2 )	b (×10 5 Pa)	ε	e (mm)	P 0 (×10 5 Pa)	V 1 (m 3 )	t c (s) (Eq. (23))
	1	1 0 -16	0	0 . 2	5 0	6	3 ×10 -3	150
	2	1 0 -18	0	0 . 1	1 0	8	3 ×10 -4	225
	3	1 0 -19	13.1	0.05	10	9	3×10 -4	276
	4	1 0 -20	0	0.05	5	10	10 -4	2250
	5	1 . 5 5 ×10 -21	30	0.05	5	10	10 -4	2250

  t ≤ t f (see Sec. II C 1 for this time constraint), p 1 evolves linearly, as shown in Figure

	1.2	t c =2250 s	t f =18520 s
	(Pa.s ) -1	Numerical simulation
	dp 1(t)/dt	Eq. (15) t>t c			Case #4
					Case #5
	1.E+00	1.E+02	1.E+04	1.E+06
		t (s)

FIG. 4. (Color online) Evolutions of dp 1 (t)/dt in each of the five cases of Table I. The characteristic time t c (see Sec. II B 2) and time constraint t f (see Sec. II C 1) are materialised by the vertical dashed lines.

and

TABLE II .

 II Characteristics of the rock samples R1 and R2.

	Sample	k l (×10 -17 m 2 )	σ k l (%)	b (×10 5 Pa)	σ b (%)	ε	e (mm)	S (×10 -3 m 2 )
	R1	3.8	0.05	1.16	0.42	0.055	60.3	1.134
	R2	9.8	0.01	0.82	0.10	0.127	60.2	1.138

  1 + δp 1 ) V 1 Noticing that p 1 (t f ) = P 1 + δp 1 while δp 1

				e
				+	ε S p(x, t f )dx,	(B1)
				0
	While we shall overestimate δp 1 as
	δp 1 <	k l S μV 1 e	(P m + b) (P 0 -P 1 ) t f =	kS μV 1 e	P m (P 0 -P 1 ) t f ,
							(B2)
	the expression p(x, t) can be extracted form Eqs. (16) and (17)
	to give
	p(x, t f ) = (P 0 +b) 2 -2	p 1 (t f )+b P 1 +b	(P m + b) (P 0 -P 1 ) e	x	1/2	.
							(B3)
							P 1 +b

ACKNOWLEDGMENTS

The financial support from EDF is gratefully acknowledged.

We thank Sandra PROFICE for raising interesting questions that helped clarifying some developments in this paper.

APPENDIX A: NUMERICAL SCHEME FOR THE 1D UNSTEADY GAS FLOW

In this Appendix, we detail the numerical procedure used to solve the unsteady gas flow through the medium. The initial boundary value problem given by Eqs. ( 6)-( 10) is first re-written using the pseudo-potential φ(x, t)

These equations, where we have denoted α = ε μ k l , β = μ V 0 k l S , and δ = μ V 1 k l S , are discretized using a finite difference scheme which is second order in space and a first order explicit Euler scheme in time. With the notation φ n i = φ ((i -1) x, (n -1) t), this scheme is given by

where m is the number of space discretization nodes. Stability of the overall algorithm is subjected to a criterion on the time step due to the explicit character of the time scheme. This criterion is a classical one for a diffusion-like equation and is such that t < min i ( ε μ x 2 2 k l (φ i ) 1/2 ).


Vassili Pustovalov

Hung Duong

Denis Pham

Kouame

Deep Unfolding RPCA for High Resolution Flow Estimation

Keywords: Blind Deconvolution (BD), Blood Flow, Clutter Separation, Robust Principal Component Analysis (RPCA), Power Doppler, High Resolution

Numerous techniques have been proposed to produce high precision blood flow or contrast-enhanced ultrasound estimates from fast ultrasound sequences. Among them, robust principal component analysis (RPCA)-based methods are known as superior to most state-of-the-art techniques. In particular, these techniques may include a deconvolution step which allows for further improvement of the resolution of the estimated blood flow images. However, they rely on many hyperparameters that have to be manually adjusted to obtain the optimal solution. To overcome this limitation, we propose a new deep unfolding neural network based on the DRPCA iterative algorithm, which enables the reconstruction of high-resolution and high-sensitivity blood flow components. Compared to other state-of-the-art methods, the proposed algorithm showed interesting performances, in terms of PSNR and SSIM on simulation data.

I. INTRODUCTION

In ultrafast ultrasound imaging and contrast-enhanced ultrasound imaging, the problem of reconstructing the blood flow and removing clutter signals from the raw data has received great interest from the community. For example, iterative RPCA-based techniques consist in formulating a mathematical model of the data and then finding a solution by optimizing the problem related to this model [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF], [START_REF] Michailovich | Iterative reconstruction of medical ultrasound images using spectrally constrained phase updates[END_REF]. This allows the use of regularizations, for example, blood sparsity or tissue low rank, in order to leverage prior knowledge for an efficient blood flow reconstruction. To further improve the reconstruction, a deconvolution step was recently introduced in the RPCA framework [START_REF] Shen | High-resolution and high-sensitivity blood flow estimation using optimization approaches with application to vascularization imaging[END_REF], [START_REF] Pham | Joint Blind Deconvolution and Robust Principal Component Analysis for Blood Flow Estimation in Medical Ultrasound Imaging[END_REF]. This approach called deconvolved RPCA (DPRCA) involves recovering the tissue reflectivity function from the observed data using a measured or estimated point spread function (PSF).

As mentioned above, all these methods are based on iterative optimization algorithms to find the solution such as the alternating direction method of multipliers (ADMM), Backward-Forward, Douglas-Rachford, iterative shrinkage-thresholding algorithm (ISTA), etc [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF] [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF] [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF]. However, the clinical applicability of these algorithms is limited due to the manual tuning of the hyperparameters and the number of iterations required to achieve the desired result. Alternatively, deep convolutional neural networks (CNNs) have been applied to ultrasound imaging and have shown higher performance than the optimization methods [START_REF] Van Sloun | Deep Learning in Ultrasound Imaging[END_REF]. However, these methods act as block-box engines which limits their interpretability.

To deal with all these issues, the unfolding (or unrolling) approaches are increasingly used as they provide a framework for expressing iterative algorithms as neural networks [START_REF] Yang | ADMM-CSNet: A Deep Learning Approach for Image Compressive Sensing[END_REF] - [START_REF] Juan | LADMM-Net: An unrolled deep network for spectral image fusion from compressive data[END_REF]. More precisely, each iteration step of an iterative algorithm is converted into a single layer of a network, set of which forms a deep CNN. Many unfolding algorithms have been proposed, among which the CORONA algorithm based on an unfolding of the ISTA, enables to efficiently separate the micro-bubbles blood signals from the tissue signals [START_REF] Solomon | Deep Unfolded Robust PCA With Application to Clutter Suppression in Ultrasound[END_REF].

In this work, we derive the DRPCA algorithm by interpreting the iterations of its ADMM steps as an unfolding CNN that we call as DRPCANet in the sequel. However, unlike CORONA [START_REF] Solomon | Deep Unfolded Robust PCA With Application to Clutter Suppression in Ultrasound[END_REF], the PSF of the system is explicitly included in the deep unfolding to achieve a high-resolution flow estimation. The remainder of the paper is structured as follows. After having recalled the DRPCA algorithm and its multiple-step ADMM in Section II, we show how to unfold this algorithm to build an unfolded CNN in Section III. Finally, we compare the proposed approach with the CORONA and the iterative DRPCA in Section V.

II. DRPCA ALGORITHM

A. RPCA model

The RPCA problem consists in recovering a low-rank matrix and a sparse matrix from a noisy observation [START_REF] Candes | Robust principal component analysis?[END_REF]. In our case, we consider that the observation matrix corresponds to the acquired ultrafast ultrasound (US) sequence, which we denote S. The first step consists in rearranging the 3D ultrafast US acquisition into a 2D matrix, called Casorati matrix. This matrix corresponds to a spatio-temporal representation of the acquisition, in which each column is a vectorized frame of the US sequence. Denoting B the sparse matrix, T the low-rank matrix and N the noise matrix associated with the observation, the observation model is:

S = B + T + N. (1)
In the reconstruction of blood flow problem, it is assumed that the tissue component is motionless. Then, due to spacetime invariance, in its Casorati form, all the columns of the tissue matrix are correlated. Therefore, the tissue matrix is assumed to have a low rank and thus corresponds to T . Regarding the blood, since the number of blood vessels is much lower than that of tissues, the blood is assumed to be sparse and as such, corresponds to the B matrix. Thus, the solutions to the RPCA problem correspond to the tissue and blood flow matrices [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF], [START_REF] Candes | Robust principal component analysis?[END_REF].

B. A new algorithm for the DRPCA model

To improve the resolution of the reconstructed blood flow, we incorporated a convolution in the RPCA model [START_REF] Shen | High-resolution and high-sensitivity blood flow estimation using optimization approaches with application to vascularization imaging[END_REF], [START_REF] Pham | Joint Blind Deconvolution and Robust Principal Component Analysis for Blood Flow Estimation in Medical Ultrasound Imaging[END_REF]. More precisely, we considered B = Hx where H denotes the system PSF and x the high-resolution blood flow to be estimated. For computational efficiency, H is assumed to be block circulant with circulant blocks [START_REF] Hansen | Deblurring images: matrices, spectra, and filtering[END_REF].

In [START_REF] Shen | High-resolution and high-sensitivity blood flow estimation using optimization approaches with application to vascularization imaging[END_REF], [START_REF] Pham | Joint Blind Deconvolution and Robust Principal Component Analysis for Blood Flow Estimation in Medical Ultrasound Imaging[END_REF], we proposed a DRPCA algorithm to recover the deconvolved high-resolution blood flow based on the above model. However, this algorithm requires strong conditions for its convergence. In this paper, we propose a new version of the DRPCA algorithm which enables an easier unfolding later. Before going into the details of this algorithm, we recall the general ADMM framework [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF].

Let f and g be convex, closed, and proper functions, ADMM can be used to solve the following minimization problem: minimize f (x) + g(ẑ)

subject to Ax + B ẑ = c. (2)
By minimizing the augmented Lagrangian:

L(x, ẑ, γ) = f (x) + g(ẑ) + ρ 2 Ax + ẑ -c + γ ρ 2 2 + cst . (3)
Which general solution is:

xk+1 = arg min x f (x k ) + ρ 2 Ax k + B ẑk -c + γ k ρ 2 2 , (4)
ẑk+1 = arg min ẑ g(ẑ k ) + ρ 2 Ax k + B ẑk -c + γ k ρ 2 2 , (5)
γ k+1 = γ k + (Ax k+1 + B ẑk+1 -c). (6)
To separate the high-resolution blood x and the tissue T from S, we use ADMM to solve the following problem:

[x, T ] = arg min x,T 1 2 ∥S -Hx -y∥ 2 F + λ ∥x∥ 1 + µ ∥y∥ * such that: x = Z, y = T, (7)
where Z and y are auxiliary variables. λ, µ are two positive hyperparameters balancing the trade-off between blood sparsity and tissue low rank. ∥ . ∥ 2 F is the Frobenius norm, ∥ . ∥ 1 is the l 1 norm, and ∥ . ∥ * is the nuclear norm. The nuclear norm corresponds to the l 1 norm of the vector of singular values of the considered matrix, thus it favors the appearance of zero singular values in the matrix T , and by extension, it favors a low rank. The l 1 norm favors the sparsity of x.

In order to match (2) with our problem in [START_REF] Van Luong | A Deep-Unfolded Reference-Based RPCA Network For Video Foreground-Background Separation[END_REF] we define:

M 1 = H I ; M 2 = I 0 ; M 3 = 0 I ; x = x y ; ẑ = Z T ; γ = γ 1 γ 2 ;
Our equality constraint is x -ẑ = 0 and we can choose:

f (x) = 1 2 ∥S -M 1 x∥ 2 F , g(ẑ) = λ ∥M 2 ẑ∥ 1 + µ ∥M 3 ẑ∥ * .
In this way, it is possible to separate the equations ( 4) to (6) in x and T (or their associated intermediate variables Z and y). And the solution of ( 7) is obtained by repeating the following iterative step:

y k+1 = (1 + ρ) -1 [S -Hx k + ρT k -γ k 2 ], (8)
x k+1 = (H T H + ρI) -1 [H T (S -y k ) + ρZ k -γ k 1 ], (9)
T k+1 = SV T µ ρ (y k+1 + γ k 2 /ρ), Z k+1 = Sof t λ ρ (x k+1 + γ k 1 /ρ), γ k+1 1 = γ k 1 + ρ(y k+1 -T k+1 ), γ k+1 2 = γ k 2 + ρ(x k+1 -Z k+1
), where Sof t λ µ (x) is the soft-thresholding function and SV T is the singular value thresholding function. The variable ρ is a positive penalty parameter that has an influence on the convergence speed, it is set equal to 1 in our case. In the following, we denote by P 1 , P 2 and P 3 , the three H dependent terms in ( 8) and ( 9):

P 1 = -H; P 2 = (H T H + ρI) -1 ; P 3 = (H T H + ρI) -1 H T ;
We thus obtained the iterative solution of the new DRPCA algorithm. The variables x, ẑ and γ are initialized as zero. After N iterations, the estimated solution of ( 7) is [x N , T N ].

III. UNFOLDING DRPCA

The iterative solution is unfolded into a multiple-layer neural network where each layer corresponds to an iteration of the DRPCA algorithm. To perform this conversion, we transform the matrix multiplications by P 1 , P 2 and P 3 into convolutional layers. Thus, in the unfolded algorithm, we replace P 1 , P 2 and P 3 with 2D convolutional kernels. Their weights are learned independently for each layer in a supervised manner. By doing so, the neural network learns to perform the deconvolution of the blood component without the knowledge of the PSF.

In addition, the regularization hyperparameters λ and µ are converted into learnable parameters of the network, they are also learned independently for each layer in a supervised manner. The soft-thresholding and the singular value thresholding operations are performed identically to the iterative DRPCA. The rest of the unfolding is straightforward as only elementary operations remain.

Thereby, we obtain a single layer derived from an iteration of the DRPCA algorithm. The diagram of this layer is given in Fig. 1. The overall neural network, referred to as DRPCANet, is obtained by concatenating several of these layers. During the training, we initialize the sparsity parameter λ as zero and the low-rank parameter µ equal to 5. The actual thresholding values for the soft-thresholding operations is σ(λ) = 1/(1 + exp(-λ)) which is the sigmoid of λ. We use the sigmoid function because we observe that it increases stability during training, especially when the hyperparameters are badly initialized.

IV. SIMULATION DATA

In order to make a fair comparison with CORONA [START_REF] Solomon | Deep Unfolded Robust PCA With Application to Clutter Suppression in Ultrasound[END_REF], we used their code to generate the US image sequences with slight changes to match our model. Compared to the original data, we doubled the tissue amplitude, changed the size of the PSF to be odd, and used only the real part of the data.

We can summarize the data simulation as follows. The ground truth of the blood component is represented by contrast agents (CAs) identified as bubbles. The number of contrast agents, their initial positions and their trajectories are random. The magnitude of the CAs are drawn from a normal distribution. The ground truth of the tissue is created from a mixture of Gaussians. The result is then post-processed and filtered to add a low motion. Finally, the sum of the tissue and CAs is convolved by the PSF and a Gaussian noise is added with a signal-to-noise ratio of 30dB.

We generated a total of 400 image-sequences of size 128×128 pixels and each composed of 50 frames. In Fig. 2 (a) and Fig. 3 (a), we show the raw data. In the next image we show the ground truth of the CAs (Fig. 2 (b) and Fig. 3 (b)). Fig. 2 depicts the Bmode a single frame of the US sequence while Fig. 3 shows the Power Doppler flow.

V. RESULTS

In this section, we provide reconstruction results for DR-PCANet trained and tested on the simulated data. Fig. 2 and Fig. 3 We compare our approach quantitatively using peak signalto-noise ratio (PSNR) and multiscale structural similarity (SSIM3) measurements on estimated CAs (Bubbles) in Tab. I. One can see that the unfolded DRPCA network outperforms the other algorithms by a large margin. As our approach differentiates by learning the deconvolution, we also trained CORONA on deconvolved data, but the results obtained were worse than the initial ones.

VI. CONCLUSION

Combining a model for the separation problem with a datadriven approach leads to improved separation of blood and tissue signals, together with noise reduction compared to popular optimization approaches. However, unlike CORONA, the PSF of the system is explicitly included in our deep unfolding process making the results more interpretable. Moreover, for DRPCANet, up to twenty layers have proven to be sufficient whereas the iterative DRPCA method requires a hundred iterations to achieve a similar result.

Fig. 1 .

 1 Fig. 1. Single Layer DiagramDRPCANet was implemented in Python (PyTorch) using its automatic differentiation framework for training. It consists of 20 layers. We trained the network over 400 simulated US sequences (see Section IV), 300 sequences were used for training and 100 for performance assessment. The loss function was a weighted sum of the mean squared errors between the predicted [x N , T N ] elements and their ground truth. DRPCANet was trained using the Adam optimizer for 100 epochs with a learning rate of 0.01 and a batch size fixed equal to 25. The experiments were performed on a high computation server with 32 GB RAM and 4 GTX-1080-TI GPUs.During the training, we initialize the sparsity parameter λ as zero and the low-rank parameter µ equal to 5. The actual thresholding values for the soft-thresholding operations is σ(λ) = 1/(1 + exp(-λ)) which is the sigmoid of λ. We use the sigmoid function because we observe that it increases stability during training, especially when the hyperparameters are badly initialized.

  present blood flow reconstruction results of the different methods against the ground truth image. Respectively, Fig. 2 (c) and Fig. 3 (c) are the CAs and Power Doppler obtained with the state-of-the-art CORONA neural network, Fig. 2 (d) and Fig. 3 (d) are the CAs and Power Doppler obtained with the iterative DRPCA algorithm, Fig. 2 (e) and Fig. 3 (e) are the CAs and Power Doppler obtained with DRPCANet. The CORONA network performs worse than our algorithms mainly because it does not perform deconvolution. Nevertheless, CORONA manages to separate the CAs from the tissue. As for DRPCA, both the iterative and the unfolded version achieved very clean results.

Fig. 2 .

 2 Fig. 2. Bmodes of a single frame of the estimated sequences of the blood component.

Fig. 3 .

 3 Fig. 3. Power Doppler images of the estimated blood component from the simulated data.

TABLE I PSNR

 I AND SSIM COMPARISON FOR THE ESTIMATIONS OF FIG. 3.

		CORONA [3] DRPCA [1]	DRPCAnet (Ours)
	PSNR (dB)	24.41	30.32	36.17
	SSIM3	0.13	0.71	0.94


Gregoire De Broglie

Louis Morge-Rollet

Denis Le Jeune

Frédéric Le Roy

Christian Roland

Charles Canaff

Jean-Philippe Diguet

Grégoire De Broglie

Frédéric Le Roy

New Methods for Fast Detection for Embedded Cognitive Radio

I. INTRODUCTION

A. Motivation

Cognitive Radio (CR) might be seen as improvement of the efficiency of the spectrum assignment policy as described in [START_REF] Mitola | Cognitive radio: Making software radios more personal[END_REF] and [START_REF] H.-S. Chen | Spectrum sensing for OFDM systems employing pilot tones[END_REF]. This will permit more accurate and adaptive emitter by quickly changing carrier frequency, bandwidth, power emission, and modulation used. It enables opportunistic or negotiated spectrum access from a Secondary User (SU) without interfering with the licensed Primary User (PU). This paper addresses the efficiency of PU's signal detection.

We will make use of the Orthogonal Frequency Division Multiple Access (OFDMA), which is mainly used in Wi-Fi (802.11ax), WiMAX, and notably by LTE-Advanced (LTE-A) cellular radio. The LTE-A allows us to do detection and non-interfering transmission operations: it is a long range communication standard compared to Wi-Fi and it defines a minimum resource for users, called the Resource Block (RB), to make scheduling easier. The RB is a two-dimensional combination of 12 subcarriers and 7 temporal symbols which are assigned by pairs to a User Equipment (UE). A pair of RBs is called subframe. This minimum resource gives us the pace to operate, but it must be done carefully to not interfere with control and synchronization data as specified in the LTE-A standard. The Fig. 1 depicts an example of LTE-A RB band where subframes can contain or not data. In this study, we provide detector solutions for a LTE-A CR network where data are send in unused RBs as showed in Fig. 2. The detectors can be embedded in systems with limited payload and processing resources such as Unmanned Aerial Vehicle (UAV). After prior synchronization, these detectors will detect the availability of each subframe of the LTE-A band as shown in Fig. 1. They will consist of small, energyefficient, low-power and low-latency computing platforms. Those communication platforms need to be autonomous from the "host" LTE-A system, so we designed them to only use the downlink synchronization signals from the Base Transceiver Station (BTS). This will also limit the knowledge of the system information of the BTS. In order to obtain low latency system, we chose to store only one temporal symbol in memory. We do precise filtering on RB to detect the LTE-A signal in a subframe. We limited the sensing time duration to one Orthogonal Frequency Division Multiplexing (OFDM) symbol. Sensing computations have to be as fast as possible to leave as much free symbols as possible for opportunistic communications. The complete frequency and time model will be shown on Section II-B.

B. Related Work

Spectrum sensing of the frequency band occupation during CR process can be carried out with many methods. We will first study one of the simplest solution: the energy detector, a useful reference to compare with. It has low computational and implementation complexities and can be used without any knowledge of the signal's characteristics. Reference [START_REF] Wang | Optimal Threshold of Welch's Periodogram for Sensing OFDM Signals at Low SNR Levels[END_REF] depicts an optimal threshold setting algorithm to minimize the error decision probability at low Signal-to-Noise Ratio (SNR) levels for different spectrum utilizations. Reference [START_REF] Harjula | Spectrum sensing in cognitive femto base stations using welch periodogram[END_REF] presents an energy detector which uses a Welch periodogram to avoid interferences between macro and femtocells. Reference [START_REF] Karunakaran | Sensing for Spectrum Sharing in Cognitive LTE-A Cellular Networks[END_REF] describes long and short term spectrum sharing with some algorithms to detect a signal produced by the PU at the same time as the SU. It analyzes a complete LTE-A subframe and assume that the next subframe would also be free. In those conditions, an unintended transmission done by PU would be seen as interferences and detected by the CR system. Reference [START_REF] Wang | Energy detection-based spectrum sensing with constraint region in cognitive LTE systems[END_REF] describes a spectrum sensing algorithm that works with Welch's energy detection and discrete wavelet packet transform based energy detection. The detection probability P d and the false alarm probability P f are often too weak at low SNR but, as mentioned in [START_REF] Wang | Energy detection-based spectrum sensing with constraint region in cognitive LTE systems[END_REF]: a higher P d or lower P f value can be obtained with the increase of the number of samples N or SNR, but increasing the number of samples leads to a slower system, which means that the low latency constraint would not be respected.

The matched filter detector stated in [START_REF] Omer | Review of spectrum sensing techniques in Cognitive Radio networks[END_REF] makes a crosscorrelation between an unknown signal and a known copy of the expected one. This implies knowledge of the signal's parameters like modulation, code, waveform. Unfortunately, that is impossible in LTE-A in a blind acquisition mode due to the number of combinations of parameters and the usage of scrambling.

The cyclostationary detector described by [START_REF] Ali | Analysis of OFDM parameters using Cyclostationary spectrum sensing in cognitive radio[END_REF]- [START_REF] Sohn | OFDM Signal Sensing Method Based on Cyclostationary Detection[END_REF] uses inherent periodicity to detect the presence of a signal. It is a very reliable detector which looks for periodicity of statistical properties in a cyclostationary signal. But this is inapplicable to our problem because of its heavy load of data and computing latency. Reference [START_REF] Pan | Covariance matrix based spectrum sensing for OFDM based cognitive radio[END_REF] uses a covariance matrix and the maximum eigenvalue detector but these are inoperable here because the first needs too much data and the second requires high complexity calculations which are too time expensive as we will see. Finally, the correlation detector can be a suitable solution but not as described in [START_REF] Chen | Spectrum Sensing for OFDM Systems Employing Pilot Tones and Application to DVB-T OFDM[END_REF] and [START_REF] Temtam | Using OFDM pilot tones for spectrum sensing with applications to mobile WiMAX[END_REF] because, they need several OFDM symbols which contradict the "one symbol" treatment used in methods considered in this paper.

Otherwise, detection of LTE-A and OFDMA signals has been explored in many ways. Reference [START_REF] Karunakaran | Sensing Algorithms and Protocol for Simultaneous Sensing and Reception-Based Cognitive D2D Communications in LTE-A Systems[END_REF] worked on Simultaneous Sensing and Transmission (SSR) and on timing alignment in CR. Reference [START_REF] Chaudhari | Autocorrelation-Based Decentralized Sequential Detection of OFDM Signals in Cognitive Radios[END_REF] developed a detection algorithm on whole frequency bandwidth which works with or without knowledge of Cyclic Prefix (CP) length. Refer-ence [START_REF] Van De Beek | ML estimation of time and frequency offset in OFDM systems[END_REF] is about uncertainties detection with a correlation in order to estimates the arrival time and the carrier frequency at the receiver. Reference [START_REF] Harjula | Spectrum sensing in cognitive femto base stations using welch periodogram[END_REF] studied a LTE's signal detector working with Welch periodogram. Reference [START_REF] Tu | Spectrum Sensing of OFDMA Systems for Cognitive Radio Networks[END_REF] proposed a methodology to do spectrum sensing and to classify the state of the primary system depending on the signal's type (control or traffic) and interfering or not with the primary system. Reference [START_REF] Ali | Machine Learning for Improved Resource Block Detection in 4G LTE Cognitive Radio Networks[END_REF] worked on an RB detector with Maximum Likehood Estimate (MLE) techniques for low SNR conditions but not on streaming. So far, we talked about opportunistic systems, but there is also coordinated spectrum sensing algorithm like [START_REF] Xu | Report the sensing results using OFDMA in cooperative spectrum sensing[END_REF] and [START_REF] Kokabian | Cooperative spectrum sensing in LTE networks[END_REF], which are more efficient than opportunistic systems and avoid hidden node problem or similar issues. We cannot use this type of sensing algorithm here because we want to develop an autonomous system which does not require to communicate with the BTS, only listen it for synchronization.

C. Contributions

In this paper, we consider and compare three new LTE-A RB signal's detectors: the energy detector, the correlator, and one that we will call "eogration". The correlator works with the 2 CPs of the same symbol, so it is different from [START_REF] Chen | Spectrum Sensing for OFDM Systems Employing Pilot Tones and Application to DVB-T OFDM[END_REF], [START_REF] Temtam | Using OFDM pilot tones for spectrum sensing with applications to mobile WiMAX[END_REF]. The eogration uses both energy and correlation to perform detection. The first two detectors are very classic but we propose here a new low latency implementation with a Discrete Fourier Transform (DFT) that allows RB filtering. We have to pay attention to LTE-A specifications because the standard use some carriers and symbol to transmit control and synchronization data. We also analyze which of those can be useful to create an autonomous CR system according to the constraints we fixed before.

D. Paper Organization

The rest of the paper is organized as follows. In section II, the signal's model is presented. Section III details theoretical results and simulations' detection about the three detection methods for one symbol over the whole bandwidth. We expose the adjustments and the simulations' detection for making the analyze over one RB in section IV. We do simulations of detection with a more realistic signal generator in section V. In section VI, we discuss the results obtained in the previous section and then we conclude in section VII.

II. SIGNAL MODEL

A. Signal model and propagation

The typical OFDM modulator used by an eNodeB (a LTE-A Base Station (BS)) in downlink communications is as follows. First, complex data [a 0 , • • • , a N -1 ] are used as an array of frequency data. They are processed by an Inverse Discrete Fourier Transform (IDFT) to obtain time domain data. The end of the symbol is then copied at the beginning to prevent intersymbol interference (ISI) produced during propagation. This step is called CP. In the following, we assume non-dispersive channel. And, the OFDM signal s[n] is only affected by a complex Additive White Gaussian Noise (AWGN) w[n]. We model each sample of the received signal x[n] as follows:

H 0 : x[n] = w[n] H 1 : x[n] = s[n] + w[n] (1)
Data [a 0 , • • • , a N -1 ] are assumed independent, so s[n] is a linear combination of independent and identically distributed (iid) variables. With enough samples, we can use the central limit theorem to infer that s[n] is a Gaussian process approximation with independent real and imaginary parts, as in [START_REF] Van De Beek | ML estimation of time and frequency offset in OFDM systems[END_REF]. Nevertheless, the added CP creates a correlation between some samples, so x[n] is not a white process and we can benefit from it for the detection.

B. Time and frequency constraints

We use the LTE-A downlink to transmit data simultaneously to the cellular network without disturbing it. This requires us to act in compliance with the LTE-A standard, especially with control and signalization data. We need to study the time/frequency structure of a LTE-A transmission. The largest temporal unit is the frame of 10 ms. Each frame is divided in subframe of 1 ms and each subframe in 2 RBs. Each RB consists of 7 time symbols, each approximately 70 µs long and 12 frequency bins. The frequency bins are spaced by 15 kHz, so 12 bins are equal to 180 kHz. This information is summarized in Fig. 3. As a key prerequisite, we consider that we need to be synchronized in time and in frequency to the targeted LTE-A network. This can be done with the help of Primary Synchronisation Signal (PSS) and In LTE-A, physical data can be "control" or "user" data. The first are transmitted at each subframe, i.e. every two RBs. They can occupy a variable length of at most the first 3 symbols of a subframe. The user data are transmitted on the remaining symbols of the subframe.

Under the constraint of real time processing, time is short to decode subframes' size. So, we must consider the worst possible case: control data can occupy the first 3 symbols of each subframe and user data uses the rest. Furthermore, the LTE-A standard spreads the control data over the whole band without obvious link with the corresponding user data. So, we cannot deduce the distribution of user data based on the control data. We only can be sure that the presence of user data on the 4 th symbol lead to a presence to the end of the subframe: the next 10 symbols.

These properties led us to the following setup: we ignore the first 3 symbols, we register the 4 th to detect if data are transmitted, the detection is done during 5 th , and if there is no user signal in the 4 th , we communicate from 6 th to 14 th . This is also represented in Fig. 3.

III. PRELIMINARY STUDY: DETECTION OF ONE SYMBOL

A. Introduction

Statistical hypotheses in discrete temporal domain are described in [START_REF] Mitola | Cognitive radio: Making software radios more personal[END_REF]. With x[n] the signal received and then analyzed, s[n] the OFDM signal, and w[n] the complex AWGN. Furthermore, we have the following probability densities:

w[n] ∼ CN (0, σ 2 w ) and s[n] ∼ CN (0, σ 2 s ).
Let N be the symbol length, N CP the size of the CP, and rect N (n) the rectangular non-null function on -N/2, N/2 . With these notations, we can create two signals as shown on Fig. 4:

x 1 [n] = x[n] rect N [n -N/2] x 2 [n] = x[n + N CP ] rect N [n -N/2] (2)
In hardware implementations, we consider that the entire OFDM symbol x[n] on Fig. 4 is stored before the calculation. This suppose that the detector is synchronized with the eNodeB as said in II-B. We have two hypotheses that correspond to two classes. We have also three estimators and we want to select the best, i.e. the one that separates as much as possible the two classes. We need a metric for that and we choose the Fisher's ratio [START_REF] Li | Fisher Linear Discriminant Analysis[END_REF]. [START_REF] Zarka | Separation and Concentration in Deep Networks[END_REF] used it in a similar situation, and the ratio takes account of the between-class variance and the within-class variance. The ratio is defined as follows:

F = (µ H0 -µ H1 ) 2 σ 2 H0 + σ 2 H1 (3)
with µ H0 and σ 2 H0 the expected value and variance for null hypothesis (H 0 ) and, µ H1 and σ 2 H1 the expected value and variance for alternative hypothesis (H 1 ).

B. Energy estimator

We present here the energy estimator we have developed. The canonical formula of an energy estimator is the following:

N -1 n=0 |x[n]| 2 = N -1 n=0 x[n] × x * [n].
We compute here the energy of an OFDM symbol and we do this in the frequency domain to allow the filtering of the RB. This is made possible by the Parseval's identity. The OFDM's properties are retained by doing the filtering with DFT and IDFT of the same length as the OFDM modulator in the eNodeB: N samples. So, we select for this calculus one of the two signals:

x 1 [n] or x 2 [
n], they are equivalent here because they contains the same amount of information: a complete OFDM symbol of N samples. This is why we have "x 1 [n] or x 2 [n]" at the left of the Fig. 5 which describes the working principle of the energy estimator. The statistical properties of the latter are showed in Table I.

Eng[m] = N points { points {x 1 } × DFT N points {x 1 } * } (4) Eng[m] = 1 N N -1 n=0 N -1 k=0 x 1 [k]e -2jπ kn N N -1 l=0 x 1 [l]e -2jπ ln N * e 2jπ nm N (5)
Eng

[m = 0] = N -1 k=0 |x[k]| 2 (6)

C. Correlation estimator

We also developed a correlation estimator showed in Fig. 6. As in III-B, it does the compute in frequency domain to allow RB filtering (but here with 2N points and a zero padding). The correlation classic formula is

R xy [τ ] = N -1 n=0 x[n + τ ]y * [n] in temporal domain and R xy [f ] = X[f ]Y * [f ] in frequency domain.
The output of the correlation estimator about the CP because it is the only repeated part in the complete symbol x[n] in Fig. 4. The statistics of this estimator are also in Table I.

R x2x1 [m] = IDFT 2N points { DFT 2N points {x 2 } × DFT 2N points {x 1 } * } (7) R x2x1 [m] = 1 2N 2N -1 n=0 2N -1 k=0 x 2 [k]e -2jπ kn 2N 2N -1 l=0 x 1 [l]e -2jπ ln 2N * e 2jπ nm 2N (8)
R x2x1 [m = N -N CP ] = NCP-1 l=0 x[l + N ]x * [l] (9)

D. Eogration estimator

The eogration estimator is a detector which is working in frequency domain like the others to permit RB filtering. The frequency product of the circularly rotated conjugated signal x 1 and the signal x 2 as described in Fig. 7 gives us a result divided in two parts: the energy over the data part and the correlation over the CP part. The statistics of the eogration are also in Table I.

Eog[m] = IDFT N points { DFT N points {x 2 } × DFT N points {x 1 } * } (10) Eog[m] = 1 N N -1 n=0 N -1 k=0 x 2 [k]e -2jπ kn N N -1 l=0 x 1 [l]e -2jπ ln N * e 2jπ nm N ( 11
) Table I STATISTICS OF THE ESTIMATORS Eng Rx 2 x 1 Eog H 0 E N σ 2 w 0 (N -N CP )σ 2 w Var N (σ 2 w ) 2 N CP (σ 2 w ) 2 N (σ 2 w ) 2 H 1 E N (σ 2 s + σ 2 w ) N CP σ 2 s N CP σ 2 s + (N -N CP )(σ 2 s + σ 2 w ) Var N (σ 2 s + σ 2 w ) 2 N CP (σ 2 s + σ 2 w ) 2 N (σ 2 s + σ 2 w ) 2 Fisher ratio N (σ 2 s ) 2 (σ 2 w ) 2 +(σ 2 s +σ 2 w ) 2 N CP (σ 2 s ) 2 (σ 2 w ) 2 +(σ 2 s +σ 2 w ) 2 N (σ 2 s ) 2 (σ 2 w ) 2 +(σ 2 s +σ 2 w ) 2
Eog[m = N -N CP ] = NCP-1 l=0 x[l + N ]x * [l] + N -1 l=NCP |x[l]| 2 (12)
The details of all calculations are showed in appendix A. Only the eogration is described because it is made of an energy and a correlation estimators.

E. Comparison of Fisher's ratios

Table I shows the statistics of the three estimators. The Fisher's ratios of the three estimators are equal at one exception: the multiplier. Indeed, the multiplier is N CP for the correlator and N for the others. In LTE-A, N CP ≈ 0.07 × N . So the Fisher's ratio of the correlation is too small compare to the Fisher's ratios of energy and eogration. The ratios of eogration and energy are equals because both got the same variance and eogration's expected value is within a constant (-N CP σ 2 w ), the same as energy making their difference equal. We also add that the Fisher's ratio does not take care of many characteristics of probability density like its skewness. So, identical Fisher's ratios do not necessarily imply identical performances, but it is a good indicator for the first step of a performances' comparison, before detection simulations.

F. Simulations

Signal generated by an OFDM process approximates a complex Gaussian process [START_REF] Van De Beek | ML estimation of time and frequency offset in OFDM systems[END_REF] We simulate two configurations H 0 and H 1 with 10 000 trials. We make an hypothesis test by setting the false alarm rate as 5 %. For this, we set the threshold as 95 % of the value of H 0 hypothesis, which let 5 % of error. For H 1 , we set the noise power σ 2 w to 1 and modify the signal power σ 2 s accordingly to SNR value. Finally, we plot the good detection probability P( Ĥ1 |H 1 ) in Fig. 8 with a false alarm probability of P( Ĥ1 |H 0 ) = 0.05.

G. Selection of good detectors

We saw that energy and eogration estimators got greater Fisher's ratio than the correlation due to the fact that they are proportional to N and N CP respectively with N > N CP . Moreover, we can see in Fig. 8 that the correlator is far less effective than the others. So, we chose to reject the correlation estimator. Concerning energy and eogration estimators, the performances of eogration are slightly better for the simulations than energy estimator ones and it depends on CP length, as we can see, for one RB, in Fig. 14. The performance's differences between them for the whole band are shown in Fig. 9.

IV. RB DETECTOR'S APPLICATION

A. RB filtering

We seek to detect the usage of an OFDM symbol over one RB, i.e. 12 subcarriers among all. The frequency bins are all independents, so we can subtract the useless ones by setting their values to 0 in frequency domain. The Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT) used to do this use a rectangular windows function because of nature of OFDM. Indeed, the bins of the OFDM are frequency spaced by the roots of sinc function, which makes them orthogonal in pairs.

B. Detector architectures

We present here the two schemes of the energy and eogration estimators in Fig. 10 and11. The sign of the delay is negative in the block diagram and became positive in the mathematical model, because one describes the data recording whereas the other process the data already registered, so we can use the + sign and stay deterministic. We can see the complexity gap between the estimators, the eogration needs much more computations on two parallel processes to obtain results.

We consider three computing units for the analysis of the computational complexity of the two detectors: the complex multiplier, the complex adder, and the inverter. The Table II is about computational complexity of each block of the Fig. 10 and Fig. 11. We consider the computation load for a RB detection (12 subcarriers) we give the results too for a typical value N = 1024.

C. Simulations of RB detection

Simulations of detection over a RB subband in Fig. 12 show us lower detection performance. For example, the 80 % threshold is reached at -8 dB in Fig. 8 and at 0 dB in Fig. 12. This drop is perfectly understandable because we massively reduced the number of available data.

V. SIMULATIONS WITH A LTE-A SIGNAL GENERATOR

To confirm our detection performances, we made a signal generator by complying with respect to a part of the 4G and 2048, the guardband and the CP length are specified in [START_REF]Specification # 36.104[END_REF], and modulations used here are Quadrature Phase Shift-Keying (QPSK), 16-Quadrature Amplitude Modulation (QAM) and 64-QAM. Results are shown in Figure 13. They are similar to those of the earlier section, which is an information about the reliability of the model used and the results.

VI. DISCUSSION

It can be seen that the detection performances of the eogration detector may change with the length of the CP. As shown in Fig. 14, we can see a discrepancy in performance between the energy detector and the one with eogration at SNR between -20 dB to 5 dB.

These two detectors will led us to a trade-off because the energy estimator is slightly less efficient than the eogration estimator with a one third less complex implementation. So, it is an advantage for an embedded system which does not have access to a powerful computing platform. We developed in this paper three new detectors for RB analysis of LTE-A Downlink (DL) signal. They work respectively with energy, correlation and eogration. We used the Fisher's ratio and detection statistics to show their usability. We also exhibit that the eogration estimator is slightly better than the energy estimator in Fig. 8 and Fig. 14, despite a computational complexity which is one and half times greater than the complexity of energy detection. This difference will be less important in the future because the estimators do not need costly math function and the number of computing units is low for state of the art calculators. We can therefore assume that the eogration estimator could be the default choice in the future. Most of the difficulties for the implementation comes from the real-time constraint of 30 µs on the detection. The remaining of the 70 µs will be used to prepare our own data to be emitted.

APPENDIX A DETAILED COMPUTATIONS OF THE EOGRATION

We compute here the eogration estimator described in Fig. 7 (x 1 's data are circularly rotated so, the end of the signal goes to the beginning):

Eog[m] = 1 N N -1 n=0 N -1 k=0 x 2 [k]e -2jπ kn N N -1 l=0 x 1 [l]e -2jπ ln N * e 2jπ nm N (13)
Eog[m] = N -1 k=0 N -1 l=0 x 2 [k]x * 1 [l] 1 N N -1 n=0 e -2jπ ((k-l)-m)n N ( 14
For k = l -N CP , we have: (23) Now, we want to find the associated statistical properties. We begin by the energy part of the eogration:

Eog[m = N -N CP ] = N -1 k=0 N -1 l=0 x 2 [k]x * 1 [l] δ[(k -l) -(N -N CP ) + N ] (19)
= N -1 l=0 x 2 [l -N CP ]x * 1 [l] (20)
= N -1 l=NCP x 2 [l -N CP ]x * 1 [l] (21)
N -1 l=NCP |x[l]| 2 .
We have for H 0 : x[n] ∼ CN (0, σ 2 w ) and for H 1 : x[n] ∼ CN (0, σ 2 w + σ 2 s ). We know in [START_REF] Papoulis | Probability, Random Variables and Stochastic Processes[END_REF] that |x[n]| 2 become an exponential law with a parameter λ equal to 1/σ 2 w or 1/(σ 2 s + σ 2 w ) depending on the hypothesis. So, the expected value and the variance are:

H 0    E N -1 n=0 |x[n]| 2 = N σ 2 w Var N -1 n=0 |x[n]| 2 = N (σ 2 w ) 2 (24)
H 1    E N -1 n=0 |x[n]| 2 = N (σ 2 s + σ 2 w ) Var N -1 n=0 |x[n]| 2 = N (σ 2 s + σ 2 w ) 2 (25)
The calculations for the correlation part are quite different. The variance formula with a complex random variables is the following:

Var(Z) = E |Z -E[Z]| 2 (26) Var(Z) = E |Z| 2 -|E[Z]| 2 (27)
and the formula of the variance for complex and independent random variables is: (28)

Var(XY ) =
So we obtain the following statistics:

H 0    E R x2x1 = 0 Var R x2x1 = N CP (σ 2 w ) 2 (29)
H 1    E R x2x1 = N CP σ 2 s Var R x2x1 = N CP (σ 2 s + σ 2 w ) 2 (30)
Finally, since the eogration is a sum of a correlation and an energy, the statistics are:

H 0    E Eog = (N -N CP )σ 2 w Var Eog = N (σ 2 w ) 2
(31)

H 1    E Eog = N CP σ 2 s + (N -N CP )(σ 2 s + σ 2 w ) Var Eog = N (σ 2 s + σ 2 w ) 2 (32)

Fig. 1 .Fig. 2 .

 12 Fig. 1. Example of a LTE-A group of subframes

Fig. 3 .

 3 Fig. 3. Fourth generation (4G) radio frame model used in this paper

Fig. 4 .

 4 Fig. 4. signals model

Fig. 5 .

 5 Fig. 5. Working principle of energy estimator

Fig. 6 .Fig. 7 .

 67 Fig. 6. Working principle of correlation estimator

Fig. 8 .

 8 Fig. 8. Good detection performance for a false alarm rate of 0.05

  and the noise in our model is an AWGN. The signals s[n] and w[n] are generated as Gaussian processes. Then, we create the CP for s[n]. According to the typical configuration for a 4G bandwidth of 3 MHz, the symbol length used is N = 256 and N CP = 18.

Fig. 9 .

 9 Fig. 9. Difference of performance detection between energy and eogration estimator

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Block diagram of the energy estimator

Fig. 12 .

 12 Fig. 12. Estimator comparison for a CP of 18 samples, a symbol of 256 samples, and on a RB subband

Fig. 14 .

 14 Fig. 14. Difference of performance detection between energy and eogration for multiple CP length

) 1 k=0 N - 1 l=0x 2 x 2 = NCP- 1 l=0x 2

 112212 We are interested by estimating Eog[m = N -N CP ], so the mean sum of exponentials is equal to1 if k -l = N -N CP or when k -l = -N CP . So, for k = l + (N -N CP ), we have: Eog[m = N -N CP ] = N -[k]x * 1 [l] δ[(k -l) -(N -N CP )] [l + (N -N CP )]x * 1 [l] (16) [l + (N -N CP )]x * 1 [l] + N ]x * [l]

=

  obtain the following estimator:Eog[m = N -N CP ] =

2 +

 2 Var(X) Var(Y ) + Var(X)|E[Y ]| Var(Y )|E[X]| 2

Table II TABLE

 II OF COMPUTE COMPLEXITY FOR A RB OF 12 SUBCARRIERS

	Operator name	multiplier	adder	inverter
	FFT & IFFT | | 2	N/2 log 2 N 12	N log 2 N 0	0 12
	( ) *	0	0	12
	multiplier	12	0	0
	operator
	Eogration	3N/2 log 2 (N ) + 12	3N log 2 (N )	12
	Eogration	15 372	30 720	12
	(N = 1024)
	Energy	N log 2 (N ) + 12	2N log 2 (N )	12
	Energy	10 252	20 480	12
	(N = 1024)

  RBs, 15 RBs, 25 RBs, 50 RBs, 75 RBs, and 100 RBs respectively (the last two are sampled at 30.72 MHz); the FFT length changes between 128

	Energy
	Eogration
	Fig. 13. Estimator comparison for a CP of 18 samples on a 15 RBs subband
	with a more realistic model
	standard: the frequency sampling of the generator varies
	between 1.92 MHz, 3.84 MHz, 7.68 MHz, 15.36 MHz, and
	30.72 MHz which correspond to 6