Yoko Watanabe
email: watanabe@onera.fr

Aurélien Veillard
email: aurelien.veillard@onera.fr

Caroline P Carvalho

Navigation and Guidance Strategy Planning for UAV Urban Operation

This paper proposes a concept of navigation and guidance strategy planner for urban operation of a VTOL-type UAV. One of major challenges of UAV autonomous navigation in an urban environment is to deal with the risk of GPS signal occlusion. In order to address this issue, various approaches have been proposed for GPS-independent UAV navigation and guidance such as visual odometry and visual servoing control. In this context, this work supposes that different navigation and guidance modes using different set of sensors are available onboard an UAV. An idea of the proposed planner is to anticipate the navigation and guidance performance degradation (or amelioration) due to unavailability (or availability) of certain mode in the path planning task. The planning problem is formulated as a 5D (3D position + selection of navigation and guidance modes) graph search problem, where the localization and path execution uncertainties are propagated according to a model of corresponding modes for each node transition. Node transition is denied if the path execution uncertainty ellipsoid intercects with any obstacle. A minimizing cost function is defined by a volume of the path execution uncertainty corridor, as it implies minimizing path distance and execution accuracy at the same time. A deterministic graph search algorithm is applied to find a flight path with specified navigation and guidance mode transitions which minimizes the defined cost function. Simulations are performed by using path planning configurations given in an existing UAV obstacle field navigation benchmark, and the results are presented to prove the proposed navigation and guidance strategy planning concept.

I. Introduction

In recent years, there are increasing demands of unmanned aerial vehicles (UAVs) performing a mission, either civil or military, in a remote complex environment; inspection of infrastructure such as railway or pipeline, scientific observation and mapping (archeology, agriculture, etc.), disaster relief and recovery, and reconnaissance & surveillance in an adversarial environment. In order to reduce workload and eligibility requirements for UAV operators in such mission operations, an UAV onboard flight system is required to have a maximum level of autonomy to ensure its safety in a complex environment.

One of the challenges associated with UAV autonomous operation in an urban (or even rurban) environment is to make its navigation and guidance system robust to GPS signal loss/degradation due to occlusion. Onboard flight system on most outdoor UAVs highly relies on the precision of GPS localization. The INS-only navigation solution diverges very quickly due to accumulation of inertial measurement bias, and consequently, UAV loses its automatic flight capability in the absence of GPS signal. In order to address this issue, one can find intensive research work on development of alternative navigation and guidance systems which do not rely on GPS information. For example, a variety of visual odometry and SLAM (Simultaneous Localization and Mapping) techniques have been proposed to estimate UAV state by either using purely vision or by fusing vision with INS (Inertial Navigation System) and/or other available sensors ( 1-3 and many others). Visual servoing approaches can be applied to UAV guidance and control relative to its environment without using UAV absolute state. For example, UAV terrain following, corridor navigation or landing can be achieved by using vision information directly in a flight controller 4-. 6 Such systems give an important key to realize an UAV autonomous navigation in absence of GPS signals. However, their localization or guidance and control performance is often degraded compared to that of a nominal flight system using GPS information.

In this context, the work presented in this paper originates in the idea of anticipating GPS signal occlusion risk and related flight performance degradation in UAV path planning task. Methods to predict GPS availability and accuracy at a given location can be found in literature. For example, methods proposed in 7 and 8 use multipath effect risk and a skyview angle, respectively, which are calculated from a 3D model of the surrounding area. A GPS availability map obtained through such methods can be very useful in path planning to choose a safe flight path which ensures no collision risk under the UAV localization uncertainty possibly degraded due to GPS unavailability.

Path planning under localization uncertainty is not new in the robotics community. Especially, a large amount of work has been done for navigating a ground mobile robot which possesses two different navigation modes; dead-reckoning and landmark-based localization 9-. 12 These works apply a deterministic or a sampling-based search algorithm to find a minimum-distance collision risk-free path while evolving a vehicle self-localization error ellipsoid. A choice among the two navigation modes are imposed in function of landmark visibility, and is not included as a planning object. One of the authors of this paper presented a safe path planner under localization uncertainty for UAV urban operation. [START_REF] Watanabe | Safe Path Planning with Localization Uncertainty for Urban Operation of VTOl UAV[END_REF] Unlike the work listed above, the proposed planner in [START_REF] Watanabe | Safe Path Planning with Localization Uncertainty for Urban Operation of VTOl UAV[END_REF] not only plans a flight path but also a navigation strategy (i.e., a choice of localization mode to be used on each path segment) by taking into account GPS unavailable zones and the UAV self-localization performance degradation in such zones.

This paper is an extension of this previous work by adding different guidance modes available onboad an UAV. In the previous work, the UAV localization uncertainty is considered to coincide with the UAV path execution accuracy. It is reasonable when assuming an absolute guidance law, such as waypoint tracking, which regulates the UAV localization. However in reality, it is unlikely to fly an UAV in a proximity of obstacles without any relative sensing and navigation function. Therefore, in addition to the UAV localization modes, this paper supposes different guidance modes with their availabilities and considers a closed-loop path execution uncertainty rather than the localization uncertainty. As a simple example, two guidance modes are given in this paper for a VTOL-type UAV; waypoint tracking (absolute guidance) and vision-based wall-following (relative guidance). The flight path execution uncertainty is modeled differently for each mode. An interesting point here is that the path execution uncertainty does not depend on the UAV absolute localization error when applying the wall following guidance. In other words, this relative guidance mode enables to ensure UAV flight safety regarding to collision risk even with a large degradation in the UAV self-localization accuracy due to GPS signal loss. Similarly to GPS or other sensors, availability of this wall-following guidance mode is limited and it depends on the surrounding environment (i.e. existance and proximity of a wall).

This paper first states a problem of this navigation and guidance strategy planning (Section II), then provides a propagation model of UAV self-localization and path execution uncertainties for each navigation and guidance mode (Section III, IV). Section V describes how the uncertainty propagation is incorporated in path planning task. Section ?? presents simulation results of the proposed guidance and navigation strategy planner. Section ?? includes concluding remarks and future perspectives.

II. Problem Statement

An objective of our path planner is to find the shortest collision risk-free path from a start point Xstart to a goal X goal in a 3D cluttered environment. The UAV onboard system is supposed to have several different navigation (i.e., localization) modes as well as guidance modes, which give different evolution of the localization and path execution error covariance. Therefore, the planner is also required to select the best combination of the navigation and guidance modes on each segment of the planned path. The planner assumes the following a-priori knowledge :

• Environment Map : a 3D occupancy map, each of whose voxel contains a binary value (free/occupied). It can be built from 3D point cloud obtained by vision or lidar scan data.

• Obstacle classification : a label of obstacles (for example, tree or building) for each occupied voxel in the 3D occupancy map. This information will be used to identify an availability of a certain guidance mode (e.g. wall following).

• Sensor Availability Maps : 3D voxel maps of the same size as the environment map, each of whose voxel contains a binary value (available/not). It can be constructed by considering characteristics of each sensor and its environment (e.g. GPS signal occlusion, visibility of landmark, wall sensing range).

• Localization Uncertainty Models : EKF (Extended Kalman Filter) process is applied for UAV self-localization. It includes a prediction step with INS measurements and an update step with a selected sensor measurement (specified for each navigation mode). This EKF-based model is used to propagate the localization error covariance along the planned path from its given initial value Pstart.

• Path Execution Uncertainty Models : a closed-loop translational dynamics model with a selected guidance mode is used to propagate the path execution error covariance along the planned path. Each guidance mode uses different information (UAV selflocalization result, sensor measurement, etc.), and hence transfers errors in the used information into the path execution error.

• Minimum Safety Distance : minimum separation distances from obstacles in each of the UAV body axis. It should be a sum of the UAV size and safety margin.

The collision risk is identified if separation distance between the path execution error ellipsoid and occupied voxels is less than the minimum safety distance.

III. UAV Navigation Modes

This section details the UAV self-localization systems which are used to propagate the localization error covariance from one node to another in the path planner. This paper considers the following four different means of localization. An EKF estimation process is applied to each of those self-localization systems.

A. Integration of INS Acceleration

Let X, V and a be the UAV position, velocity and acceleration in the locally-fixed inertial reference frame (North-East-Up frame). Assuming the precise attitude estimation, the INS acceleration measurement is modeled by

aINS = LBI a + ∆a + wa
where LBI is a direction cosine matrix from the reference frame to the UAV body frame. ∆a is an additive acceleration bias and wa is a measurement noise. Let

x = [ X T V T ∆a T ] T (1)
be the estimation state vector. Then its dynamics can be given by

ẋ =   Ẋ V ∆a   =   V a w∆a   = F x + GaINS + Gww (2)
where w = [ w T a w T ∆a ] T is a Gaussian process noise and

F =   O I O O O -L T BI O O O   , G =   O L T BI O   , Gw =   O O -L T BI O O I  
Given an estimate xk-1 , its error covariance P k-1 at time t k-1 and the INS acceleration input at time t k , the estimation state and its error covariance are propagated to time t k according to

x1 k = Φx k-1 + GaINS k ∆t k P 1 k = ΦP k-1 Φ T + GwQG T w ∆t k (3)
where

Φ = I + F ∆t k , ∆t k = t k -t k-1
is the state transition matrix, Q is a covariance matrix of the process noise w.

B. Sensor Measurement Updates

When a measurement of one of the sensors is available, the INS-propagated estimation x1 k is corrected through EKF update process. Let

z m k = hm(x(t k )) + ν m k
be a measurement obtained by m-th sensor at time t k . Then, the update process becomes

xm k = x1 k + K m k (z m k -hm(x 1 k )) P m k = (I -K m k H m k )P 1 k (4)
where

K m k = P 1 k H m k T (H m k P 1 k H m k T + R m k ) -1 is a Kalman gain, H m k is a Jacobian matrix of the measurement model hm(x 1
k ) and R m k is a covariance matrix of the measurement noise ν m k . The measurement models for the four different sensors listed at the beginning of the section are summarized below;

• GPS : The absolute position and velocity measurement,

h2(x) = [ X T V T ] T .
• Optical flow : Suppose that a camera is fixed on a UAV at a position X B c in the UAV body frame looking nadir. Then the camera position and velocity become

Xc = X + L T BI X B c
(5)

V c = V + L T BI (ωB × X B c )
where ωB is UAV angular velocity in the body frame. When UAV pitch and roll angles are small, the optical flow field on a flat ground surface Z = Z0 can be approximated by the following affine model.

ṗ = wc rc -rc wc p -f uc + qc vc -pc = Ap + b (6)
where p is a pixel coordinates, f is a focal length in pixel, and

ω C = LCBω B =   pc qc rc   , v C = LCI V c Z0 -Zc =   uc vc wc  
Suppose that the UAV system has an image processor which estimates the affine model of the optical flow field (A, b). [START_REF] Watanabe | Simultaneous Visual Target Tracking and Navigation in a GPS-denied Environment[END_REF] Then we can extract the measurement of

h3(x) = v C (x) = LCI V + LCB(ωB × X B c ) Z0 -Z -e T 3 L T BI X B c
where Z0 is a ground elevation provided by the environment model.

• Landmark : Let X l be a position of a landmark in the reference frame, which is known a-priori to the UAV system. If this landmark is inside the field of view of the onboard camera, a landmark detection algorithm can measure its pixel coordinates;

h4(x) = p l (x) = f Z C l X C l Y C l , X C l = LCI (X l -Xc)
where the camera position Xc given by (5).

IV. UAV Guidance Modes

This section describes the UAV guidance laws which are used to propagate the path execution error covariance from one node to another in the path planner. This paper considers the following two different guidance modes. 1) Waypoint tracking : Path following by using the UAV self-localization result (Section III). This mode is available everywhere, but its path execution accuracy directly depends on the UAV localization accuracy.

2) Wall following : Path following relative to a wall by using onboard vision-based sensor measurement but not the UAV localization result. It uses the optical flow and range measurements, which give relative velocity information to the detected wall, directly in the guidance law. This mode is available in a given sensing range from voxels labeled "wall", when the path direction is close to parallel to the wall orientation.

A. Waypoint Tracking Guidance

Consider a transition from one voxel X k-1 in the 3D occupancy map to one of its 26 neighboring voxels X k . Let u be an unit vector from X k-1 to X k , and V0 is a desired cruise speed. Then, a reference trajectory for this node transit can be defined by

X * (t) = X k-1 + V0(t -t k-1 )u (7) V * (t) = V0u (8)
for t k-1 ≤ t ≤ t k with t = t k as a time at the initial node X k-1 . t k is a time at the destination node X k , and defined by

t k = t k-1 + X2 -X1 /V0 = t k-1 + ∆t k . Let e = [ (X -X * ) T (V -V * ) T ] T
be a reference trajectory tracking error. Then its dynamics can be written by

ė = O I O O e + O I (a + w) = Ae + B(a + w)
where w is a process noise whose covariance matrix is given by Q. By assuming a small deviation in the initial condition, a simple linear guidance law can be applied to follow this reference trajectory. a

= -Kp( X(t) -X * (t)) -K d ( V (t) -V * (t)) (9)
where X and V are the UAV absolute position and velocity estimation results from Section III. Then the closed-loop trajectory tracking error becomes

ė = O I -Kp -K d e + O O O Kp K d O x + O I w = Ãe + B x + Bw ( 10
)
where x is the estimation error of the state defined in (1). By integrating this dynamic equation from t k-1 to t k , the trajectory tracking error at time t k is given by e(t k ) = e Ã∆t k e(t k-1 ) +

t k t k-1 e Ã(t k -s) B x(s)ds + t k t k-1
e Ã(t k -s) Bw(s)ds (11) Finally, the trajectory tracking error covariance, denoted by E 1 (t k ), can be calculated as follows.

E 1 (t k ) = E e(t k )e T (t k ) = e Ã∆t k E(t k-1 )e ÃT ∆t k + t k t k-1 e Ã(t k -s) BP (s) BT e ÃT (t k -s) ds + t 0 e Ã(t-s) BQB T e ÃT (t-s) ds
where E(t k-1 ) is the initial trajectory tracking error covariance, P (t) is the UAV selflocalization error covariance resulting from the EKF process in Section III. Assuming a small time interval, we can approximate an evolution of the trajectory tracking error covariance from time t k-1 to t k by

E 1 k = E k-1 + ÃE k-1 + E k-1 ÃT + BP k-1 BT + BQB T ∆t k (12)
The tracking error covariance at the destination node X k is given by this E 1 k .

B. Wall Following Guidance

The second guidance mode consists in a wall following by using optical flow and range measurements, but without using the UAV absolute localization results. Suppose a vertical wall whose orientation is ψ wall in the 3D world. We consider the same node transition and the same reference trajectory as in the previous subsection, assuming that the both nodes are within the proximity sensing range from the wall. Taking a fixed point on the wall Xp as an origin, a wall-fixed frame is defined by

X W = L3(ψ wall )(X -Xp)
where L3 represents Z-axis rotation matrix. The wall distance d is defined by d = -Y W , a negative of the second component of this X W . Figure 1 illustrates this wall frame. The reference trajectory (7-8) is transformed to this wall frame as follows.

X * W (t) = X W k-1 + V0tL3(ψ wall )u = X W k-1 + V0tu W (13) V * W (t) = V0L3(ψ wall )u = V0u W (14)
According to, 4 from the optical flow measurement, one can deduce the information v C = V C /d, the UAV velocity in the onboard camera frame scaled by the wall distance. Besides, the wall distance d can be measured by a range sensor. Assuming a perfect knowledge on the camera attitude (hence a rotation matrix from the camera frame to the wall frame LW C ), the following measurement can be used in the reference trajectory tracking guidance.

Ŷ W = Y W -Ỹ W = - d V W = V W - Ṽ W = dLW C vC
where Ṽ W and Ỹ W are the measurement errors.

Ỹ W = -d (15) Ṽ W ≃ dLW C ṽC + dLW C vC (16)
The wall following guidance law can be designed in a similar manner as the waypoint tracking guidance law;

a W = -Kp( XW (t) -X * W (t)) -K d ( V W (t) -V * W (t)) (17)
But as we do not have a position information parallel to the wall, there is no position error feedback on X W and Z W axes and it results in the linear gain in a form of

Kp =   0 0 0 0 kp 0 0 0 0  
The trajectory tracking error dynamics becomes

ėW = O I -Kp -K d e W +    0 O   0 kp 0   K d    Ỹ W Ṽ W + O I w W = Ãe W + B Ỹ W Ṽ W + Bw W (18)
From the similar procedure as in Section IV.A, the trajectory tracking error covariance in the wall frame can be propagated from t k-1 to t k by

E 2W k = E W k-1 + ÃE W k-1 + E W k-1 ÃT + BR k-1 BT + BQ W B T ∆t k (19)
R k is a measurement error covariance matrix of (15-16) given by

R k ≃ H k r d 0 T 0 ROF H T k , H k = -1 0 T LW C vC dLW C
where r d and ROF are the range and optical flow sensor error covariances respectively. This covariance matrix can be easily converted from the wall frame to the inertial frame by applying a rotation matrix L3(-ψ wall ). We recall here that, unlike a case of the waypoint tracking guidance, the uncertainty in the trajectory tracking error (19) does not depend on the UAV self-localization uncertainty.

V. Path Planning

The main idea of this paper is to apply a graph search algorithm to find a safe and short path under the path execution uncertainty E which evolves differently with a combination from the four different UAV self-localization modes presented in Section III and the two different guidance modes presented in Section IV.

A. Environment Model

In this work, a 3D voxel model is used to represent the environment as a graph. Let ∆X be a resolution and X0 be a position of the most South-West-Bottom corner node of the model. Then, a position corresponding to a (i, j, k)-voxel is

X (i,j,k) = X0 +   i -1 0 0 0 j -1 0 0 0 k -1   ∆X
An occupancy of the (i, j, k)-voxel is a-priori known. Furthermore, each occupied voxel has an obstacle label which identifies if it corresponds to a building-type obstacle where the wall following guidance mode is applicable or not. This object classification can be done by using a method proposed in [START_REF] Saux | Rapid Semantic Mapping: Learn Environment Classifiers on the Fly[END_REF] for example, and walls can be detected as contours of the building-type obstacles. The 3D voxel model is given by map(X (i,j,k) ) = 2 if occupied and wall 1 if occupied but not wall 0 if free On the graph, a neighborhood of the (i, j, k)-voxel is defined by its 26 adjacent ones; (i + δi, j + δj, k + δk) where δi, δj, δk ∈ {-1, 0, 1}. In addition to the occupancy map, an availability map for each of the four navigation modes and of the two guidance modes is supposed to be created from the 3D voxel model map(X (i,j,k) ). An availability of the m-th navigation mode at the (i, j, k)-voxel is given by navmap(X (i,j,k) , m) = 1 if the m-th navigation mode is available 0 if the m-th navigation mode is unavailable

Similarly, an availability of the l-th guidance mode at the (i, j, k)-voxel is given by guimap(X (i,j,k) , l) = 1 if the l-th guidance mode is available 0 if the l-th guidance mode is unavailable

As states in Section III and IV, the navigation and guidance mode availability condition is given by sensor characteristics used in each mode. Distance to a wall < sensing range, and Angle between path direction and wall < threshold

B. Uncertainty Propagation Model

Consider a UAV displacement from a voxel X0 to its neighbor X1. Suppose the UAV localization and path execution error covariance matrices at X0 are given by P0 and E0 respectively. If navmap(X1, m) = 1, the localization error covariance can be propagated to P m 1 through the corresponding observation model in Section III. If guimap(X1, l) = 1, the path execution error covariance can be propagated to E l 1 through the corresponding model in Section IV, which may use the UAV localization error covariance P0.

C. Posability and Traversability Conditions

Let EX > 0 be the position tracking error covariance. In this paper, an error ellipsoid is defined by the conventional 2σ-error ellipsoid of EX augmented by the safety distance d saf e . Assuming the zero mean error X = X0, the UAV is posable at the voxel X0 with the path execution error covariance E0 if map(X) = 0 for ∀ X within the augmented ellipsoid of EX0. Now we consider a traversability from a voxel X0 with the path execution error covariance E0 to a voxel X1 with the updated covariance E1. The first condition of the traversability is that the UAV must be posable at the both initial and destinate positions. Then, one should check if there is a collision between the evolving error ellipsoid and occupied voxels along the edge connecting the two. For calculation simplicity, the evolution of the error ellipsoid is approximated by a truncated elliptic cone shown in Figure 2 in this work. The UAV is traversable from X0 with E0 to X1 with E1, if map(X) = 0 for all nodes within this elliptic cone.

D. Uncertainty Corridor

In many existing work of path planning under uncertainty, a minimum distance (or time, energy) path is calculated with a traversability constraint on the error ellipsoid. Some propose a cost function as a weighed sum of the distance and a localization uncertainty cost (e.g. [START_REF] Luders | Robust Sampling-based Motion Planning with Asymptotic Optimality Guarantees[END_REF] ) to integrate the uncertainty minimization in the planning objective. However, such a cost function requires an ad-hoc parametrization to decide how much priority to be put on minimum uncertainty over minimum distance. This paper inherits a concept of uncertainty corridor proposed in. [START_REF] Watanabe | Safe Path Planning with Localization Uncertainty for Urban Operation of VTOl UAV[END_REF] The uncertainty corridor is defined by a space swept by the path execution error ellipsoid E evolving along a planned path from the start to the goal. Figure 3 shows an example of the uncertainty corridor. As done in, [START_REF] Watanabe | Safe Path Planning with Localization Uncertainty for Urban Operation of VTOl UAV[END_REF] a volume of this uncertainty corridor is chosen as a cost function in the planning task. Minimizing the volume of the uncertainty corridor implies minimizing both the path length and the accumulated path execution uncertainty. The uncertainty corridor can be approximated by an ensemble of the elliptic cones (Figure 2) between two successive voxels on the path. Therefore, the cost function is given by

J = π 6 N n=2 X (in,jn,kn) -X (i n-1 ,j n-1 ,k n-1 ) • u2 maxn u3 max n-1 + u2 max n-1 u3 maxn + 2(u2 max n-1 u3 max n-1 + u2 maxn u3 maxn ) (20)
where {X1, X2, • • • , XN } constitutes a planned path.

E. Path Search Algorithms

In order to find a collision risk-free path, which minimizes a cost function defined in (20), one can apply a classical graph search algorithm such as A*. [START_REF] Hart | A Formal Basis for the Heuristic Determination of Minimum Cost Paths[END_REF] In this work, a graph is built in a 5D search space where a node corresponds to a combination of a 3D voxel of the environment map X (i,j,k) , a navigation mode m and a guidance mode l. One node has maximum 208 (26 neighbouring voxels ×4 navigation modes ×2 guidance modes) connections. A connection to a node (X (i,j,k) , m, l) is established if the m-th navigation and l-th guidance modes are available at the voxel X (i,j,k) , and the edge is "traversable" with the corresponding evolution of the path execution uncertainty. A* search algorithm can be applied with the following cost and heuristic functions.

• cScore : a transition cost from one node to another. A volume of the approximated uncertainty corridor (the truncated elliptic cone) is used.

• hScore : a heuristic function estimating a cost-to-go at the current node to the goal. An infimum limit of the volume of the uncertainty corridor from the current node to the goal, corresponding to the case of direct connection with EX = O, is used.

• gScore : a cost to reach to the current node from the start, a summation of cScores of edges on the path.

• fScore : an estimated total cost, a sum of the gScore and hScore. This is used as a priority queue in the algorithm.

Differences from the conventional minimum-distance A* search are; i) it optimizes the cost (distance + uncertainty) over selections of a path along with the navigation and guidance modes, ii) the traversability condition under the uncertainty is imposed to confirm a connection between two nodes, iii) it outputs not only a sequence of nodes (X path ) but also a sequence of the navigation and guidance modes (M path , L path ) associated with each path segment.

F. Post Processing

A resulting path X path is represented as a sequence of adjacent voxels in the 3D occupancy map. In order to execute the planned path, it needs to be converted to a more simplified flight plan which can be given by a sequence of sparse waypoints. In order to eliminate intermediate voxels on a given path X path , Theta*-like process is performed on it. Theta* algorithm is an any-angle path planner proposed in. [START_REF] Nash | Theta*: Any-Angle Path Planning on Grids[END_REF] It adds to A* algorithm an extra process of checking the connectivity from a parent of the current node to a neighbored node before checking a direct connectivity from the current node to the neighbor. If the neighbored node can be connected to the parent, it establishes this connection by omitting the current node. The similar process is applied to omit the intermediate points on the path, by checking the traversability and the availability of navigation and guidance mode specified at each node.

VI. Simulation Results

The concept of the navigation and guidance strategy planner proposed in this paper has been validated through simulations with obstacle configurations defined in the UAV obstacle field navigation benchmark. [START_REF] Mettler | Benchmarking of Obstacle Field Navigation Algorithms for Autonomous Helicopters[END_REF] This benchmark is open to public, and obstacle models and baseline trajectory data are available online. It contains six different simple obstacle configurations, and a terrain model of San Diego downtown for real urban environment. This paper uses four of the simple obstacle configurations.

A. Simulation Settings

This subsection describes simulation and path planning settings used in our simumlation validation.

• Path Planning Configurations: Four different configurations shown in Figure 4 The minimum safety distance is set as 8m in horizontal axes and 6m in vertical. The nominal UAV speed is given by 3.5m/s. For simplification, the estimation of the accleration bias ∆a is omitted. The initial UAV self-localization error covariance is given by

Pstart = (2) 2 I (m 2 ) O O (0.4) 2 I ((m/sec) 2 )
The inital path execution error covariance Estart coincides with Pstart. It means that the map contains 101×101× 14 voxels. All obstacles are considered as "building" and their contours corresnpond to "wall" to which the wall-following guidance mode is applicable.

• Navigation and guidance mode availability map: For each voxel on the 3D occupancy map, availability of each navigation and guidance mode is calculated according to the conditions in Table 1. 90 • is used for the sky-view angle threshold for the GPS availability. The onboard camera is supposed to have 60 • field-of-view for each axis.

As the tested environment has parfectly flat ground surface, the elevation variance threshold for the optical flow-based navigation mode is set to 0. For the wall-following guidance mode, 15m wall sensing range and ±30 • angle threshold are used.

• Parameters in the navigation modes: The following covariance matrices for the process noise and the measurement noise are used in the EKF process. The EKF prediction and update process is iterated with a sampling time 0.1 sec, except the GPS update process with 0.5 sec, on each path segment.

• Parameters in the guidance modes: The linear control gain is set as kp = 0.1 and k d = 1.4 kp, which gives a damping ratio ξ = 0.7 and a natural frequency ωn = √ 0.1 in the closed-loop system. The measurement noise covariance of the optical flow is the same as that in the optical flow-based navigation mode, and that of the distance measurement is (0.1) 2 . The process noise is set to zero.

B. Results

The following four cases are tested in simulations. • TEST2: One of the solutions to navigate an UAV in a partially GPS-denied environment is to prepare an alternative GPS-free navigation system. In this second simulation set, the optical flow-based navigation system is added to the TEST1 case.

• TEST3: In the third set, instead of adding the vision-based navigation system, the vision-based wall-following guidance system is added to the TEST1 case.

• TEST4: In the last simulation set, all the four navigation modes and the two guidance modes are supposed to be available onbaord.

VII. Conclusion

This paper proposed an UAV navigation and guidance strategy planner which takes into account an evolution of the path execution uncertainty, in function of availability and choice of different UAV self-localization and guidance modes, in order to ensure the flight safety. In our previous paper, [START_REF] Watanabe | Safe Path Planning with Localization Uncertainty for Urban Operation of VTOl UAV[END_REF] we have shown that the UAV navigation strategy planner can be effective in planning a safer flight trajectory under GPS signal occlusion risk. This paper extended this work by considering path execution uncertainty, rather than localization uncertainty, with different guidance modes. While some guidance mode rely on the UAV localization result, the others do not. The concept of this navigation and guidance strategy planning is validated through simulations.

The work of this paper applied a classical graph search algorithm to find the optimal strategy and flight plan to reach a given goal. However, this approach has some limitations. For instance, the path execution uncertainty is used for the traversability check and the cost calculation. But it is not really treated as a probability distribution of the node transition (i.e. UAV displacement). That is, with such path execution uncertainty, the node transition should include more than one possibilities with corresponding probabilities for an applied guidance law. Hence, the planner needs to deal with this "uncertain" node transition associated with the path execution uncertainty. Moreover, we would like to include an "observation" process in the strategy plannig. During a mission execution, the UAV onboars system is capable to detect a current availability of each of the navigation and guidance mode (e.g. GPS signal reception, landmark or wall detection). Therefore, the planner must adapt a strategy in function of this observation result. In order to treate those two aspects, the authors are currently working on formulating the poblem in POMDP (Partially Observable Markov Decision Process) framework [START_REF] Smallwood | The Optimal Control of Partially Observable Markov Processes over a Finite Horizon[END_REF] . [START_REF] Kaelbling | Planning and Acting in Partially Observable Stochastic Domains[END_REF] POMDP handles decisional problems where the result of actions applied to the system are uncertain, and where the observations of discrete states are partial or imprecise. Finally, our future perspective includes development of an online embeddable strategy planner. As our problem dimension is larger than that of conventional flight path planners, it is a big challenge for us to develop an online planner which runs onboard an UAV during the mission execution flight in order to replan a strategy upon new environment and UAV state perception, and upon new observation of sensor availability.

Figure 1 :

 1 Figure 1: Wall Frame

Figure 2 :

 2 Figure 2: Elliptic Cone Approximation of the Error Ellipsoid Evolution between Two Voxels

  are used. Every obstacles have 20m of height above a ground surface. Start and goal points are separated by 100m, both at the altitude 10m above ground.

Figure 4 : Obstacle Configurations 19 •

 419 Figure 4: Obstacle Configurations 19

  [START_REF] Sanfourche | Minimum Uncertainty Robot Path Planning using a POMDP Approach[END_REF] I , RGP S = diag ([ (0.2) 2 (0.2) 2 (0.5) 2 (0.1) 2 (0.1) 2 (0.2) 2 ]) , ROF = diag ([ (0.02) 2 (0.02) 2 (0.04) 2 ]) , RLM = (0.03) 2 I

•

  TEST1: Most outdoor UAVs rely on the GPS/INS self-localization and use waypoints to navigate. If GPS signal is lost, they have no choice but to use the INS-only localization mode. In the first simulation set, we assume only the INS-only and GPS/INS navigation modes and the waypoint tracking guidance mode available onboard an UAV.

Figure 5 :

 5 Figure 5: TEST1: Strategy Planning Results with Nominal System

Figure 6 :

 6 Figure 6: TEST2: Strategy Planning Results with Vision-based Navigation System

Table 1 :

 1 Table 1 summarizes the availability condition for each mode. GPS availability condition is simplified with a sky-view angle threshold in this work. However, it is desirable to use a more precise GPS availability model such as those found in 7 and. 8 Availability Condition of Each Navigation and Guidance Modes

		Modes	Availability Condition
	Navigation	1 INS-only 2 GPS 3 Optical flow Elevation variance of a surface appeared on the image < threshold Everwhere Sky-view angle > threshold 4 Landmark Landmark is within the onboard camera field-of-view
	Guid.	1 Waypoint 2 Wall	Everywhere


Karine Zidane
email: karine.zidane@isae.fr

Jérôme Lacan
email: jerome.lacan@isae.fr

Mathieu Gineste
email: mathieu.gineste@thalesaleniaspace.com

Caroline Bes
email: caroline.bes@cnes.fr

Camille Bui
email: camille.hc.bui@gmail.com

Enhancement of MARSALA Random Access with Coding Schemes, Power Distributions and Maximum Ratio Combining

Several random access (RA) techniques have been proposed recently for the satellite return link. The main objective of these techniques is to resolve packets collisions in order to enhance the limited throughput of traditional RA schemes. In this context, Multi-Replica Decoding using Correlation based Localisation (MARSALA) has been introduced and has shown good performance with DVB-RCS2 coding scheme and equi-powered transmissions. However, it has been shown in the literature that alternative coding schemes and packets power distributions can have a positive impact on RA performance. Therefore, in this paper, we investigate the behaviour of MARSALA with various coding schemes and various packet power distributions, then we propose a configuration for optimal performance. This paper also introduces the enhancement of MARSALA RA scheme by adding MRC to optimize replicas combination and study the impact on the throughput. We compare two different MRC techniques and we evaluate, via simulations, the gain achieved using MRC with different coding schemes and unbalanced packets. The simulation results demonstrate that the proposed enhancements to MARSALA show substantial performance gain, i.e. throughput achieved for a target Packet Loss Ratio (PLR).

I. INTRODUCTION

The use of Random Access (RA) techniques over the satellite return link has gained a huge research interest recently. With the large spread of low-cost satellite terminals and Machine to Machine (M2M) applications, researchers have started to investigate the performance of RA techniques on the satellite return link in order to handle larger networks with low latencies. A number of RA use cases has been presented in the recent standard DVB-RCS2 [1]. In particular, the combination of RA and dedicated access is of interest for communications with long silent periods and very short packets transmissions. However, synchronous RA methods like Slotted Aloha (SA) [START_REF] Roberts | ALOHA Packet System with and without Slots and Capture[END_REF] and Diversity Slotted Aloha (DSA) [START_REF] Choudhury | Diversity ALOHA-A Random Access Scheme for Satellite Communications[END_REF] do not perform well on a satellite return link, due to destructive packets collisions. Although DSA is a variant of SA with multiple replicas per packet, it performs slightly better than SA but only for small channel loads.

In 2007, Contention Resolution Diversity Slotted Aloha (CRDSA) [START_REF] Casini | Contention Resolution Diversity Slotted ALOHA (CRDSA) : An Enhanced Random Access Scheme for Satellite Access Packet Networks[END_REF] has been introduced as the first synchronous RA technique to combine packets replication with Successive Interference Cancellation (SIC). In CRDSA and its variants [START_REF] Del Rio | A high-performance MAC protocol for Consumer Broadband Satellite Systems[END_REF], [START_REF] Gaudenzi | Advances in Random Access Protocols for Satellite Networks[END_REF], each terminal transmits multiple replicas of the same packet on different timeslots of the frame. Each replica contains pointers to the locations of the other replicas of the same packet. The receiver stores the frame then reads the timeslots in a sequential order. On each timeslot, if a packet is decoded successfully (due to the absence of collisions or the capture effect), the receiver locates the other replicas then removes the packet and its replicas from the frame. This SIC process allows to cancel the interference contribution of the successfully decoded packets and their replicas, and possibly recover additional packets. CRDSA significantly enhances the throughput and the Packet Loss Ratio (PLR) compared to SA and DSA, and it has been included in the DVB-RCS2 standard.

Irregular Repetition Slotted Aloha (IRSA) [START_REF] Liva | Graph-Based Analysis and Optimization of Contention Resolution Diversity Slotted ALOHA[END_REF] follows the same concept as CRDSA, except that the terminals can transmit an irregular number of packet replicas. IRSA enables to enhance the throughput by computing an optimal probability distribution of the number of replicas to be transmitted by each terminal. However, the performance of IRSA for a PLR lower than 10 -3 is less important than CRDSA with a fixed number of 3 or 4 replicas. Other recent RA techniques have been proposed in the literature like Coded Slotted Aloha (CSA) [START_REF] Chiani | The Marriage between Random Access and Codes on Graphs: Coded Slotted Aloha[END_REF] and Multi-Slot Coded Aloha (MuSCA) [START_REF] Bui | An Enhanced Multiple Random Access Scheme for Satellite Communications[END_REF]. CSA and MuSCA also use the SIC principle, but they are both based on packets coding and fragmentation rather than packets replication. In CSA, each packet is divided into k fragments, then the fragments are encoded with an erasure code (n,k) and the resulting n coded fragments are transmitted on randomly chosen timeslots of the frame. As for MuSCA, each packet is encoded with a robust Forward Error Correction (FEC) code of rate R, then it is divided into N b fragments. A signalling field containing the location of the other fragments is coded separately and added to each packet. A variant of MuSCA with an irregular fragmentation rate has also been proposed in [START_REF] Bui | Multi-slot coded ALOHA with Irregular Degree Distribution[END_REF]. MuSCA allows to significantly enhance the throughput compared to the existing synchronous RA methods, however it requires additional signalling overhead.

The challenge to compromise on low signalling overhead and enhanced throughput has motivated the proposition of a new decoding and demodulation technique for CRDSA called Multi Replica Decoding using Correlation Based Localisation (MARSALA) [START_REF] Bui | A Multi-Replica Decoding Technique for Contention Resolution Diversity Slotted Aloha[END_REF]. MARSALA introduces a technique for replicas localisation and decoding, in order to recover nondecoded packets in CRDSA. The transmission scheme is the same as in CRDSA, the only modifications to take into account are at the receiver side. To localise the replicas of a given packet, MARSALA performs correlation of the signal received on one timeslot with the signals received on other timeslots of the frame. When correlation peaks are detected, the replicas are combined to obtain a higher Signal to Noise plus Interference Ratio (SNIR). Moreover, a method to correct the timing offsets and phase shifts between the replicas, and ensure coherent replicas combination has been presented in [START_REF] Zidane | Estimation of Timing Offsets and Phase Shifts Between Packet Replicas in MARSALA Random Access[END_REF].

MARSALA has shown good performance in a scenario of equi-powered packets using the DVB-RCS2 turbo code for linear modulation. However, it has been shown in previous studies [START_REF] Del Rio | A high-performance MAC protocol for Consumer Broadband Satellite Systems[END_REF], [START_REF] Gaudenzi | Advances in Random Access Protocols for Satellite Networks[END_REF], [START_REF] Collard | On the Optimum Packet Power Distribution for Spread Aloha Packet Detectors With Iterative Successive Interference Cancelation[END_REF], that packets power unbalance can significantly enhance the SIC performance. In addition, it has been proven that the choice of the coding scheme has also an important impact on triggering more SIC iterations. For all these reasons, these enhancements will be studied for MARSALA in this paper. Overall, three main contributions are proposed in this paper, having the objective to enhance the performance of MARSALA RA scheme:

• Adding Maximum Ratio Combining (MRC) to MARSALA. Two techniques for applying MRC with MARSALA are proposed and evaluated via simulations.

• Evaluating MARSALA with packets power unbalance, as well as presenting an adequate packets power distribution to further enhance the performance.

• Study of MARSALA with various coding schemes.

II. SYSTEM MODEL

We consider a frame structure over one frequency carrier for the return link RA channel. The frame of length T f is divided into N s timeslots. The transmission scheme for each terminal is the same as defined in CRDSA. We denote by λ the total number of users (i.e. terminals) sharing the transmission on one frame duration. Each users sends N b replicas of the same packet on randomly chosen timeslots. We suppose that only one packet (with its replicas) per user is sent on the duration of one frame. Before transmission, each packet is encoded with a turbo code of rate R and modulated with Quadrature Phase Shift Keying (QPSK). Finally, signalling fields for the purpose of channel estimation are added.

We consider a scenario with fixed terminals communicating with a geostationary satellite. Therefore, we assume a stationary channel model with the following parameters:

(1) carrier frequency and timing offsets randomly distributed among different users, but constant among the replicas of a same user over the duration on one frame, (2) constant replicas amplitude over the frame duration and (3) random phase shifts among the replicas of a same packet. However, the phase shift is considered constant over the duration of one timeslot and phase noise fluctuations on a timeslot duration are neglected.

The processing scheme at the receiver side is illustrated in Fig. 1. First, the frame is stored, then each timeslot is processed in order to recover the received packets. In a first Fig. 1: Processing scheme for each frame at the receiver side, with combination of CRDSA and MARSALA.

step, the receiver attempts to recover the packets by using CRDSA. Each recovered packet is removed from the frame (after channel estimation) along with its replicas localised with the decoded pointers. The next step is to re-scan the frame and repeat the process of packets recovery with SIC until all the packets are decoded successfully. However, if CRDSA is blocked, i.e. the frame is scanned iteratively with no additional packets retrieved due to strong collisions, MARSALA is applied.

As described in [START_REF] Bui | A Multi-Replica Decoding Technique for Contention Resolution Diversity Slotted Aloha[END_REF], [START_REF] Zidane | Estimation of Timing Offsets and Phase Shifts Between Packet Replicas in MARSALA Random Access[END_REF], MARSALA chooses one timeslot as a reference slot T S ref and correlates its signal with the other slots of the frame. The correlation peaks are used to identify the replicas of the packets present on T S ref . Then the timing offsets and phase shifts between the replicas are corrected as detailed in [START_REF] Zidane | Estimation of Timing Offsets and Phase Shifts Between Packet Replicas in MARSALA Random Access[END_REF], and coherent replicas combination is performed. With the resulting SNIR and depending on the encoding scheme, the packet has a higher successful decoding probability. Thus, if the packet is decoded successfully thanks to this SNIR gain, it is removed from the frame. Then CRDSA scheme is applied again in the next iteration. It is worth clarifying that, although frequency, timing and phase shifts between replicas were taken into account in [START_REF] Zidane | Estimation of Timing Offsets and Phase Shifts Between Packet Replicas in MARSALA Random Access[END_REF], the fluctuating phase noise was not. According to the DVB-RCS2 guidelines [START_REF]Second Generation DVB Interactive Satellite System (DVB-RCS2)[END_REF], the phase noise is lower than -16 dBc/Hz (with dBc/Hz denoting noise power relative to the carrier contained in a 1 Hz bandwidth centered at different offsets from the carrier). Given that those fluctuations are relatively low compared to packet collisions distortions, the results presented in this work assume a negligible impact of phase noise, and open perspectives for evaluating it in future work.

III. MARSALA SCHEME WITH MRC

In previous work, MARSALA with Equal Gain Combining (EGC) has been presented. In this section, we propose to add Maximum Ratio Combining (MRC) [START_REF] Lo | Maximum Ratio Transmission[END_REF], [START_REF] Jasper | Diversity Reception Communication System with Maximum Ratio Combining Method[END_REF] to MARSALA in order to effectively use the information from all the received packet replicas. As a matter of fact, MRC has been widely used in diversity reception communication systems such as Multiple-Input Multiple-Output (MIMO) and Single-Input Multiple-Output (SIMO). Given that MARSALA is also a diversity-based transmission/reception method, we propose two MRC techniques for MARSALA: MRC based on packet SNIR knowledge and MRC based on received power per timeslot.

At this point, it is important to note that a method involving MRC of replica packets chunks has been proposed for an asynchronous RA scheme in [START_REF] Clazzer | Enhanced Contention Resolution Aloha -ECRA[END_REF], [START_REF] Clazzer | Exploiting Combination Techniques in Random Access MAC Protocols: Enhanced Contention Resolution ALOHA[END_REF]. However, the MRC technique used requires the knowledge of a symbol-by-symbol SNIR and this information is not easily retrieved especially in high channel traffic regimes.

In order to describe the MRC scheme for MARSALA, let us detail the signal received on the reference timeslot T S ref as shown below

y 1 (t) = s 1 (t)e j2π∆f1t+φ1 + n 1 (t) + i 1 (t) , (1)
where s 1 is the signal corresponding to the first replica of a given packet on T S ref , n 1 and i 1 denote respectively the noise and interference terms on T S ref . ∆f 1 , φ 1 are respectively the frequency offset and phase shift relative to the first replica. We suppose that the frequency offset is the same for all replicas, however the phase shift varies randomly from one replica to another. Supposing that the k th replicas of s ref are localised using the correlation procedures, we express the signal containing the k th replica as follows

y k (t) = s 1 (t -∆τ k )e j2π∆f1t+φ k + n k (t) + i k (t) , (2)
with φ k denoting the phase shift of the k th replica and ∆τ k representing the timing offset between the replicas. n k , i k are the noise and interference terms on the timeslot containing the k th replica. For sake of simplicity, we suppose that timing offsets and phase shifts between the replicas are perfectly compensated. Then, each signal is multiplied by a weighting factor α. The resulting combined signal is shown below:

y sum (t) = N b k=1 α k s 1 (t)e j2π∆f1t+φ1 + N b k=1 α k (n k (t) + i k (t)) . (3)
The equivalent SNIR obtained with MARSALA is

SN IR eq = N b k=1 α k 2 P s1 N b k=1 α 2 k (N k + I k ) = N b k=1 α k 2 N b k=1 α 2 k SN IR -1 k , (4)
where P s1 is the power of s 1 and N + I is the noise plus interference power term. If EGC is applied, α k = 1, otherwise α k is defined according to the MRC technique used.

A. MRC based on packet SNIR knowledge

According to (4), the optimal value of SN IR eq can be obtained with α k = SN IR k . Then SN IR eq would be equal to:

SN IR eq,max = N b k=1 SN IR k . (5)
Therefore, obtaining the maximum gain of MRC requires the knowledge of the received SN IR for each replica on each timeslot. A number of SNIR estimation techniques as well as the impact of estimation errors on MRC, can be found in the literature [START_REF] Lo | Maximum Ratio Transmission[END_REF], [START_REF] Jasper | Diversity Reception Communication System with Maximum Ratio Combining Method[END_REF]. It has been demonstrated in [START_REF] Jasper | Diversity Reception Communication System with Maximum Ratio Combining Method[END_REF] that data-aided channel estimation only slightly degrades MRC performance. Therefore, in this paper, as a first contribution to analyse the impact of MRC on MARSALA, we will consider perfect SNIR knowledge per timeslot.

B. MRC based on received power per timeslot

In case we do not have proper SNIR estimation for each replica on each timeslot, we propose to use MRC based on received power per timeslot. This technique requires the condition that all the replicas of a given packet are equipowered, however the interference packets can have different power levels. We also consider that in a frame duration, the attenuation is constant; so the useful signal received powers are equal for each set of packet replicas. The concept is the following: given that replicas of a same packet are equipowered, then we can deduce the interference level on each timeslot by measuring P k , the total received power on each timeslot. In other words, among the N b signals on the timeslots containing replicas of a given packet, the signal having the highest received power, contains the highest interference level. Therefore, once the interference level is known for each replica with the power measuring method, we can use this criteria to choose the MRC weighting factor as α k = (P k ) -1 .

C. Evaluation of MARSALA with MRC via simulations

Following the two MRC techniques described above, we evaluate the performance gain of MARSALA with MRC compared to EGC, in terms of throughput and PLR. In the simulations, CRDSA and MARSALA are combined according to the system model described in Section II. For MARSALA, the impact of imperfect replicas combination on the performance is taken into account, following the model defined in [START_REF] Zidane | Estimation of Timing Offsets and Phase Shifts Between Packet Replicas in MARSALA Random Access[END_REF]. Also, the effect of residual channel estimation errors caused by imperfect interference cancellation [START_REF] Zidane | Improved Channel Estimation for Interference Cancellation in Random Access Methods for Satellite Communications[END_REF], [START_REF] Zidane | Effect of Residual Channel Estimation Errors in Random Access Methods for Satellite Communications[END_REF] is taken into consideration. The simulations environment is provided by a satellite communications simulator developed by Thales Alenia Space and CNES. In order to only evaluate the MRC gain metric, we consider in the following that all the packets are equi-powered and QPSK modulation with DVB-RCS2 turbo code for linear modulation of rate 1/3 are used. To facilitate the recognition of several MARSALA and CRDSA versions, we denote by MARSALA-2 and CRDSA-2, the MARSALA and CRDSA systems where each terminal transmits 2 replicas of a given packet. The same notation is taken for MARSALA-3 and CRDSA-3. Fig. 2 and Fig. 3 show the normalised throughput and PLR in function of the normalised load on the frame, obtained respectively with MARSALA-2 and MARSALA-3. The normalized load (G) is expressed in bits per symbol and computed as shown below:

G = R * log 2 (M ) * λ N s , (6)
with R being the code rate and M the modulation order. The normalized throughput (T) is given by

T = G (1 -P LR(G)) , (7)
where PLR(G) is the probability that a packet is not decoded for a given G and a given SNIR.

We can observe that the performance of MARSALA with MRC based on SNIR knowledge and MARSALA with MRC based on received power per timeslot, is nearly the same. Therefore, the complexity of SNIR estimation procedures in very low SNIR regimes, can be replaced with less complex operations of received power measurement per timeslot.

Based on the results in Fig. 2 and Fig. 3, let us discuss the impact of MRC techniques on the performance of MARSALA. We can see that the performance enhancement is more visible for MARSALA-3 than MARSALA-2. On one hand we notice that the throughput of MARSALA-2 at 7 dB is increased by 20% but it is still nearly the same at 10 dB. Yet, the PLR plots do not show important enhancement for MARSALA-2. On the other hand, we can observe that MARSALA-3 throughput is increased by 20% at 3 dB and 7 dB. At 10 dB, the throughput of MARSALA-3 with MRC is enhanced by 12% and it has a maximum value of 1.2 bits/symbol. The PLR at 10 -3 is also enhanced. In addition, we can notice that the throughput obtained with MARSALA-3 and MRC at 3 dB is approximately equal to the throughput obtained at 10 dB without MRC. However, the performance in terms of PLR is less significant. In the rest of the paper, MRC will be taken into account for MARSALA in the simulations scenarios.

IV. MARSALA WITH PACKETS POWER UNBALANCE

It has been shown in [START_REF] Del Rio | A high-performance MAC protocol for Consumer Broadband Satellite Systems[END_REF], [START_REF] Gaudenzi | Advances in Random Access Protocols for Satellite Networks[END_REF] that received packets power unbalance between different users has a positive impact on the performance of CRDSA. In fact, power unbalance enables CRDSA to decode packets in collision thanks to the capture effect [START_REF] Whitehouse | Exploiting the Capture Effect for Collision Detection and Recovery[END_REF]. In other words, the strongest packets are decoded first with a higher successful decoding probability, then they are removed with SIC iterations. Thus, the weaker packets are less interfered and have a higher successful decoding probability as well.

In this section, we apply the same concept to MARSALA, given that it is also a diversity based RA method with SIC. First, in order to give a comparison with CRDSA, we study the same packets power distribution as described in [START_REF] Del Rio | A high-performance MAC protocol for Consumer Broadband Satellite Systems[END_REF], [START_REF] Gaudenzi | Advances in Random Access Protocols for Satellite Networks[END_REF]: lognormal distribution. Then, we analyse three other probability distribution functions in order to further enhance the performance.

A. Lognormal packets power distribution

Fig. 4 shows the throughput and PLR of MARSALA-2 and MARSALA-3, with lognormally distributed packets power. We suppose that all the users transmit their packets at an (E s /N 0 ) ratio equal to 10 dB. Then, at the receiver side, the received packets are attenuated following a lognormal Probability Distribution Function (PDF) of parameters µ = 0 dB and σ = 0 dB (equi-powered packets case) or σ = 2 dB.

As expected, we can observe that the performance enhancement with σ = 2 dB is significant compared to the equipowered packets case. With lognormal PDF, the throughput is increased by around 50% in MARSALA-2 and MARSALA-3.

The PLR enhancement at 10 -3 is only visible for MARSALA-3 with G = 1.6 bits/symbol.

B. Proposed packets power distributions for MARSALA

Given the important positive impact of packets power unbalance on the performance of RA methods employing SIC, recent research has proposed to apply power control techniques on the terminals and to derive an optimal PDF for the packets transmission power. In [START_REF] Collard | On the Optimum Packet Power Distribution for Spread Aloha Packet Detectors With Iterative Successive Interference Cancelation[END_REF], an optimal packets power distribution has been derived for spread Aloha packet detectors with iterative SIC. The authors have found that a PDF uniform in dB is the optimal solution for their scheme. Also, as stated in [START_REF] Del Rio | A high-performance MAC protocol for Consumer Broadband Satellite Systems[END_REF], ongoing work in the European Space Agency (ESA) is studying the optimal packets power distribution for CRDSA. As for MARSALA, an analytical study shall be done in close future. However, in this work, we consider three PDFs for packets power, as shown in Fig. 5: • The uniform distribution (in dB), with E s /N 0 varying in the interval [0 -13] dB.

• The half normal distribution (in linear scale), with parameters µ = 1 and σ = 7 (in linear scale) and E s /N 0 varying in [0-13] dB. We pick the value σ = 7 in order to have a large variance and to approach the uniform distribution, but still have lower probabilities at the rightmost edges of the PDF. • The reversed half normal distribution (in linear scale), with µ = 1 and σ = 7 (in linear scale).

Fig. 6 shows the performance of MARSALA-3 with the packets power distributions described above compared to the equipowered packets case as well as the lognormal distribution with parameters µ = 0 dB, σ = 2 dB and E s /N 0 = 10 dB. Based on the results obtained, throughput values obtained for a PLR < 10 -3 . We can notice that, among the three proposed PDFs, the packets power following the half normal distribution, presents a maximum gain (63.63%) compared to the equi-powered packets case. The advantage of using this distribution is that a high probability of users is given a low E s /N 0 ratio (between 0 dB and 4 dB).

V. PERFORMANCE COMPARISON OF MARSALA WITH VARIOUS CODING SCHEMES

Lately, MARSALA performance has been evaluated only with DVB-RCS2 turbo code for linear modulation. Yet, it has been shown in previous research [START_REF] Del Rio | A high-performance MAC protocol for Consumer Broadband Satellite Systems[END_REF] that the coding scheme behaviour in the Packet Error Rate (PER) region between 0.1 and 1, has an impact on RA schemes using SIC. For this reason, the authors in [START_REF] Del Rio | A high-performance MAC protocol for Consumer Broadband Satellite Systems[END_REF] have proven that the choice of 3GPP code for CRDSA instead of DVB-RCS2, enables to trigger more SIC iterations and achieve a better performance.

Seen that MARSALA is also a RA method that uses SIC, the choice of the coding scheme will also affect its performance. Therefore, in this section we will compare the throughput and PLR of MARSALA with three encoding schemes:

• DVB-RCS 2 turbo FEC code for linear modulation, which is a 16-states double binary Circular Recursive Systematic Convolutional (CRSC) code. In the simulations, we consider the reference waveform id-13 as defined in the DVB-RCS2 standard (payload burst length = 1476 symbols, with QPSK modulation and code rate R = 1/3).

• 3GPP turbo code [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA)[END_REF], with a burst length equal to 150 bits.

• Consultative Committee for Space Data Systems (CCSDS) [START_REF] Ccsds | TM Synchronization and channel coding: Recommended Standard[END_REF] turbo code, of rate 1/3 constructed from information block lengths of 456 bits. The CCSDS turbo code used in the simulations is provided by the Coded Modulation Library (CML) [START_REF] Valenti | Coded Modulation Library[END_REF]. Fig. 7 shows the theoretical PER curve for each of the 3 coding schemes. We can remark that for the PER region [0.1 -1], the 3GPP turbo code has the best performance. The impact of the coding schemes on the throughput and PLR of MARSALA-3 with equi-powered packets is shown in Fig. 8. As expected, MARSALA-3 with 3GPP and CCSDS achieves largely better performance than DVB-RCS2, and the best results are obtained with 3GPP. The throughput is enhanced by 54% compared to the DVB-RCS2 coding scheme and it is almost doubled compared to CRDSA. In addition, the PLR floor is significantly lower.

VI. CONCLUSION

MARSALA is a diversity based RA method employing SIC, therefore its performance is highly affected by the following design metrics: the method used for replicas combination, the packets power unbalance and the coding scheme behaviour. In this paper, we have proposed to add MRC to MARSALA. Two MRC techniques have been described and they have both shown a performance gain of around 20% in terms of throughput and PLR. We have also concluded that, similarly to CRDSA, packets power unbalance enables MARSALA to achieve a better performance. Thus, we have shown that the throughput can reach 1.8 bits/symbol with the half normal packets power distribution. In addition, we have evaluated MARSALA with various coding schemes, and the maximum throughput has been obtained with 3GPP turbo code. In conclusion, combining MARSALA with MRC, half normal packets power distribution and 3GPP turbo code can present an optimal design for an optimal performance, by maximizing the throughput and minimizing the PLR.

Fig. 2 :

 2 Fig. 2: Comparison between MARSALA with EGC, with MRC based on packet SNIR (MRC-SNIR) and with MRC based on received power per timeslot (MRC-P). N b = 2 replicas. QPSK modulation, DVB-RCS2 turbo code R = 1/3, E s /N 0 = 10 dB. Equi-Powered packets. (a) Throughput. (b) PLR.

Fig. 3 :

 3 Fig. 3: Comparison between MARSALA with and without MRC techniques. N b = 3 replicas. QPSK modulation, DVB-RCS2 turbo code R = 1/3, E s /N 0 = 10 dB. Equi-Powered packets. (a) Throughput. (b) PLR.

Fig. 4 :

 4 Fig. 4: Comparison between MARSALA-2 and MARSALA-3 with equi-powered packets and lognormally distributed packets power. QPSK modulation, DVB-RCS2 turbo code R = 1/3, E s /N 0 = 10 dB. (a) Throughput. (b) PLR.

Fig. 5 :

 5 Fig. 5: PDFs proposed for packets power in MARSALA-3.

Fig. 6 :

 6 Fig. 6: MARSALA-3 performance with various proposed packets power distributions. QPSK modulation, DVB-RCS2 turbo code R = 1/3. (a) Throughput. (b) PLR.

Fig. 7 :

 7 Fig. 7: PER vs. E s /N 0 for 3GPP, DBB-RCS 2 and CCSDS turbo codes. QPSK modulation with code rate R = 1/3.

Fig. 8 :

 8 Fig. 8: Comparison between MARSALA-3 with DVB-RCS2, 3GPP and CCSDS turbo codes. QPSK modulation, code rate R = 1/3. Equi-powered packets. (a) Throughput. (b) PLR.

TABLE I :

 I Comparison of MARSALA-3 maximum throughput (in bits/symbol) and performance gain at PLR < 10 -3 , with various packets power distributions.

  Table I summarizes the maximum

		0.16			Uniform Half Normal
		0.14		Half Normal Reversed
		0.12
		0.1
	PDF	0.08
		0.06
		0.04
		0.02
		0	0	2	4	6	8	10	12
						E s /N 0 (dB)


Student Member, IEEE Jean-Marc Belloir
email: jean-marc.belloir@isae.fr

Jean-Baptiste Lincelles

Alice Pelamatti
email: alice.pelamatti@airbus.com

Student Member, IEEE Clémentine Durnez

Member, IEEE Vincent Goiffon
email: vincent.goiffon@isae.fr

Member, IEEE Cédric Virmontois
email: cedric.virmontois@cnes.fr

Senior Member, IEEE Philippe Paillet
email: philippe.paillet@cea.fr

Member, IEEE Pierre Magnan
email: pierre.magnan@isae.fr

Member, IEEEx Olivier Gilard
email: olivier.gilard@cnes.fr

Editor A Bermak

J.-M Belloir

J.-B Lincelles

Dark Current Blooming in Pinned Photodiode CMOS Image Sensors

This paper demonstrates the existence of dark current blooming in pinned photodiode (PPD) CMOS image sensors (CISs) with the support of both experimental measurements and TCAD simulations. It is usually assumed that blooming can appear only under illumination, when the charge collected by a pixel exceeds the full well capacity (FWC) (i.e., when the photodiode becomes forward biased). In this paper, it is shown that blooming can also appear in the dark by dark current leakage from hot pixels in reverse bias (i.e., below the FWC). The dark current blooming is observed to propagate up to nine pixels away in the experimental images and can impact hundreds of pixels around each hot pixel. Hence, it can be a major image quality issue for the state-of-the-art PPD CISs used in dark current limited applications, such as low-light optical imaging and should be considered in the dark current subtraction process. This paper also demonstrates that one of the key parameter for dark current optimization, the transfer gate bias during integration, has to be carefully chosen depending on the application because the optimum bias for dark current reduction leads to the largest dark current blooming effects. Index Terms-Active pixel sensor, blooming, dark current, CMOS image sensor (CIS), diffusion current, full well capacity (FWC), pinned photodiode (PPD), p-n junction, reverse bias, thermionic emission.

I. INTRODUCTION

.

These enhancements allow them to be considered in an increasing number of low-light imaging applications. In the past years, significant efforts have been made to better understand the electrical properties and functioning of the PPD, especially regarding the pinning voltage, the PPD capacitance, the full well capacity (FWC), the transfer gate (TG) functioning, and the charge transfer efficiency [START_REF] Krymski | Estimates for scaling of pinned photodiodes[END_REF]- [START_REF] Sarkar | Temperature effects on feedforward voltage in standard CMOS pinned photodiodes[END_REF].

Nowadays, the low-light performance of the state-of-theart PPD CIS remains mainly limited by "hot" pixels, which have a higher dark current than the mean dark current of the pixel array (and thus, a lower dynamic range and sensitivity). For most of them, the higher dark current is due to the presence of one or several Shockley-Read-Hall recombinationgeneration (SRH R-G) centers [START_REF] Sah | Carrier generation and recombination in P-N junctions and P-N junction characteristics[END_REF] in the depletion region of the photodiode, which thermally generate electron-hole pairs when pnn 2 i < 0 (where p and n are the hole and electron concentrations and n i is the intrinsic free carrier concentration. This condition is verified only when the p-n junction of the PPD is reverse biased, i.e., when the charge stored in the PPD is lower than the equilibrium FWC (EFWC), which is charge stored in the PPD at equilibrium (i.e., for zero bias) [START_REF] Pelamatti | Temperature dependence and dynamic behavior of full well capacity in pinned photodiode CMOS image sensors[END_REF]. The highest reverse bias of the PPD is called the pinning voltage [START_REF] Pelamatti | Temperature dependence and dynamic behavior of full well capacity in pinned photodiode CMOS image sensors[END_REF] and is reached when the PPD is completely empty.

On the other hand, blooming effects in CIS are usually expected to appear only in forward bias, i.e., when the charge stored in the PPD is higher than the EFWC [START_REF] Pelamatti | Temperature dependence and dynamic behavior of full well capacity in pinned photodiode CMOS image sensors[END_REF]. Indeed, it is only in forward bias that the net current through a p-n junction is positive, i.e., that more electrons flow from the n side to the p side than from the p side to the n side [START_REF] Sze | Semiconductor Devices: Physics and Technology[END_REF]. In that case, there are more electrons flowing from the PPD to the epitaxy than from the epitaxy to the PPD, which means that electrons are leaking from the PPD to the epitaxy, and eventually, they are collected by the surrounding pixels. The PPD can become forward biased (collected charge higher than the EFWC) under illumination, but also with effects, such as electroluminescence or impact ionization by hot carrier injection [START_REF] Maestre | Hot carriers effects and electroluminescence in the CMOS photodiode active pixel sensors[END_REF], which act, respectively, as external sources of photons and excess minority carriers. Similar to light, these effects can overfill the impacted pixels and produce light-like blooming in the dark. However, in the dark and without these effects, the charge stored in the PPD cannot exceed the EFWC (and the photodiode cannot become forward biased), because the generation rate of SRH R-G centers cancels at equilibrium ( pnn 2 i = 0) [START_REF] Sah | Carrier generation and recombination in P-N junctions and P-N junction characteristics[END_REF]. For these reasons, it is usually assumed that hot pixels cannot produce dark current blooming (i.e., that the dark current generated by hot pixels cannot leak toward neighboring pixels).

In this paper, it is shown that dark current blooming is possible in PPD CIS contrary to what is usually believed. Indeed, it is shown that blooming can exist for small reverse biases, allowing the dark current generated in hot pixels to leak toward surrounding pixels. The dark current blooming is observed to propagate up to nine pixels away from hot pixels, potentially impacting hundreds of pixels around each hot pixel. Therefore, the dark current blooming can have major effects on low-light images and should be considered in the dark current subtraction process.

II. DEVICES UNDER TEST

Two custom four-transistor PPD (4T-PPD) CIS is tested in this paper. These devices are fabricated in a commercially available 0.18-µm technology dedicated to imaging and are made of 256 × 256 pixel arrays with different pixel pitches (4.5 and 7 µm). Fig. 1 shows the schematic cross section of a typical 4T-PPD pixel. The PPD can be seen as a spherical p-n junction formed by an n-type photodiode encapsulated in a grounded (biased at 0 V) p-type silicon. When light impinges on the pixel, the photo-generated electrons are collected by the p-n junction and stored in the n-type photodiode. At the end of the exposure time, the electrons are transferred to the floating diffusion by applying a positive TG bias (V TG = 3.3 V here), which empties the photodiode. The floating diffusion allows converting the signal from the charge domain to the voltage domain and the output signal is a voltage proportional to the collected charge. The PPD electrostatic potential (V PPD ) corresponds to the reverse bias of the p-n junction formed by the photodiode and the surrounding p-type silicon (which is grounded): it is positive in reverse bias (when the collected charge is lower than the EFWC) and negative in forward bias.

III. EXPERIMENTAL RESULTS

A. Long Integration Time Images

First of all, Fig. 2 shows frame crops acquired with both CIS in the dark and with long integration (exposure) times. Blooming effects are clearly visible around all the hot pixels of both CIS, and the magnified areas show that many pixels can be impacted around the hottest pixels. Hence, blooming effects seem possible in the dark around all regular SRH R-G hot pixels in a 4T-PPD CIS.

Fig. 3 shows the integrated charge of each pixel of the 4.5-µm pixel pitch CIS depending on the integration time in the dark. The pixels, which have a dark current (measured before saturation) higher than twice the mean dark current of the CIS, are considered as hot pixels, as opposed to the other pixels which will be called "cold" pixels. For very long integration times (>30 s), both hot and cold pixels reach a similar final saturation level, which is the EFWC. Indeed, after a very long time spent in the dark, all the photodiodes have returned to equilibrium by integrating their dark current. For integration times between 10 and 20 s, an intermediate saturation level lower than the EFWC is observed for hot pixels and persists until the cold pixels reach a similar charge level. At the intermediate level, the PPD of a hot pixel is reverse biased (which corresponds to a positive PPD potential), because the charge is lower than the EFWC. Thus, the PPD has a positive dark current I dark , which is mainly due to the generation current of SRH R-G centers for a hot pixel, and which can be simply expressed as [START_REF] Sze | Semiconductor Devices: Physics and Technology[END_REF]

I dark = I 0,dark 1 -exp - qV PPD 2kT forV ppd > 0. (1)
In ( 1), V PPD is the PPD potential and I 0,dark is the maximum dark current of the pixel reached at high reverse biases (qV PPD 2 kT). The idea that the integrated charge in the hot pixels is constant in the intermediate level suggests that the dark current (positive under reverse bias) is compensated by another current leaking out of the PPD. This hypothesis is supported by the idea that the blooming effects observed on the dark frames (Fig. 2) appear for integration times between 10 and 20 s, which corresponds to the situation where hot pixels are at the intermediate saturation level.

Fig. 4 shows the integrated charge and the dark current as a function of integration time for a 7-µm pixel pitch CIS hot pixel, as well as the mean charge and the mean dark current of its surrounding pixels depending on their distance from the hot pixel. The first neighbors correspond to pixels adjacent to the hot pixel, the second neighbors correspond to pixels adjacent to the first neighbors, and so on. The normal dark current of the CIS (mean dark current of the cold pixels) is 550 e-/s at T = 60 °C, whereas it is about 50 000 e-/s for the hot pixel studied here. As soon as the hot pixel reaches the intermediate level [about 20 kilo electrons (ke-) after 0.4 s, Fig. 4(a)], a large dark current rise is visible in the first neighbors [from 550 to more than 2000 e-/s, Fig. 4(b)] and their charge increases more rapidly than the normal charge curve of cold pixels [Fig. 4(a)]. A slight dark current rise is also visible in the second and third neighbors, which suggests that the dark current of the hot pixel is leaking in the epitaxy and can diffuse several pixels away in the epitaxy before being collected by surrounding pixels. After about 10 s, the first neighbors also reach the intermediate level [18 ke-, Fig. 4(a)] and a dark current rise is observed in the second neighbor. Successive dark current augmentations are visible in the surrounding pixels at increasing distances from the hot pixel [Fig. 4(b)], which shows that the dark current leaks from pixels that reach the intermediate level and propagates toward further pixels by diffusing in the epitaxy. The dark current blooming is observed up to nine pixels away from the hot pixel [Fig. 4(b)], at which point all the cold pixels of the CIS reach the intermediate saturation level [at about 30 s, Fig. 4(a)]. After that, the dark current blooming stops and both hot and cold pixels reach the EFWC [∼21 ke-after about 40 s, Fig. 4(a)]. In conclusion, the analysis of the long integration time images in the two PPD CIS suggests that blooming is possible in reverse bias, and that the dark current blooming from hot pixels is responsible for the observed blooming effects on the images.

B. Effect of Temperature

Fig. 5 shows the average integrated charge of all the hot pixels of the 4.5-µm pixel pitch CIS as a function of integration time in the dark and at different temperatures. The EFWC (final saturation level) increases with temperature, which is expected due to the increase of the pinning voltage [START_REF] Pelamatti | Temperature dependence and dynamic behavior of full well capacity in pinned photodiode CMOS image sensors[END_REF]. On the contrary, the intermediate saturation level reduces with temperature, which means that the PPD potential at this level increases with temperature (i.e., higher reverse bias). In this level, it has been assumed that the dark current of the PPD [generation current given by (1)] is balanced by a leakage current toward the surrounding pixels. This leakage current can be explained by an electron diffusion current within the epitaxy. Indeed, the diffusion current I diff of the p-n junction formed by the PPD and the epitaxy can be expressed as [START_REF] Sze | Semiconductor Devices: Physics and Technology[END_REF]

I diff = I 0,diff exp - qV PPD kT . ( 2
)
In ( 2), I 0,diff is the diffusion current at equilibrium. This equation shows that the diffusion current is nonnegligible at small reverse bias (it is permitted by thermionic emission over the p-n junction potential barrier [START_REF] Han | A charge transfer model for CMOS image sensors[END_REF]). Once electrons have been emitted from the hot pixel PPD to the epitaxy, they can diffuse toward surrounding pixels. In the low doped p-type epitaxy (typically 10 15 cm -3 ), the electron diffusion length [START_REF] Sze | Semiconductor Devices: Physics and Technology[END_REF] (typically 100 µm for 10 15 cm -3 ) is much larger than the distance between two PPD (which is roughly equal to the pixel pitch). Hence, the diffusion current can be approximated by the "short diode" model [START_REF] Sze | Semiconductor Devices: Physics and Technology[END_REF], which corresponds to a diffusion limited thermionic emission [START_REF] Biljanovic | Thermionic emission process in carrier transport in PN homojunctions[END_REF].

In that case, the recombination within the epitaxy is negligible and all the electrons emitted by the hot pixel are collected by the surrounding pixels (the diffusion current of the hot pixel photodiode equals the sum of the drift currents of the surrounding pixels photodiodes). Additionally, the amount of electrons emitted by the surrounding pixels can be considered negligible, because their charge is much lower than the intermediate saturation level. Hence, their reverse bias is higher and their diffusion current is negligible from (2). Thus, there are only two electron currents in the hot pixel PPD: the generation (dark) current filling the PPD and the limited thermionic emission emptying the PPD. These two currents compensate for a given PPD potential V PPDeq obtained by equating ( 1) and ( 2)

V PPDeq = 2kT q ln 1 2 4I 0,diff + I 0,dark I 0,dark + 1 2 . ( 3
)
The PPD potential at the intermediate saturation level V PPDeq increases with temperature, because I 0,diff rises more rapidly with temperature than I 0,dark . Indeed, the diffusion current (which is proportional to n 2 i ) varies more rapidly than the generation current (which is proportional to n i ), which explains why the intermediate saturation level decreases with temperature. Additionally, it can be seen from ( 3) that the PPD potential is lower if I 0,dark is higher, which explains why the intermediate saturation level is higher for the hot pixel than for cold pixels, as shown in Fig. 4. Indeed, for hot pixels, the generation current is higher and requires a larger diffusion current (thus a lower reverse bias) to be compensated.

TABLE I EXPERIMENTAL VARIATION OF THE INTERMEDIATE SATURATION LEVEL AND THE EFWC WITH TEMPERATURE

Table I presents the experimental values for the intermediate saturation level and the EFWC extracted from Fig. 5 but also for lower temperatures. Between 25 °C and 65 °C, the intermediate saturation level decreases with temperature as explained earlier. On the other hand, below 25 °C, the intermediate saturation level increases with temperature (i.e., follows the variation of the EFWC). This can be explained by the idea that the reverse bias at which the two currents compensate V PPDeq remains similar at low temperatures. For each temperature, V PPDeq can be estimated from the intermediate saturation level and the EFWC. Indeed, if the PPD capacitance is assumed not to change with PPD potential, then the PPD potential V PPD and the collected charge Q(V PPD ) are linked with the EFWC and the pinning voltage V pin by the equation

Q(V PPD ) ≈ EFWC 1 - V PPD V pin . ( 4
)
Therefore, if Q(V PPD ) is the intermediate saturation level, V PPD corresponds to V PPDeq . For the 4.5 µm-pixel pitch CIS, V pin is about 0.5 V; the estimated values for V PPDeq are reported in Table I. Below 25 °C, V PPDeq remains low (about 30 mV) and does not vary with temperature, which explains why the intermediate saturation level follows the variation of the EFWC (increases with temperature). In conclusion, the analysis of the effect of temperature supports the idea that the blooming effects observed happen in reverse bias conditions and that they correspond to dark current blooming.

IV. TCAD SIMULATIONS

A 2-D TCAD structure comprising three 7-µm pixel pitch 4T-PPD pixels similar to the pixel represented in Fig. 1 was designed to validate the existence of the dark current blooming under simulation. The doping concentration of the elementary pixel in this structure is represented in Fig. 6. First of all, a quasi-stationary simulation is performed with the middle pixel in the integration mode (TG OFF, V TG = -0.2 V) and the two surrounding pixels in the reset mode (TG ON, V TG = 3.3 V). In this state, the middle pixel (which represents a hot pixel) should reach the EFWC while the surrounding pixels (which represent cold pixels) are still empty. Fig. 7(a) shows that a generation current still exists in the PPD of the middle pixel, in agreement with the slightly positive PPD potential (about 30 mV here). In Fig. 7(b), it can be seen that an electron current is diffusing from the middle pixel toward the surrounding pixels through the epitaxy. The total diffusion current is equal to the generation current in the middle pixel in Fig. 7(a), hence, the dark current of the hot pixel is leaking to the surrounding pixels. The PPD potential of the hot pixel is slightly positive (reverse bias) and the integrated charge in the PPD (8900 e-/µm in the direction perpendicular to Figs. 6 and7) is also lower than the EFWC (9500 e-per µm); this confirms that the dark current blooming occurs in reverse bias at an intermediate saturation level. This simulation was performed at different temperatures in order to study the evolution of the intermediate saturation level and of the EFWC with temperature. As reported in Table II, the EFWC always increases with temperature as experimentally observed (Table I). On the other hand, the intermediate saturation level increases with temperature at low temperature (like experimental results) and remains similar at higher temperatures, showing that it does not follow the variation of the EFWC at high temperature. The PPD potential (reverse bias at the intermediate saturation level) increases with temperature as observed experimentally.

Fig. 8 shows the integrated charge of the middle pixel (hot pixel) and the surrounding pixels during a transient simulation. The simulation starts from a pseudoequilibrium state with all the pixels in integration mode, i.e., at EFWC (9500 e-). At 0 s, the surrounding pixels are put in reset mode and the middle pixel is kept in integration mode. Because of the dark current blooming, the charge of the middle pixel reduces from the EFWC (9500 e-and V PPD = 0 V) to the intermediate level (8900 e-and V PPD > 0 V), which corresponds to the charge of the middle pixel in the quasi-stationary simulation of Fig. 7. At 30 s, the surrounding pixels are put back in integration mode. The dark current blooming persists (the middle pixel stays at the intermediate saturation level) until the surrounding pixels also reach this level (∼50 s); then, all the pixels go back to the EFWC. In conclusion, the TCAD simulations agree with the experimental results and show that the blooming effects can be explained by the dark current blooming at small reverse bias.

V. DARK CURRENT BLOOMING COUNTERMEASURE

The effect of the TG bias (V TG ) during integration on the dark current blooming has been investigated (Fig. 9). A slight negative bias (-0.2 V in this paper) is usually chosen for low-light imaging because it ensures a very low mean dark current of the CIS. Indeed, the oxide interface under the TG is mainly accumulated with holes in that case, which neutralize the interface generation centers [START_REF] Pelamatti | Comparison of pinning voltage estimation methods in pinned photodiode CMOS image sensors[END_REF]. However, the potential barrier under the TG is higher than the potential barrier between the PPD and the epitaxy, because the substrate is grounded (0 V). Therefore, the dark current of the PPD is rather emitted toward the epitaxy than under the TG by thermionic emission, which leads to the dark current blooming (Fig. 9). If the TG is biased positively (+0.2 V in Fig. 9), which is a common light induced blooming countermeasure, the potential barrier under the TG becomes lower than toward the epitaxy and the dark current will rather diffuse under the TG (and will be collected by the floating diffusion) than toward the epitaxy. In that case, it can be seen in Fig. 9 that the dark current blooming is supressed. However, the mean dark current of the CIS is higher due to the depleted oxide interface under the TG.

VI. CONCLUSION

The existence of the dark current blooming in 4T-PPD CIS has been demonstrated based on both experimental measurements and TCAD simulations. All the results are in agreement with the analytical theory of the p-n junction and suggest that blooming is possible for small reverse bias. Hence, the dark current generated by SRH R-G hot pixels can leak toward surrounding pixels through the epitaxy. This effect, permitted by thermionic emission over the p-n junction potential barrier followed by the diffusion of electrons in the epitaxy, can produce blooming effects in hundreds of pixels around each hot pixel. Hence, it seems critical to take the dark current blooming into account when performing the dark signal subtraction in long integration time images, e.g., for low-light imaging. Increasing the TG bias during integration suppresses the dark current blooming but strongly increases the overall dark current of the CIS, hence this setting should be carefully chosen depending on the application.

Fig. 1 .

 1 Fig. 1. Schematic cross section of a 4T-PPD pixel. The readout circuit contains three transistors (not represented here) used to sample the output voltage of the floating diffusion.

Fig. 2 .

 2 Fig. 2. Dark frames (in electrons) in (a) 4.5-µm pixel pitch CIS and (b)7-µm pixel pitch CIS with exposure times of, respectively, 13 and 25 s and at temperatures of, respectively, 50 °C and 60 °C.

Fig. 3 .

 3 Fig. 3. Integrated charge as a function of integration time in the dark for all the pixels of the 4.5-µm pixel pitch CIS at 50 °C.

Fig. 4 .

 4 Fig. 4. (a) Charge and (b) dark current of a hot pixel and its surrounding pixels as a function of integration time in the dark at T = 60 °C in the 7-µm pixel pitch CIS.

Fig. 5 .

 5 Fig. 5. Average charge curve of all the hot pixels of the 4.5-µm pixel pitch CIS in the dark at different temperatures. For each curve, the integration time is normalized with respect to the transition between the intermediate saturation level and the EFWC for clarity purposes.

Fig. 6 .Fig. 7 .

 67 Fig. 6. 2-D TCAD doping structure of the 7-µm pixel pitch 4T-PPD pixel. This structure is repeated three times to form the complete simulated structure.

Fig. 8 .

 8 Fig. 8. Charge in middle and surrounding pixels during the transient simulation in the dark (T = 300 K).

Fig. 9 .

 9 Fig. 9. Dark frames of an identical region of the 7-µm pixel pitch CIS, in identical experimental conditions except the TG bias during integration: -0.2 V on the left and +0.2 V on the right.

TABLE II SIMULATED

 II VARIATION OF THE INTERMEDIATE SATURATION LEVEL AND THE EFWC WITH TEMPERATURE