PX4 Firmware
PX4 Autopilot Software http://px4.io
terrain_estimator.cpp
Go to the documentation of this file.
1 /****************************************************************************
2  *
3  * Copyright (c) 2015 Roman Bapst. All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright
10  * notice, this list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright
12  * notice, this list of conditions and the following disclaimer in
13  * the documentation and/or other materials provided with the
14  * distribution.
15  * 3. Neither the name PX4 nor the names of its contributors may be
16  * used to endorse or promote products derived from this software
17  * without specific prior written permission.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
22  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
23  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
24  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
25  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
26  * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
27  * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
28  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
29  * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  *
32  ****************************************************************************/
33 
34 /**
35  * @file terrain_estimator.cpp
36  * A terrain estimation kalman filter.
37  */
38 
39 #include "terrain_estimator.h"
40 #include <lib/ecl/geo/geo.h>
41 #include <px4_platform_common/defines.h>
42 
43 #define DISTANCE_TIMEOUT 100000 // time in usec after which laser is considered dead
44 
46  _distance_last(0.0f),
47  _terrain_valid(false),
48  _time_last_distance(0),
49  _time_last_gps(0)
50 {
51  _x.zero();
52  _u_z = 0.0f;
53  _P.setIdentity();
54 }
55 
57 {
58  return (distance < 40.0f && distance > 0.00001f);
59 }
60 
62  const struct sensor_combined_s *sensor,
63  const struct distance_sensor_s *distance)
64 {
65  matrix::Dcmf R_att = matrix::Quatf(attitude->q);
66  matrix::Vector3f a{sensor->accelerometer_m_s2[0], sensor->accelerometer_m_s2[1], sensor->accelerometer_m_s2[2]};
68  u = R_att * a;
69  _u_z = u(2) + CONSTANTS_ONE_G; // compensate for gravity
70 
71  // dynamics matrix
73  A.setZero();
74  A(0, 1) = 1;
75  A(1, 2) = 1;
76 
77  // input matrix
79  B.setZero();
80  B(1, 0) = 1;
81 
82  // input noise variance
83  float R = 0.135f;
84 
85  // process noise convariance
87  Q(0, 0) = 0;
88  Q(1, 1) = 0;
89 
90  // do prediction
91  matrix::Vector<float, n_x> dx = (A * _x) * dt;
92  dx(1) += B(1, 0) * _u_z * dt;
93 
94  // propagate state and covariance matrix
95  _x += dx;
96  _P += (A * _P + _P * A.transpose() +
97  B * R * B.transpose() + Q) * dt;
98 }
99 
101  const struct distance_sensor_s *distance,
102  const struct vehicle_attitude_s *attitude)
103 {
104  // terrain estimate is invalid if we have range sensor timeout
105  if (time_ref - distance->timestamp > DISTANCE_TIMEOUT) {
106  _terrain_valid = false;
107  }
108 
109  if (distance->timestamp > _time_last_distance) {
110  matrix::Quatf q(attitude->q);
111  matrix::Eulerf euler(q);
112  float d = distance->current_distance;
113 
115  C(0, 0) = -1; // measured altitude,
116 
117  float R = 0.009f;
118 
120  y(0) = d * cosf(euler.phi()) * cosf(euler.theta());
121 
122  // residual
123  matrix::Matrix<float, 1, 1> S_I = (C * _P * C.transpose());
124  S_I(0, 0) += R;
125  S_I = matrix::inv<float, 1> (S_I);
126  matrix::Vector<float, 1> r = y - C * _x;
127 
129 
130  // some sort of outlayer rejection
131  if (fabsf(distance->current_distance - _distance_last) < 1.0f) {
132  _x += K * r;
133  _P -= K * C * _P;
134  }
135 
136  // if the current and the last range measurement are bad then we consider the terrain
137  // estimate to be invalid
139  _terrain_valid = false;
140 
141  } else {
142  _terrain_valid = true;
143  }
144 
145  _time_last_distance = distance->timestamp;
146  _distance_last = distance->current_distance;
147  }
148 
149  if (gps->timestamp > _time_last_gps && gps->fix_type >= 3) {
151  C(0, 1) = 1;
152 
153  float R = 0.056f;
154 
156  y(0) = gps->vel_d_m_s;
157 
158  // residual
159  matrix::Matrix<float, 1, 1> S_I = (C * _P * C.transpose());
160  S_I(0, 0) += R;
161  S_I = matrix::inv<float, 1>(S_I);
162  matrix::Vector<float, 1> r = y - C * _x;
163 
165  _x += K * r;
166  _P -= K * C * _P;
167 
168  _time_last_gps = gps->timestamp;
169  }
170 
171  // reinitialise filter if we find bad data
172  bool reinit = false;
173 
174  for (int i = 0; i < n_x; i++) {
175  if (!PX4_ISFINITE(_x(i))) {
176  reinit = true;
177  }
178  }
179 
180  for (int i = 0; i < n_x; i++) {
181  for (int j = 0; j < n_x; j++) {
182  if (!PX4_ISFINITE(_P(i, j))) {
183  reinit = true;
184  }
185  }
186  }
187 
188  if (reinit) {
189  _x.zero();
190  _P.setZero();
191  _P(0, 0) = _P(1, 1) = _P(2, 2) = 0.1f;
192  }
193 
194 }
Type theta() const
Definition: Euler.hpp:132
bool is_distance_valid(float distance)
matrix::Vector< float, n_x > _x
Definition of geo / math functions to perform geodesic calculations.
#define DISTANCE_TIMEOUT
float accelerometer_m_s2[3]
Type phi() const
Definition: Euler.hpp:128
Quaternion class.
Definition: Dcm.hpp:24
Matrix< Type, N, M > transpose() const
Definition: Matrix.hpp:353
static constexpr float CONSTANTS_ONE_G
Definition: geo.h:51
uint64_t _time_last_distance
void setIdentity()
Definition: Matrix.hpp:447
matrix::Matrix< float, 3, 3 > _P
void measurement_update(uint64_t time_ref, const struct vehicle_gps_position_s *gps, const struct distance_sensor_s *distance, const struct vehicle_attitude_s *attitude)
Vector< float, 6 > f(float t, const Matrix< float, 6, 1 > &, const Matrix< float, 3, 1 > &)
Definition: integration.cpp:8
Euler angles class.
Definition: AxisAngle.hpp:18
void zero()
Definition: Matrix.hpp:421
Quaternion< float > Quatf
Definition: Quaternion.hpp:544
float dt
Definition: px4io.c:73
void setZero()
Definition: Matrix.hpp:416
void predict(float dt, const struct vehicle_attitude_s *attitude, const struct sensor_combined_s *sensor, const struct distance_sensor_s *distance)