A Haptic-Based Proximity Sensing System for Buried Object in Granular Material

The prototype of GRAINS is simple and compact, as explicated in Fig. 1. It is composed of three parts, i.e., the 6-DoF industrial robotic manipulator (UR5, Universal Robots, Denmark), the 6-axis F/T sensor (Force Torque Sensor FT 300, Robotiq, Canada), and a 3D-printed rod with $1$ cm diameter and $14$ cm in length. As proof of concept, the robot arm UR5 used here is only to automatically control the probe motion and easily acquire its position $\mathbf{x}_{t}$. In the future version, the probe can be driven by a motor or other mechanism and even installed on a mobile platform using other positioning devices. Similarly, the current F/T sensor could be replaced with a low-cost tactile sensor in the future, as long as the measurements $\mathcal{T}=\{(\mathbf{x}_{t},\mathbf{f}_{t})\}$ could be obtained simultaneously.

Due to the fan-shaped distribution (from top-view) of the failure wedge zone ahead of the probe, only the buried object directly in front of the probe’s motion can compress GMs, resulting in a jamming state. So, the objects located on two sides of the movement cannot be detected by the probe because they do not enter the failure wedge zone. To enlarge the perception area, we design a novel trajectory to control the probe searching objects in all directions. We call this motion as sprial trajecotry, combining the linear advancement and circular vibration, where the linear motion ensures forward exploration and circular vibration rotates the failure wedge zone. As illustrated in Fig. 3, the spiral trajectory is defined by three parameters, that is, the circular radius (CR), advance velocity (AV), and motion velocity (MV). Here CR (unit: m) refers to the radius of circular movement, and AV (unit: m) indicates the forward distance of the probe after one circular motion. From Fig. 3, it can be found that CR and AV fully determine the path of spiral movement, and MV is a dimensionless variable from UR5 controller, ranging from 0 to 1, presenting the ratio of its end-effector speed.

The path of the spiral trajectory can be formulated as follows. Given start and end positions $\mathbf{x}_{s}$ and $\mathbf{x}_{e}$, the vector from the starting point to goal can be expressed as

The unit vector of $\mathbf{u}$ and its angle could be calculated by $\hat{\mathbf{u}}={\mathbf{u}}/{||\mathbf{u}||_{2}}$ and $\theta_{u}=\arctan(y_{u}/x_{u})$. By definition of the spiral trajectory, the center point $\textbf{x}_{c}$ of the circular motion will move along the direction of $\hat{\mathbf{u}}$. At the same time, the probe will rotate around this moving center. Let’s consider one circular motion in the spiral trajectory for simplicity, as depicted in the enlarged view in Fig. 3. In this case, the probe will advance the distance of AV, which is also the distance of the center $\textbf{x}_{c}$ after one circular motion. Without loss of generality, we suppose the probe has been rotated $\theta\in[0,2\pi]$ around the center. Then the distance of the center that has been moved along $\hat{\mathbf{u}}$ can be calculated by

Thus, the current position of the center is

Then, the current position of the probe $\mathbf{x}\triangleq[x,y]^{\intercal}$ along the spiral trajectory can be formulated by

where $\Theta$ is the angle of vector $\mathbf{v}$ and can be determined according to the geometrical relationship from Fig. 3 as:

If we uniformly discretize $[0,2\pi]$ into $H$ segments, then for each angle, i.e.,

we can generate the pointwise positions $\mathbf{x}_{i}$ along the spiral trajecotry by substituting (9) into (5)-(8), yielding

From (10), it can be found that the path of spiral trajectory only relies on CR and AV. Given the path, we can control the motion velocity of the probe along this path by tuning the dimensionless variable MV.

The goal of the spiral trajectory we design is to enlarge the proximity perception area via rotating the failure wedge zone during the exploration in GMs. The following experiment validates our design. Fig. 4 demonstrates the snapshots of the failure wedge zone ahead of the probe along linear and spiral trajectories, respectively. Note that here we employ the event camera [7], a dynamic vision sensor with high sensitivity to dynamic motions, to record the probe motion along two trajectories. In simple terms, the event camera is capable of capturing the pixels of moving objects, but it cannot detect stationary objects. Consequently, as illustrated in Fig. 4-(a), the region highlighted in front of the probe (outlined by the dashed curve) exhibits a fan-shaped zone, that is, the failure wedge zone. From Fig. 4-(b), the failure wedge zone rotates due to the circular motion involved in spiral trajectory, enabling the probe to detect surrounding objects around $360^{\circ}$.

According to Sec. 4, we know our proximity sensor perceives the nearby object in advance via the variations in the haptic feedback. However, due to the complicated properties of GMs, the analytical interaction force model can not be obtained, compared to the well-studied contact models for rigid objects and liquids. Prior works primarily focus on the interaction force value under different simplifications, such as RFT model [18, 15] or a pre-defined threshold [14], to deal with the complex interaction signals.

In this work, we solve this problem by examining the force pattern rather than focusing on some discrete force values.

Note that $M$ in (11) should be smaller than $L$ in (1). Let $\mathbf{t}=[1,\cdots,M]\in\mathbb{R}^{M}$. We can assume that the force magnitude $f_{t}$ is related to every $t$ as follows

where $p(t)$ is a Gaussian Process (GP) specified by a mean function and a covariance function:

and $\epsilon_{t}$ is independently and identically distributed Gaussian noise:

If defining $\mathbf{f}=\{f_{1},\cdots,f_{M}\}$ and $\epsilon=\{\epsilon_{1},\cdots,\epsilon_{M}\}$, then (12) can be reformulated equivalently as

Given the force pattern $\mathcal{F}$, we can predict the real-value output $f_{*}$ for a new input index $t_{*}$ by the Gaussian regression process (GPR) [22]. Here $\mathcal{F}$ is deemed as the training set. Then, given a set of test inputs $\mathbf{t}_{*}=[M+1,\cdots,M+M_{*}]\in\mathbb{R}^{M_{*}}$, the posterior distribution $p(\mathbf{t}_{*})$ conditioned on $\mathcal{F}$ is a Gaussian distribution expressed as

where $I$ is the identity matrix, and the matrix $K(\mathbf{t},\mathbf{t}_{*})$ is an $M\times M_{*}$ matrix containing covariances evaluated at all pairs of training and test data, and similarly for other matrices $K(\mathbf{t},\mathbf{t})$, $K(\mathbf{t}_{*},\mathbf{t})$, and $K(\mathbf{t}_{*},\mathbf{t}_{*})$.

In this work, due to the spiral trajectory, the force pattern $\mathcal{F}$ exhibits a periodic pattern in the time domain, which can be modeled using the Exp-Sine-Squared (ESS) kernel (aka periodic kernel) [5]:

where $t_{i},t_{j}\in\mathbf{t}$, and $\sigma_{p}^{2}$ is the overall variance, and $l$ is the length scale, and $T$ refers to the periodicity of the kernel. To account for the high randomness from GMs, the white noise kernel [5] is used to estimate the noise of the measured forces:

where $\sigma_{w}^{2}$ is the noise variance and $I$ is the identity matrix. The sum-kernel

is employed as the prior for the covariance function in (13).

Finally, the force magnitude distributions $\mathbf{f}_{*}$ corresponding to $\mathbf{t}_{*}$ can be predicted by mean and covariance functions by evaluating (16). For example, the predicted force magnitude $f_{*}$ at the time $t_{*}\in\mathbf{t}_{*}$ can be expressed by the GP distribution, denoted as:

Therefore, by learning a period of historical force pattern $\mathcal{F}$, we can shortly generate the predicted force pattern distribution $\mathcal{F}_{*}=\{(t_{*},f_{*})|t_{*}\in\mathbf{t}_{*},f_{*}\sim\mathcal{GP}(\mu_{*},\sigma_{*})\}$ over the near future via GPR. Since the $\mathbf{t}$, $\mathbf{t_{*}}$ and $\mathbf{f}$ are all one-dimensional (1D) vectors, the computational burden from GPR is quite small, enabling the proposed method work in a real-time manner.

To find the function $\xi$ in Problem 3.1, we can determine the jamming state by monitoring the anomaly occurring in the force pattern, which is evaluated by the z-score between real force measurement $f_{t}$ and its prediction $f_{*}$ at time $t_{*}$:

When there is no object around the probe, the z-scores of new observations are expected to stay within a high confidence interval (CI). Normally, for the Gaussian distribution, the z-score for 99% CI is 2.576. In contrast, a sudden and sharp increase in z-scores indicates the anomaly in the force pattern, probably resulting from a subterranean object in the vicinity of the probe. We online monitor the changes in the z-score, and once a certain threshold $\bar{z}$ is exceeded, our system could release a proximity warning, followed by a stop command sent to the robot arm. Thus, we can mathematically express our haptic-based proximity sensing system as the solution of Problem 3.1 as follows:

It is obvious that the selection of the z-score threshold $\bar{z}$ in (22) will be vital since it directly determines our tolerance for outliers and, in turn, when to stop moving to avoid the potential collision. In the following, we will introduce how to adaptively select system parameters (including $\bar{z}$) for different particles through the proposed parameter calibration experiments.

To adapt to different granular scenarios, we have developed an automated parameter calibration process, which can autonomously select the optimal parameter set $\mathcal{D}$ for the current GM, including $\mathcal{D}=\{\text{MV}^{*},T^{*},\bar{z}\}$. Recall that the variable MV is a dimensionless parameter, ranging from $0$ to $1$, to control the speed of the end-effector of the robot arm, as introduced in Sec. 5.2. The hyperparameter $T$ used in (17) represents the periodicity of the force pattern $\mathcal{F}$. Even if $T$ will be tuned during the GPR learning, an optimal value $T^{*}$ as the initial value for the online learning will reduce the computational cost. Most importantly, the periodicity $T$ only relies on the parameters of spiral trajectory, i.e., CR, AV, MV, and be determined by the following equation:

where $V_{max}$ is the max speed of the end-effector and $f_{s}$ is the sampling frequency of F/T sensor. In this work, $V_{max}=0.08968$ m/s and $f_{s}=62.5$ Hz according to the setup used. Also, $l_{p}$ refers to the length of the spiral path for one circular motion, as in Fig. 3. It can be calculated by the sum of segment-wise paths from discretized path points, i.e.,

Given CR and AV, we can obtain the $l_{p}$ by substituting $\mathbf{x}_{i},\mathbf{x}_{i+1}$ from (10) into (24). Then the optimal periodicity $T^{*}$ can be determined in advance by submitting $l_{p}$ and the given MV^∗ into (23). In addition, from our preliminary experiments, we find that different types of GMs display varying z-score thresholds to identify the jamming state. To increase credibility and reduce the ratio of false positive cases in various GMs, the GM-specific parameter $\bar{z}$ should be calibrated in advance as well.

To do so, we take samples from current granules and store them in a container. Our GRAINS then controls the probe raking in this container along a spiral trajectory with given CR and AV, as well as MV_i from the set $\{$MV${}_{i}\}$. For each MV_i, GRAINS calculates the periodicity $T_{i}$ by (23). Then we can formulate an ESS kernel by substituting $T_{i}$ in (17). By using the sum-kernel (19) from (17) and (18), we can initialize a GP model from (13), denoted by $GP_{i}$. Note that other parameters in (17) and (18) will be given fixed values, and they will be tuned by GPR learning. At each MV_i, the GRAINS controls the probe to explore $20$ cm in the container. After that, GRAINS successively divides the online measurements $\mathcal{T}$ into several force patterns and feeds the $k$-th force pattern $\mathcal{F}_{k}$ into $GP_{i}$. Based on the predicted mean and standard deviation from $GP_{i}$, GRAINS computes the z-scores of the $(k+1)$-th force pattern $\mathcal{F}_{k+1}$. This process is repeated to obtain z-scores for all force patterns except the first one. GRAINS then calculates the Root Mean Squared Error (RMSE) of z-scores at each MV_i, denoted as $\text{RMSE}_{i}$. Finally, GRAINS determines the optimal MV^∗ corresponding to the minimum RMSE values. The periodicity prior $T^{*}$ at MV^∗ is then chosen as the initial hyperparameter for (17). In addition, since no object is located under GMs, the maximum z-score at MV^∗ during this process is considered as the threshold $\bar{z}$ for (22).

So, it can be seen that the above calibration process is all automated, and optimal (hyper)parameters in $\mathcal{D}$ are all determined autonomously by the proposed system GRAINS according to the type of granules, without human factors.

First of all, we demonstrate the near-field perception of the proposed proximity sensing system GRAINS in sands, as shown in Fig. 5. Here the sand is composed of silica with diameters ranging from $0.5$ to $1.6$ mm.

To determine the optimal (hyper)parameters for sands, we set the speed MV${}_{i}=\{0.2,0.3,\cdots,0.7\}$ and generate the spiral path by CR $=0.02$ m and AV $=0.01$ m. These settings are based on our experimental observations. If MV is too small, low search efficiency is not acceptable. Also, MV cannot be too large; otherwise, it is too late for an emergency stop. Large CR and AV values may cause the exploration to exceed the existing sandbox space. After the automatic calibration, the resulting RMSE values corresponding to MV_i are revealed in Table 1. Thus, we obtain the optimal (hyper)parameters as $\mathcal{D}=\{\text{MV}^{*}:0.2,T^{*}:439,\bar{z}:3.9\}$. After that, we validate the proximity sensing performance in sands, where a wooden cylinder is buried, as depicted in Fig. 5-(f). From the depth assumption in Sec. 3, the penetration depth of the probe in the experiment is fixed at $4$ cm, and the underground objects are located at a relatively shallow depth, less than $4$ cm, to ensure the overlap in the depth direction between the probe and the buried objects. Before the real experiment, the F/T sensor will be zeroed out to prevent drift. During the experiment, the force sensor will use a low-pass filter to filter the original data to reduce the influence of noise. Then GRAINS monitors the filtered force measurements $f_{t}$ in a real-time manner and employs the sliding window to clip the force pattern $\mathcal{F}$ with $M=2000$ and predict the force pattern distribution $\mathcal{F}_{*}$ over the next $M_{*}=1000$ by (16) with pre-calibrated optimal parameters $\mathcal{D}$. The z-score at $t_{*}\in\mathbf{t}_{*}$ is constantly computed from (21) once the force measurement $f_{t}|_{t=t_{*}}$ is available.

The results of proximity sensing from GRAINS are illustrated in Fig. 5-(f). The measurements $\mathcal{T}$ during the spiral trajectory can shown in Fig. 5-(b). The window sliding operation is shown in Fig. 5-(c). A window of $M=2000$ iterations is sliding along $\mathcal{T}$ to form the $\mathcal{F}$ (green curves in Fig. 5-(d)), and z-scores $\mathbf{z}$ of the next $M_{*}=1000$ iterations are calculated accordingly, as illustrated in the insets of Fig. 5-(d), by the real measurements $f_{t}|_{t\in\mathbf{t}_{*}}$ (pink curves in Fig. 5-(d)) and predicted distribution $\mathcal{F}_{*}$ (blue dashed curve and blue shaded area in Fig. 5-(d)). It can be seen that when the probe is far away from the object, e.g., $0$-th eps and $3$-th eps in Fig. 5-(d), the predicted force pattern is consistent with the real measurements, as $f_{t}|_{t=t_{*}}$ is within CI $99\%$, or equivalently $z_{i}<\bar{z}=3.9,\ \forall z_{i}\in\mathbf{z}$. So, by definition of (22), GRAINS safely outputs $\xi(\mathcal{T})=0$. When the probe approaches the object, the failure wedge zone contacts the object before the probe as expected, and granule jamming occurs. Thus, the force pattern induced by the jamming state starts to diverge from the high CI area with the growth of z-scores. Once a z-score value $z_{L}$ exceeds the predetermined threshold $\bar{z}$ at time $t_{L}$ during the $6$-th eps as shown in Fig. 5-(d), GRAINS identifies the presence of subsurface stuff nearby, and outputs $\xi(\mathcal{T}|_{t=t_{L}})=1$. Therefore GRAINS sends a cease command to stop the robot arm at the current position $\mathbf{x}_{t}|_{t=t_{L}}$ to avoid the potential collision accordingly. From Fig. 5-(f), the GRAINS perceives the presence of near objects beneath sands at a distance of $2.1$ cm in advance. From (3), we can say the proximity sensing range for this case is $\zeta(\xi)=2.1$ cm.

We compare our method with the baseline approach, which uses the pre-defined threshold for the force value, as in [14]. In this case, we set the force bar as $15$ N, as shown in the dashed line in Fig. 5-(a), considering the maximum static friction force from still particles that the probe should overcome. From Fig. 5-(e), we can observe that the probe fails to detect the proximity of underground objects and breaks into several pieces finally. In addition, if the force threshold is lower than $15$ N, then the probe could stuck at the initial stage since this value should be larger than the maximum static friction force, as depicted in Fig. 5-(a). Therefore, the baseline method requires great care to determine the force threshold. And, for different granules, these human efforts need to be repeated. However, our method can autonomously update parameters to adapt to different GMs, as evaluated in the following.

We test our GRAINS in other granules, as shown in Fig. 6-(b), including cassia seed (about $3.5$ mm in length and $2.5$ mm in width), cat litter (about $3\sim 4$ mm in diameter), and soybean (about $6.5$ mm in diameter). Similar to the calibration process for the sand, the optimal parameters for these GMs are determined by GRAINS autonomously and presented in Table 2. Then we conduct proximity perception in each GM for $20$ times and resulting proximity sensing ranges $\zeta(\xi)$ defined by (3) are recorded. From Fig. 6-(a), we can find that the median value of $\zeta(\xi)$ in cat litter ($\sim 4.5$ cm) is similar to that in the sand ($\sim 4.2$ cm) but with a smaller dispersion. This discrepancy can potentially be attributed to the notably larger particulate size of cat litter, which leads to increased instability in the conduction of force chains compared to sand. Since a rough granule surface is beneficial for forming a highly compressed state, the smoother surface of cassia seed requires the probe to get closer to objects to sense force anomaly signals from granular jamming, resulting in a shorter $\zeta(\xi)$ ($\sim 2.7$ cm), even if cassia seed has a similar grain size to cat litter. Soybean has the shortest $\zeta(\xi)$ ($\sim 1.3$ cm) due to its largest particle size and smoothest surface.

In addition, we test the effects of non-optimal parameters on $\zeta(\xi)$. As shown in Fig. 6-(a) about the cat litter, the probe stops noticeably closer to the object with a non-optimal MV$=0.5$ than that with the optimal MV${}^{*}=0.3$, i.e., $0.6$ cm and $4.5$ cm, respectively. Furthermore, several outliers can be identified in the case of MV $0.5$. This highlights the necessity and importance of parameter calibration in GRAINS.