{'fraction_non_alphanumeric': 0.039941812136325854, 'fraction_numerical': 0.011596009975062344, 'mean_word_length': 4.310306775546237, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Contact Surface Estimation via Haptic Perception

Legged systems need to optimize contact force in order to maintain contacts. For this, the controller needs to have the knowledge of the surface geometry and how slippery the terrain is. We can use a vision system to realize the terrain, but the accuracy of the vision system degrades in harsh weather, and it cannot visualize the terrain if it is covered with water or grass. Also, the degree of friction cannot be directly visualized. In this paper, we propose an online method to estimate the surface information via haptic exploration. We also introduce a probabilistic criterion to measure the quality of the estimation. The method is validated on both simulation and a real robot platform.


I. INTRODUCTION
Legged robots have many practical applications in both industry and the public sector; particularly, in places that are too dangerous for humans, such as mining inspections in underground tunnels or the aftermath of an earthquake. The need for legged systems instead of wheeled systems is apparent since the passages are normally uneven and difficult to maneuver.
Yet, vision systems may degrade in harsh conditions, such as in dusty and foggy environments. Moreover, vision systems are unable to detect/realize terrain under cover (e.g., A LiDAR system may not be able to distinguish a water surface or high grass from a solid contact.) The challenge here is how can the robot maneuver over an unknown environment when the perception system fails?
Model-based control has shown promising results in the robotic community, ranging from bi-manual manipulation [1], quadruped locomotion [2], to whole-body humanoid control [3] [4]. This control mechanism is particularly important for legged locomotion, where the robots need to maintain appropriate contact forces at the point-of-contacts to avoid slippage. However, those methods require the terrain friction and surface geometry to be provided a priori, since this information is needed for contact force optimization [5] [6]. Additionally, the surface properties are also required for footstep planning over rough terrain [7] [8].
The surface geometry can be gathered via vision system [9] [10]; on the other hand, the accuracy can be disturbed in harsh weather conditions. Besides, the surface friction can not be directly visualized. Without sensing skills or prior information about the environment, this framework could only be applied successfully in laboratory settings.  One potential solution is to utilize a single foot as a haptic sensor while standing on the other three. This concept was motivated by how humans explore the environment. For example, when entering a dark room for the first time, we use our hand to touch our surroundings to avoid bumping into a wall, before making a move. The robot can use a single foot to touch the surface near the next desired footstep position and estimate the terrain information.
Previous work on estimation of surface normal takes the learning by demonstration approach, where the users move the robot on the contact surface to sense the constraint, and then estimate the constraint offline [11] [12][13] [14]. However, it is not practical to provide demonstrations to all environments for the robots where the constraints vary. Instead of learning from human demonstration, in [15], the robot performs a set of trajectories to collect the motion data. Nevertheless, the data collection took 2 minutes in this work, which is also not practical in real robotic applications.
We can also find work related to the online estimation of terrain information in the literature of quadruped locomotion. [16] estimates the height of the support surface using proprioception and LiDAR, and [17] proposed a method to find the surface normal under the assumption that terrain is consistent for all standing legs. In contrast, our work relaxes such assumptions.
Additionally, the estimation of the friction coefficient has been carried out in many works. While most measures the contact force using force-torque sensor [18] [19], no sensor is required in our work. On top of that, previous work only attempted to find the friction coefficient on a flat surface.
In this paper, we propose a method to estimate the surface properties when the robot loses information about the terrain. The contribution of this paper includes (1) an efficient online method that estimates the surface normal, (2) a haptic sensing strategy that senses the environment by the robot itself, and (3) a probabilistic criterion that evaluates the result of the estimation, so the robot can stop exploring when it is confident about the estimation. These methods have been  Figure 1).

II. BACKGROUND
First, we define the problem and clarify what quantity we aim to estimate in §II-A and briefly review the Bayesian inference method in §II-B. Finally, we highlight our contribution in §II-C.

A. Problem Definition
Let q,q,q ∈ R d+6 be the generalized positions, velocities, and accelerations of a floating-based n degree-of-freedom robotics system, the rigid body dynamics with k contacts can be described as where M ∈ R (d+6)×(d+6) is the inertia matrix, h ∈ R (d+6) is the Coriolis, centrifugal, and gravitational force, B ∈ R (d+6)×d is the selection matrix, τ ∈ R d is the control torque, J c,i ∈ R 3×(d+6) is the constraint Jacobian that relates the joint space to the i th stance foot, and λ c,i ∈ R 3 is the reaction force at the i th contact. For model-based controllers, such as [2][6], the goal is to find the control torques for the desired tasks while satisfying a set of kinematics and dynamics constraints. A typical constraint that avoids robots from slippage is the friction cone constraints [21]. Given a unit vectorn normal to the contact surface, we can divide the contact force λ c into three components such that λ n =n λ c is normal to the contact surface, and λ x , λ y are tangential to the contact surface. By Coulomb's Law, the magnitude of the tangential force should not exceed the friction coefficient µ times the magnitude of the normal force to avoid slipping where µ is the friction coefficient which depends on the material of the surface. An example can be seen in Figure 2. The surface normal n determines the direction of the friction cone, and the friction coefficient µ determines the shape of the friction cone. Geometrically, the controller needs to find the torque τ such that the resulting contact force λ c lies with the friction cone.
Ideally, we need to provide then and µ to the controller before the optimization of control torques. Instead of providing this information to the controller a priori, the goal of this paper is to explore and estimate these values by the robot itself.

B. Bayesian inference
Given a random variable X with unknown distribution and a set of observations D, Bayesian inference [22] is a method to infer the probability of the unknown variables X. The prior distribution P(X) is the distribution of the parameters before the observations. The likelihood P(D|X) is the distribution of the observations, conditional on the variable. The posterior probability P(X|D) is the conditional probability after some observations, which is determined by Bayes' rule If the prior and posterior have the same distribution, the prior is the conjugate prior of the likelihood distribution.
In this work, we would like to know the probability of our estimationn and µ without knowing the true distribution of these two variables. In §III-C and §IV-B, we will use Bayesian inference with conjugate prior to approximate the probability of our estimation.

C. Contributions
In this paper, our goal is to estimaten and µ with minimum prior information. The contributions of this paper include the following: • Previous work on estimation ofn takes the offline learning approach [11][12] [13], while we introduced an efficient online. • In [15], the robot explores the environment by sliding on the contact surface. This work took 2 minutes to finish one estimation since it was hard to slide on an unknown surface. Here, the robot explores by touching and finishes in a realistic time frame. • Instead of using force/torque sensor to estimate µ [18] [19], our method does not require the F/T sensor, which tends to be noisy. • While previous work cannot quantify how well the estimation is, we proposed a probabilistic criterion that evaluates the results of the estimation, so the robot may stop exploring when it is confident about the estimation.

III. SURFACE NORMAL ESTIMATION
In this section, we describe our method for surface normal estimation. When the robot loses its visual accuracy, we will collect data of the surface via haptic exploration. Namely, we use a single foot as a haptic sensor while standing on the other three.

A. Data Collection via Haptic Exploration
Let p,ṗ,p ∈ R 3 denotes the position, velocities, and accelerations of the feet. The joint space and end-effector space are related through the Jacobianṗ = J(q)q and p = J(q)q +Jq. Let p * ∈ R 3 be the desired foot position for the current swing leg according to some predefined footstep planner. Before making that contact and shift the center-of-mass, the swing foot tries to explore the area around p * to sense the area. An example can be found in Figure 3, the left front leg is the swing leg, and the red point is the desired foot position p * for this leg.
The robot explores the area around p * by making K contacts with that area according to where α defines the magnitude of δ i (i.e., how far it is from the p * ) and θ k defines the angle between [1, 0] and δ k . In an ideal world, only K = 3 points are needed to define the normal vector. In reality, more contact points will be needed due to sensory noise.
Example Let p * = [0, 0, 0] be the desired position for the n/ext footstep. We use (3) to set the exploration points. Let α = 1 and θ k = [0, 90, 180, 270], the robot will try to make contact at p * . For each k exploration, when the robot detects collision between the swing foot and the contact surface, the contact positions of the foot c ∈ R 3 are recorded. (Note that, c may not be identical with p, since there might be an early collision and/or tracking error) A set of K data points will be collected P = [c 1 , c 2 , . . . , c k ]. An example of exploration can be found in Figure 4.
In the above example, the magnitude of the exploration is set to α = 1. In practice, a reasonable choice of α is proportional to the size of the foot. Recall that, the reason for gathering the terrain properties is to optimize the contact force for the stance legs, and the only information that matters is the terrain underneath each foot. This simplifies the complexity of dealing with highly uneven terrain.

B. Estimation of surface normal
After the exploration, we have a set of K points that v k = c i − c j ∀ i = j is a vector that connects two contact points (see Figure 5).
Remark D can be interpreted as a set of vectors that lie in the contact surface S. This also means, all vectors in D are orthogonal to the surface normal., i.e., v kn = 0, ∀v k ∈ D. Therefore, we seek to solve the homogeneous linear equations D n = 0 .
Let D = U S V be the singular value decomposition of D such that U is the matrix of left singular vectors, S is a diagonal matrix of singular values, and V is a matrix of the right singular vectors. Let A ∈ R 3×3 = DD such that Remark In (5), we aim to make A a Hermitaian matrix (i.e., A ≡ A ) such that the eigenvalues S S ∈ R are always real, and it is always possible to find an orthogonal basis. The surface normaln can be found by taking the columns of V that corresponds to the smallest singular value.
The rank of the constraint can be seen from the rank of S S: • rank(S S) = 1: The observations D spans a plane.
There is a 1-dimensional constraint and the robot can stand on a surface. This observation means that A has a vanishing singular value, and any linear combination of the corresponding right-singular vectors is a valid solution. Since we are looking for a unit vector, the solution is the right-singular vector.
• rank(S S) = 2: There is a 2-dimensional constraint and the robot can stand on a line. If this is the case, we know the desired footstep position p * is not a good choice for the robot since the foot is likely to be stuck. This is a special case and out of the scope of this work.

C. Stopping Criterion
In this section, we introduce a criterion to evaluate the confidence of the estimation. Specifically, we would like to know how well our estimated contact normaln is without the knowledge of the true contact normaln * . We will use the conjugate prior [22] to represent the unknown mean and variance of the ground truth, and then infer the probability of our estimation. 1) Likelihood: Given a set of observations D; we would like to know what is the likelihood of this observation given the current estimate of contact normal. We assume that the inner product of v andn has normal distribution v n ∼ N (µ D , σ 2 D ) where µ D is the observation mean, and σ 2 D is the observation variance. The likelihood of the observations given the estimation P(D|n) can be inferred from the probability density function If the estimation is correct, then (4) is satisfied, so we can assume µ D = 0. Let = 1 σ 2 D denotes the precision of the observations. The above equation can be re-written as P(D|n, ) = 1 (2π) n/2 n/2 exp − 2 k (v kn ) 2 2) Prior distribution: The conjugate prior of a normal distribution with unknown mean and variance can be represented by Normal-gamma distributionn ∼ N G(ñ, κ, α, β) whereñ, κ, α, β are the unknown hyper-parameters of the distribution. The probability ofn given is defined by 3) Initial values: Initially, we assume the robot is going to walk over flat terrain, so the initial guess of the contact normal is set toñ 0 = [0, 0, 1] The parameters of the prior distribution are set to α 0 = 0.1, β 0 = 0.1, κ 0 = 1.

4) Posterior distribution:
The likelihood of the estimation n k given the observed data samples D, or posterior distribution, can be written as P(n k |D) = P(D|n k )P(n k ). By the property of Normal-gamma distribution, if the prior is a Normal-gamma distributionn k ∼ N G(ñ k , κ k , α k , β k ), the posterior can be proven to be a Normal-gamma distribution n k+1 ∼ N G(ñ k+1 , κ k+1 , α k+1 , β k+1 ).
Remark By multiplying (7) and (8) and simplifying the expression, the posterior distribution can be described by another Normal-Gamma distribution. For each iteration, we As soon as the cumulative density is high enough, i.e., P(n k |D) ≥ Λ , we can stop the exploration.

IV. SURFACE FRICTION ESTIMATION
After learning the surface normal, we can decompose the contact force into normal and tangential parts. In this section, we describe a method to estimate the friction coefficient µ.

A. Surface Friction Estimation
Since we use the swing leg to perform the exploration and standing with the other three legs, the force exerted by the robot is the only force that acts on the contact surface. Therefore, the action force (e.g., the force applied by the foot) is equal to the reaction force (e.g., the ground reaction force) but the opposite direction; specifically, where f w ∈ R 3 is the force exerted by the robot at the swing foot. We use f ⊥ and f to denote the normal and tangential components of f w , and f w = f ⊥ + f (see Figure 6). During the last haptic exploration, we collected a set of contact points c 1 , c 2 , ..., c K . Since all these points are on the contact surface, we want to move the feet around these points. Namely, we generate the desired trajectory by interpolating the positions between all c i , and this trajectory will provide the desired foot position for the impedance controller Λp + Dṗ + Kp (11) wherep = p * w − p is the error between c i and the foot position, and p * w denotes the desired foot position during exploration. To ensure that the motion is moving on the contact surface, we project the motion onto the contact surface by Remark Sincen is a unit vector, the pseudo inverse is equal to its transpose. Therefore, the projection is defined by I 3 −n n instead of I 3 −n †n From (2), we know that the foot cannot slide on the contact surface if µλ z ≥ λ 2 x + λ 2 y , and the vice versa. Substituting (10) into (2) indicates that −µf ⊥ ≥ ||f || → foot does not move −µf ⊥ < ||f || → foot may move (13) Therefore, we can vary the values of µ ∈ [0.2, 1.0] to get the normal force During the exploration, we start with a large friction coefficient. If there is a slippage, the friction coefficient will be reduced. Finally, the total force at the swing foot will be f w = f ⊥ + f , and the torque of the swing leg is

B. Stopping Criteria
During the exploration, we would like to know how likely our estimated friction coefficient µ is smaller or equal to the true friction value µ * , so we can stop exploration if we are confident about our estimation. For this, we will use the Bernoulli distribution with Beta distribution as the conjugate prior to [23] model the probability distribution of our estimation. 1) Likelihood distribution: During the exploration, the swing foot is either sliding or non-sliding on the contact surface. We represent the outcome of the k th step as a boolean variable z k ∈ {0, 1} where z k = 1 denotes noslippage and z k = 0 denotes slippage. Since we don't have the true friction coefficient, we cannot use (13) to determine slippage; instead, the slippage is determined by how much the desired end-effector velocity v * = 0 is different from the observed end-effector velocity v on the contact surface: Since the observation is binary, we represent the likelihood as Bernoulli distribution z k |µ k ∼ Bernoulli(P(µ k )), and the probability density is defined as Namely, if there is no slippage, the probability that the observation should not slip is the probability that the current friction coefficient setting is sufficient.
2) Prior distribution: The Beta distribution is the conjugate prior of the Bernoulli distribution. Although we do not have the true distribution of µ, the prior distribution can be described as where B denotes beta distribution and a k , b k > 0 are the unknown shape parameters of the distribution. The probability of the prediction P(µ k ) can be estimated by the probability density function of beta distribution.
3) Initial Guess: The initial guess of the friction coefficient is µ 0 = 1.0, since it is larger than friction coefficients in most scenarios. These shape parameters are set to a 0 = 1, b 0 = 1 at the beginning of the experiment. This is equivalent to a uniform distribution on the interval of 0 and 1 where all potential solutions have the same probability.

4) Posterior distribution:
If the prior is of the form of a Beta distribution, then the posterior will be of the same distribution [23]. Therefore, the posterior distribution µ k+1 , given samples z k also has a beta distribution Remark For each iteration, we only need to update the hyper-parameters of the distribution to evaluate the estimation.
Then, we update the prior by substituting in the posterior and repeat. This process continues until we have enough confidence in the estimate. Specifically, we want to ensure the estimated friction coefficient is smaller than the true friction coefficient. The estimation is completed if where P(µ * > µ k ) = 1 − P(µ * < µ k ) ∈ [0, 1], µ is a predefined threshold function, and P(µ * ≤ µ k ) is the cumulative distribution function of beta distribution.

V. EVALUATION
In this section, we describe how we validate our proposed method, in both simulation and on a real robotic platform.

A. Learning surface normal
In the first experiment, we would like to demonstrate the idea in simulation. First, we made different slopes in simulated scenarios, where the true inclination of the surface can be controlled. We tested on the slope with inclination = On average, the mean and standard deviation of the prediction accuracy (n −n * ) 2 is (5.61 × 10 −4 ) ± (6.63 × 10 −4 ). Figure 7 shows an example of the resulting friction cone using the estimated surface normal. We can see that the resulting friction cone is perpendicular to the contact surface, and the resulting contact force falls within this friction cone.

B. Learning surface frictions
In this experiment, we would like to test the accuracy of friction estimation in simulated scenarios. We generated terrains with different friction coefficients µ * ∈ [0.4, 0.8] and used the method in §IV to estimate the friction coefficient.
Overall, the mean and standard deviation of accuracy is (3.08 ± 1.81) × 10 −3 over 5 experiments. Figure 8 shows the resulting friction cones, from the left to right are µ * = 0.4, 0.6, 0.8. We can see that the shape of the friction cone changes as the friction coefficient gets bigger. Figure 9 shows the progression of probability distribution over time. In this example, the true friction coefficient is µ * = 0.5. In Figure 9(a), the x-axis is the friction coefficient, the y-axis is the probability density. The blue curve is the final estimation (at step k), and the red lines are the estimation at formal steps k − 10, k − 20, k − 30. We can see that the posterior distribution of the estimate slowly converges to the unknown mean, and the variance decreases as we receive more observations. Figure 9(b) shows the cumulative distribution of the same example. In this figure, the x-axis is the iteration, the yaxis is the cumulative probability, and the red curve is the cumulative distribution over iterations. The blue line is the cut-off point where the cumulative distribution is good enough, indicates that the robot can stop exploration. Finally, we conduct experiments using ANYmal [20], a torque-controlled quadruped robot made by ANYbotics 1 . The robot weights approximately 35 kg and has 12 joints actuated by Series Elastic Actuators (SEAs). Currently, the soft real-time control cycle is 2.5 ms. The control software is developed based on Robot Operating System (ROS).
The terrain consisted of wedges with random inclination (see Figure 10). Each wedge was about 50 × 50 cm. The robot used the proposed methods to estimate surface normal and friction before moving forward.
We do not have the ground truth of the surface normal, and the ground truth is approximated by measuring the dimensions of the wedges. The accuracy of predicting the surface normal is 8.02 × 10 −3 ± 1.01 × 10 −2 . The friction coefficient is also not well defined, especially, it depends on the materials that interact. Therefore, we cannot evaluate the accuracy of the friction coefficient. The evaluation is based on whether the robot can maintain contact with the terrain after the estimation. Please see the supplementary video https://youtu.be/SnafS_5361g for the result on the real robotic platform.

VI. CONCLUSION
This paper proposes a method for estimating surface information via haptic exploration. We use one of the legs to act as a haptic sensor for exploring the environment before walking. While collecting data, we estimate the surface normal and friction coefficient of the terrain, so the controller can optimize the contact force accordingly. Probabilistic criteria are introduced to evaluate the quality of the estimate. The method was validated on both a simulated environment and the quadruped robot ANYmal.
Future work will focus on incorporating vision data with the haptic exploration (e.g., explore the area where vision data are too noisy to evaluate) and applying this technique for footsteps planning. In addition, one issue we notice during the experiment is that the accuracy of contact detection (i.e., determining whether the robot is in contact or not) affects the accuracy of data collection. We will incorporate a better contact estimation [24][25] with our haptic exploration.
ACKNOWLEDGEMENT This work was funded by the European Commission Horizon 2020 Work Programme: THING ICT-2017-1 780883,