Paris Oikonomou

Athanasios Dometios

Mehdi Khamassi

Costas S Tzafestas

Task Driven Skill Learning in a Soft-Robotic Arm

published or not. The documents may come

L'archive ouverte

Task Driven Skill Learning in a Soft-Robotic Arm

Paris Oikonomou1 , Athanasios Dometios 1 , Mehdi Khamassi 1,2 and Costas S. Tzafestas 1 Abstract-In this paper we introduce a novel technique that aims to dynamically control a two-module bio-inspired softrobotic arm in order to qualitatively reproduce a path defined by sparse way-points. The main idea behind this work is based on the assumption that a complex trajectory may be derived as a combination of a discrete set of parameterizable simple movements, as suggested by Movement Primitive (MP) theory. Capitalising on recent advances in this field, the proposed controller uses a Probabilistic MP (ProMP) model which initially creates an abstract mapping in the primitive-level between the task and the actuation space, and subsequently guides the movement's composition by exploiting its unique properties -conditioning and blending. At the same time, a learning-based adaptive controller updates the composition parameters by estimating the inverse kinematics of the robot, while an auxiliary process through replanning ensures that the trajectory complies with the new estimation. The learning architecture is evaluated on both a simulation model, and a real soft-robotic arm. The research findings show that the proposed methodology constitutes a novel approach that successfully manages to simplify the trajectory control task for robots of complex dynamics when high-precision is not required.

Index Terms-Robot Learning, Probabilistic Movement Primitives, Reinforcement learning, Soft Robotics

I. INTRODUCTION

Well organized health care systems are increasing the life expectancy of modern societies, according to World Health Organization's research on health and aging [START_REF] Organization | World report on ageing and health[END_REF]. Personal care (showering or bathing) is included among the first Activities of Daily Living (ADL), which incommode an elderly's life [START_REF] Millán-Calenti | Prevalence of functional disability in activities of daily living (adl), instrumental activities of daily living (iadl) and associated factors, as predictors of morbidity and mortality[END_REF]. The present work capitalises on the outcomes of the I-SUPPORT project (EU H2020 grant agreement no. 643666), which contributed to robotics research effort for daily assistance through the development of an innovative, modular robotic system integrating a soft-robotic arm (Fig. 1), that aims at supporting frail older adults to safely and independently complete various physically and cognitively demanding bathing tasks, such as properly washing their back and their lower limbs.

In this context, bathing tasks should be executed in a human-friendly way in order to increase the comfort of an elderly user. Thus, proper washing motions for each task should be learned through demonstrations by health care experts. Expert demonstration in the task space might raise some requirements for each task, in terms of execution time and motion complexity. However, decomposition into simpler primitive motions is necessary for a robotic device for technical reasons. The fusion of such primitive actions with different parameters (e.g. duration, amplitude etc.) can reproduce more delicate and human-friendly actions. Dynamic Movement Primitives (DMPs) proposed in [START_REF]Dynamic movement primitives-a framework for motor control in humans and humanoid robotics[END_REF] are nonlinear attractor systems which are mostly used for trajectory control since they provide stability and scalability properties. A more recent work described in [START_REF] Paraschos | Probabilistic movement primitives[END_REF] and [START_REF] Paraschos | Using probabilistic movement primitives in robotics[END_REF] presents the concept of Probabilistic Movement Primitives (ProMPs) which acts as a probabilistic framework with enhanced properties, such as conditioning or blending, that allows the exploitation of the properties of trajectory distributions for representing and learning primitives actions.

In addition, an interactive bathing application, which involves direct human-robot physical contact, is way more demanding in terms of safety than other assistive tasks. The execution of such tasks by a rigid robotic manipulator is considered risky, hence the research attention has shifted to soft robots [START_REF] Laschi | Soft robotics: Technologies and systems pushing the boundaries of robot abilities[END_REF] with inherent or structural compliance, which gives them the ability to actively interact with the environment and undergo large deformations. Many continuum manipulators have already been presented with tendon [START_REF] Renda | A general mechanical model for tendondriven continuum manipulators[END_REF] or pneumatic actuation [START_REF] Grzesiak | The bionic handling assistant: a success story of additive manufacturing[END_REF] or a combination of those [START_REF] Ansari | Towards the development of a soft manipulator as an assistive robot for personal care of elderly people[END_REF]. The more complex the mechanical and actuation structure of the robot is, the more sophisticated kinematic analysis and control scheme is required. Although, analytic kinematic models based on constant curvature assumption have been established [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF], powerful control strategies for continuum manipulators are still being developed [START_REF] Thuruthel | Control strategies for soft robotic manipulators: A survey[END_REF].

Recent research work [START_REF] Kapadia | A model-based sliding mode controller for extensible continuum robots[END_REF] specifically addressed the problems related to the use of continuum robots to perform dynamic control through model-based approaches. A similar approach is presented in [START_REF] Falkenhahn | Dynamic control of the bionic handling assistant[END_REF], where model-based controllers are developed based on suitable models using a combination of feedforward control and decoupled PD-controllers, applied to a pneumatically actuated manipulator. A different approach is proposed in [START_REF] Thuruthel | Learning dynamic models for open loop predictive control of soft robotic manipulators[END_REF], based on the design of openloop predictive controllers and machine learning dynamic models directly from the actuation to the task space. Based on a different set of techniques, the work presented in [START_REF] Godage | Dynamics for variable length multisection continuum arms[END_REF] applies novel spatial dynamics to variable length multisection continuum arms assuming circular arc deformation of continuum sections without torsion. A relevant approach is presented in [START_REF] Braganza | A neural network controller for continuum robots[END_REF] where the authors use a feedforward neural network component to compensate for dynamic uncertainties.

In the literature, the number of control approaches that are applied on the same soft-robot, allowing for a fair performance comparison, is limited and only two of them focus on dynamic control. In our previous work [START_REF] Oikonomou | Periodic movement learning in a soft-robotic arm[END_REF], a model-free neurodynamic control scheme is proposed in combination with Central Pattern Generators (CPGs) assigned with the task of generating cyclic motion patterns of desired features by the end-effector of the soft-robotic arm. Despite the complexity of the achieved task, this implementation is strictly limited to produce periodic movements in the task space. In this work, we focus on a more robust implementation that extends and broadens the capabilities of the approach proposed in [START_REF] Oikonomou | Periodic movement learning in a soft-robotic arm[END_REF], enabling the soft-arm to follow a more generic path.

Another control architecture evaluated on the same hardware is proposed in [START_REF] Thuruthel | Model-based reinforcement learning for closed-loop dynamic control of soft robotic manipulators[END_REF], where a closed-loop predictive controller was implemented with a model-based policy learning algorithm and trained using trajectory optimization and supervised learning. However, such control schemes require high computational resources, big amount of data, and a lot of iterations in order to successfully reproduce a single trajectory, since they focus on accuracy at every time step, while convergence during training is time-consuming. As a consequence, it is hard to generalize their performance to unseen trajectories, thus failing to ensure robustness. Besides, in most applications that use soft-robots, high accuracy in the trajectory level is not a principal requirement, while a qualitative reproduction with similar motion properties of the demonstrated trajectory is preferred, which constitutes the main focus of this work. Additionally, in [START_REF] Thuruthel | Model-based reinforcement learning for closed-loop dynamic control of soft robotic manipulators[END_REF], achieving some key motion execution properties, such as the capability to control the completion of the trajectory execution within a desired duration, also known as temporal scalability, is not evident.

In this paper, we present an architecture that aims to dynamically control a two-module bio-inspired soft-robotic arm. Our approach focuses on the implementation of a learning-based control scheme inspired by motor control theory, that has the ability to qualitatively reproduce a path defined by sparse way-points. This simplification of the trajectory control task is proved to be useful for robots of complex unmodeled dynamics. In particular, our approach combines: (a) Movement Primitive (MP) theory which provides a mathematical model capable of building a mapping in the trajectory-level between task and actuation spaces, while also producing complex movements as the result of the combination of primitives that are formed by simple demonstrations, and (b) a Reinforcement Learning (RL) algorithm assigned with the task of approaching the unknown inverse kinematics of the soft-robot in a multi-task learning manner as it is required by the nature of the task. In parallel, an extension of the ProMP's properties is designed that provides the ability to replan the route of the trajectory during the course of the execution. The efficiency of the complete methodology is evaluated in both a simulation model and a real soft-robot. The set of preliminary experimental results presented in this paper highlights the ability of the proposed framework to dynamically control a complex robotic system with unknown dynamics.

The following section introduces the reader to the problem statement before proceeding to the description of the control architecture in Section III. Experimental evaluation results are presented in Section IV, while concluding remarks together with indicative future research directions are provided in Section V.

II. PROBLEM STATEMENT

A soft-robot like the one examined in the present work (Fig. 1), is not often assigned with the task of path following, especially when high-precision is required; the mechanical properties of its design is mostly exploited in tasks where safety must be ensured, such as those involving humanrobot interaction or manipulation of fragile materials. The task examined in this paper does not differ from this class of applications, hence our focus lies on the qualitative reproduction of primitive actions defined in a subspace of a soft robot's workspace. The desired motion's properties are derived by analysis and decomposition of human expert's demonstrations which are usually conducted in the task space. In some approaches, kineasthetic teaching [START_REF] Gomez-Gonzalez | Adaptation and robust learning of probabilistic movement primitives[END_REF] of the robot by a user is employed, in order to overcome difficulties related to the demonstration's mapping to the actuation space. However, in the present application where kineasthetic teaching is impossible due to the mechanical structure of the soft-robot, generation from demonstrations should be redefined to cover the required subspace.

The mapping between the task and joint space might be provided by an inverse kinematic model, like in [START_REF] Paraschos | Using probabilistic movement primitives in robotics[END_REF] where a model is extracted based on the known geometry of a rigid manipulator. However, in cases where the design of such a mathematical model is difficult, because of the complexity or the stochasticity introduced by the robot's dynamics, different methods are recommended. In recent years, common approaches that cope with the difficulty to map the task-space to the actuation-level, involves the design of a learning-based controller that focuses on model learning [START_REF] Nguyen-Tuong | Model learning for robot control: a survey[END_REF], [START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: a survey[END_REF]. The methods proposed in these papers, either use supervised-learning techniques training a model offline without being capable of adjusting online its parameters and hence its behavior, or focus on extracting a dynamic model for the robot, which constitutes a different problem since the system usually is highly dependent on previous states.

As already stated in [START_REF] Oikonomou | Periodic movement learning in a soft-robotic arm[END_REF], the non-linearity as well as the stochasticity of a soft robot, introduced by its mechanical structure, and soft-properties do not encourage the use of a classical control scheme and the design of a fixed mathematical model that predicts its behavior, thus the design of a learning-based controller is the only way forward. At the same time, two of the main requirements that must be satisfied by the chosen class of learning algorithms are the generalizability and the online adaptability to any potential change of the robot's dynamics -which constitute a usual phenomenon in bio-inspired systems. Under these conditions, a Reinforcement Learning (RL) control scheme [START_REF] Sutton | Reinforcement Learning: An Introduction, ser. Adaptive Computation and Machine Learning series[END_REF] seems to be suitable for approximating the inverse static model under an unsupervised framework.

III. CONTROL ARCHITECTURE

A. Preliminaries

In this section, we briefly describe the simulation model, while we also present the experimental setup with the real robot where some preliminary results were derived in order to evaluate the performance of the proposed methodology.

1) Constant Curvature Approach for Simulation Model: The present simulation model is based on a piece-wise constant curvature (PCC) approximation that models the forward kinematics of Festo's "Bionic Handling Assistant" [START_REF] Rolf | Constant curvature continuum kinematics as fast approximate model for the bionic handling assistant[END_REF]. The motion of each physical module constituting the robot is driven by three radially symmetric arranged independent actuators, which change the configuration of the module after modifying their length, resulting in extension, contraction, or omnidirectional bending.

2) Experimental Setup: The performance of the proposed learning controller was evaluated experimentally on a twomodule soft manipulator (Fig. 1), which is described in [START_REF] Ansari | Towards the development of a soft manipulator as an assistive robot for personal care of elderly people[END_REF], [START_REF] Manti | Soft assistive robot for personal care of elderly people[END_REF]. This robotic module has been adapted from the soft-arm which was developed by the Biorobotics Institute (Pisa, Italy) in the frame of the EU project I-SUPPORT. Each module comprising the robot is made up of hybrid actuation, including three radially symmetric tendons driven by three motors, in combination with pneumatic chambers, whose actuation is considered to be fixed in this work. Therefore, the real robot actuation is based on six inputs at the motor control level. Regarding position feedback, a 3D magnetic tracker (3D Guidance trakSTAR Class 1 Type B by Ascension Technology Corp.) was used, whose 6-DoF electromagnetic probe is attached at the end of the manipulator, providing its position and orientation.

B. Trajectory extraction from demonstrations

One of the main ideas behind the ProMP theory introduced in [START_REF] Paraschos | Probabilistic movement primitives[END_REF] that is exploited in this work, is based on its capability to form a single primitive from demonstrations that are not identical but they accomplish similar tasks, due to its property to model sequence of states as a trajectory distribution under a probabilistic framework. When experimenting with soft-robots, this seems to be quite convenient since the stochasticity introduced by their soft-properties does not guarantee reproduction of the same result even when the same actuation is applied. In addition, the ability to reach novel targets by conditioning and blending multiple primitives could cope with the requirement to pass through some pre-determined way-points at specific timestamps during the execution of the trajectory.

In some robotic applications, the design of a dynamic control scheme based on classical optimal control theory is not feasible due to the absence of a fixed mathematical model that approximates robot's dynamics -the same phenomenon is observed in the present case. Movement Primitive (MP) theory could provide a straightforward solution that focuses on building a mapping between task and actuation spaces in a trajectory-level, in which a demonstration in the first corresponds to a sequence of actions in the motors, and viceversa. Exploiting the idea of ProMPs where a single primitive consists of multiple similar demonstrations performing the same task, a correspondence between primitives in both spaces -task and actuation -is formed. Following this strategy, a task-driven skill learning process is implemented in the sense that the composition of the requested trajectory is performed at a high-level using the properties provided by ProMP theory -conditioning, combining and blending primitives. Subsequently, the composition's parameters as extracted in the previous step are transferred unchanged in the actuation space, where the corresponding primitives are activated accordingly, resulting in the generation of a trajectory by the robot's End-Effector (EE) which is similar to the desired one.

The training process of ProMP model and more specifically the extraction of demonstrations is crucial for the proposed methodology's performance. At the beginning, a subspace within the workspace is defined as a region of interest for the requested trajectories, where the demonstrations should lie in. Subsequently, each motor is assigned with the task to start moving from a random point defined by a mean and a standard deviation, towards another end-point defined accordingly as shown in Fig. 2. It is evident that both extreme points should lie within the subspace. Randomness between demonstrations is enhanced by forcing the motors to pass through some median random points, lying between the two extremes, while smooth transition between consecutive points is guaranteed with linear interpolation. Therefore, a demonstration is derived as the coordinated activation and motion of all motors simultaneously. The same process is followed for the extraction of various demonstrations, ensuring that all of them results in movements of similar directions in the task-space. Before proceeding to training, the resulting trajectories are grouped into classes, under a similarity criterion. In the present case, it is assumed that all demonstrations generated by a distribution around the same extreme points in the actuation-space, belong to the same primitive. During the last step, the estimation of each class' parameters takes place so that each primitive, defined as trajectory distribution, represents the variations of all corresponding demonstrations.

C. Learning-based controller for inverse static model

In the frame of this work, a variation of Continuous Actor-Critic Learning Automaton (CACLA) proposed in [START_REF] Van Hasselt | Reinforcement learning in continuous action spaces[END_REF], [START_REF] Van Hasselt | Reinforcement learning in continuous state and action spaces[END_REF] has been designed, due to its ability to handle continuous variables in both state and action spaces. The goal of the present application is to execute a trajectory whose only requirement -apart from the evident that must be similar to the paths obtained by demonstrations -is that the robot's EE passes through the consecutive way-points defined by the user in the task-space. Before proceeding to the description of the CACLA's variation, it should be clarified that each way-point P cond [n] constitutes a target, for which the learning-based controller is assigned to the task to find the corresponding actuation a (n) ; n is the id of the conditioning point in set P cond . This fact implies that, during the execution of the trajectory the target is changed right after the EE passes through a way-point at the corresponding timestamp, so the actuation of the succeeding way-point is requested. In this way, the conditioning points are assigned periodically as targets, assuming that the each trajectory is executed more than once until the proximity to the way-points satisfies a convergence criterion.

To address the constant change of controller's target, a modified version of CACLA is designed, named Continuous Actor-Critic Learning Automaton for Multi-Task Learning (CACLA-MTL) which differs from the original in the sense that it keeps track of the current target's state in order to Algorithm 1: CACLA-MTL(n, n + , a (n) Compute & store error g

, P cond [n + ]) Def __init__(n, N, P cond [n]): k = 0 θ V k = RestoreCriticT heta() θ Ac k = RestoreActorT heta() k n = 0, for all n = 1, .., N g (n) kn = 0, for all n = 1, .., N a (n) = Ac k (P cond [n]) return a (n) end Def __call__(n, n + , a (n) , P cond [n + ]): k = k + 1 k n = k n + 1 //
(n)
kn (Eq. 9) Compute reward r

(n) k (Eq. 10) Compute TD-error δ k (Eq. 6) Update Critic θ V k (Eq. 4) if δ k > 0 then
Update Actor θ Ac k (Eq. 5) end // For target P cond [n + ] :: Sample a (n+) (Eq. 7) return a (n+) end recall it when the same target is assigned to the controller. At the same time, the actuation of the succeeding target is computed after updating the actor-critic model for the current one. A main requirement that must be satisfied is that the multiple parallel tasks must concern the same plant since all targets share the same actor-critic model θ V , θ Ac . The steps of the proposed CACLA-MTL written in pseudo-code is presented in Alg. 1.

Note that, CACLA-MTL presented in this paper refers and is applied only to the conditioning points of the ProMP model. Thus in the rest of the current paragraph, subscript cond is omitted from all variables. Also, for reasons of simplicity, the following applies:

V k (s (n) kn ) = V (s (n) kn , θ V k ) Ac k (s (n) kn ) = Ac(s (n) kn , θ Ac k ) (1)
where subscripts k denotes the global time-step, while k n is the time-step at which the update of the conditioning point n takes place.

Here, the state-space of the conditioning point n of the ProMP model at time k n was chosen to be composed of 6 continuous features in the task-space level; 3 of them represent the desired EE's position of the conditioning point n, while the remaining indicate the corresponding current EE's position:

s (n) kn = p desired , p current (n) kn (2)
where p = {p x , p y , p z } denotes the EE's position along the x,y,z-axis, respectively. As already implied, in this application the desired position in the task-space coincides to the way-point P cond [n].

At this point, it should be noted that the critic model is dependent of all the variables of the state-space as presented above, while the actor model uses only the first componentthe p desired , indicating the desired EE's position; we expect the proposed learning algorithm to compute directly the actuation that corresponds to the desired EE's position as a unitary movement, without depending on its current position, and/or learning a sequence of actions to reach the target.

Besides, the conditioning point n of the ProMP model in the actuation level is composed of the following 6 continuous variables, each one of which refers to an actuator of the robot:

a (n) =
where a ij for i = {1, 2} and j = {1, 2, 3} represents the actuation of motor j that belongs to module i of the robot.

To deal with the continuity of both the state and action spaces, radial basis functions with equally distributed Gaussian functions are used as linear function approximators (FA). This means that a parameter vector θ V is assigned to the critic FA, while a similar θ Ac is used for the actor FA. As in the original CACLA algorithm, here the parameters of each FA are updated based on the following equations, with the only difference being the introduction of superscript (n), which indicates that a variable refers to the conditioning point n, as shown below:

θ V k = θ V k-1 + α V δ k ∇ θ V V k-1 (s (n) kn-1 ) (4)
θ Ac k = θ Ac k-1 +α Ac a (n) -Ac k-1 (s (n) kn-1 ) ∇ θ Ac Ac k-1 (s (n) kn-1 ) (5)
with

δ k = r (n) k + γV k-1 (s (n) kn ) -V k-1 (s (n) kn-1 ) (6)
where s

(n)
kn is the state of conditioning point n at time-step k n , and V k (s (n) kn ) represents its state value function. α V is a learning rate, while γ is a discount factor. Regarding the actor, α Ac is a learning rate, while Ac k denotes the action that the actor's FA outputs at time-step k. a (n) is the action that is sampled from a Gaussian distribution with mean Ac k , and corresponds to the conditioning point n. Using this kind of exploration method, the policy is defined by:

π k s (n) kn , a (n) = 1 √ 2πσ exp   - a (n) -Ac k (s (n) kn ) 2 2σ 2    (7
) where σ denotes the standard deviation of the Gaussian exploration through a logistic function:

σ(g) = h 1 + e -c(g+b) (8)
where h, b and c are positive constants determining its maximum value, the bias on the horizontal axis, and its width, respectively. σ is also dependent on the variable g, which is the last observed error for conditioning point n in the task-space, and is defined as the corresponding Euclidean distance between the desired and the current EE's position:

g (n) kn = p (n) kn,desired -p (n) kn,current (9)
It is evident that when the error is large, implying a considerable deviation from the target, the Gaussian's width is large, thus enhancing the exploration, while in the opposite case a small error results in an action that is close to the actor's output with high probability.

The reward function at time-step k has been chosen to be a heuristic that represents how close to, or how far from the target position of the conditioning point n, the EE has been moved in relation to its previous state s (n) kn-1 , and is given as follows: r

(n) k = g (n) kn-1 -g (n) kn (10)
It can be easily seen that the instant reward could be either positive or negative, in case the EE at the conditioning point n moved towards or away from the desired position, respectively.

D. Re-planning in the ProMP-level

As stated earlier in this section, the proposed learningbased controller CACLA-MTL is assigned with the task to provide the ProMP model with an estimation of the actuation that corresponds to the desired way-point in the task-space, as the latter is defined by the user. The update of the inverse statics' estimation takes place during the execution of the trajectory by the soft-robot, implying that a new estimation of the succeeding actuation is available right after the sequence of actuations passes through a conditioning point. On the other hand, the execution of the trajectory is preceded by the weight formation of the ProMP framework which ensures that the sequence of actions taking place in the motorlevel passes through the estimated conditioning points. The last argument implies that the trajectory extracted by the ProMP model is invariant during the execution, under normal circumstances.

Re-planning in the ProMP-level facilitates the exploitation of the new estimation of the anchor points, as they extracted by the updated CACLA-MTL model, performing necessary changes online rather than after the end of the execution. In the present work, claiming that the trajectory re-plans its route, means that the power of the dominant instance of blended primitives is gradually decreased over the course of the execution, and is replaced by another instance that involves the updated estimation of the succeeding actuation. Fig. 3

depicts an illustrative example of transitions between instances of blended primitives.

At this point it should be clarified that as implied previously, the CACLA-MTL controller updates the actuation's estimation model right after the sequence of actions passes through a conditioning point in the actuation-level, indicating that the controller is updated as many times as the number of way-points during the course of a single execution. As a consequence, the number of activation functions as well as the number of instances of blended primitives that operate in parallel, coincide with the number of conditioning points.

In the rest of this section it is assumed that N groups of primitives u i are introduced, which are of identical structure, meaning that each one contains a copy of all primitives u i,j , trained using the process described in Section III-B.

The intuition behind the process of re-planning in the ProMP-level is quite similar to the one that is followed when applying the blending property in the ProMP theory as described in [START_REF] Paraschos | Probabilistic movement primitives[END_REF], [START_REF] Paraschos | Using probabilistic movement primitives in robotics[END_REF], that it could be seen as a blending of blended primitives. Indeed, a time-varying activation function b(t) ensures stable and smooth transition between the trajectories as it changes gradually from 0 to 1 and viceversa in order to activate or disable, respectively, a sequence of actions. Concretely, it only differs in the aspect that the first is computed as the product of the exponentiation of every single primitive u i,j to the power of the corresponding activation function a j (t), while the latter is the sum of the product between every instance of the blended primitives u i and the corresponding activation function b i (t), as indicated in the following equation:

v(t) = N i=1 b i (t)   M j=1 u i,j (t) aj (t)   (11)
where the actuation at time t is computed. Another small but not minor difference between the two procedures is that there is a time-interval around the critical area of each conditioning point where the activation function b(t) remains fixed in order to facilitate the evaluation of the controller's performance. The example in Fig. 3 depicts the evolution of activation function b(t) as well as the replanned trajectory of v(t).

E. Complete Methodology

The ProMP training process is followed by the execution of the proposed methodology, presented step-by-step in Alg. 2, resulting in the generation of trajectories by the robot's EE. Few elements of Alg. 2 are given below: At the beginning, the user should define the set P cond of conditioning points in the task-space and the corresponding timestamps T cond where the first takes place. Afterwards,

if t == T cond [n] then n + = n (mod N ) + 1 a (n+) cond = CACLA-MTL(n, n + , v(t), P cond [n + ]) u n+,Q[n+] = UpdatePrim a (n+) cond n = n + end wait() t = getCurrentT ime() -t start end end
each element P cond [n] is classified to the primitive Q n whose parameters -weights' mean and standard deviation -differ the least, using Mahalanobis distance as heuristic. Regarding the next steps, note that the activation functions are formed depending only on T cond , while they remain invariant as long as P cond and T cond are too. Subsequently, the algorithm runs constantly until a convergence criterion is satisfied regarding the estimation of conditioning points. The estimation update a (n+) cond of conditioning point n + is performed after the actuation passes through the precedent point a (n) cond when its evaluation takes place. Later when conditioning point a (n+) cond is available, reformation of u n+,Q[n+] 's weights is performed in order to re-plan the route of the actuation's sequence and pass through the updated estimation.

IV. EXPERIMENTAL EVALUATION

In this work, more than one algorithms are designed and proposed in order to fulfil the requirements that have been initially set. Therefore, the evaluation should focus on their capability to perform as an entity in a collaborative way rather than on the efficiency of each algorithm individually. A set of preliminary experimental results is presented demonstrating the capabilities of the proposed approach and the novelty of this work. The proposed dynamic control scheme is trained, while its performance is evaluated on both the simulation and the real robot, which are presented in Section III-A. Before proceeding to the execution of the proposed methodology, demonstrations are extracted and grouped into classes following the process described in Section III-B, capturing not only the generated paths in the task-space illustrated in Fig. 4b, but also the sequence of actions resulting in each path.

As a consequence, a dataset D f wd of EE's points with the corresponding actuation is formed that is exploited as follows: a Radial Basis Function Neural Network (RBF-NN) [START_REF] Haykin | Neural networks and learning machines[END_REF] is trained using the extracted information in order to form a model that approximates the forward kinematics of the robot. Such a model facilitates the learning process since it gives the possibility to train the proposed controller offline, providing a good initial estimate of the model before the online execution. This technique is recommended when the convergence is proved to be time-consuming due to the absence of a fixed mathematical model that approximates the robot's behavior. Note that, the same process is followed for the simulation even though the latter is executed in the software level. The evaluation of the RBF-NN on D f wd 's observations (error in simulation: [3.7+ -1.6]mm, error in real softarm: [6.3+ -3.4]mm) ensures that it approximates adequately the forward kinematics.

The training of the CACLA-MTL is then performed offline for both configurations -simulation and real robot -using the RBF-NN model. It should be noted that at the beginning this process takes place independently of the ProMP, since the controller is randomly initialized and potentially results in actuations that are out of the predefined subspace of demonstrations. Therefore, even though all algorithms should be evaluated as an entity, CACLA-MTL's performance is clearly independent of the ProMP. At every epoch, N points are selected randomly in the task-space considered as targets constituting the set P cond of conditioning points, and CACLA-MTL is executed continuously until a criterion is satisfied; all conditioning points of P cond simultaneously result in error that is under a predetermined threshold, or the maximum number of iterations is exceeded. Considering the distance between the RBFs forming the linear FA (30mm in critic, and 15mm in actor), the evaluation of CACLA-MTL on D f wd 's observations (error in simulation: [17.3+ -9.4]mm, error in real softarm: [22.2+ -12.8]mm) proves that the controller sufficiently models the inverse kinematics.

The training of each MP as described in [START_REF] Paraschos | Using probabilistic movement primitives in robotics[END_REF] using the corresponding set of demonstrations is followed by the assessment of the complete methodology. At this step, the learned models are used by the robot in order to perform some desired tasks, and are evaluated in terms of repeatability and skill transfer. As proposed in Alg. 2 the user defines a set of conditioning points P cond in the task-space along with the corresponding timestamps T cond and the algorithm is executed, resulting in the generation of the trajectory only once. As implied in Alg. 1 during the first execution the component of Gaussian exploration is removed for all P cond [n] and the pure estimation of the learned model is used without performing online correction. Figures 4a and4c depict three trajectories in the task-space produced in simula-tion and the real robot, along with the corresponding desired conditioning points. From Table I we see that the mean error in each conditioning point for consecutive execution (10 times) of three different trajectories is small, demonstrating the robustness of our algorithm. These results show that the dynamic control scheme proposed in this paper is capable of producing and generalizing in complex movements out of simple demonstrations on a soft-robotic arm.

A more illustrative presentation of the real robot's performance is given in the video accompanying the manuscript.

V. CONCLUSION AND DISCUSSION

In the frames of this work, apart from the application of the existing ProMP theory on a soft-robotic arm, two new methods are proposed: (i) a variation of CACLA for multi-task learning named CACLA-MTK and (ii) a technique for replanning in the ProMP-level. Both implementations are used as auxiliary algorithms to the ProMP in order to facilitate the trajectory generation in the task-space, providing online correction and adaptation capabilities. The results show that the proposed architecture satisfies the requirements that have been set from the beginning regarding the ability of the algorithm to qualitatively reproduce a trajectory defined by sparse way-points. In future work, we plan to improve the RL controller, accelerating the convergence in order to learn faster when using the real robot. One of the main drawbacks of the methodology is that the requested trajectories are limited to be similar with the derived demonstrations, e.g. in terms of trajectory direction/flow, which is planned to be improved. Eventually, the action-set could also be extended to include the pneumatic actuation, offering the ability to physically interact with the environment, handling external loads and applying forces.

Fig. 1 :

 1 Fig. 1: The two module Soft-Robotic manipulator developed during the i-Support project, (a) installed in a clinical environment and (b) a replicated version in the lab.

Fig. 2 :

 2 Fig. 2: Execution of a right to left demonstration motion generated with coordinated activation and motion of all motors simultaneously. Various demonstrations are generated in order to cover the desired subspace of the robot's workspace.

Fig. 3 :

 3 Fig. 3: (a) Stars (*) indicate way-points, and circles (o) the timeinterval where bi(t) remains fixed. (b) The activation of Red-Green-Blue trajectories is indicated by the corresponding bi(t) in (a), while the Cyan is the resulted trajectory generated by replanning.

Algorithm 2 :

 2 MethodologySet P cond and corresponding T cond N : number of conditioning points M : number of primitives n : ID of current conditioning point n + : ID of subsequent conditioning point u i : group of primitives i u i,j : primitive j of group iQ = ClassifyToPrimitives(P cond ) a(t) = ActivationFunA(T cond ) b(t) = ActivationFunB(T cond ) Set n = 1 a (n) cond = CACLA-MTL-init(n, N, P cond [n]) u n,Q[n] = UpdatePrim a(n) cond while Not Converged do t start = getCurrentT ime() t = 0 while Not the End of Trajectory do Compute & execute actuation v(t) (Eq. 11)

Fig. 4 :

 4 Fig. 4: Experimental results presented both for simulated and real robot. (a) The execution of three different trajectories defined by their conditioning points in time are shown for simulation environment with magenta, black and cyan, respectively. (b) The generated demonstrations on the real robot are grouped in different ProMPs marked with RGB colors. These demonstrations cover the whole subspace of interest in which the robot will operate. (c) The execution of three different trajectories, defined by their conditioning points in time, are shown for the real robot with magenta, black and cyan, respectively. The proposed approach is able to generalize in novel motions with the subspace of interest.

  For action a (n) :: Observe & store current state for point n s

	(n)
	kn

  {a 11 , a 12 , a 13 , a 21 , a 22 , a 23 }(n)

TABLE I :

 I Mean error values (in cm) for 10 consecutive executions of experimental trajectories described with their Dominant ProMP (Fig.4c) with real Soft-arm. The errors are measured in each conditioning point (CP) in the respective time

	Trajectory	Dominant CP1 CP2 CP3 CP4 CP5
	ID	ProMP
	1 (Magenta)	11232	3.1	1.9	2.5	3.3	4.3
	2 (Black)	32221	2.9	3.2	2.3	1.2	3.9
	3 (Cyan)	31133	2.0	2.4	2.6	1.0	4.0

All authors are with the School of Electrical and Computer Engineering, National Technical University of Athens, Greece.

Mehdi Khamassi is also with Sorbonne Université, CNRS, Institute of Intelligent Systems and Robotics, Paris, France.Email: oikonpar@mail.ntua.gr

ACKNOWLEDGMENT

This research has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 101017054 (project: SoftGrip). This work has also been partially supported by the French Centre National de la Recherche Scientifique (CNRS) PICS international scheme no. 279521.


Haiyan Liu

Huaming Guo
email: hmguo@cugb.edu.cn

Olivier Pourret

Zhen Wang

Zhanxue Sun

Weimin Zhang

Maohan Liu

Distribution of rare earth elements in sediments of the North China Plain: A probe of sedimentation process

Keywords: Lanthanides, Geochemistry, Sedimentation, Fractionation, Critical Zone

Rare earth elements (REE) are powerful tracers for our understanding of sedimentary provenance and depositional processes. This study investigated geochemical characteristics of REE in 226 sediment samples collected from piedmont, central and littoral plains of the North China Plain. Results showed that total REE concentrations exhibited a decreasing trend as piedmont > central > littoral, and 82%, 70% and 67% samples had total REE concentrations higher than the value of Upper Continental Crust. The REE concentrations were controlled by sedimentary provenance and sediment physicochemical properties. Ternary diagrams of sandstone classification and Al 2 O 3 -CaO+Na 2 O-K 2 O reflected that sediments originated from greywackes and they experienced an incipient to intermediate chemical weathering. Average Chemical Index of Alteration values were 68, 51 and 55 for sediments from piedmont, central and littoral sediments, respectively. Sediment North American Shale Composite (NASC)-normalized patterns were characterized by light REE (LREE) enrichment over heavy REE (HREE) and coherent negative Eu anomalies

(Eu/Eu* ranging between 0.57 and 1.00). Value of (La/Yb) NASC ranged from 0.86 to 1.15 (average 1.23), and from 1.03 to 1.49, and from 0.91 to 1.49 in three zones, indicating REE fractionation from piedmont to littoral occurred during sedimentation processes. A positive correlation between (La/Yb) NASC and total REE concentrations was observed in piedmont and central sediment samples, implying that LREE were preferentially enriched over HREE with an accumulation of REE. This positive correlation was not found in littoral sediment samples, where a general decreasing trend along depth for REE concentrations existed. Results of this study shows 3 implications for using REE as a tracer in provenance studies from a basin scale.

Introduction

Rare earth elements are the lanthanides with atomic numbers ranging from 57 (La, lanthanum) to 71 (Lu, lutetium) [START_REF] Henderson | Rare earth element geochemistry[END_REF] as well as yttrium (Y, atomic number 39) and scandium (Sc, atomic number 21;[START_REF] Tyler | Rare earth elements in soil and plant systems-a review[END_REF]. The geochemical signatures of rare earth elements and yttrium (thereafter denoted REE) in sediments and sedimentary rocks are powerful tracers to delineate the geochemical processes.

These include magmatism, chemical weathering, diagenesis, erosion and water-rock interactions [START_REF] Mclennan | Rare earth elements in sedimentary rocks: Influence of provenance and sedimentary processes[END_REF]. As a result, the concentrations and fractionations of REE in sediments and soils have received considerable attention in the past few decades.

Numerous investigations have been devoted to studying geochemical behaviors and mobility of REE in sediments and soils, including lake, river and marine sediments [START_REF] Pourret | Continental shelves as potential resource of rare earth elements[END_REF][START_REF] Xu | From source to sink: Rare-earth elements trace the legacy of sulfuric dredge spoils on estuarine sediments[END_REF][START_REF] Wang | Discrimination of rare earth element geochemistry and co-occurrence in sediment from Poyang Lake, the largest freshwater lake in China[END_REF][START_REF] Consani | Baseline of distribution and origin of Rare Earth Elements in marine sediment of the coastal area of the Eastern Gulf of Tigullio (Ligurian Sea, North-West Italy)[END_REF], mining surroundings [START_REF] Hu | Physiological and biochemical effects of rare earth elements on plants and their agricultural significance: a review[END_REF][START_REF] Wang | Geochemical fractions of rare earth elements in soil around a mine tailing in Baotou[END_REF], wetlands [START_REF] Davranche | Rare earth element patterns: a tool for identifying trace metal sources during wetland soil reduction[END_REF][START_REF] Cheng | Vertical distribution of rare earth elements in a wetland soil core from the Sanjiang Plain in China[END_REF], agricultural areas [START_REF] Silva | Trace element fluxes in sediments of an environmentally impacted[END_REF][START_REF] Han | Geochemistry of rare earth elements in soils under different land uses in a typical karst area, Guizhou Province, Southwest China[END_REF], alluvial deposits [START_REF] Pédrot | Upper soil horizons control the rare earth element patterns in shallow groundwater[END_REF][START_REF] Xie | Impact of sedimentary provenance and weathering on arsenic distribution in aquifers of the Datong basin, China: Constraints from elemental geochemistry[END_REF], unconsolidated clay sediments (Guo et al., 2010;[START_REF] Sá Paye | Distribution and fractionation of the rare earth elements in Brazilian soils[END_REF] and sedimentary rocks [START_REF] Mclennan | Relationships between the trace element composition of sedimentary rocks and upper continental crust[END_REF][START_REF] Šmuc | Geochemical characteristics of rare earth elements (REEs) in the paddy soil and rice (Oryza sativa L.) system of Kočani Field, Republic of Macedonia[END_REF]). It's recognized that REE concentrations and distributions in sediments depend on sedimentary provenance and sediment physicochemical properties. Significant correlations between REE concentrations and contents of clay and metals (i.e. Fe, Mn, Al) were observed by [START_REF] Mihajlovic | Rare earth elements in German soils-A review[END_REF]; while pH value and cation exchange capacity were shown not important factors in their study. Grain size impacted REE concentrations in sediment as well. This is indicated by the results of [START_REF] Kimoto | Applicability of rare earth element oxides as a sediment tracer for coarse-textured soils[END_REF], showing that decreasing grain size generally increased REE concentrations due to that clay minerals could host REE [START_REF] Laveuf | A review on the potentiality of rare earth elements to trace pedogenetic processes[END_REF]. The impact of clay, silty and sand fractions on solid REE was recently investigated in different soil profiles [START_REF] Mihajlovic | Rare earth elements in German soils-A review[END_REF]. Presence of Fe/Mn oxides was another factor influencing REE in sediments/soils, and more REE can be accumulated by amorphous Fe/Mn oxides as compared to crystalline ones [START_REF] Yan | Sequential leachates of multiple grain size fractions from a clay-rich till, Saskatchewan, Canada: implications for controls on the rare earth element geochemistry of porewaters in an aquitard[END_REF][START_REF] Compton | Rare earth element behavior in soils and salt pan sediments of a semi-arid granitic terrain in the Western Cape, South Africa[END_REF].

Constituents of REE in sediments and soils for a large part result from natural processes. A global distribution pattern of REE in the soils/sediments was recently evaluated by [START_REF] Adeel | Cryptic footprints of rare earth elements on natural resources and living organisms[END_REF], which was achieved using spatial visualization of REE concentration distributions in soils/sediments derived from all over the world.

Results obtained from the study showed that total REE concentrations of agriculture soils/sediments ranged from 83 to 9840 mg/kg and concentrations of light REE (LREE) were higher than heavy REE (HREE). This indicates that, in addition to the natural processes, anthropogenic inputs (i.e. fertilizer) contribute to the REE concentrations in sediments/soils in some individual regions. Indeed, fertilization with REE-bearing fertilizers and infiltration of residuals was suggested to be the main mechanism of agriculture-sourced REE entering into soils and sediments [START_REF] Hu | Physiological and biochemical effects of rare earth elements on plants and their agricultural significance: a review[END_REF]. However, mechanisms for REE mobilization in different soils/sediments and the potential risks of long-term exposure to REE-abundant environment are not adequately understood.

Average REE concentration in China's soil was established to be 177 mg/kg [START_REF] Liang | Environmental biogeochemical behaviors of rare earth elements in soil-plant systems[END_REF], being comparable to the value (189 mg/kg) of the earth's crust [START_REF] Wei | A survey on the background contents of 15 rare earth elements in Chinese soil[END_REF]. China is one of the countries that early widely applied REE in various sectors [START_REF] Pang | Application of rare-earth elements in the agriculture of China and its environmental behavior in soil[END_REF][START_REF] Hu | Rare earth elements in soils[END_REF]; such as agriculture, the REE have been used as a prompter for crops since 1990s. This has led to high REE pollution levels (i.e. La >80 mg/kg; Nd >85 mg/kg; Sm >45 mg/kg) [START_REF] Adeel | Cryptic footprints of rare earth elements on natural resources and living organisms[END_REF]and references therein). Indeed, sediment and soil contaminations caused by elevated REE have been reported in various regions across China, including in Poyang Lake [START_REF] Wang | Discrimination of rare earth element geochemistry and co-occurrence in sediment from Poyang Lake, the largest freshwater lake in China[END_REF], Baotou City [START_REF] Li | Effects of the accumulation of the rare earth elements on soil macrofauna community[END_REF], Tibetan Plateau [START_REF] Wu | Distribution, pollution, and ecological risks of rare earth elements in soil of the northeastern Qinghai-tibet plateau[END_REF], Jiaozuo Bay [START_REF] Zhang | Rare earth elements in surface sediments of a marine coast under heavy anthropogenic influence: The Bohai Bay, China[END_REF], Sanjiang Plain [START_REF] Cheng | Vertical distribution of rare earth elements in a wetland soil core from the Sanjiang Plain in China[END_REF] and Liaodong Bay [START_REF] Lin | Distribution of rare earth elements in the estuarine and coastal sediments of the daliao river system, China[END_REF]. The North China Plain (NCP) is one of the largest sedimentary basin in Asia [START_REF] Xing | Groundwater hydrochemical characteristics and processes along flow paths in the North China Plain[END_REF], and is the provincial politics, economy and culture center of China. The plain is characterized by large population (> 0.35 billion) and highly-developed agricultural industry, which plays an important role in the nation's grain manufacture [START_REF] Chen | Quaternary Geology in Hebei[END_REF][START_REF] Kendy | Groundwater recharge from irrigated cropland in the North China Plain: case study of Luancheng County, Hebei Province, 1949-2000[END_REF]. The densely-populated and industrialized areas were shown to be accompanied by increased discharge of REE (i.e. Gd and Sm) into environments [START_REF] Kulaksız | Anthropogenic dissolved and colloid/nanoparticle-bound samarium, lanthanum and gadolinium in the Rhine River and the impending destruction of the natural rare earth element distribution in rivers[END_REF][START_REF] Hatje | Increases in anthropogenic gadolinium anomalies and rare earth element concentrations in San Francisco Bay over a twenty-year record[END_REF], due to a large amount of REE applied in industry, medicine as well as agriculture.

Studies concerning REE concentrations in sediments of the NCP have been performed [START_REF] Wei | Rare earth elements geochemistry and its implication for sediment provenance in the Beijing plain[END_REF][START_REF] Liu | Geochemical behaviors of rare earth elements in groundwater along a flow path in the North China Plain[END_REF][START_REF] Zhang | Rare earth elements in surface sediments of a marine coast under heavy anthropogenic influence: The Bohai Bay, China[END_REF]Liu et al., 2019a, b), but a full of investigation of mobility, distributions and behaviors of REE in sediments is still lacking. This would hinder a systematical evaluation of REE levels as well as the possible anthropogenic impact, which is vital for establishing a baseline of distribution and background of REE concentration in the contaminated sediments and for administrating regulation guidelines [START_REF] Alfaro | Rare-earth-element geochemistry in soils developed in different geological settings of Cuba[END_REF][START_REF] Consani | Baseline of distribution and origin of Rare Earth Elements in marine sediment of the coastal area of the Eastern Gulf of Tigullio (Ligurian Sea, North-West Italy)[END_REF].

Therefore, the objectives of this study are to: (i) investigate REE concentration and fractionation characteristics in sediments collected from three zones (The zonation was previously specified by [START_REF] Chen | Isotopes and sustainability of ground water resources, North China Plain[END_REF]), (ii) characterize the relationships between REE distribution and sediment properties, (iii) assess the feasibility of using REE as a probe of sedimentation process.

Material and methods

Study area

The North China Plain (NCP) is located in the eastern China with longitude of 112°30'~119°30'E, and latitude of 34°46'~40°25'N. The terrains are mostly plains with two mountainous regions (Yanshan and Tanghang mountains) bounding to the north and south, and an oceanic region (Bohai Bay) to the east [START_REF] Chen | Groundwater in Hebei[END_REF], respectively. The plain has a total area of approximately 31×10 4 km 2 and population of about 0.35 billion inhabitants. The NCP has a semi-arid humid climate with a temperature range of -23 °C~42.7 °C (in history). Average annual rainfall ranges from 500 mm to 600 mm, mostly (>74%) dominating in summer. Average annual evaporation is 1.5 to 3 times of annual rainfall.

The NCP is one of largest sedimentary basins in Asia [START_REF] Xing | Groundwater hydrochemical characteristics and processes along flow paths in the North China Plain[END_REF].

According to [START_REF] Chen | Groundwater in Hebei[END_REF], the NCP's geological strata are comprised of four stratigraphic areas: Yinshan-Nuruerhu, Yanshan, Taihangshan and Hebei, with ages ranging from the oldest Archaeozoic to the newest Holocene. Formations in study area have been declining since the Cenozoic due to neotectonics. Geographically, the plain generally decreases from the north to the east with a slope of 0.3‰~0.6‰, and thus it's divided into three zones: piedmont alluvial-proluvial plain, central alluvial-lacustrine plain and littoral plain.

Surface sediments in the NCP mainly are Cenozoic deposits with ages from Paleogene to Holocene. Lower Tertiary stratum is composed of red mudstone which is interbeded by red and gay sandstones. Upper Tertiary strata are characterized by consolidated to semi-consolidated alluvial-diluvial deposits, including siltstone, mudstone, sandstone and pebbly sandstone. Quaternary sediments are divided into four groups from top to bottom: Holocene (Q4), Upper (Q 3 m), Middle (Q 2 c) and Lower (Q 1 n) Pleistocene [START_REF] Chen | Groundwater in Hebei[END_REF]. The Holocene sediments group has a thickness of 5 m to 40 m, and mainly consists of yellow and gray sandy loam, fine sand, clay and sandy gravel. The Upper Pleistocene is alluvial, diluvial and eolian losses and yellow loam. The Middle Pleistocene stratums are widely distributed in intermountain basin and foothills, and consist of red and brown clay, loam and sandy loam, and red mud pebble layers. The Lower Pleistocene sediments group has a buried depth of 350 m to 500 m with a thickness of 80 m to 100m, and it's mainly loam, sand loam and lenticular sand layers.

Piedmont alluvial-proluvial sediments are suggested to derive from Taihang and Yanshan mountains by weathering, transport and sedimentation process. Central alluvial-lacustrine sediments are deposited from materials carried by Yellow River and Zhangweinan River running from SW to NE since middle Pleistocene [START_REF] Shao | Geomorphical structures and evolution of the Huanghuaihai plain in China[END_REF]. The littoral plain forms with interactions between terrene and ocean, and sedimentary deposits come from terrestrial rivers draining into Bohai Bay.

Sediment sampling

Sediment samples were taken from three boreholes (Y09, S30 and H02) drilled in three different zones (Fig. 1) of the NCP in the July of 2013. The core sediments were sampled with an interval generally ranging between 1 m and 2 m. Slurries and muds attached on the surface of the core sediments were scrapped off to avoid possible cross contaminations. The samples were parceled with silver paper and were stored in clean plastic bags which were immediately filled with pure N 2 gas for a further preservation. All the sealed core sediments were placed in a shaded environment with temperature below 4 °C during sampling. Finally, the samples were transported to laboratory and preserved in a freezer before analysis.

Sediment analysis

All samples used for major and minor constituents analyses were dried under room temperature in the lab, and were ground to grain size <74 μm with an agate mortar for analysis. The X-ray fluorescence (XRF) (ARL Advant X) technique was used for measurement of various oxides (SiO 2 , Fe 2 O 3 , Al 2 O 3 , MgO, CaO, Na 2 O, K 2 O, MnO, P 2 O 5 and TiO 2 ) following the national standard method (GB/T 14506).

Sediment samples were prepared with a glass flux sheet method. Specifically, 0.7 g (±0.001) sample was first placed in a 25 mL porcelain crucible; followed by addition of 5.2 g (±0.001) anhydrous lithium tetraborate, 0.4 g (±0.001) lithium fluoride and 0.3 g (±0.01) ammonium nitrate, the mixture was evenly stirred. Then the mixture was transferred to a platinum (95%)-gold (5%) alloyed crucible for fusing 10 to 15 mins under temperature 1150 to 1250 °C. Finally, fusant was made into glass beads using a fusion machine. All glass beads were preserved in a clean and dry environment before analysis. Reference materials GSS (shown below) were prepared with the same method for subsequent calibration of XRF. The sample intensity (I S ) was calculated by the difference between peak (I p ) and background (I b ) values, as shown with Eq. (2-1).

I S = I p -I b (2-1)
Rare earth elements were quantified by employing ICP-MS (AG 7500) with the digested samples. The samples were digested following the method previously described in [START_REF] Liu | Geochemical behaviors of rare earth elements in groundwater along a flow path in the North China Plain[END_REF]. Briefly, 9 mL concentrated HCl (37% v/v), 2 mL HNO 3 (65% v/v) and 9 mL HF (40% v/v) were used for digesting 0.5 g dried samples for 24 h under 108 °C condition. Then the mixture was dried and dissolved again with purified HNO 3 (2% v/v). The obtained solutions were used for REE analysis. To check the stability of the equipment system and the analytical accuracy, reference soil sample including GSS1, GSS2, GSS3 and GSS5, obtained from the Center of China National Standard Reference Material, were ran simultaneously. The recommended and test values of REE and metal oxides in standard references were shown in supplementary materials (Table S1). Furthermore, the parallel sample was performed every ten samples for further monitoring of the analytical precisions. The results show that elements mostly had accuracies (RSD, δ1) ranging between 2% and 5%. The standard materials had accuracy generally better than 3%, except for P 2 O 5 , which had an accuracy better than 8% in some samples. The analytical precisions of REE in all samples were generally better than 5%.

REE normalization

Rare earth element concentrations in all sediment samples were normalized with average REE concentrations of the North American Shale Composite (NASC). The corresponding normalization REE values were taken from [START_REF] Taylor | The Continental Crust: Its Composition and Evolution[END_REF]. Fractionation between LREE and HREE was quantified using fractionation where the superscript * indicates the geogenic background.

Calculation of mass transfer coefficient

To evaluate REE mobility, the mass transfer coefficient (τ) was calculated by the Eq. (2-5). element i. This study chooses Ti as the reference because it has been used previously [START_REF] Wu | Distribution, pollution, and ecological risks of rare earth elements in soil of the northeastern Qinghai-tibet plateau[END_REF].

τ i, REE = (C REE /C i ) sample / (C REE /C i ) NASC -1 (2-5)

Results and discussion

Sediment geochemistry and provenance

Sediment lithology in three boreholes were mainly clay, silt and silty clay (Fig. 2). Fine and course sands were generally interbedded with silt and silty clay at different sampling depths. This was well observed in the central (S30) and littoral (H02) borehole bar graphs (Figs. 2b andc). In piedmont borehole, a gravel layer occurred at depth between 38 m and 61 m (Fig. 2a). Gravel was not found in boreholes S30 and HH02. Overall, fine grains prevailed in central and littoral sediment samples. Colors of the fine grains changed in three zones as well. Pale yellow and pale brownish-yellow clay and silty were mainly distributed in piedmont sediments. For central borehole sediments, the color changed from pale brownish-yellow to brown and tan vertically. For littoral sediments, taupe and tawny silty and silty clay were mainly in upper and middle parts, and pale brownish-yellow and pale brown in lower part (Fig. 2c).

The component of sediment samples was dominated by SiO 2 ( Changes of sediment lithology were controlled by the sedimentary provenances and processes [START_REF] Bhuiyan | Provenance, tectonics, and source weathering of modern fluvial sediments of the Brahmaputra-Jamuna River, Bangladesh: inference from geochemistry[END_REF][START_REF] Blake | Application of REE geochemical signatures for Mesozoic sediment provenance to the Gettysburg Basin, Pennsylvania[END_REF]. This would account for that sediments were deposited alternately with clay, silty clay and fine sand (Fig. 2).

Terrigenous clastic deposition and land-ocean interactions led to the formation of littoral sedimentation. Weathering and erosion of the bedrocks around the piedmont plain and subsequent transport and sedimentation contributed to the NCP sediments.

Ternary diagram of sandstone classification [START_REF] Blatt | Origin of Sedimentary Rocks (2nd edition)[END_REF] showed that the majority of sediment samples were located in the greywacke domain, and three central samples in the arkose domain indicating that sediments originally mostly were greywacke (Fig. 3a). The dominant SiO 2 contents and fine-grained components reflected that the sedimentary materials had been undergone intensive denudations and a long-distance transport before sink. Indeed, the ternary diagram of

Al 2 O 3 -CaO+Na 2 O-K 2 O (A-CN-K)
showed that the weathering trend of the sediment samples were mostly along the axis of (CaO+Na 2 O) (Fig. 3b). This probably indicated that sediment samples from the NCP were in an early weathering stage where minerals such as plagioclase was incongruently decomposed in association with leaching of Na and Ca generating products of illite, kaolinite and smectite [START_REF] Nesbitt | Chemical processes affecting alkalis and alkaline earths during continental weathering[END_REF]. The offset to the uppermost of K 2 O axis was attributed to the addition of terrigenous detrital in the later sedimentary stage (Cox et a1., 1995).

To quantify the weathering intensity in different samples, the chemical index of alteration (CIA) values (initially was valid for granitic rocks) were calculated following the method proposed by [START_REF] Nesbitt | Early Proterozoic climates and plate motions inferred from major chemistry of lutites[END_REF] [START_REF] Panahi | Behavior of major and trace elements (including REE) during Paleoproterozoic pedogenesis and diagenetic alteration of an Arehean granite near Ville Marie, Québee, Canada[END_REF]].

[CIA = (Al 2 O 3 /(Al 2 O 3 + K 2 O + Na 2 O + CaO*))×100; CaO* = [CaO -(10/3×P 2 O 5 )
Chemical index of alteration can be an index reflecting sample in a weathering trend from fresh status to intensely-altered status. Granites and granodiorites usually have CIA values ranging from 45 to 50, and fresh basalts have values between 30 and 45.

Muscovite has a value of 75, and illite, montmorillonites and beidellites give values from 75 to 85 [START_REF] Nesbitt | Early Proterozoic climates and plate motions inferred from major chemistry of lutites[END_REF]. Results showed that CIA value of piedmont sediment samples ranged from 55 to 78 with an average value of 68, and central and littoral samples ranged from 37 to 70 (average 51), and from 43 to 70 (average 55), respectively. This indicates that sediment samples were incipiently (CIA: 50 to 60) to intermediately (CIA: 60 to 80) chemically weathered, according to the classification standards suggested by [START_REF] Fedo | Unraveling the effects of potassium metasomatism in sedimentary rocks and paleosols with implications for paleoweathering conditions and provenance[END_REF]. It's also observed that piedmont sediment samples generally had higher CIA values as compared to those of central and littoral sediments, which means a stronger weathering may have occurred in the piedmont sediment samples. These differences probably were associated with the hydrogeological conditions such as REE residence time, outcrop lithology and aquifer replenishment conditions [START_REF] Condie | Another look at rare earth elements in shales[END_REF][START_REF] Blake | Application of REE geochemical signatures for Mesozoic sediment provenance to the Gettysburg Basin, Pennsylvania[END_REF]. The results were consistent with piedmont samples being located in the upper part of the A-CN-K diagram, while central and littoral sediment samples below the piedmont samples (Fig. 3b). The role of chemical weathering on sedimentary processes was studied previously based on sediment geochemical compositions (LaMaskin et a1., 2008;[START_REF] Dostal | Geochemistry of low-grade ctastic rocks in the Acatl An Complex of southern Mexico: Evidence for local provenance in felsie-intermediate igneous rocks[END_REF][START_REF] Xie | Impact of sedimentary provenance and weathering on arsenic distribution in aquifers of the Datong basin, China: Constraints from elemental geochemistry[END_REF]. Combined to studies of the residence time of REE [START_REF] Condie | Another look at rare earth elements in shales[END_REF] [START_REF] Lan | Distribution characteristics of rare earth elements in surface sediment and their provenance discrimination in the eastern Bohai and northern Yellow Seas[END_REF]. Others showed that core sediment ΣREE concentrations were slightly higher than the value of surface sediment (Liu et al., 2019b).

Sediment ΣREE concentrations varied along depth. In the piedmont borehole, ΣREE concentrations were relatively constant within 11 m below land surface. A small decrease in ΣREE concentrations was observed at depth from 11 m to 16 m before they increased to >300 mg/kg at 25 m. After a new decrease to about 63 m, ΣREE concentrations fluctuated for the remaining sampling depths (Fig. 2a). ΣREE concentrations in central sediments were highly variable along borehole, although a decreasing trend occurred from land surface to 30 m. Vertical changes of ΣREE concentrations were mainly related to sediment physicochemical properties [START_REF] Kimoto | Applicability of rare earth element oxides as a sediment tracer for coarse-textured soils[END_REF][START_REF] Mihajlovic | Rare earth elements in German soils-A review[END_REF]. Sediments with different grain sizes, such as clay, silt and silty clay, had different capacities in hosting REE. The major minerals in the aquifer sediments of study area were shown to be quartz, K-feldspar, hornblende and calcite accounting for 5%-14% in the piedmont [START_REF] Chen | Quaternary Geology in Hebei[END_REF].

Clay minerals (i.e. kaolinite, illite and montmorillonite) were more commonly distributed in the central and littoral plains as compared to the piedmont plain [START_REF] Chen | Isotopes and sustainability of ground water resources, North China Plain[END_REF]. In the piedmont, where fine sand prevailed (i.e. from 38 m to 62 m), ΣREE concentrations progressively decreased (Fig. 2a). Low ΣREE concentrations of central sediments were found at depth of approximately 30 m, 65 m to 68 m, 107 m, 151 m to 154 m, where fine sands were sandwiched with silt and silty clay (Fig. 2b).

Silt sand and fine sand more commonly occurred at deeper depth, leading to a decreasing trend of ΣREE concentrations (Fig. 2c). This demonstrated that presence of clay was one of factors controlling ΣREE concentrations. Lanthanides were thought to be hosted by clay component containing in the sediments [START_REF] Yang | The rare earth element compositions of the Changjiang (Yangtze) and Huanghe (Yellow) river sediments[END_REF]. The disparities of clay content along depth for a large part accounted for REE variations along depth. This was supported by the vertical variations of Si/Al values, which were proxies for grain size classification (i.e. clay, silt, and silty clay) (Fig. S2).

Statistical analyses showed that ΣREE concentrations were negatively correlated to SiO 2 and were significantly positively correlated to Fe 2 O 3 and Al 2 O 3 (except for littoral sediment samples) ( [START_REF] Ohta | Rare earth element partitioning between Fe oxyhydroxide precipitates and aqueous NaCl solutions doped with NaHCO 3 : Determinations of rare earth element complexation constants with carbonate ions[END_REF][START_REF] Pourret | Rare earth element sorption onto hydrous manganese oxide: A modeling study[END_REF][START_REF] Liu | Rare earth elements sorption to iron oxyhydroxide: Model development and application to groundwater[END_REF][START_REF] Mihajlovic | Rare earth elements in soil profiles of various ecosystems across Germany[END_REF]. These findings are in line with the results documented by [START_REF] Wang | The surface sediment type and their rare earth element characteristics from the continental shelf of the northern South China Sea[END_REF], showing that fine fractions with abundant clay minerals and Fe/Al/Mn-containing oxides tended to accumulate REE, while coarse fractions with less clay minerals and more quartz led to less REE contents. However, further information on mineralogical compositions of sediment samples is needed for better constraining the impact of mineralogy on REE concentrations from a quantitative perspective. In all, sediment REE concentrations were essentially regulated by their genetic properties and were revealed by their vertical variations.

REE fractionation patterns

All sediment samples had relatively flat NASC-normalized REE patterns with a gentle enrichment of LREE over HREE (Fig. 4). The degree of LREE enrichment was quantified by (La/Yb) NASC values, which ranged between 0.86 and 1.15 with an average of 1.23 in piedmont sediments. Central and littoral sediment samples had (La/Yb) NASC ranges of 1.03 to 1.14 (average 1.21) and 0.91 to 1.49 (average 1.09), respectively (Table 1). Sediments mostly showed negative Ce anomalies with Ce/Ce* ranging between 0.69 and 1.17 (average 0.91), and between 0.85 and 0.98 (average 0.90), and between 0.80 and 0.98 (average 0.89) in piedmont, central and littoral plains, respectively. Negative Eu anomalies commonly occurred in investigated sediments with Eu/Eu* ranges of 0.82 to 1.07 (average 0.91), 0.83 to 1.31 (average 0.95), and 0.86 to 1.34 (average 0.95) in three zones, respectively. Correlation analysis indicated that ΣREE concentrations generally increased with an increase of (La/Yb) NASC in piedmont and central sediments, while this trend was not observed in littoral sediment samples (Fig. 5a). These results reflected REE were nonuniformly accumulated by sediments during their transport. The correlation between (La/Yb) NASC and CIA values showed that the larger (La/Yb) NASC values occurred in samples containing higher ΣREE concentrations and generally having higher CIA values (Fig. 5b). Thus, LREE enrichments could result from the chemical weathering of silicates in crustal materials, which led to higher levels of LREE in sediments [START_REF] Oliveira | Rare earth element patterns in lake sediments as studied by neutron activation analysis[END_REF][START_REF] Lin | Distribution of rare earth elements in the estuarine and coastal sediments of the daliao river system, China[END_REF]. Larger CIA values were obtained from piedmont and central zones as compared to the values in littoral zones (shown above).

The different domains of littoral sediment samples from their piedmont and central counterparts might be related to impact of ocean and/or human activities. However, more studies are needed to verify this speculation.

Sediment samples mostly had negative Eu anomalies (Fig. S2). Average Eu/Eu* anomalies were quite similar in three different zones. This indicated that sedimentary processes, which evolved weathering of bedrocks around plains and migration and deposition in the flat plains, did not fractionate Eu to a larger extend with respect to its neighboring Sm and Gd. Early investigations demonstrated that Eu compositions and distributions in continental crust were fundamentally controlled by shallow, intracrustal differentiation involving Eu-bearing minerals (i.e. plagioclase) [START_REF] Mclennan | Relationships between the trace element composition of sedimentary rocks and upper continental crust[END_REF], which was shown to possess 4% to 23% (wt) in piedmont sediments [START_REF] Liu | Geochemical behaviors of rare earth elements in groundwater along a flow path in the North China Plain[END_REF]. However, Eu/Eu* values were negatively correlated to ΣREE concentrations (Fig. 5d), suggesting that evolution of Eu anomalies was related to changes of redox conditions as well during sedimentation. Preferential mobilization of Eu(II) over other trivalent REE leads to a segregation of Eu from lanthanide series under reducing conditions [START_REF] Lee | Rare earth elements as indicators of groundwater environment changes in a fractured rock system: evidence from fracture-filling calcite[END_REF], and hence higher ΣREE concentrations occur in lower-Eu/Eu* samples, as shown in Figure 5d. This mechanism has been proposed for using Eu anomalies as an indicator for investigating REE origins as well as mobility of redox-sensitive metal(oid) (Guo et al., 2010;[START_REF] Davranche | Rare earth element patterns: a tool for identifying trace metal sources during wetland soil reduction[END_REF].

Another possibility resulting in negative correlation between Eu anomalies and ΣREE concentrations would be related to minerals contained in sediments. Such as plagioclase, which tends to preferentially accumulate Eu(II) over trivalent REE by means of isomorphism with Ca 2+ , Na + and K + [START_REF] Banks | Distribution of rare earth elements in crystalline bedrock groundwaters: Oslo and Bergen regions, Norway[END_REF].

Cerium anomalies (Ce/Ce*) were relatively constant in central and littoral sediments despite variable ΣREE concentrations, while they were distributed scattered as a function of ΣREE concentrations in piedmont sediment samples (Fig. S2). The fluctuation of redox conditions could account for this phenomenon due to fact that Ce is readily to be precipitated as insoluble Ce(IV) via oxidative scavenging, and when conditions evolve into reducing, the precipitated Ce will be released into solution by reductive dissolutions [START_REF] Bau | Oxidative scavenging of cerium on hydrous Fe oxide: Evidence from the distribution of rare earth elements and yttrium between Fe oxides and Mn oxides in hydrogenetic ferromanganese crusts[END_REF][START_REF] Steinmann | Controls on transport and fractionation of the rare earth elements in stream water of a mixed basaltic-granitic catchment basin (Massif Central, France)[END_REF].

Cerium cycle was believed to be biologically mediated [START_REF] Moffett | Microbially mediated cerium oxidation in seawater[END_REF][START_REF] Tanaka | A specific Ce oxidation process during sorption of rare earth elements on biogenic Mn oxide produced by Acremonium sp. Strain KR21-2[END_REF] and involved metal (i.e. Fe, Mn) oxides [START_REF] Yu | Cerium sequestration and accumulation in fractured crystalline bedrock: The role of Mn-Fe (hydr-) oxides and clay minerals[END_REF]. Results of our recent study showed that Ce precipitation from oxic aqueous solution was more likely to occur in the piedmont aquifers where low-Fe/Mn-concentration groundwater prevailed (Liu et al., 2019a), though no pronounced positive Ce anomalies were commonly observed in the piedmont aquifer sediments. Thus, both positive and negative Ce anomalies, and no Ce anomalies can be found in solid phases in nature, as the results shown in present study.

Implication for REE signatures

Statistical analyses showed that individual REE concentrations were significantly correlated with correlation coefficient (R 2 ) all better than 0.8. [START_REF] Verplanck | Aqueous Stability of Gadolinium in Surface Waters Receiving Sewage Treatment Plant Effluent, Boulder Creek, Colorado[END_REF][START_REF] Hatje | Increases in anthropogenic gadolinium anomalies and rare earth element concentrations in San Francisco Bay over a twenty-year record[END_REF], agricultural activities [START_REF] Alfaro | Rare-earth-element geochemistry in soils developed in different geological settings of Cuba[END_REF] and atmospheric depositions [START_REF] Tyler | Rare earth elements in soil and plant systems-a review[END_REF], although the investigated sediments were unlikely to suffer from these anthropogenic pollutions due to that most of them were not formed nowadays. Positive REE anomalies will show up in REE normalized patterns when individual elemental concentration is greater than its background value, such as Gd anomalies, which have been reported worldwide [START_REF] Ebrahimi | Gadolinium as an emerging microcontaminant in water resources: threats and opportunities[END_REF]and references therein). In summary, sedimentary REE patterns conveys important information regarding the composition of the continental crust and anthropogenic impact.

Conclusions

The 226 sediment samples were collected from three different zones (i.e.

Figure 1 Figure 2 Figure 3 Figure 4 Figure 5

 12345 Figure 1 Study area and sediments sampling locations (boreholes)

  Fig. 1

  Fig. 3

  The τ i, REE values are an indicator of depletion or enrichment of REE after normalization against reference element i in standard materials (NASC for this study). Positive τ i, REE values mean enrichment of REE (total or individual) in the sediments (s) with respect to NASC; negative values indicate depletion and zero indicate REE are immobile. Elements include Fe, Al, Mn, Ti, and Ca can be used as reference

Table 1

 1 Iron oxide (Fe 2 O 3 ) was the second abundant metal oxides, with a range from 1.86 wt% to 9.86 wt% and an average of 5.46 wt% in piedmont sediments. The average content (wt%) of Fe 2 O 3 in central and littoral sediments was 4.45 wt% and 4.06 wt%, respectively, exhibiting slightly lower than that of piedmont sediments. Other metal oxides including Na 2 O, K 2 O, CaO, MgO, P 2 O 5 and TiO 2 showed relatively lower content in most of the measured sediment samples as compared to SiO 2 and Al 2 O 3 (Table1). Average value (wt%) of loss on ignition (LOI) was 3.88, 8.56 and 8.02 in piedmont, central and littoral zones with ranges of1.13 to7.36, 3.12 to 19.91 and 2.5 to 17.56, respectively.

	). For

Table 2 )

 2 , reflecting that REE mainly adsorbed onto/incorporated into Fe and Al-containing clay minerals. Normalization of REE to Ti showed that τ i, REE values mostly changed slightly, and generally were positive for piedmont sediments and negative for central and littoral sediment. This suggested that dilution of weathered residuals (i.e. SiO 2 ) did not significantly modify REE concentrations in sediment samples. Moreover, trends of SiO 2 /Al 2 O 3 values were generally consistent with those of CIA, suggesting that weathering of silicates like feldspar played a role in this process (Fig.2). In other words, in addition to clay content, Fe and Al-bearing oxides could be another factor influencing REE concentrations in investigated sediment samples. Therefore, mixture of clay minerals and Fe-Mn oxides accounted for normalized REE values. Metal oxides containing Fe, Mn and Al provide important binding sites for sorption of REE

Table 1

 1 Concentrations of REE and major oxides and REE fractionation indices in the sediment samples from three zones

	indicated that sediments mostly were derived from greywacke, as suggested by
	triangular diagram of sandstone classification. All sediment samples had been
	undergone incipient to intermediate weathering with average CIA of 68, 51 and 55 in
	three zones, respectively. Total REE concentrations ranged from 190.24 to 353.87
	mg/kg (average 228.62 mg/kg), and from 88.85 to 290.80 mg/kg (average 182.88
	mg/kg), and from 87.90 to 237.61 mg/kg (average 176.69 mg/kg) in piedmont, central
	and littoral plains, exhibiting a decreasing trend as piedmont > central > littoral.
	Sediment physicochemical properties such as mineral composites and metal oxide
	abundance controlled REE concentrations, as observed along depth in three sediment
	sampling profiles. All sediment samples had coherent NASC-normalized REE
	patterns that were characterized by light REE (LREE) enrichment over heavy REE
	(HREE) with (La/Yb) NASC ranging between 0.86 and 1.49. Negative Ce and Eu
	anomalies generally occurred in the investigated sediments with Ce/Ce* and Eu/Eu*
	ranges of 0.69 to 1.17 and 0.82 to 1.34, respectively. The present study provides
	important dataset for future comparative studies in this region.

piedmont, central and littoral) of the north China plain to investigate the concentrations and fractionations of rare earth element (REE). Sediment geochemistry

Table 2

 2 Pearson's correlations between REE concentrations and major oxides, CIA

Acknowledgments

This investigation has been funded by National Natural Science Foundation of China (Nos. 41902243, 41222020 and 41172224), National Natural Science Foundation of Jiangxi Province, China (20202BABL211018), the program of China Geology Survey (No. 12120113103700), and the Fok Ying-Tung Education Foundation, China (Grant No. 131017). East China University of Technology Research Foundation for Advanced Talents (Nos. DHBK2019094 and SHT201901).