基于openmv的卡尔曼滤波(简化版)

一些函数的定义 


import sensor, image, time,utime
from pyb import RTC


class Vector:
    def __init__(self, lis):
        self._values = lis  # 在vector中设置私有变量values存储数组数据

    # 返回 dim 维零向量
    @classmethod
    def zero(cls, dim):
        return cls([0] * dim)

    # 返回向量的模
    def norm(self):
        return math.sqrt(sum(e ** 2 for e in self))

    # 返回向量的单位向量
    def normalize(self):
        if self.norm() < 1e-8:
            raise ZeroDivisionError("Normalize error! norm is zero.")
        return 1 / self.norm() * Vector(self._values)

    # 向量相加,返回结果向量
    def __add__(self, other):
        assert len(self) == len(other), "error in adding,length of vectors must be same"
        return Vector([a + b for a, b in zip(self, other)])

    # 向量相减,返回结果向量
    def __sub__(self, other):
        assert len(self) == len(other), "error in subbing,length of vectors must be same"
        return Vector([a - b for a, b in zip(self, other)])

    # 向量点乘,返回结果向量
    def dot(self, other):
        assert len(self) == len(other), \
            "error in doting,length of vectors must be same."
        return sum(a * b for a, b in zip(self, other))

    # 向量左乘标量,返回结果向量
    def __mul__(self, k):
        return Vector([a * k for a in self])

    # 向量右乘标量,返回结果向量
    def __rmul__(self, k):
        return Vector([a * k for a in self])

    # 向量数量除法,返回结果向量
    def __truediv__(self, k):
        return 1 / k * self

    # 向量取正
    def __pos__(self):
        return 1 * self

    # 向量取负
    def __neg__(self):
        return -1 * self

    # 取出向量的第index个元素,调用时直接vec[]
    def __getitem__(self, index):
        return self._values[index]

    # 返回向量长度,调用时直接len(vec)
    def __len__(self):
        return len(self._values)

    # 返回向量Vector(...)
    def __repr__(self):  # __repr__和__str__都在调用类时自动执行其中一个,倾向位置在最后一个
        return "Vector({})".format(self._values)

    # 返回向量(...),调用时直接print(vec)
    def __str__(self):
        return "({})".format(",".join(str(e) for e in self._values))


class Matrix:

    def __init__(self, list2d):
        self._values = [row[:] for row in list2d]

    @classmethod
    def zero(cls, r, c):
        """返回一个r行c列的零矩阵"""
        return cls([[0] * c for _ in range(r)])

    def T(self):
        """返回矩阵的转置矩阵"""
        return Matrix([[e for e in self.col_vector(i)]
                       for i in range(self.col_num())])

    def __add__(self, another):
        """返回两个矩阵的加法结果"""
        assert self.shape() == another.shape(), \
            "Error in adding. Shape of matrix must be same."
        return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))]
                       for i in range(self.row_num())])

    def __sub__(self, another):
        """返回两个矩阵的减法结果"""
        assert self.shape() == another.shape(), \
            "Error in subtracting. Shape of matrix must be same."
        return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))]
                       for i in range(self.row_num())])

    def dot(self, another):
        """返回矩阵乘法的结果"""
        if isinstance(another, Vector):
            # 矩阵和向量的乘法
            assert self.col_num() == len(another), \
                "Error in Matrix-Vector Multiplication."
            return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())])

        if isinstance(another, Matrix):
            # 矩阵和矩阵的乘法
            assert self.col_num() == another.row_num(), \
                "Error in Matrix-Matrix Multiplication."
            return Matrix([[self.row_vector(i).dot(another.col_vector(j)) for j in range(another.col_num())]
                           for i in range(self.row_num())])

    def __mul__(self, k):
        """返回矩阵的数量乘结果: self * k"""
        return Matrix([[e * k for e in self.row_vector(i)]
                       for i in range(self.row_num())])

    def __rmul__(self, k):
        """返回矩阵的数量乘结果: k * self"""
        return self * k

    def __truediv__(self, k):
        """返回数量除法的结果矩阵:self / k"""
        return (1 / k) * self

    def __pos__(self):
        """返回矩阵取正的结果"""
        return 1 * self

    def __neg__(self):
        """返回矩阵取负的结果"""
        return -1 * self

    def row_vector(self, index):
        """返回矩阵的第index个行向量"""
        return Vector(self._values[index])

    def col_vector(self, index):
        """返回矩阵的第index个列向量"""
        return Vector([row[index] for row in self._values])

    def __getitem__(self, pos):
        """返回矩阵pos位置的元素"""
        r, c = pos
        return self._values[r][c]

    def ni_matrix(self):
        """返回矩阵的逆(二维矩阵)"""
        a = self._values[0][0]
        b = self._values[0][1]
        c = self._values[1][0]
        d = self._values[1][1]
        ni1 = Matrix(([d, -b], [-c, a]))
        k = 1 / (a * d - b * c)
        return Matrix([[e * k for e in ni1.row_vector(i)]
                       for i in range(ni1.row_num())])

    def size(self):
        """返回矩阵的元素个数"""
        r, c = self.shape()
        return r * c

    def row_num(self):
        """返回矩阵的行数"""
        return self.shape()[0]

    __len__ = row_num

    def col_num(self):
        """返回矩阵的列数"""
        return self.shape()[1]

    def shape(self):
        """返回矩阵的形状: (行数, 列数)"""
        return len(self._values), len(self._values[0])

    def __repr__(self):
        return "Matrix({})".format(self._values)

    __str__ = __repr__

 初始化


# ________________________________Apriltags相关系数___________________________________________

f_x = (2.8 / 3.984) * 160  # find_apriltags defaults to this if not set
f_y = (2.8 / 2.952) * 120  # find_apriltags defaults to this if not set
c_x = 160 * 0.5  # find_apriltags defaults to this if not set (the image.w * 0.5)
c_y = 120 * 0.5  # find_apriltags defaults to this if not set (the image.h * 0.5)
K_x = 23
# ________________________________卡尔曼滤波器模型___________________________________________
##预测
# pre_Pk=(A.dot(Pk_1)).dot(A.T())+Q
# pre_Xkhat = A.dot(pre_Xk_1hat)
# pre_Xk_1hat = pre_Xkhat

##校正
# Xk_hat = pre_Xk + Kk.dot(Zk - pre_Xk)
# Kk=pre_Pk.dot((pre_Pk_1+R).ni_matrix())
# Pk = (I -Kk).dot(pre_Pk)
# ________________________________初始化噪声、误差协方差矩阵及相关变量___________________________________________

q1 = 1
q2 = 1
r1 = 1
r2 = 1
p1 = 1
p2 = 1

T = 0
A = Matrix([[1, T], [0, 1]])
H = Matrix([[1, 0], [0, 1]])
pre_xpk = Matrix([[0, 0], [0, 0]])
est_xpk_1 = Matrix([[0, 0], [0, 0]])
pre_xk = Matrix([[0], [0]])
est_xk = Matrix([[0], [0]])
est_xk_1 = Matrix([[0], [1]])
est_xpk = Matrix([[p1, 0], [0, p2]])

pre_ypk = Matrix([[0, 0], [0, 0]])
est_ypk_1 = Matrix([[0, 0], [0, 0]])
pre_yk = Matrix([[0], [0]])
est_yk = Matrix([[0], [0]])
est_yk_1 = Matrix([[0], [1]])
est_ypk = Matrix([[p1, 0], [0, p2]])  # 初始值

#K = 0
tag_cx_1 = 0
tag_cy_1 = 0
X_est = 0
XV_est = 0
Y_est = 0
YV_est = 0
MEA_XV = 0
MEA_cx = 0

MEA_YV = 0
MEA_cy = 0
rect2 = 0
rect3 = 0
#pre_V = 0
#pre_YV = 0

R = Matrix([[1, 0], [0, 1]])  # 测量噪声协方差矩阵
Q = Matrix([[5, 0], [0, 1]])  # 过程噪声协方差矩阵  第一个参数代表位置,第二个参数代表速度0
P = Matrix([[p1, 0], [0, p2]])  # 误差协方差矩阵
I = Matrix([[1, 0], [0, 1]])  # 单位矩阵

def kalman_x_filter(cx, vx):
    global A, T, K, pre_xk, pre_xpk, est_xk_1, est_xpk_1, est_xpk, est_xk

    # 预测阶段
    pre_xk = A.dot(est_xk_1)  # k时刻状态(x)预测值
    pre_xpk = A.dot(est_xpk_1).dot(A.T()) + Q

    # 更新阶段
    # 卡尔曼增益
    x_Zk = Matrix([[cx], [vx]])  # 观测矩阵
    x_Kk = pre_xpk.dot(H.T()).dot((H.dot(pre_xpk).dot(H.T()) + R).ni_matrix())  # 观测矩阵H可视为1
#    x_Kk = pre_xpk.dot((pre_xpk + R).ni_matrix())

    # 状态更新
    est_xk = pre_xk + x_Kk.dot(x_Zk - H.dot(pre_xk))
    est_xk_1 = est_xk

    # 协方差更新
    est_xpk = (I - x_Kk).dot(pre_xpk)
    est_xpk_1 = est_xpk
    return est_xk.__getitem__([0, 0]), est_xk.__getitem__([1, 0])


def kalman_y_filter(cy, vy):
    global A, T, K, pre_yk, pre_ypk, est_yk_1, est_ypk_1, est_ypk, est_yk

    pre_yk = A.dot(est_yk_1)
    pre_ypk = A.dot(est_ypk_1).dot(A.T()) + Q

    y_Zk = Matrix([[cy], [vy]])
    y_Kk = pre_ypk.dot(H.T()).dot((H.dot(pre_ypk).dot(H.T()) + R).ni_matrix())
#    y_Kk = pre_ypk.dot((pre_ypk + R).ni_matrix())

    est_yk = pre_yk + y_Kk.dot(y_Zk - H.dot(pre_yk))
    est_yk_1 = est_yk

    est_ypk = (I - y_Kk).dot(pre_ypk)
    est_ypk_1 = est_ypk
    return est_yk.__getitem__([0, 0]), est_yk.__getitem__([1, 0])

主函数

#________________________________主函数___________________________________________

sensor.reset()
sensor.set_pixformat(sensor.RGB565)
sensor.set_framesize(sensor.QQVGA)
sensor.skip_frames(time = 2000)
sensor.set_auto_gain(False)
clock = time.clock()
rtc = RTC()
rtc.datetime((0, 0, 0, 0, 0, 0, 0, 0))
date1 = [0, 0, 0, 0, 0, 0, 0, 0]

while(True):
    clock.tick()
    img = sensor.snapshot()
    find_tag = img.find_apriltags()
    if find_tag:
        date2 = rtc.datetime()
        for tag in img.find_apriltags(fx=f_x, fy=f_y, cx=c_x, cy=c_y):

            T = 1/255*(date1[7] - date2[7]) + date1[6] - date2[6]  # 计算公式
            date1[7] = date2[7]
            date1[6] = date2[6]

            MEA_XV = (tag.cx() - tag_cx_1) / T
            MEA_VY = (tag.cy() - tag_cy_1) / T

            MEA_cx=tag.cx()
            MEA_cy=tag.cy()
            tag_cx_1 = tag.cx()
            tag_cy_1 = tag.cy()

            A = Matrix([[1,T],[0,1]])
            X_est, XV_est = kalman_x_filter(MEA_cx, MEA_XV)
            Y_est, VY_est = kalman_y_filter(MEA_cy, MEA_VY)

            # img.draw_rectangle(int(tag.cx()-(tag.rect()[2])/2),int(tag.cy()-(tag.rect()[3])/2),tag.rect()[2],tag.rect()[3],color = (255, 0, 0))
            img.draw_cross(tag.cx(), tag.cy(), color = (255, 0, 0))
            img.draw_rectangle(tag.rect(),color = (255, 0, 0))

            img.draw_rectangle(int(X_est-(tag.rect()[2])/2),int(Y_est-(tag.rect()[3])/2),tag.rect()[2],tag.rect()[3],color = (0, 0, 255))
            img.draw_cross(int(X_est), int(Y_est), color = (0, 0, 255))
            rect2, rect3=tag.rect()[2],tag.rect()[3]

    else:  # 脱离目标采用卡尔曼滤波器进行预测
#        A = Matrix([[1,T],[0,1]])
#        a=1.6
        # X_est, XV_est=kalman_x_filter(X_est,a*MEA_V)
        # Y_est, YV_est=kalman_y_filter(Y_est,a*MEA_YV)
        X_est, XV_est=kalman_x_filter(X_est,MEA_XV)
        Y_est, YV_est=kalman_y_filter(Y_est,MEA_YV)
        img.draw_rectangle(int(X_est-rect2/2),int(Y_est-rect3/2),rect2,rect3,color = (0, 0, 255))
        img.draw_cross(int(X_est), int(Y_est), color = (0, 0, 255))


  • 5
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值