The present invention relates to the field of positioning of equipment and, in particular, discloses a system for determining the spatial pose of swing loading equipment utilised in mining operations such as electric mining shovels.
Various solutions to the problem of determining the position and orientation of mobile equipment units have been previously proposed. Solutions invariably take advantage of a variety of localization sensors, including some of those given above.
Duddek et al. (1992) discloses a method of determining the position and orientation of the end of a excavator bucket utilising GPS sensors and a receiver in the vicinity of the bucket wheel.
Kalafut et al. (2002) proposes a system by which the position and heading of a machine can be determined through the use of a single positioning sensor. Readings are taken from the positioning sensor over time, and a motion profile is generated to estimate the heading of the machine. This approach is particularly applicable to machines that are commonly in motion, and have well-defined dynamic characteristics. In a mining application, haul trucks are a good candidate for this type of approach, so long as they are in motion.
Another example of a single-sensor positioning system is that proposed by Sahm et al. (1995) which uses a single sensor, capable of collecting (x, y, z)-position measurements connected to the boom of a mining shovel. If the shovel's undercarriage is assumed to be stationary during a dig cycle, then a set of points can be measured over time to generate the plane in which the sensor exists. This estimate, along with the current measurement of position from the sensor, can be used to estimate the current position of the shovel bucket.
The method of localizing from an estimated plane is further explored by Dizchavez (2001). Two GPS antennas are mounted on the machine house at known locations of equal elevation. During the operation of the machine, rotation of the house can be measured, and using calculations based on standard-deviation analysis, an estimate of the plane in which the two antennas lie is formed. From this plane, and the current position and orientation of the sensors within the plane, another part of the machine can be localized given a kinematic model and appropriate joint position information.
It is desirable to provide an improved method and apparatus for determining the spatial pose of mining equipment or the like.
Any discussion of the prior art throughout the specification should in no way be considered as an admission that such prior art is widely known or forms part of common general knowledge in the field.
Unless the context clearly requires otherwise, throughout the description and the claims, the words “comprise”, “comprising”, and the like are to be construed in an inclusive sense as opposed to an exclusive or exhaustive sense; that is to say, in the sense of “including, but not limited to”.
In accordance with a first aspect of the present invention, there is provided a method of determining the global pose of a mining shovel, the method including the step of applying a multi stage calculation, including: (a) as a first stage computing the location of the mining shovel carbody (c-frame) relative to a local geodetic frame (g-frame) utilising a global positioning system, an inclinometer, and a swing axis resolver; (b) as a second stage computing a house pose (h-frame) relative to the c-frame using a global positioning system, an axis inertial sensor and a swing axis resolver. (c) as a third stage computing a bucket pose (b-frame) relative to the h-frame using crowd and hoist axis resolvers.
The steps (a) and (b) are preferably carried out utilising an extended Kalman filter. The step (a) can be carried out utilising an iterative routine until convergence. The inclinometer can be a twin axis inclinometer. The inertial sensor can be a six axis inertial sensor. The first portion of the shovel can comprise the machine house.
In accordance with a further aspect of the present invention, there is provided a method determining the global pose of electric mining shovels as a three stage calculation process which: (a) at a first stage computes the location of the carbody (c-frame) relative to a local geodetic frame (g-frame) utilising a global positioning system, an dual axis inclinometer, and a swing axis resolver until convergence; (b) at a second stage computes the house pose (h-frame) relative to the c-frame using a global positioning system, a six axis inertial sensor (three rate gyroscopes and three linear accelerations) and a swing axis resolver; (c) at a third stage computes the bucket pose (b-frame) relative to the h-frame using crowd and hoist axis resolvers.
In accordance with another aspect of the present invention, there is provided a method for determining the global spatial pose of a mining shovel, the method comprising the steps of (a) designating a first Earth-Centred-Earth-Fixed (ECEF) frame or e-frame of reference; (b) designating a local geodetic coordinate frame, denoted a g-frame, in the vicinity of the mining shovel, defined as a set of Cartesian coordinate axes in the e-frame; (c) designating a set of Cartesian coordinate axes, denoted a c-frame, in the close vicinity to the carbody or under-carriage of the mining shovel; (d) determining the location of the c-frame within the g-frame; (e) designating a set of Cartesian coordinate axes, denoted a h-frame, in the vicinity of the machine house of the mining shovel; (1) determining the location of the h-frame within the c-frame; (g) designation a set of Cartesian coordinate axes, denoted the b-frame, fixed to the close vicinity of the shovel handle and dipper assembly; and (h) determining the location of the b-frame within the h-frame.
Preferred forms of the present invention will now be described with reference to the accompanying drawings in which:
As illustrated in
As shown in
The preferred embodiment exploits such operational characteristics of mining shovels 1 to address the problem of determining the pose of the shovel.
The preferred embodiment also exploits the combinations of several available complimentary sensor measurements, including
The preferred embodiment presents a formulation of a recursive algorithm based on the extended Kalman filter that determines the global shovel pose using combinations of these measurements.
Knowing shovel pose in real-time is useful for several purposes, which include
The prior art solutions ignore the estimation theory that can be adapted to the estimation problem. Specifically, it is possible to formulate the problem as a state estimation exercise, in which the relative positions and orientations of the g-, c-, h- and b-frames can be expressed as the states of a dynamic system and knowledge of the causal relationship (the ‘process model’) between the measured operator command references and the resulting motion of the machine is used to propagate current knowledge of shovel pose (in the form of a probability distribution) forward in time to facilitate fusion with a combination of measurements from the sensors identified previously.
The geometry relevant to the problem including the various coordinate frames are shown in
The pose of the shovel bucket is computed in two stages.
The objective in the first stage is to compute the location of the c-frame relative to the h-frame from the following measurements
The objective in the second stage is, having found the location of the c-frame relative to the g-frame, to compute the location of the h-frame relative to the c-frame using the following measurements
The objective in the third stage is to compute the location of the b-frame relative to the h-frame using the following measurements
The calculation process to determine shovel pose is at follows
Underpinning this staged calculation process is the idea that the measurements used at Stage 1 can provide rich information about the low frequency motions of the machine sufficient to accurately determine the position of the c-frame relative to the g-frame. During the normal to-and-fro motion associated with normal production, in addition to large scale swing motion, there is the potential for the machine house to undertake smaller amplitude rocking motions, particularly during digging. The sensor measurements used at Stage 2 aim towards accurate determination of these motions. In this sense, the Stage 2 filter aims towards a higher bandwidth of estimation.
The methodology for the calculation of Stages 1 and 2 is the extended Kalman filter (EKF), Gelb (1974), The EKF requires a system model of the form
{dot over (x)}=f(x,u,t)+w,w˜N(0,Q)
z
k
=h(xk,uk,k)+vk,vk˜N(0,R) (1)
where f(x,u,t) is a vector-valued function describing the dynamics of the system that is used to propagate the current estimate of state and state covariance forward in time based on measurement of the operator command reference, so it can be combined with newly obtained measurement data. The vector-valued function h(xk,uk,k) expresses the measurements in terms of the state vector x and inputs u.
The EKF requires linearization of f(x,u,t) and h(xk,uk,k) about the estimated state trajectory {circumflex over (x)} and the conversion of the linearized continuous dynamics to a discrete time form. It is desirable to use the following notation:
where Δt is the measurement update rate and
The vectors w and v in Equations 1 are termed the process and measurement noise and are assumed to be generated by zero mean, Gaussian processes with covariances Q and R respectively.
The computational scheme of the EKF involves the following five steps (Gelb, 1974):
1. State estimate propagation
xk+1−=F{circumflex over (x)}k+ (7)
2. State covariance propagation
P
k+1
−
=FP
k
+
+P
k
+
F
T
+Q (8)
3. Calculation of the Kalman gain
K
k+1
=P
k+1
−
∇h
k+1
T(∇hk+1Pk+1−∇Vhk+1T+R)−1 (9)
4. State estimate update
{circumflex over (x)}
k+1
+
={circumflex over (x)}
k+1
−
+K
k+1(zk+1−h({circumflex over (x)}k+1−,k)) (10)
5. Error covariance update at time step k
P
k+1
+=(I−Kk+1∇hk+1)Pk+1− (11)
Equations 7 to 11 define the EKF algorithm which provides the best linear state estimator for a non-linear system measured by the minimum mean squared error. The superscripts ‘−’ and ‘+’ in Equations 7 to 11 indicate evaluation of quantities before and after a measurement has been made.
The dynamic model used for propagation of shovel pose in Stages 1 and 2 includes, as a common element, a causal model relating operator joystick reference to swing motions within the vector-valued function f (x,u,t). The preferred embodiment of this model for Centurion enables P&H shovels is given below.
P&H Centurion enabled electric mining shovels use an ABB DCS/DCF600
Multi-Drive controller to regulate motor speed, armature current and field current in each of the three DC motors. The controller is made up of four integral components; a PID or PI motor speed control loop, an armature current saturation limiter, a PI current control loop and an EMF-field current regulator.
The swing drive uses a combination of torque control and bang-bang speed control, where by the swing joystick position generates a piecewise speed reference and an armature current saturation limit. A schematic of the swing drive model is shown in
Modelling the shovel's drives effectively requires a means for incorporating the non-linear saturation effects seen in the motor armature currents. To include these effects into the prediction models a sinusoidal input describing function is used. The describing function, which will be abbreviated to DF, was developed primarily for the study of limit cycles in non-linear dynamic systems, see Gelb and Vander Velde (1968), and Graham and McRuer (1961). The basic idea of the describing function approach is to replace each non-linear element in a dynamic system with a quasi-linear descriptor or describing function equivalent gain whose magnitude is a function of the input amplitude.
Armature current saturation is modelled using a DF equivalent gain for sinusoidal saturation.
where b, is the first or fundamental term of the output's Fourier series and A is the amplitude of the input.
Implementing the DF into the prediction models requires the equivalent gain to be evaluated at each time step. If the input to the current saturation block is larger than the saturation limit, the saturated output is calculated by multiplying input magnitude by the equivalent gain.
The assumption that the higher order terms of the saturated output's Fourier transform are not dominant can be supported by the notion that the shovel's drive dynamics act as low pass filters and the fundamental frequency of the output is much less attenuated in passing through the system than the output's higher order harmonics.
The drive prediction models are presented as continuous, linear state space systems with the form
{dot over (x)}=Ax+Bu.
The input vector u, contains the reference motor speeds ωsd generated from the joystick signals, the static torque load on the motor due to gravitational effects Ts and a coulomb friction disturbance input fs. The state vector x, for each model, contains swing armature current Is, the swing motor speed ωs, the swing motor position θs, and the integrals of the error in the speed and current controllers, ∫eω
The swing drive model also contains the swing reference armature current prior to the saturation limit Is
The describing function gains Gs appears as an element in drive system and input matrices which is recalculated at each time step. The input to the current limiter is the swing reference armature current state Is
The swing armature current saturation limits can be determined from the swing joystick position.
It is to be noted that to be effective the Stage 1 and Stage 2 models should contain so-called ‘shaping stages’ that accommodate bias in sensors and account for physical artefacts such as transmission backlash.
It is to be noted that the Stage 2 model can also be implemented by the use of so-called Ornstein-Uhlenbeck stochastic process whose parameters may be determined from subsequent autocorrelation analysis.
Given the position (xa,ya,za) of the nth GPS antenna in the h-frame, Eqn. 12 expresses the GPS measurement made in the g-frame in terms of the position of the shovel carbody (xc,yc,zc) and the direction cosine matrices Rc2g and Rh2c, describing the 3D rotations between the c- and g-, and the h- and c-frames, respectively. These matrices can be calculated in a number of ways, e.g. Euler angles or quaternions. The parameters describing these matrices are states of the estimator.
z
gps(xc,yc,zc)T+Rc2gRh2c(xa,ya,za)T (12)
GPS measurements are based around an approximation of the Earth's surface in the form of a bi-axial ellipsoid. The dimensions of this ellipsoid are defined by one of several standards or datums.
p
0
=N
0 cos φ0
r
b
=r
G
×r
a (16)
p
n
local
=T·(pnglobal−r0G). (18)
The pitch and roll inclinations of the machine house, as measured by the inclinometer, are determined as the angle between unit vectors in the xh- and yh-axes and the (xg, yg) plane. Expressing the rotation of the c-frame relative to the g-frame as the Euler angles (φc,θc,φc), using the standard Z-X-Z rotation convention, the pitch (α) and roll (β) angles are
α=sin−1(sin θc sin φc) (19)
β=sin−1(sin θc cos φc) (20)
Raw inclinometer measurements may be improved by factoring away accelerations of the sensor.
The acceleration of the IMU will be measured as the global acceleration of its location in the machine house, (xi,yi,zi), rotated to be aligned with the orthogonal sensor axes. Given the position of the shovel carbody (xc,yc,zc) and the direction cosine matrices Rc2g and Rh2c, describing the 3D rotations between the c- and g-, and the h- and c-frames, respectively, the acceleration measurements are
zacc=Rg2cRc2hai (21)
Where the acceleration of the IMU in the g-frame is found from
The g-frame is assumed to be non-accelerating and non-rotating.
The angular velocity of the IMU will be measured as the global angular velocity of the machine house, rotated to be aligned with the orthogonal sensor axes. Expressing the rotation of the h-frame relative to the c-frame as the Euler angles (φh,θh,φh), using the standard Z-X-Z rotation convention, the measured angular velocities about the RPY axes are
The Stage 3 calculations are dependent on the machine in question. In the preferred embodiment the calculations have been carried out for a P&H-class mining shovels
The c-frame is denoted Ocxcyczc; the h-frame is denote Ohxhyhzh and is embedded in the machine house; the Omxmymzm-m-frame in the saddle; and the Obxbybzb b-frame in the dipper. The x- and z-axes of all the body-fixed frames are in the sagittal plane of the machine house, that is the plane parallel to the plane of projection shown in
Four-by-four homogeneous transformation matrices can be used to describe the relationship between frames. We denote the matrix describing the transform from Oixiyizi to Oixiyizi by Di→j and note that the action of this matrix maps (homogeneous) points in the j-frame to the i-frame. For instance, if p is a point whose location is known in Obxbybzb, fixed in the bucket, then the coordinates of that point in frame Ocxcycz, can be found from.
p′=D0→3p.
The structure of Di→j is
where Ri→j is a 3×3 rotation matrix and ti→j is a 3-dimensional translation vector. Four-by-four homogeneous transformation matrices commute according to
Di→k=Di→jDj→k.
The origin of the c-frame is located at the interface between the upper surface of the tracks and the underside surface of the machine house. The zc-axis collinear with the swing axis. The xc-axis points in the direction of forward travel of the crawler tracks, and the yc-axis completes a right-handed trihedral coordinate frame.
The origin Oh of frame Ohxhyhzh is coincident with Oc and zh is collinear with zc. When θ1=0, frames Ocxcyczc and Ohxhyhzh coincide. A positive angle θ1 corresponds to an anti-clockwise rotation of the machine house relative to the tracks when viewed from above. The homogeneous transformation matrix Dc→h is given by
Frame Omxmymzm is fixed to the saddle with O2 at the center of rotation of the saddle. When θ2 is equal to 0, the coordinate directions of Omxmymzm are parallel to those of Ohxhyhzh. The displacement matrix describing the rigid body displacement from Frame h to Frame n is given by
where design parameters l1 and φ1 are as shown in
The origin Ob of Obxbybzb is located as follows. The saddle angle (θ2) is set equal to 90 degrees so that the handle is horizontal. The handle is then displaced such that the hoist rope falls vertically (θ5=90 degrees and θ6=0 degrees). The origin Ob is located at the intersection of the pitch-line of the handle-rack and the hoist rope; zb is set to be collinear with the axis of the hoist rope; xb is set parallel to the pitch-line of the handle-rack. Note that axis xm is orthogonal to axis xb. The displacement matrix describing the rigid body displacement Dm→b is given by
The multiplication of the above Eqns gives:
The swing angle, θ1, the pivot angle, θ2, and the crowd extension, d3, parameterize the displacement and rotation of the body fixed frames relative to the world frame. These configuration variables can be grouped as follows:
θ=(θ1,θ2,d3)T.
The displacements of the swing motor θs, crowd motor θc, and hoist motor θh can be similarly grouped:
ψ=(θs,θc,θh)T.
The values of θ determine ψ and vice versa. These mappings are not bijective. However, within the physical working range of these variables their correspondence is one-to-one. Note that the specification of either θ or ψ determines the inclination of the hoist rope, a seventh variable, labeled θ5 in
To build up the constraint equations, start by noting the dependent coordinates θs and θc are related θ1 and d3 by transmission ratios leading to
where Gs is the transmission ratio of the swing drive and Gc is the transmission ratio of the crowd drive.
The constraint equations relating θ2 and θ5 to θc and θh can be developed using the vector loop shown in
Summing vectors in the vector loop of
0=ζ(θ,ψ,θ5)=z1+z2+z3+z4+z5−z6−z7,
which can be expanded to give,
0=ζ(θ,ψ,θ5)=l1eiφ
where the variables li, di, θi and φi are as defined in
Eqn. 24 can be written as,
ζ(θ,ψ,θ5)=l1eiφ
To remove d5 from this equation it is convenient to first introduce the variable dh. As shown in
where Gh is the hoist transmission ratio. The above expression and Equation 25 can be used to relate d5 and dh,
The last term on the right hand side of Eqn. 27 accounts for the angle of wrap of the hoist rope around the sheave. Substituting Eqn. 27 into Eqn. 26 gives
Taking the real and imaginary components of Eqn. 28 relates θh and θ5 to the generalized coordinates,
Note that use has been made of the following trigonometric relationships in arriving at Eqn. 29,
Concatenating Eqns. 23a and 29 gives
The kinematic tracking problem is to determine the values of θ and θ5 given ψ or to determine ψ and θ5 given ψ. We call the first problem forward kinematic tracking and the second inverse kinematic tracking. For notational convenience in distinguishing between these two problems, we write
θF=(θ1,θ2,θ3,θ5)
ψF=(θs,θc,θh)
when working in the domain of the forward kinematic tracking problem and
θl=(θ1,θ2,θ3)
ψl=(θs,θc,θh,θ5)
The difference here is in with the grouping of θ5. Both problems amount mathematically to solving the non-linear constraint equations. We chose to do this iteratively using a multi-variable Newton's method. The Jacobian matrices developed in expressing a Newton's method solution are used to refer motor inertias to the configuration variables and to solve the statics problem.
Applying Taylor's series expansion to Eqn 30
The objective is to find a valid configuration, i.e. Γ(θF+ΔθF,ψF+ΔψF)=0. It follows that
This leads to the iteration equation
The algorithm below, gives an algorithm for kinematic tracking for P&H-class shovels. The algorithm takes the current motor positions, ψFk and uses Eqns. 31 and 32 to find the new values of θFk consistent with the constraint equations. For reliable convergence the algorithm requires a good initial values θF0. In practice this can be achieved by initializing from a well defined configuration such as provided in
It can be seen that the preferred embodiments provide an accurate method for maintaining a close approximation of the shovel position at all times.
Although the invention has been described with reference to specific examples it will be appreciated by those skilled in the art that the invention may be embodied in many other forms.
Number | Date | Country | Kind |
---|---|---|---|
2008900081 | Jan 2008 | AU | national |
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/AU09/00019 | 1/7/2009 | WO | 00 | 7/8/2010 |