CN109141475A - 一种dvl辅助sins鲁棒行进间初始对准方法 - Google Patents

一种dvl辅助sins鲁棒行进间初始对准方法 Download PDF

Info

Publication number
CN109141475A
CN109141475A CN201811075323.1A CN201811075323A CN109141475A CN 109141475 A CN109141475 A CN 109141475A CN 201811075323 A CN201811075323 A CN 201811075323A CN 109141475 A CN109141475 A CN 109141475A
Authority
CN
China
Prior art keywords
indicate
vector
indicates
moment
navigation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201811075323.1A
Other languages
English (en)
Other versions
CN109141475B (zh
Inventor
徐祥
徐大诚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Weishi dark blue Technology Co., Ltd
Original Assignee
Suzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Suzhou University filed Critical Suzhou University
Priority to CN201811075323.1A priority Critical patent/CN109141475B/zh
Publication of CN109141475A publication Critical patent/CN109141475A/zh
Application granted granted Critical
Publication of CN109141475B publication Critical patent/CN109141475B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
    • G01C25/005—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass initial alignment, calibration or starting-up of inertial devices

Landscapes

  • Engineering & Computer Science (AREA)
  • Manufacturing & Machinery (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种DVL辅助SINS鲁棒行进间初始对准方法，包括以下步骤，步骤(1)：获取传感器实时数据；步骤(2)：建立矢量视运动参数方程；步骤(3)：在建立参数方程的基础上，利用自适应鲁棒卡尔曼滤波估计参数矩阵；步骤(4)：建立基于估计参数的重构矢量；步骤(5)：利用最优基四元数姿态确定方法，并计算确定姿态与真实姿态之间的误差角；步骤(6)：初始对准过程运行时间为M，获取实时数据的时刻为k，若k＝M，则输出初始对准结果，完成初始对准过程，若k<M，表示初始对准过程未完成，则重复上述步骤(1)至步骤(5)，直至初始对准过程结束，该方法解决了当系统量测噪声出现异常噪声时，对准精度差、对准过程发散的问题。

Description

一种DVL辅助SINS鲁棒行进间初始对准方法
技术领域：
本发明涉及捷联惯性导航系统初始对准领域，具体而言涉及一种DVL辅助SINS鲁棒行进间初始对准方法。
背景技术：
捷联惯性导航系统初始对准技术是系统正常导航定位的关键技术之一，采用DVL辅助惯性导航系统进行行进间初始对准具有对准精度高、可靠性好、自主性强等优点。当前，初始对准可以分为粗对准和精对准两个过程，其中粗对准主要实现粗略的姿态估计；精对准则是在粗对准的基础上进行姿态精估计。众多学者都对捷联惯导系统初始对准技术进行了深入的研究。在精对准方面，通过引入鲁棒卡尔曼滤波技术可以实现对准过程的鲁棒化，提高系统的稳定性。但在粗对准方面，当前的研究热点均是基于最优基姿态确定的解析方法，难以通过鲁棒滤波技术对其进行直接改进。这也使得粗对准过程易受外部辅助信息噪声干扰，造成对准结果稳定性差的缺点，而且传统方法无法进行鲁棒化粗对准。
发明内容：
本发明的目的是为了提高初始对准过程的稳定性及抗干扰性，提出了一种DVL辅助SINS鲁棒行进间初始对准方法，该方法能够实现粗对准过程的鲁棒化，提高系统对准的稳定性。
本发明的技术方案具体如下：
一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，包括以下几个步骤：
步骤(1)：获取传感器实时数据；
步骤(2)：建立矢量视运动参数方程；
步骤(3)：在建立参数方程的基础上，利用自适应鲁棒卡尔曼滤波估计参数矩阵；
步骤(4)：建立基于估计参数的重构矢量；
步骤(5)：利用最优基四元数姿态确定方法，并计算确定姿态与真实姿态之间的误差角；
步骤(6)：初始对准过程运行时间为M，获取实时数据的时刻为k，k＝M若k＝M，则输出初始对准结果，完成初始对准过程，若k＜M，表示初始对准过程未完成，则重复上述步骤(1) 至步骤(5)，直至初始对准过程结束。
进一步的，步骤(1)中所述传感器实时数据包括陀螺仪数据和加速度计数据。
进一步的，步骤(2)中建立矢量视运动参数方程的方式为，
步骤(2.1)：定义解算所需的参考坐标系，
b-载体坐标系，表示捷联惯性导航系统三轴正交坐标系，其x轴、y轴和z轴分别指向载体的右-前-上；
n-导航坐标系，表示载体所在位置的地理坐标系，其三轴分别指向当地东向、北向和天向；
e-地球坐标系，表示原点在地心，x轴为地心指向本初子午线与赤道交点，z轴为地心指向北极点，y轴与x轴和z轴构成右手坐标系；
i-惯性坐标系，表示惯性空间非旋转坐标系；
b0-初始载体坐标系，表示惯导系统开机运行时刻的载体坐标系，并在整个对准过程中相对于惯性空间保持静止；
n0-初始导航坐标系，表示惯导系统开机运行时刻的导航坐标系，并在整个对准过程中相对于惯性空间保持静止；
e0-初始地球坐标系，表示惯导系统开机运行时刻的地球坐标系，并在整个对准过程中相对于惯性空间保持静止；
步骤(2.2)：矢量视运动建模，
由坐标变换关系可知：
式中：表示导航系速度的微分；表示载体系相对于导航系变化的方向余弦矩阵；表示方向余弦矩阵的微分；vb表示载体系速度；表示载体系速度微分；
由方向余弦矩阵计算关系可知：
式中：表示方向余弦矩阵的微分；表示载体系相对于导航系变化的方向余弦矩阵；表示载体系相对于导航系的转动角速度在载体系上的投影；×表示矢量叉乘运算；
由惯性导航系统比力方程可知：
式中：表示导航系速度的微分；表示载体系相对于导航系变化的方向余弦矩阵；fb表示比力；表示地球系相对于惯性系的旋转角速度在导航系的投影；表示导航系相对于地球系的转动角速度在导航系的投影；×表示矢量叉乘运算；vn表示导航系速度；gn表示重力加速度在导航系下的投影；
通过推导可知：
式中：表示载体系相对于导航系变化的方向余弦矩阵；表示载体系相对于导航系的转动角速度在载体系上的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；表示初始导航系相对于导航系的方向余弦矩阵；表示初始载体系相对于初始导航系的方向余弦矩阵；表示载体系相对于初始载体系的方向余弦矩阵；fb表示比力；表示地球系相对于惯性系的旋转角速度在导航系的投影；表示导航系相对于地球系的转动角速度在导航系的投影；×表示矢量叉乘运算；vn表示导航系速度；gn表示重力加速度在导航系下的投影；
对上式进行积分运算可得：
式中：表示载体系相对于初始载体系的方向余弦矩阵；表示地球系相对于惯性系的旋转角速度在载体系的投影；表示导航系相对于地球系的转动角速度在载体系的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；fb表示比力；表示初始载体系相对于初始导航系的方向余弦矩阵；表示导航系相对于初始导航系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
定义矢量构成：
式中：β表示观测矢量；α表示参考矢量；表示载体系相对于初始载体系的方向余弦矩阵；表示地球系相对于惯性系的旋转角速度在载体系的投影；表示导航系相对于地球系的转动角速度在载体系的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；fb表示比力；表示导航系相对于初始导航系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
构造矢量观测器：
式中：β表示观测矢量；α表示参考矢量；表示初始载体系相对于初始导航系的方向余弦矩阵；
对参考矢量和观测矢量进行积分运算并表示为离散形式：
式中：βk表示k时刻的观测矢量；表示k时刻的载体系相对于初始载体系的方向余弦矩阵； vb(k)表示k时刻的载体系速度；vb(0)表示初始时刻的载体系速度；β′k表示k时刻的中间变量；β′k-1表示k-1时刻的中间变量；表示k-1时刻载体系相对于初始载体系的方向余弦矩阵；表示k-1时刻导航系相对于k-1时刻载体系的方向余弦矩阵；Δvn(k-1)表示k-1时刻一个 DVL采样周期内速度增量的中间变量；Δvb(k-1)表示k-1时刻一个DVL采样周期内比力增量；αk表示k时刻的参考矢量；αk-1表示k-1时刻的参考矢量；Δtd表示DVL采样周期；I3表示三维单位矩阵；表示导航系相对于惯性系的旋转角速度在导航系的投影；×表示矢量叉乘运算；gn表示重力加速度在导航系下的投影；
由重力矢量视运动可知：
式中，α表示参考矢量；表示初始地球系相对于初始导航系的方向余弦矩阵；表示地球系相对于初始地球系的方向余弦矩阵；表示导航系相对于地球系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
对上式进行积分运算并做处理可得：
α＝ΦΓ(t)
式中：Γ(t)表示随时间变化的矢量；Φ表示参数矩阵；且可以表示为：
式中：Φ表示参数矩阵；g表示地球重力加速度；ωie表示地球自转角速率；L表示当地纬度；t表示时间；
观测矢量的真值可以表示为：
β＝ΞΓ(t)
式中：表示初始导航系相对于初始载体系的方向余弦矩阵，Φ表示参数矩阵；β表示观测矢量；Γ(t)表示随时间变化的矢量；
当代入DVL量测，并考虑量测矢量含有噪声时，可以表示为：
式中：表示k时刻的量测观测矢量；Ξ表示常值参数矩阵；Γ(k)表示k时刻的离散矢量；表示k时刻载体系相对于初始载体系的方向余弦矩阵；δvb表示量测噪声；
利用参数估计理论进行参数建模可得：
式中：Ξz，k表示k时刻参数矩阵第三行分量；wz，k-1表示过程噪声；表示k时刻观测矢量的第三个元素；Γ(k)表示k时刻时间矢量；表示k时刻量测噪声的第三个元素。
进一步的，步骤(3)中在建立参数方程的基础上，利用自适应鲁棒卡尔曼滤波估计参数矩阵的方式为，
步骤(3.1)：自适应鲁棒卡尔曼滤波时间更新：
Pz，k|k-1＝Pz，k-1|k-1+Qz，k-1|k-1
式中：表示k时刻的参数矩阵第三行分量的一步预测；表示k-1时刻的参数矩阵第三行分量的后验估计；Pz，k|k-1表示k时刻误差协方差一步预测；Pz，k-1|k-1表示k-1时刻误差协方差后验估计；Qz，k-1|k-1表示k-1时刻过程噪声；
步骤(3.2)：鲁棒滤波：
由HuberM估计理论可知，鲁棒滤波的最小化代价函数为：
式中：表示k时刻后验估计参数矩阵第三行分量；argmin表示最小化操作；Jk()表示代价函数；Ξz，k表示k时刻参数矩阵第三行分量的变量形式；ρ()表示核函数；ζz，k表示变量；
变量ζz，k可以计算为：
式中：Rz，k|k表示量测噪声；表示k时刻观测矢量的第三个元素；Ξz，k表示k时刻参数矩阵第三行分量；Γ(k)表示k时刻时间矢量；
核函数可以表示为：
式中：ρ()表示核函数；ζz，k表示变量；γ表示高斯阈值；
由此可得权值函数为：
式中：表示权值函数；sgn()表示符号函数；γ表高斯阈值；ζz，k表示变量；
卡尔曼滤波残差可以表示为：
式中：表示计算残差；表示k时刻观测矢量的第三个元素；表示k时刻的参数矩阵第三行分量的一步预测；Γ(k)表示k时刻时间矢量；
定义中间变量：
式中：表示定义的中间变量；Rz，k|k表示量测噪声；表示计算残差；
重构鲁棒化矢量：
式中：表示鲁棒化观测矢量；表示k时刻的参数矩阵第三行分量的一步预测；Γ(k)表示k时刻时间矢量；表示由变量决定的权值；表示计算残差；
步骤(3.3)：自适应鲁棒卡尔曼滤波量测更新：
利用上面鲁棒化矢量，可以得到卡尔曼滤波的量测更新为：
Gz，k＝Pz，k|k-1Γ(k)((Γ(k))TPz，k|k-1Γ(k)+Rz，k|k)
Pz，k|k＝Pz，k-1|k-1+Qz，k-1|k-1
式中：Gz，k表示卡尔曼滤波增益；Pz，k|k-1表示k时刻误差协方差一步预测；Γ(k)表示k时刻时间矢量；Rz，k|k表示量测噪声；表示k时刻后验估计参数矩阵第三行分量；表示k时刻的参数矩阵第三行分量的一步预测；表示鲁棒化观测矢量；Pz，k|k表示k时刻误差协方差后验估计。
进一步的，步骤(4)中建立基于估计参数的重构矢量的方式为，
利用自适应鲁棒卡尔曼滤波方程进行矢量重构可得：
式中：表示k时刻重构矢量；表示k时刻后验估计参数矩阵；Γ(k)表示k时刻时间矢量。
进一步的，步骤(5)中计算确定姿态与真实姿态之间的误差角的方式为，
由最优基四元数姿态确定方法可知，姿态K矩阵可以表示为：
式中：Kk表示k时刻的K矩阵；Kk-1表示k-1时刻的K矩阵；表示k时刻重构矢量；αk表示k时刻参考矢量。
进一步的，步骤(1)中，陀螺仪量测常值漂移误差为εb＝[0.02 0.02 0.02]T°/h，陀螺仪量测随机游走误差为输出频率为200Hz；加速度计量测常值漂移误差为陀螺仪量测随机游走误差为输出频率为200Hz。
进一步的，步骤(3)中，高斯阈值设定为γ＝1.345，卡尔曼滤波参数初值为 Pj，0|0＝(105m/s)2I4，量测噪声和过程噪声协方差设定为Rj，k|k＝(0.1m/s)2，Qj，k-1|k-1＝ (10-3m/s)2I4。
进一步的，步骤(6)中，M＝600。
借由上述方案，本发明至少具有以下优点：
(1)本发明采用矢量视运动建模与分析方法，实现了矢量视运动的量化表示，通过参数识别剔除异常噪声；
(2)本发明采用自适应鲁棒滤波方法，对参数矩阵进行了最有估计，同时弱化了异常噪声的影响；
(3)本发明采用矢量重构方法，有效的消除了异常噪声，实现了矢量的降噪过程，提高了对准结果的稳定性。
上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，并可依照说明书的内容予以实施，以下以本发明的较佳实施例并配合附图详细说明如后。
附图说明：
图1是DVL辅助SiNS鲁棒行进间初始对准流程图；
图2是载体运动仿真曲线图；
图3是矢量视运动重构曲线图；
图4是纵摇角误差曲线图；
图5是横摇角误差曲线图；
图6是航向角误差曲线图；
具体实施方式：
下面结合附图和实施举例对本发明作进一步的详细说明：
本实施例将本发明提出的一种DVL辅助SINS鲁棒行进间初始对准方法(在附图中标示为“新方法”)通过Matlab仿真软件进行仿真验证，从而证明对准过程的鲁棒化。仿真硬件环境均为Intel (R)Core(TM)T9600 CPU 2.80GHz，4G RAM，Windows 7操作系统。如图2和图3所示，为行进间对准过程载体运动曲线图及矢量视运动重构曲线图，图2中表示姿态的图中虚线代表横摇角，实线代表纵摇角。从图中可以看出，采用矢量重构技术可以有效的消除异常噪声干扰，矢量视运动可以较好的恢复。图4、图5和图6为DVL辅助SINS行进间初始对准误差图，从图中可以看出，采用矢量重构技术之后，对准结果有效的抑制了量测异常噪声的干扰，航向角对准误差可以在150s 左右达到1°的对准精度，而传统方法则受到外部异常噪声的干扰造成对准不稳定现象。
本发明是一种DVL辅助SINS鲁棒行进间初始对准方法，算法流程如图1所示，包括以下几个步骤：
步骤1：获取传感器实时数据，所述传感器实时数据包括陀螺仪数据和加速度计数据；
步骤2：建立矢量视运动参数方程，定义解算所需的参考坐标系如下：
b-载体坐标系，表示捷联惯性导航系统三轴正交坐标系，其x轴、y轴和z轴分别指向载体的右-前-上；
n-导航坐标系，表示载体所在位置的地理坐标系，其三轴分别指向当地东向、北向和天向；
e-地球坐标系，表示原点在地心，x轴为地心指向本初子午线与赤道交点，z轴为地心指向北极点，y轴与x轴和z轴构成右手坐标系；
i-惯性坐标系，表示惯性空间非旋转坐标系；
b0-初始载体坐标系，表示惯导系统开机运行时刻的载体坐标系，并在整个对准过程中相对于惯性空间保持静止；
n0-初始导航坐标系，表示惯导系统开机运行时刻的导航坐标系，并在整个对准过程中相对于惯性空间保持静止；
e0一初始地球坐标系，表示惯导系统开机运行时刻的地球坐标系，并在整个对准过程中相对于惯性空间保持静止；
由坐标变换关系可知：
式中：表示导航系速度的微分；表示载体系相对于导航系变化的方向余弦矩阵；表示方向余弦矩阵的微分；vb表示载体系速度；表示载体系速度微分；
由方向余弦矩阵计算关系可知：
式中：表示方向余弦矩阵的微分；表示载体系相对于导航系变化的方向余弦矩阵；表示载体系相对于导航系的转动角速度在载体系上的投影；×表示矢量叉乘运算；
由惯性导航系统比力方程可知：
式中：表示导航系速度的微分；表示载体系相对于导航系变化的方向余弦矩阵；fb表示比力；表示地球系相对于惯性系的旋转角速度在导航系的投影；表示导航系相对于地球系的转动角速度在导航系的投影；×表示矢量叉乘运算；vn表示导航系速度；gn表示重力加速度在导航系下的投影；
由上面推导可知：
式中：表示载体系相对于导航系变化的方向余弦矩阵；表示载体系相对于导航系的转动角速度在载体系上的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；表示初始导航系相对于导航系的方向余弦矩阵；表示初始载体系相对于初始导航系的方向余弦矩阵；表示载体系相对于初始载体系的方向余弦矩阵；fb表示比力；表示地球系相对于惯性系的旋转角速度在导航系的投影；表示导航系相对于地球系的转动角速度在导航系的投影；×表示矢量叉乘运算；vn表示导航系速度；gn表示重力加速度在导航系下的投影；
对上式进行积分运算可得：
式中：表示载体系相对于初始载体系的方向余弦矩阵；表示地球系相对于惯性系的旋转角速度在载体系的投影；表示导航系相对于地球系的转动角速度在载体系的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；fb表示比力；表示初始载体系相对于初始导航系的方向余弦矩阵；表示导航系相对于初始导航系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
定义矢量构成：
式中：β表示观测矢量；α表示参考矢量；表示载体系相对于初始载体系的方向余弦矩阵；表示地球系相对于惯性系的旋转角速度在载体系的投影；表示导航系相对于地球系的转动角速度在载体系的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；fb表示比力；表示导航系相对于初始导航系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
所以可以构造矢量观测器：
式中：β表示观测矢量；α表示参考矢量；表示初始载体系相对于初始导航系的方向余弦矩阵；
对参考矢量和观测矢量进行积分运算并表示为离散形式：
式中：βk表示k时刻的观测矢量；表示k时刻的载体系相对于初始载体系的方向余弦矩阵； vb(k)表示k时刻的载体系速度；vb(0)表示初始时刻的载体系速度；β′k表示k时刻的中间变量；β′k-1表示k-1时刻的中间变量；表示k-1时刻载体系相对于初始载体系的方向余弦矩阵；表示k-1时刻导航系相对于k-1时刻载体系的方向余弦矩阵；Δvn(k-1)表示k-1时刻一个 DVL采样周期内速度增量的中间变量；Δvb(k-1)表示k-1时刻一个DVL采样周期内比力增量；αk表示k时刻的参考矢量；αk-1表示k-1时刻的参考矢量；Δtd表示DVL采样周期；I3表示三维单位矩阵；表示导航系相对于惯性系的旋转角速度在导航系的投影；×表示矢量叉乘运算；gn表示重力加速度在导航系下的投影；
由重力矢量视运动可知
式中，α表示参考矢量；表示初始地球系相对于初始导航系的方向余弦矩阵；表示地球系相对于初始地球系的方向余弦矩阵；表示导航系相对于地球系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
对上式进行积分运算并做处理可得：
α＝ΦΓ(t)
式中：Γ(t)表示随时间变化的矢量；Φ表示参数矩阵；且可以表示为：
式中：Φ表示参数矩阵；g表示地球重力加速度；ωie表示地球自转角速率；L表示当地纬度；t表示时间；
因此，观测矢量的真值可以表示为：
β＝ΞΓ(t)
式中：表示初始导航系相对于初始载体系的方向余弦矩阵，Φ表示参数矩阵；β表示观测矢量；Γ(t)表示随时间变化的矢量；
当代入DVL量测，并考虑量测矢量含有噪声时，可以表示为：
式中：表示k时刻的量测观测矢量；Ξ表示常值参数矩阵；Γ(k)表示k时刻的离散矢量；表示k时刻载体系相对于初始载体系的方向余弦矩阵；δvb表示量测噪声；
利用参数估计理论进行参数建模可得：
式中：Ξz，k表示k时刻参数矩阵第三行分量；wz，k-1表示过程噪声；表示k时刻观测矢量的第三个元素；Γ(k)表示k时刻时间矢量；表示k时刻量测噪声的第三个元素；
步骤3：在建立参数方程的基础上，利用自适应鲁棒卡尔曼滤波估计参数矩阵；
(1)自适应鲁棒卡尔曼滤波时间更新
Pz，k|k-1＝Pz，k-1|k-1+Qz，k-1|k-1
式中：定示k时刻的参数矩阵第三行分量的一步预测；表示k-1时刻的参数矩阵第三行分量的后验估计；Pz，k|k-1表示k时刻误差协方差一步预测；Pz，k-1|k-1表示k-1时刻误差协方差后验估计；Qz，k-1|k-1表示k-1时刻过程噪声；
(2)鲁棒滤波：
由HuberM估计理论可知，鲁棒滤波的最小化代价函数为：
式中：表示k时刻后验估计参数矩阵第三行分量；argmin表示最小化操作；Jk()表示代价函数；Ξz，k表示k时刻参数矩阵第三行分量的变量形式；ρ()表示核函数；ζz，k表示变量；
变量ζz，k可以计算为：
式中：Rz，k|k表示量测噪声；表示k时刻观测矢量的第三个元素；Ξz，k表示k时刻参数矩阵第三行分量；Γ(k)表示k时刻时间矢量；
核函数可以表示为：
式中：ρ()表示核函数；ζz，k表示变量；γ表示高斯阈值；
由此可得权值函数为：
式中：表示权值函数；sgn()表示符号函数；γ表高斯阈值；ζz，k表示变量；
卡尔曼滤波残差可以表示为：
式中：表示计算残差；表示k时刻观测矢量的第三个元素；表示k时刻的参数矩阵第三行分量的一步预测；Γ(k)表示k时刻时间矢量；
定义中间变量：
式中：表示定义的中间变量；Rz，k|k表示量测噪声；表示计算残差；
因此重构鲁棒化矢量：
式中：表示鲁棒化观测矢量；表示k时刻的参数矩阵第三行分量的一步预测；Γ(k)表示k时刻时间矢量；表示由变量决定的权值；表示计算残差；
(3)自适应鲁棒卡尔曼滤波量测更新
利用上面鲁棒化矢量，可以得到卡尔曼滤波的量测更新为：
Gz，k＝Pz，k|k-1Γ(k)((Γ(k))TPz，k|k-1Γ(k)+Rz，k|k)
Pz，k|k＝Pz，k-1|k-1+Qz，k-1|k-1
式中：Gz，k表示卡尔曼滤波增益；Pz，k|k-1表示k时刻误差协方差一步预测；Γ(k)表示k时刻时间矢量；Rz，k|k表示量测噪声；表示k时刻后验估计参数矩阵第三行分量；表示k时刻的参数矩阵第三行分量的一步预测；表示鲁棒化观测矢量；Pz，k|k表示k时刻误差协方差后验估计；
步骤4：建立基于估计参数的重构矢量；
利用上面自适应鲁棒卡尔曼滤波方程进行矢量重构可得：
式中：表示k时刻重构矢量；表示k时刻后验估计参数矩阵；Γ(k)表示k时刻时间矢量；
步骤5：利用最优基四元数姿态确定方法，并计算确定姿态与真实姿态之间的误差角；
由最优基四元数姿态确定方法可知，姿态K矩阵可以表示为：
式中：Kk表示k时刻的K矩阵；Kk-1表示k-1时刻的K矩阵；表示k时刻重构矢量；αk表示k时刻参考矢量；
步骤(6)：初始对准过程运行时间为M，获取实时数据的时刻为k，k＝M若k＝M，则输出初始对准结果，完成初始对准过程，若k＜M，表示初始对准过程未完成，则重复上述步骤(1) 至步骤(5)，直至初始对准过程结束。
MATLAB仿真实验，在以下的仿真条件下，对该方法进行仿真实验：
陀螺仪量测常值漂移误差为εb＝[0.02 0.02 0.02]T°/h，陀螺仪量测随机游走误差为输出频率为200Hz；加速度计量测常值漂移误差为陀螺仪量测随机游走误差为输出频率为200Hz。高斯阈值设定为γ＝1.345，卡尔曼滤波参数初值为Pj，0|0＝(105m/s)2I4，量测噪声和过程噪声协方差设定为Rj，k|k＝(0.1m/s)2，Qj，k-1|k-1＝(10-3m/s)2I4。步骤(6)中，M＝600。
该DVL辅助SINS鲁棒行进间初始对准方法提高对矢量视运动进行建模分析，采用自适应鲁棒滤波方法对量测异常噪声进行弱化，提高初始对准过程的稳定性。为克服传统方法无法进行鲁棒化粗对准的问题，采用视运动矢量分析与建模技术，通过对参数矩阵的鲁棒化估计，可以实现粗对准过程的鲁棒化，提高系统对准的稳定性。
以上所述仅是本发明的优选实施方式，并不用于限制本发明，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (9)

1.一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，包括以下几个步骤：
步骤(1)：获取传感器实时数据；
步骤(2)：建立矢量视运动参数方程；
步骤(3)：在建立参数方程的基础上，利用自适应鲁棒卡尔曼滤波估计参数矩阵；
步骤(4)：建立基于估计参数的重构矢量；
步骤(5)：利用最优基四元数姿态确定方法，并计算确定姿态与真实姿态之间的误差角；
步骤(6)：初始对准过程运行时间为M，获取实时数据的时刻为k，k＝M若k＝M，则输出初始对准结果，完成初始对准过程，若k&lt;M，表示初始对准过程未完成，则重复上述步骤(1)至步骤(5)，直至初始对准过程结束。
2.根据权利要求1所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(1)中所述传感器实时数据包括陀螺仪数据和加速度计数据。
3.根据权利要求1所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(2)中建立矢量视运动参数方程的方式为，
步骤(2.1)：定义解算所需的参考坐标系，
b—载体坐标系，表示捷联惯性导航系统三轴正交坐标系，其x轴、y轴和z轴分别指向载体的右-前-上；
n—导航坐标系，表示载体所在位置的地理坐标系，其三轴分别指向当地东向、北向和天向；
e—地球坐标系，表示原点在地心，x轴为地心指向本初子午线与赤道交点，z轴为地心指向北极点，y轴与x轴和z轴构成右手坐标系；
i—惯性坐标系，表示惯性空间非旋转坐标系；
b0—初始载体坐标系，表示惯导系统开机运行时刻的载体坐标系，并在整个对准过程中相对于惯性空间保持静止；
n0—初始导航坐标系，表示惯导系统开机运行时刻的导航坐标系，并在整个对准过程中相对于惯性空间保持静止；
e0—初始地球坐标系，表示惯导系统开机运行时刻的地球坐标系，并在整个对准过程中相对于惯性空间保持静止；
步骤(2.2)：矢量视运动建模，
由坐标变换关系可知：
式中：表示导航系速度的微分；表示载体系相对于导航系变化的方向余弦矩阵；表示方向余弦矩阵的微分；vb表示载体系速度；表示载体系速度微分；
由方向余弦矩阵计算关系可知：
式中：表示方向余弦矩阵的微分；表示载体系相对于导航系变化的方向余弦矩阵；表示载体系相对于导航系的转动角速度在载体系上的投影；×表示矢量叉乘运算；
由惯性导航系统比力方程可知：
式中：表示导航系速度的微分；表示载体系相对于导航系变化的方向余弦矩阵；fb表示比力；表示地球系相对于惯性系的旋转角速度在导航系的投影；表示导航系相对于地球系的转动角速度在导航系的投影；×表示矢量叉乘运算；vn表示导航系速度；gn表示重力加速度在导航系下的投影；
通过推导可知：
式中：表示载体系相对于导航系变化的方向余弦矩阵；表示载体系相对于导航系的转动角速度在载体系上的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；表示初始导航系相对于导航系的方向余弦矩阵；表示初始载体系相对于初始导航系的方向余弦矩阵；表示载体系相对于初始载体系的方向余弦矩阵；fb表示比力；表示地球系相对于惯性系的旋转角速度在导航系的投影；表示导航系相对于地球系的转动角速度在导航系的投影；×表示矢量叉乘运算；vn表示导航系速度；gn表示重力加速度在导航系下的投影；
对上式进行积分运算可得：
式中：表示载体系相对于初始载体系的方向余弦矩阵；表示地球系相对于惯性系的旋转角速度在载体系的投影；表示导航系相对于地球系的转动角速度在载体系的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；fb表示比力；表示初始载体系相对于初始导航系的方向余弦矩阵；表示导航系相对于初始导航系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
定义矢量构成：
式中：β表示观测矢量；α表示参考矢量；表示载体系相对于初始载体系的方向余弦矩阵；表示地球系相对于惯性系的旋转角速度在载体系的投影；表示导航系相对于地球系的转动角速度在载体系的投影；×表示矢量叉乘运算；vb表示载体系速度；表示载体系速度微分；fb表示比力；表示导航系相对于初始导航系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
构造矢量观测器：
式中：β表示观测矢量；α表示参考矢量；表示初始载体系相对于初始导航系的方向余弦矩阵；
对参考矢量和观测矢量进行积分运算并表示为离散形式：
式中：βk表示k时刻的观测矢量；表示k时刻的载体系相对于初始载体系的方向余弦矩阵；vb(k)表示k时刻的载体系速度；vb(0)表示初始时刻的载体系速度；β′k表示k时刻的中间变量；β′k-1表示k-1时刻的中间变量；表示k-1时刻载体系相对于初始载体系的方向余弦矩阵；表示k-1时刻导航系相对于k-1时刻载体系的方向余弦矩阵；Δvn(k-1)表示k-1时刻一个DVL采样周期内速度增量的中间变量；Δvb(k-1)表示k-1时刻一个DVL采样周期内比力增量；αk表示k时刻的参考矢量；αk-1表示k-1时刻的参考矢量；Δtd表示DVL采样周期；I3表示三维单位矩阵；表示导航系相对于惯性系的旋转角速度在导航系的投影；×表示矢量叉乘运算；gn表示重力加速度在导航系下的投影；
由重力矢量视运动可知：
式中，α表示参考矢量；表示初始地球系相对于初始导航系的方向余弦矩阵；表示地球系相对于初始地球系的方向余弦矩阵；表示导航系相对于地球系的方向余弦矩阵；gn表示重力加速度在导航系下的投影；
对上式进行积分运算并做处理可得：
α＝ΦΓ(t)
式中：Γ(t)表示随时间变化的矢量；Φ表示参数矩阵；且可以表示为：
式中：Φ表示参数矩阵；g表示地球重力加速度；ωie表示地球自转角速率；L表示当地纬度；t表示时间；
观测矢量的真值可以表示为：
β＝ΞΓ(t)
式中： 表示初始导航系相对于初始载体系的方向余弦矩阵，Φ表示参数矩阵；β表示观测矢量；Γ(t)表示随时间变化的矢量；
当代入DVL量测，并考虑量测矢量含有噪声时，可以表示为：
式中：表示k时刻的量测观测矢量；Ξ表示常值参数矩阵；Γ(k)表示k时刻的离散矢量；表示k时刻载体系相对于初始载体系的方向余弦矩阵；δvb表示量测噪声；
利用参数估计理论进行参数建模可得：
式中：Ξz,k表示k时刻参数矩阵第三行分量；wz,k-1表示过程噪声；表示k时刻观测矢量的第三个元素；Γ(k)表示k时刻时间矢量；表示k时刻量测噪声的第三个元素。
4.根据权利要求3所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(3)中在建立参数方程的基础上，利用自适应鲁棒卡尔曼滤波估计参数矩阵的方式为，
步骤(3.1)：自适应鲁棒卡尔曼滤波时间更新：
Pz,k|k-1＝Pz,k-1|k-1+Qz,k-1|k-1
式中：表示k时刻的参数矩阵第三行分量的一步预测；表示k-1时刻的参数矩阵第三行分量的后验估计；Pz,k|k-1表示k时刻误差协方差一步预测；Pz,k-1|k-1表示k-1时刻误差协方差后验估计；Qz,k-1|k-1表示k-1时刻过程噪声；
步骤(3.2)：鲁棒滤波：
由HuberM估计理论可知，鲁棒滤波的最小化代价函数为：
式中：表示k时刻后验估计参数矩阵第三行分量；argmin表示最小化操作；Jk()表示代价函数；Ξz,k表示k时刻参数矩阵第三行分量的变量形式；ρ()表示核函数；ζz,k表示变量；
变量ζz,k可以计算为：
式中：Rz,k|k表示量测噪声；表示k时刻观测矢量的第三个元素；Ξz,k表示k时刻参数矩阵第三行分量；Γ(k)表示k时刻时间矢量；
核函数可以表示为：
式中：ρ()表示核函数；ζz,k表示变量；γ表示高斯阈值；
由此可得权值函数为：
式中：表示权值函数；sgn()表示符号函数；γ表高斯阈值；ζz,k表示变量；
卡尔曼滤波残差可以表示为：
式中：表示计算残差；表示k时刻观测矢量的第三个元素；表示k时刻的参数矩阵第三行分量的一步预测；Γ(k)表示k时刻时间矢量；
定义中间变量：
式中：表示定义的中间变量；Rz,k|k表示量测噪声；表示计算残差；
重构鲁棒化矢量：
式中：表示鲁棒化观测矢量；表示k时刻的参数矩阵第三行分量的一步预测；Γ(k)表示k时刻时间矢量；表示由变量决定的权值；表示计算残差；
步骤(3.3)：自适应鲁棒卡尔曼滤波量测更新：
利用上面鲁棒化矢量，可以得到卡尔曼滤波的量测更新为：
Pz,k|k＝Pz,k-1|k-1+Qz,k-1|k-1
式中：Gz,k表示卡尔曼滤波增益；Pz,k|k-1表示k时刻误差协方差一步预测；Γ(k)表示k时刻时间矢量；Rz,k|k表示量测噪声；表示k时刻后验估计参数矩阵第三行分量；表示k时刻的参数矩阵第三行分量的一步预测；表示鲁棒化观测矢量；Pz,k|k表示k时刻误差协方差后验估计。
5.根据权利要求4所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(4)中建立基于估计参数的重构矢量的方式为，
利用自适应鲁棒卡尔曼滤波方程进行矢量重构可得：
式中：表示k时刻重构矢量；表示k时刻后验估计参数矩阵；Γ(k)表示k时刻时间矢量。
6.根据权利要求5所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(5)中计算确定姿态与真实姿态之间的误差角的方式为，
由最优基四元数姿态确定方法可知，姿态K矩阵可以表示为：
式中：Kk表示k时刻的K矩阵；Kk-1表示k-1时刻的K矩阵；表示k时刻重构矢量；αk表示k时刻参考矢量。
7.根据权利要求1所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(1)中，陀螺仪量测常值漂移误差为εb＝[0.02 0.02 0.02]To/h，陀螺仪量测随机游走误差为输出频率为200Hz；加速度计量测常值漂移误差为陀螺仪量测随机游走误差为输出频率为200Hz。
8.根据权利要求4所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(3)中，高斯阈值设定为γ＝1.345，卡尔曼滤波参数初值为Pj,0|0＝(105m/s)2I4，量测噪声和过程噪声协方差设定为Rj,k|k＝(0.1m/s)2，Qj,k-1|k-1＝(10-3m/s)2I4。
9.根据权利要求1所述的一种DVL辅助SINS鲁棒行进间初始对准方法，其特征在于，步骤(6)中，M＝600。
CN201811075323.1A 2018-09-14 2018-09-14 一种dvl辅助sins鲁棒行进间初始对准方法 Active CN109141475B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811075323.1A CN109141475B (zh) 2018-09-14 2018-09-14 一种dvl辅助sins鲁棒行进间初始对准方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811075323.1A CN109141475B (zh) 2018-09-14 2018-09-14 一种dvl辅助sins鲁棒行进间初始对准方法

Publications (2)

Publication Number Publication Date
CN109141475A true CN109141475A (zh) 2019-01-04
CN109141475B CN109141475B (zh) 2020-06-05

Family

ID=64825497

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811075323.1A Active CN109141475B (zh) 2018-09-14 2018-09-14 一种dvl辅助sins鲁棒行进间初始对准方法

Country Status (1)

Country Link
CN (1) CN109141475B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110108301A (zh) * 2019-05-14 2019-08-09 苏州大学 模值检测动基座鲁棒对准方法
CN110987018A (zh) * 2019-12-19 2020-04-10 苏州大学 比力微分的位置法dvl误差标定方法及系统
CN111854747A (zh) * 2020-08-25 2020-10-30 东南大学 一种载体大机动情况下的dvl辅助sins粗对准方法
CN111912427A (zh) * 2019-05-10 2020-11-10 中国人民解放军火箭军工程大学 一种多普勒雷达辅助捷联惯导运动基座对准方法及系统
CN112525218A (zh) * 2020-11-23 2021-03-19 哈尔滨工程大学 一种ins/dvl组合导航系统鲁棒智能协同校准方法
CN113670335A (zh) * 2021-08-18 2021-11-19 河海大学 基于dvl辅助和向量截断化k矩阵的水下载体初始对准方法
CN114485723A (zh) * 2021-02-08 2022-05-13 北京理工大学 一种自适应鲁棒矩阵卡尔曼滤波的高旋体空中对准方法
CN116295511A (zh) * 2022-12-16 2023-06-23 南京安透可智能系统有限公司 一种用于管道潜航机器人的鲁棒初始对准方法及系统
CN116499492A (zh) * 2022-12-05 2023-07-28 华中光电技术研究所(中国船舶集团有限公司第七一七研究所) 一种基于dvl辅助的匀速直航下捷联罗经粗对准方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103278163A (zh) * 2013-05-24 2013-09-04 哈尔滨工程大学 一种基于非线性模型的sins/dvl组合导航方法
CN105043415A (zh) * 2015-07-13 2015-11-11 北京工业大学 基于四元数模型的惯性系自对准方法
CN105806363A (zh) * 2015-11-16 2016-07-27 东南大学 基于srqkf的sins/dvl水下大失准角对准方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103278163A (zh) * 2013-05-24 2013-09-04 哈尔滨工程大学 一种基于非线性模型的sins/dvl组合导航方法
CN105043415A (zh) * 2015-07-13 2015-11-11 北京工业大学 基于四元数模型的惯性系自对准方法
CN105806363A (zh) * 2015-11-16 2016-07-27 东南大学 基于srqkf的sins/dvl水下大失准角对准方法

Non-Patent Citations (8)

* Cited by examiner, † Cited by third party
Title
KANGHUA TANG等: "A Novel INS and Doppler Sensors Calibration Method for Long Range Underwater Vehicle Navigation", 《SENSORS》 *
WU Y等: "Velocity/Position Integration Formula Part I:Application to In-Flight Coarse Alignment", 《IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS》 *
YI ZHANG: "An approach of DVL-AIDED SDINS alignment for in-motion vessel", 《OPTIK》 *
安德玺等: "一种基于滤波参数在线辨识的鲁棒自适应滤波器", 《自动化学报》 *
张朝飞等: "抗扰动的捷联惯导系统回溯参数辨识对准法", 《中国惯性技术学报》 *
徐祥等: "一种改良Kalman滤波参数辨识粗对准方法", 《中国惯性技术学报》 *
童余德等: "基于MATLAB的卡尔曼滤波法参数辨识与仿真", 《航电技术》 *
罗莉: "水下航行器捷联惯导系统粗对准方法研究", 《中国优秀硕士学位论文全文数据库工程科技II辑》 *

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111912427B (zh) * 2019-05-10 2022-03-01 中国人民解放军火箭军工程大学 一种多普勒雷达辅助捷联惯导运动基座对准方法及系统
CN111912427A (zh) * 2019-05-10 2020-11-10 中国人民解放军火箭军工程大学 一种多普勒雷达辅助捷联惯导运动基座对准方法及系统
CN110108301A (zh) * 2019-05-14 2019-08-09 苏州大学 模值检测动基座鲁棒对准方法
CN110987018B (zh) * 2019-12-19 2023-11-24 苏州大学 比力微分的位置法dvl误差标定方法及系统
CN110987018A (zh) * 2019-12-19 2020-04-10 苏州大学 比力微分的位置法dvl误差标定方法及系统
CN111854747B (zh) * 2020-08-25 2022-08-12 东南大学 一种载体大机动情况下的dvl辅助sins粗对准方法
CN111854747A (zh) * 2020-08-25 2020-10-30 东南大学 一种载体大机动情况下的dvl辅助sins粗对准方法
CN112525218A (zh) * 2020-11-23 2021-03-19 哈尔滨工程大学 一种ins/dvl组合导航系统鲁棒智能协同校准方法
CN114485723A (zh) * 2021-02-08 2022-05-13 北京理工大学 一种自适应鲁棒矩阵卡尔曼滤波的高旋体空中对准方法
CN114485723B (zh) * 2021-02-08 2024-02-27 北京理工大学 一种自适应鲁棒矩阵卡尔曼滤波的高旋体空中对准方法
CN113670335A (zh) * 2021-08-18 2021-11-19 河海大学 基于dvl辅助和向量截断化k矩阵的水下载体初始对准方法
CN113670335B (zh) * 2021-08-18 2023-10-24 河海大学 基于dvl辅助和向量截断化k矩阵的水下载体初始对准方法
CN116499492A (zh) * 2022-12-05 2023-07-28 华中光电技术研究所(中国船舶集团有限公司第七一七研究所) 一种基于dvl辅助的匀速直航下捷联罗经粗对准方法
CN116295511A (zh) * 2022-12-16 2023-06-23 南京安透可智能系统有限公司 一种用于管道潜航机器人的鲁棒初始对准方法及系统
CN116295511B (zh) * 2022-12-16 2024-04-02 南京安透可智能系统有限公司 一种用于管道潜航机器人的鲁棒初始对准方法及系统

Also Published As

Publication number Publication date
CN109141475B (zh) 2020-06-05

Similar Documents

Publication Publication Date Title
CN109141475A (zh) 一种dvl辅助sins鲁棒行进间初始对准方法
CN110207697B (zh) 基于角加速度计/陀螺/加速度计的惯性导航解算方法
EP3819592B1 (en) Positioning device, positioning method, and program
CN105180937B (zh) 一种mems‑imu初始对准方法
Li et al. Optimization-based INS in-motion alignment approach for underwater vehicles
CN106342284B (zh) 一种飞行载体姿态确定方法
CN104655131B (zh) 基于istssrckf的惯性导航初始对准方法
CN108680186B (zh) 基于重力仪平台的捷联式惯导系统非线性初始对准方法
CN103743414B (zh) 一种里程计辅助车载捷联惯导系统行进间初始对准方法
CN109596144B (zh) Gnss位置辅助sins行进间初始对准方法
CN105737823B (zh) 一种基于五阶ckf的gps/sins/cns组合导航方法
CN109163735B (zh) 一种晃动基座正向-正向回溯初始对准方法
Stančić et al. The integration of strap-down INS and GPS based on adaptive error damping
CN109931955B (zh) 基于状态相关李群滤波的捷联惯性导航系统初始对准方法
CN107289930A (zh) 基于mems惯性测量单元的纯惯性车辆导航方法
CN107390247A (zh) 一种导航方法、系统及导航终端
Sun et al. Adaptive sensor data fusion in motion capture
CN103017787A (zh) 适用于摇摆晃动基座的初始对准方法
CN112432642A (zh) 一种重力灯塔与惯性导航融合定位方法及系统
CN112902956A (zh) 一种手持式gnss/mems-ins接收机航向初值获取方法、电子设备、存储介质
CN113340298A (zh) 一种惯导和双天线gnss外参标定方法
Sun et al. Coarse alignment based on IMU rotational motion for surface ship
Liu et al. Attitude determination for MAVs using a Kalman filter
CN109443378A (zh) 速度辅助行进间回溯初始对准方法
CN110095117A (zh) 一种无陀螺惯性量测系统与gps组合的导航方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20211111

Address after: 100192 1217b, 12th floor, building a 1, Qinghe Jiayuan East District, Haidian District, Beijing

Patentee after: Beijing Weishi dark blue Technology Co., Ltd

Address before: 215000 No. 8, Jixue Road, Xiangcheng District, Suzhou City, Jiangsu Province

Patentee before: Suzhou University

TR01 Transfer of patent right