CN109000782A - 一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法 - Google Patents

一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法 Download PDF

Info

Publication number
CN109000782A
CN109000782A CN201811128365.7A CN201811128365A CN109000782A CN 109000782 A CN109000782 A CN 109000782A CN 201811128365 A CN201811128365 A CN 201811128365A CN 109000782 A CN109000782 A CN 109000782A
Authority
CN
China
Prior art keywords
parameter
kalman filtering
equation
state
error
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.)
Pending
Application number
CN201811128365.7A
Other languages
English (en)
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.)
Harbin Engineering University
Original Assignee
Harbin Engineering 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 Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN201811128365.7A priority Critical patent/CN109000782A/zh
Publication of CN109000782A publication Critical patent/CN109000782A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01HMEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H9/00Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means
    • G01H9/004Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means using fibre optic sensors

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Filters That Use Time-Delay Elements (AREA)

Abstract

本发明属于相位生成载波解调技术领域,公开了一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法,随机选取五个数据点,代入椭圆方程中,构建出一个五元代数方程组;选取椭圆参数估计值的偏差作为状态向量;将待拟合的数据点依次利用卡尔曼滤波算法,不断的更新状态向量和状态协方差矩阵;新更新的协方差矩阵作为下一时刻的误差协方差矩阵;将更新获得的新参数的偏差值与初始的参数估计值进行相加,得到椭圆参数新的估计结果;当更新过程中相邻两次估计的结果变化均小于给定误差10‑8时认为达到收敛状态,停止迭代;估计值作为最优估计值,即最优的椭圆参数。本发明相对幅度和谐波抑制比有了很大的改善,可以有效校正非线性误差,提高解调的精度。

Description

一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法
技术领域
本发明属于相位生成载波解调技术领域,尤其涉及一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法。
背景技术
相位生成载波解调(PGC)技术是一种广泛应用于干涉型光纤传感器的零差解调法。PGC主要有反正切以及交叉微分相乘(DCM)算法。反正切解调方法较DCM而言,因其受光源光强波动影响小、受滤波器通带纹波影响小、线性度好、动态范围大等优点,已经作为实用的解调技术应用于光纤水听器阵列解调。由于反正切解调在干涉信号经过乘载波以及低通滤波之后的处理是一种非线性的算法过程。非线性误差成为影响PGC解调性能的关键因素。为了提高测量精度,尽可能降低非线性误差带来的测量误差,进行系统的非线性校正是系统中必不可少的重要工作。PGC解调技术对相位的解调精度会因为非线性误差的存在而严重变差。如果能够对误差进行校正,解调的精度会得到很大的提升。传统的方法中对于低频部分激光强度带来的松弛振荡噪声所产生的寄生调幅项,将泵浦功率的改善与光电负反馈技术结合起来共同控制电路的输出,保证噪声频率的一半等于调制频率半倍频的奇数倍,从激光器端抑制光强噪声滤除混叠噪声,同时光电负反馈的引入,弥补了多次修改泵浦功率的不足;基于数字解调器的反正切方法,利用数字解调处理的方法,与传统的微分和交叉乘法(DCM)相比具有测量范围扩大、电路相对简单、操作性强、噪声的性能更好等优点;采用3×2耦合器可以得到两路相位差固定的信号,再利用微分交叉相乘的方法,直接对信号进行调制,可以有效消除误差;进行基于PGC-DCM改进的相位载波调制技术,构造出DCM技术解调后的多余项,然后与DCM解调结果相除消除误差项。随着技术的发展更新,上述方法在精度上已经不能满足人们的需要。现在已经提出更高精度的蒙特卡洛法、增益调整法、周期误差补偿法,并将非线性误差校正技术中已有的神经网络、椭圆拟合法等方法应用于PGC解调过程中非线性误差的校正。蒙特卡洛法通过概率统计确定非线性误差大小,是一种比较新型的非线性误差估计法。增益调整法通过增益调整法模拟调整检测器调整增益,再利用周期补偿法降低周期误差,使得非线性误差降低到0.01nm的量级。以往单频干涉仪在解调瞬时相位时多会用到反三角函数,反三角函数由于具有多值映射的特性,求解过程中会包含分段函数。神经网络通常是光滑连续的函数,有利于系统性能的综合分析,在测量干涉仪瞬时相位时达到较高的精度。
综上所述,现有技术存在的问题是:易受孤立点和噪声的影响,少量孤立点就可以使结果有较大变化,稳健性比较差而当待拟合数据集中在较短的弧段时,不同偏差大小的噪声也会使拟合结果发生改变。且当给的样本点数比较大的时候计算量过大,需要时间和空间代价都过于庞大,并不适合进行实时的椭圆拟合。
发明内容
针对现有技术存在的问题,本发明提供了基于卡尔曼滤波的椭圆拟合非线性误差校正方法。
本发明借鉴拓展卡尔曼滤波方法,将非线性系统泰勒展开并略去二阶及以上项得到一个近似的线性系统,然后应用卡尔曼滤波完成对目标的滤波估计。本发明是一种基于卡尔曼滤波的椭圆拟合校正方法,所述基于卡尔曼滤波的椭圆拟合校正方法技术方案如下:
步骤一,随机选取五个数据点(xi,yi),i=1~5代入椭圆方程ax2+2bxy+(1-a)y2+2dx+2ey+f=0中,构建出一个五元代数方程组,方程组的解为椭圆的五个参数(a,b,d,e,f),作为椭圆参数初始值;
步骤二,选取椭圆参数估计值的偏差作为状态向量,初始值为[0,0,0,0,0]T,由于初始参数和真实参数差距可能性很大,所以选取无限大的对角矩阵作为状态协方差矩阵的初值;
步骤三,将待拟合的数据点依次利用卡尔曼滤波算法,不断的更新状态向量和状态协方差矩阵,卡尔曼滤波算法流程及状态更新方法会在下面描述;
步骤四,每一时刻状态的更新都由前一时刻的估计和这一时刻的观测数据计算得到,每一步都将新更新的参数偏差作为下一时刻参数偏差的预测值,将新更新的协方差矩阵作为下一时刻的误差协方差矩阵;
步骤五,将更新获得的新参数的偏差值与初始的参数估计值进行相加,得到椭圆参数新的估计结果;
步骤六,不断的重复执行步骤三~步骤五,每次的初始偏差为[0,0,0,0,0]T,初始参数为上次输出的最优偏差与上次的初始参数之和,误差协矩阵为上依次循环最终输出的误差协方差矩阵;当更新过程中相邻两次估计的结果变化均小于给定误差10-8时认为达到收敛状态,停止迭代;估计值作为最优估计值,即最优的椭圆参数。
所述基于卡尔曼滤波的椭圆拟合校正方法的椭圆及一般二次曲线的描述为:
f(x,y)=ax2+2bxy+cy2+2dx+2ey+f=0;
其中对于椭圆方程满足a+c>0,对椭圆方程归一化满足条件:
a+c=1;
椭圆参数用一个五维向量来表示:
X=(a,b,d,e,f)';
估计用到的观测值为系统的待拟合点表示:
Ni=(xi,yi)';
观测数据与真实数据的偏差为:
vi=Yi-Ni
其中:vi为高斯噪声,其对应的噪声的协方差矩阵为:
cov[vi]=Ri
进一步,所述基于卡尔曼滤波的椭圆拟合校正方法的椭圆方程是一个非线性的约束,待测参数与待估计点代入作为观测方程:
F(Xi,Yi)=0;
F(X,Y)为代入椭圆真实参数时的值,偏差为0;当给出估计的参数M和观测的数据N,将待拟合的数据点代入估计的椭圆方程中得到的偏差为F(M,N)≠0,存在一定的偏差。
进一步,所述基于卡尔曼滤波的椭圆拟合校正方法的卡尔曼滤波对信号进行线性化处理:对F(Xi,Yi)=0在椭圆参数Xi的估计值Mi和Yi的观测数据Ni处利用泰勒级数对非线性函数进行线性化处理,得到:
F(Mi,Ni)表示将估计参数Mi和观测值Ni代入式,得到的偏离椭圆的距离,Mi=[aa,bb,dd,ee,ff]'对于系统-F(Mi,Ni)表示如下:
-F(Mi,Ni)=-(aaxi 2+2bbxiyi+ccyi 2+2ddxi+2eeyi+ff);
其中通过求导得:
正常卡尔曼滤波的观测方程为Z=HX+v,上式中令ΔX=X-M则得到:
偏差观测值:
Z=-F(Mi,Ni);
观测矩阵:
观测噪声:
观测噪声期望和协方差矩阵为:
进一步,所述基于卡尔曼滤波的椭圆拟合校正方法的将每次更新后的参数代入椭圆方程,与真实椭圆的偏差作为系统观测值,对应标准卡尔曼滤波的公式,具体步骤位:
步骤一,系统建模:
(1)状态方程:
ΔX(k)=A×ΔX(k-1)+w(k-1);
系统不包含状态方程,视作A=1,w(k-1)=0;
(2)观测方程:
Z(k)=-F(Mi,Ni)=-(aaxi 2+2bbxiyi+ccyi 2+2ddxi+2eeyi+ff);
步骤二,初始化:
(1)状态变量的初始化:
ΔX(0)=[0,0,0,0,0]T
(2)状态协方差的初始化:
选取[0,0,0,0,0]T作为状态向量的初值;选择初值为1的无限大的对角矩阵作为状态协方差矩阵的初值;
步骤三,预测过程:
(1)状态预测:
Xekf(k|k-1)=A×Xekf(k-1);
其中状态转移矩阵A=1,上式为:
Xekf(k|k-1)=Xekf(k-1);
直接将上一时刻的估计值作为这一时刻的一步预测值;
(2)状态协方差预测矩阵:
Pekf(k|k-1)=A×Pekf(k-1)×A'+Q(k-1);
其中A=1,w=0,Q是过程噪声w的协方差矩阵,Q=0,等效为:
Pekf(k|k-1)=Pekf(k-1);
步骤四,更新过程:
(1)增益矩阵的更新方程:
Kg(k)=Pekf(k|k-1)×H'/(H×Pekf(k|k-1)×H'+R);
(2)状态估计更新方程:
Xekf(k)=Xekf(k|k-1)+Kg(k)×(Z(k)-H×Xekf(k|k-1));
(3)误差协方差更新方程:
Pekf(k)=(I-Kg(k)×H)×Pekf(k|k-1);
再将k时刻的状态估计值和误差协方差估计值送入下一时刻,进入不断迭代的过程。
本发明的另一目的在于提供一种利用所述基于卡尔曼滤波的椭圆拟合校正方法的光纤水听器。
本发明基于卡尔曼滤波的椭圆拟合的校正算法,利用卡尔曼滤波法计算出各个校正参数。由于本发明采用了迭代的方法,使得存储量少,简单易行,便于实现,且极大的加快了运算速度,可以用在各种对实时性要求较高的领域。同时,进行系统估计时,不仅考虑了动态测量值,还考虑系统内部的动态变化,将预测值与测量值做权重和,使得测量误差平方和的期望最小。实验结果(图6、7、8)表明:较直接解调和最小二乘方法,卡尔曼滤波的同频相对幅度和谐波抑制比有了很大的改善。本发明可以有效校正非线性误差,提高解调的精度。
附图说明
图1是本发明实施例提供的基于卡尔曼滤波的椭圆拟合校正方法流程图。
图2是本发明实施例提供的校正误差的整体实现流程图。
图3是本发明实施例提供的伴生调幅指数m在卡尔曼滤波算法校正下对系统指标的影响示意图;
图中:(a)RAB随伴生调幅指数变化
(b)HSR随伴生调幅指数变化
图4是本发明实施例提供的载波幅度偏差程度在卡尔曼滤波算法校正下对系统指标的影响示意图;
图中:(a)RAE随Gp/Hp的变化
(b)HSR随Gp/Hp的变化
图5是本发明实施例提供的载波相移在卡尔曼滤波算法下对系统指标的影响示意图;
图中:(a)RAE随检波相移的变化
(b)HSR随检波相移的变化
图6是本发明实施例提供的伴生调幅指数下最小二乘与卡尔曼滤波的对比示意图;
图中:(a)RAE随伴生调幅指数的变化
(b)HSR随伴生调幅指数的变化
(c)相位噪声随伴生调幅指数的变化
图7是本发明实施例提供的两路载波幅度不等时最小二乘与卡尔曼滤波的对比示意图;
图中:(a)RAE随Gp/Hp的变化
(b)HSR随Gp/Hp的变化
(c)相位噪声随Gp/Hp的变化
图8是本发明实施例提供的检波相移下最小二乘与卡尔曼滤波的对比示意图;
图中:(a)RAE随检波相移的变化
(b)HSR随检波相移的变化
(c)相位噪声随检波相移的变化
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
本发明构建了非线性误差校正问题的观测模型,通过直接解调、最小二乘算法与卡尔曼滤波之间的对比,验证了本发明在PGC解调过程中对非线性误差校正的有效性。
下面结合附图对本发明的应用原理作详细的描述。
如图1所示,本发明实施例提供的基于卡尔曼滤波的椭圆拟合校正方法包括以下步骤:
S101:在大量待拟合数据点中随机选取五个数据点,代入椭圆方程中,构建出一个五元代数方程组,方程组的解即为椭圆的五个参数,将这组参数作为椭圆参数初始值;
S102:选取椭圆参数估计值的偏差作为状态向量,初始值为[0,0,0,0,0]T,选取无限大的对角矩阵作为状态协方差矩阵的初值;
S103:将待拟合的数据点依次利用卡尔曼滤波算法,不断的更新状态向量和状态协方差矩阵;
S104:每一时刻状态的更新都由前一时刻的估计和这一时刻的观测数据计算得到,每一步都将新更新的参数偏差作为下一时刻参数偏差的预测值,将新更新的协方差矩阵作为下一时刻的误差协方差矩阵;
S105:将更新获得的新参数的偏差值与初始的参数估计值进行相加,得到椭圆参数新的估计结果;
S106:不断的重复执行S105~S105,每次的初始偏差为[0,0,0,0,0]T,初始参数为上次输出的最优偏差与上次的初始参数之和,误差协矩阵为上依次循环最终输出的误差协方差矩阵;当更新过程中相邻两次估计的结果变化均小于给定误差10-8时认为达到收敛状态,停止迭代;此时的估计值作为最优估计值,即最优的椭圆参数。
下面结合附图对本发明的应用原理作进一步的描述。
1.非线性误差校正算法
1.1PGC解调信号预处理
在实际情况中,如果直接调制干涉信号,在光源频率发生变化的基础上同时叠加一个变化的光功率项,同时本地载波与相位载波可能存在相移,使得干涉信号形式变为:
(1+mcos(ω0t))即为光强波动产生的寄生调幅项,m是伴生调幅指数。
同时两路载波幅度GP≠HP不相等,这样解调输出是两路信号变成非正交的,其形式表示为:
式中:h,k为信号的直流分量大小,a,b为信号的交流幅度大小,为待解调的相位值,δ为信号间的相位差。
非线性误差的存在导致输出两路信号中交流幅度大小不等、直流偏移变化量也不相等,同时两路信号之间也会产生相位差。实验过程中选取调制深度固定为2.36,不考虑其对系统的影响,则误差来源主要是伴生调幅指数m、两路载波幅值GP与HP的偏差程度以及本地载波与相位载波的相位差这些因素的存在会使得解调出来的相位精度变差。想要提高测量精度,必须将非正交的信号变为正交的两路信号,再通过反正切解调出相位信号。
不存在非线性误差时,利用两路解调信号绘制的图形是一个圆,当非线性误差存在时两路信号绘制的图形为椭圆形,并且在该椭圆方程中,h,k表示椭圆的中心位置,a,b为椭圆长轴短轴大小,δ为两组曲线的固定相位差。非线性误差的存在使圆发生了一些变换,h,k不再为0,a,b不再为1。与传统方法不同,本发明在于对解调之后非正交信号进行正交化处理,具体实现过程如图2所示。在非线性校正部分如果能够求出椭圆的五个未知参数,便可转化为理想情况下圆的方程。
根据非线性误差存在时解调信号的一般表达式,要校正其中存在的非线性误差,首先需要对信号进行预处理:
对公式(2)进行整理可以得到:
对式(3)进行平方和运算,进一步整理可得:
与标准的椭圆方程对比:
x2+Bxy+Cy2+Dx+Ey+F=0 (5)
可以得到椭圆参数h,k,a,b,δ与标准椭圆参数之间的对应关系:
利用式(3)便可求得两路正交的信号。然后利用反正切方法,即可求解出被测相位的大小。
根据上面的推导,选择椭圆匹配算法[11]完成非线性误差的校正,这样的重点转化为求解椭圆的参数B,C,D,E,F。椭圆拟合校正是一种非常实用的统计算法,它利用某种准则消除畸变的参数,算法设计简单、编程容易,结果稳定可靠已经被广泛应用于很多场合。
1.2基于最小二乘的校正算法
估计椭圆参数的方法很多,最小二乘的原理是用方程观测值与真实值差值的平方和最小来建立残差方程,通过偏导数值为0转化为正规方程,进一步对未知的椭圆模型参数B,C,D,E,F进行估计。
将实测的干涉信号I解调输出的两路非正交信号(Ixi,Iyi)代入椭圆方程表达式中,对应最小二乘的残差表达式为:
vi=Ixi 2+BIxiIyi+CIyi 2+DIxi+EIyi+F (7)
vi=0表示满足椭圆方程,偏差为0。当估计参数代入椭圆方程后会产生一定的偏差,定义为残差vi≠0。
令L=[v1,v2,···,vi]T,M=[Ix1 2,Ix2 2,···,Ixn 2]T,待估计的参数为N=[B,C,D,E,F]T,同时令:
则残差方程的矩阵形式为:
L=M+A×N (9)
利用最小二乘原理,等精度测量时选取目标函数为:
通过J对B,C,D,E,F分别求偏导数,保证残差的平方和最小:
得到正规方程组,并令结果为0,满足下面的公式:
AT×L=0结合式(9)可得:
N=(ATA)-1ATM (13)
借助矩阵进行求解,得到如下矩阵形式:
求解上面的矩阵方程,便可以得到关于椭圆参数的线性方程组,从而解出椭圆参数B,C,D,E,F,进而由上面公式(6)推导出h,k,a,b,δ。达到非线性误差的校正目的,完成解调。
1.3基于卡尔曼滤波的校正算法
卡尔曼滤波可以用来进行参数的估计,使用之前需要构建系统模型,系统模型建立的越准确最终估计值就越精确。椭圆及一般二次曲线的描述如下式所示:
f(x,y)=ax2+2bxy+cy2+2dx+2ey+f=0 (15)
其中对于椭圆方程必须满足a+c>0。对椭圆方程归一化需满足条件:
a+c=1;
表达式归一化之后可以避免求解得到双曲线,同时表达式中没有利用常规的f=1进行椭圆中心化,因为f=1时会导致椭圆每个参数都很小,容易引起计算误差,同时f=1时不能包含经过原点的椭圆。将椭圆参数用一个五维向量来表示:
X=(a,b,d,e,f)' (16)
对这个向量进行估计。估计用到的观测值为系统的待拟合点,用如下方法表示:
Ni=(xi,yi)' (17)
观测数据与真实数据的偏差为:
vi=Yi-Ni (18)
其中:vi为高斯噪声,其对应的噪声的协方差矩阵为:
cov[vi]=Ri (19)
对于每个待拟合点来说,椭圆方程就是一个非线性的约束,将待测参数与待估计点代入作为观测方程:
F(Xi,Yi)=0 (20)
F(X,Y)为代入椭圆真实参数时的值,偏差为0。当给出估计的参数M和观测的数据N,将待拟合的数据点代入估计的椭圆方程中得到的偏差为F(M,N)≠0,存在一定的偏差。
卡尔曼滤波应用于非线性系统时需要对信号进行线性化处理。
对F(Xi,Yi)=0在椭圆参数Xi的估计值Mi和Yi的观测数据Ni处利用泰勒级数对非线性函数进行线性化处理,得到:
即:
F(Mi,Ni)可以表示将估计参数Mi和观测值Ni代入式,得到的偏离椭圆的距离,假设Mi=[aa,bb,dd,ee,ff]'对于系统而言-F(Mi,Ni)还可以表示如下:
-F(Mi,Ni)=-(aaxi 2+2bbxiyi+ccyi 2+2ddxi+2eeyi+ff) (23)
其中通过求导可得:
正常卡尔曼滤波的观测方程为Z=HX+v,上式中令ΔX=X-M则可得到:
偏差观测值:
Z=-F(Mi,Ni) (26)
观测矩阵:
观测噪声:
观测噪声期望和协方差矩阵为:
利用卡尔曼滤波进行校正是以椭圆估计的参数偏差ΔX为状态向量,对偏差进行预测然后更新修正,不断迭代求出最优的偏差值,估计的最优偏差加上初始的参数便可得到最优的估计参数。系统中没有状态方程所以不包含状态转移矩阵。估计过程中直接把上一时刻的更新值作为下一时刻的预测值。观测方程利用上述的方法进行构建。将每次更新后的参数代入椭圆方程,与真实椭圆的偏差作为系统观测值,对应标准卡尔曼滤波的公式,具体步骤可以写成如下的形式:
1、系统建模:
(1)状态方程:
ΔX(k)=A×ΔX(k-1)+w(k-1) (30)
系统不包含状态方程,可以视作A=1,w(k-1)=0;
(2)观测方程:
Z(k)=-F(Mi,Ni)=-(aaxi 2+2bbxiyi+ccyi 2+2ddxi+2eeyi+ff) (31)
2、初始化:
(1)状态变量的初始化:
ΔX(0)=[0,0,0,0,0]T (32)
(2)状态协方差的初始化:
选取的参数偏差初值对于系统影响很小,都会在后期不断更新修正,为了方便,选取[0,0,0,0,0]T作为状态向量的初值;同时初始参数与真实参数之间的偏差可能性很大,选择初值为1的无限大的对角矩阵作为状态协方差矩阵的初值;
3、预测过程:
(1)状态预测:
Xekf(k|k-1)=A×Xekf(k-1);
其中状态转移矩阵A=1,上式即为:
Xekf(k|k-1)=Xekf(k-1) (34)
卡尔曼滤波的过程中,系统状态保持不变,没有状态方程,直接将上一时刻的估计值作为这一时刻的一步预测值;
(2)状态协方差预测矩阵:
Pekf(k|k-1)=A×Pekf(k-1)×A'+Q(k-1);
其中A=1,w=0,Q是过程噪声w的协方差矩阵,所以Q=0,等效为:
Pekf(k|k-1)=Pekf(k-1) (35)
4、更新过程:
(1)增益矩阵的更新方程:
Kg(k)=Pekf(k|k-1)×H'/(H×Pekf(k|k-1)×H'+R) (36)
(2)状态估计更新方程:
Xekf(k)=Xekf(k|k-1)+Kg(k)×(Z(k)-H×Xekf(k|k-1)) (37)
(3)误差协方差更新方程:
Pekf(k)=(I-Kg(k)×H)×Pekf(k|k-1) (38)
再将k时刻的状态估计值和误差协方差估计值送入下一时刻,进入不断迭代的过程。
综合各方面因素,卡尔曼滤波整体实现的算法步骤可概括如下:
(1)在大量待拟合数据点中随机选取五个数据点,代入椭圆方程中,构建出一个五元代数方程组,方程组的解即为椭圆的五个参数,将这组参数作为椭圆参数初始值;
(2)选取椭圆参数估计值的偏差作为状态向量,初始值为[0,0,0,0,0]T,初始偏差的可能性很多,因此选取无限大的对角矩阵作为状态协方差矩阵的初值;
(3)将待拟合的数据点依次利用卡尔曼滤波算法,不断的更新状态向量和状态协方差矩阵;
(4)每一时刻状态的更新都由前一时刻的估计和这一时刻的观测数据计算得到。每一步都将新更新的参数偏差作为下一时刻参数偏差的预测值,将新更新的协方差矩阵作为下一时刻的误差协方差矩阵;
(5)将更新获得的新参数的偏差值与初始的参数估计值进行相加,得到椭圆参数新的估计结果;
(6)不断的重复执行步骤(3)~(5),每次的初始偏差为[0,0,0,0,0]T,初始参数为上次输出的最优偏差与上次的初始参数之和,误差协矩阵为上依次循环最终输出的误差协方差矩阵。当更新过程中相邻两次估计的结果变化均小于给定误差10-8时认为达到收敛状态,停止迭代。将此时的估计值作为最优估计值,即最优的椭圆参数。
利用上述卡尔曼滤波算法得到椭圆参数之后,便可以得到两路正交的信号,完成非线性误差的校正。
下面结合实验对本发明的应用效果作详细的描述。
1.实验部分
1.1卡尔曼滤波与直接解调
为了验证算法的可行性,与最小二乘算法相同,对影响正交性的三个因素:伴生调幅指数m、两路载波幅值偏差程度和本地载波与相位载波的相位差分别制定以下三个实验。
1、伴生调幅指数m≤0.3的范围内各项指标能否达到要求,并与不加算法直接解调时进行对比,验证卡尔曼滤波算法在伴生调幅指数m增大时对非线性误差的校正效果。
仿真中,信号幅度DP=3,直流偏移AP=0,随机噪声为0,伴生调幅指数m=0:0.05:0.3,两路载波幅度相等GP=HP,本地载波与相位载波的相位差仿真结果如图3所示。
从图3中可以看出,可以看出,同一伴生调幅指数m时,卡尔曼滤波的同频相对幅度和谐波抑制比均优于直接解调。随着伴生调幅指数m的增大,直接解调的同频相对幅度误差逐渐增大,谐波抑制逐渐减小,两个指标变化范围较大,m较大时失真程度较大。利用卡尔曼滤波校正之后,随着伴生调幅指数的增大,指标保持稳定状态。卡尔曼滤波算法对伴生调幅带来的误差可以进行很好的校正。
2、两路载波幅度不同时各项指标能否达到要求,验证卡尔曼滤波算法在两路载波幅度差增大时对非线性误差的校正效果。
仿真中,信号幅度DP=3,伴生调幅指数m=0,随机噪声为0,直流偏移AP=0,两路载波幅度偏差GP/HP=1:10,本地载波与相位载波的相位差仿真结果图如图4所示:
从图4中可以看出,同一载波幅值差距时卡尔曼滤波的同频相对幅度和谐波抑制比均优于直接解调,随着两路载波幅值偏差程度不断增大,直接解调的同频相对幅度误差逐渐增大,谐波抑制逐渐减小,两个指标的变化范围较大,幅度偏差较大时失真程度较大。利用卡尔曼滤波校正之后,随着伴生调幅指数的增大,两个指标基本保持稳定状态。相位噪声两种方法解调输出的相差不大,卡尔曼滤波算法校正了非线性误差。
3、本地载波与相位载波的相位差对系统指标的影响,过程中选取相位差为0~π/4,验证卡尔曼滤波算法在不同相位差时的校正效果。
仿真中,信号幅度DP=3,直流偏移AP=0,伴生调幅指数m=0,两路载波幅度相等GP=HP,本地载波与相位载波的相位差仿真结果如图5所示;从图5可以看出,同一载波相移时卡尔曼滤波的同频相对幅度和谐波抑制比均优于直接解调。随着载波相移的增大,直接解调的同频相对幅度误差逐渐增大,谐波抑制逐渐减小,两个指标变化范围较大,本地载波与相位载波的相位差较大时不能满足系统指标要求。利用卡尔曼滤波校正之后,随着伴生调幅指数的增大,三个指标保持稳定状态。卡尔曼滤波算法校正了载波相移引起的非线性误差。
1.2卡尔曼滤波与最小二乘
1、伴生调幅指数m≤0.3的范围内,最小二乘算法和卡尔曼滤波两种方法校正下各项指标的变化情况。
仿真中,信号幅度DP=3,直流偏移AP=0,随机噪声为0,伴生调幅指数m=0:0.05:0.3,两路载波幅度相等GP=HP,本地载波与相位载波的相位差仿真如图6所示:
由图6可以看出,最小二乘的同频相对幅度和谐波抑制比在伴生调幅指数变化时始终保持稳定状态,卡尔曼滤波有波动,最小二乘的稳定性优于卡尔曼滤波。相位噪声两种方法相差不大。
2、载波幅度偏差程度不同时,最小二乘算法和卡尔曼滤波两种方法校正下各项指标的变化情况。
仿真中,信号幅度DP=3,伴生调幅指数m=0,随机噪声为0,直流偏移AP=0,两路载波幅度不等GP/HP=1:10,本地载波与相位载波的相位差仿真图为图7所示。
由图7可以看出,最小二乘法整体来看更加稳定,对误差抵抗性能较好,相同的载波幅值偏差时,卡尔曼滤波的同频相对幅度和谐波抑制比优于最小二乘,卡尔曼滤波算法会使校正效果更优。
3、本地载波与相位载波的相位差变化时,最小二乘算法和卡尔曼滤波两种方法校正下各项指标的变化情况。
仿真中,信号幅度DP=3,直流偏移AP=0,伴生调幅指数m=0.1,两路载波幅度相等GP=HP,本地载波与相位载波的相位差两种方法的对比如图8所示:
从图8中可以看出,最小二乘算法的同频相对幅度和谐波抑制比卡尔曼滤波稳定,当检波相移逐渐增大时,卡尔曼滤波的同频相对幅度和谐波抑制比指标逐渐变好,优于最小二乘。
1.3对比结果分析
从上面的实验分析,最小二乘算法较卡尔曼滤波明显的优点是稳定性好,卡尔曼滤波较最小二乘在载波幅度和检波相移发生变化时同频相对幅度和谐波抑制比较好,具体根据实际系统的影响因素选择不同的算法。
对于最小二乘拟合算法而言,待拟合数据点较少时易受偏差较大点的影响,当待测数据非常多时,会使误差较大的测量点对拟合曲线的精度影响减小,提高数据处理的精确度。本次实验选取了400k个数据点进行拟合,利用计算机解决矩阵问题,实现起来比较简单,易于编程,结果稳定。同时最小二乘对变量的约束少,灵活性很强可以使用于多种场合。但是拟合只有保证在大量数据点时才能准确,增加了运算量,而且伴随着数据之间的运算可能在系统中引入本不存在的误差。另外,最小二乘法直接令残差平方和作为目标函数,令其最小然后求解得到参数,是一次性的估计过程。
卡尔曼滤波与最小二乘法相比,是一种迭代法,先给定一个估计值,在若干次迭代更新之后逼近最优估计。进行系统状态估计时不仅考虑动态测量值,还考虑系统内部的动态变化,将预测值与测量值做权重和,使测量误差平方和的期望值最小。本质上是一个递推反馈算法,重视的是预测值与观测值的最优融合问题。当更新过程中时间固定的话,卡尔曼滤波相当于经典的最小二乘法。算法过程中涉及到的中间量比较少只需要存储上一时刻的值,占用内存少,数据处理效率高,是一种实时性很强的校准算法。卡尔曼滤波的缺点在于具有很多的限制条件,需要事先知道信号和噪声的统计特性,这在实际应用中是很难满足的,实验过程中一般选择的是高斯随机噪声;构建的系统模型的精确度也会直接影响校正效果,在实际情况中,系统模型往往是变化的,当模型发生变化时容易导致卡尔曼滤波器的滤波发散,准确性变差。
本发明提出基于卡尔曼滤波的椭圆拟合的校正算法,利用卡尔曼滤波法计算出各个校正参数。实验结果表明:较直接解调和最小二乘方法,卡尔曼滤波的同频相对幅度和谐波抑制比有了很大的改善。本发明可以有效校正非线性误差,提高解调的精度。
以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于卡尔曼滤波的椭圆拟合校正方法,其特征在于,所述基于卡尔曼滤波的椭圆拟合校正方法包括:
步骤一,随机选取五个数据点,代入椭圆方程中,构建出一个五元代数方程组,方程组的解为椭圆的五个参数,作为椭圆参数初始值;
步骤二,选取椭圆参数估计值的偏差作为状态向量,初始值为[0,0,0,0,0]T,选取无限大的对角矩阵作为状态协方差矩阵的初值;
步骤三,将待拟合的数据点依次利用卡尔曼滤波算法,不断的更新状态向量和状态协方差矩阵;
步骤四,每一时刻状态的更新都由前一时刻的估计和这一时刻的观测数据计算得到,每一步都将新更新的参数偏差作为下一时刻参数偏差的预测值,将新更新的协方差矩阵作为下一时刻的误差协方差矩阵;
步骤五,将更新获得的新参数的偏差值与初始的参数估计值进行相加,得到椭圆参数新的估计结果;
步骤六,不断的重复执行步骤三~步骤五,每次的初始偏差为[0,0,0,0,0]T,初始参数为上次输出的最优偏差与上次的初始参数之和,误差协矩阵为上依次循环最终输出的误差协方差矩阵;当更新过程中相邻两次估计的结果变化均小于给定误差10-8时认为达到收敛状态,停止迭代;估计值作为最优估计值,即最优的椭圆参数。
2.如权利要求1所述的基于卡尔曼滤波的椭圆拟合校正方法,其特征在于,所述基于卡尔曼滤波的椭圆拟合校正方法的椭圆及一般二次曲线的描述为:
f(x,y)=ax2+2bxy+cy2+2dx+2ey+f=0;
其中对于椭圆方程满足a+c>0,对椭圆方程归一化满足条件:
a+c=1;
椭圆参数用一个五维向量来表示:
X=(a,b,d,e,f)';
估计用到的观测值为系统的待拟合点表示:
Ni=(xi,yi)';
观测数据与真实数据的偏差为:
vi=Yi-Ni
其中:vi为高斯噪声,其对应的噪声的协方差矩阵为:
cov[vi]=Ri
3.如权利要求2所述的基于卡尔曼滤波的椭圆拟合校正方法,其特征在于,所述基于卡尔曼滤波的椭圆拟合校正方法的椭圆方程是一个非线性的约束,待测参数与待估计点代入作为观测方程:
F(Xi,Yi)=0;
F(X,Y)为代入椭圆真实参数时的值,偏差为0;当给出估计的参数M和观测的数据N,将待拟合的数据点代入估计的椭圆方程中得到的偏差为F(M,N)≠0,存在一定的偏差。
4.如权利要求1所述的基于卡尔曼滤波的椭圆拟合校正方法,其特征在于,所述基于卡尔曼滤波的椭圆拟合校正方法的卡尔曼滤波对信号进行线性化处理:对F(Xi,Yi)=0在椭圆参数Xi的估计值Mi和Yi的观测数据Ni处利用泰勒级数对非线性函数进行线性化处理,得到:
F(Xi,Yi)=F(Mi,Ni)+[▽XF(Mi,Ni)'(Xi-Ni)]+[▽YF(Mi,Ni)'(Yi-Ni)]=0;
-F(Mi,Ni)=▽XF(Mi,Ni)'(Xi-Ni)]+[▽YF(Mi,Ni)'(Yi-Ni);
F(Mi,Ni)表示将估计参数Mi和观测值Ni代入式,得到的偏离椭圆的距离,Mi=[aa,bb,dd,ee,ff]'对于系统-F(Mi,Ni)表示如下:
-F(Mi,Ni)=-(aaxi 2+2bbxiyi+ccyi 2+2ddxi+2eeyi+ff);
其中通过求导得:
XF(Mi,Ni)=(xi 2-yi 2,2xiyi,2xi,2yi,1)';
YF(Mi,Ni)=2(aaxi+bbyi+dd,bbxi+ccyi+ee)';
正常卡尔曼滤波的观测方程为Z=HX+v,上式中令ΔX=X-M则得到:
偏差观测值:
Z=-F(Mi,Ni);
观测矩阵:
H=▽XF(Mi,Ni);
观测噪声:
vi=▽YF(Mi,Ni)'(Yi-Ni);
观测噪声期望和协方差矩阵为:
E[vi]=▽YF(Mi,Ni)'E[(Yi-Ni)]=0;
var[vi]=[▽YF(Mi,Ni)]R[▽YF(Mi,Ni)]'
=4[(aaxi+bbyi+dd)2+(bbxi+ccyi+ee)2]∑2I。
5.如权利要求1所述的基于卡尔曼滤波的椭圆拟合校正方法,其特征在于,所述基于卡尔曼滤波的椭圆拟合校正方法的将每次更新后的参数代入椭圆方程,与真实椭圆的偏差作为系统观测值,对应标准卡尔曼滤波的公式,具体步骤位:
步骤一,系统建模:
(1)状态方程:
ΔX(k)=A×ΔX(k-1)+w(k-1);
系统不包含状态方程,视作A=1,w(k-1)=0;
(2)观测方程:
Z(k)=-F(Mi,Ni)=-(aaxi 2+2bbxiyi+ccyi 2+2ddxi+2eeyi+ff);
步骤二,初始化:
(1)状态变量的初始化:
ΔX(0)=[0,0,0,0,0]T
(2)状态协方差的初始化:
选取[0,0,0,0,0]T作为状态向量的初值;选择初值为1的无限大的对角矩阵作为状态协方差矩阵的初值;
步骤三,预测过程:
(1)状态预测:
Xekf(k|k-1)=A×Xekf(k-1);
其中状态转移矩阵A=1,上式为:
Xekf(k|k-1)=Xekf(k-1);
直接将上一时刻的估计值作为这一时刻的一步预测值;
(2)状态协方差预测矩阵:
Pekf(k|k-1)=A×Pekf(k-1)×A'+Q(k-1);
其中A=1,w=0,Q是过程噪声w的协方差矩阵,Q=0,等效为:
Pekf(k|k-1)=Pekf(k-1);
步骤四,更新过程:
(1)增益矩阵的更新方程:
Kg(k)=Pekf(k|k-1)×H'/(H×Pekf(k|k-1)×H'+R);
(2)状态估计更新方程:
Xekf(k)=Xekf(k|k-1)+Kg(k)×(Z(k)-H×Xekf(k|k-1));
(3)误差协方差更新方程:
Pekf(k)=(I-Kg(k)×H)×Pekf(k|k-1);
再将k时刻的状态估计值和误差协方差估计值送入下一时刻,进入不断迭代的过程。
6.一种利用权利要求1~5任意一项所述基于卡尔曼滤波的椭圆拟合非线性误差校正方法的光纤水听器。
CN201811128365.7A 2018-09-27 2018-09-27 一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法 Pending CN109000782A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811128365.7A CN109000782A (zh) 2018-09-27 2018-09-27 一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811128365.7A CN109000782A (zh) 2018-09-27 2018-09-27 一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法

Publications (1)

Publication Number Publication Date
CN109000782A true CN109000782A (zh) 2018-12-14

Family

ID=64589430

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811128365.7A Pending CN109000782A (zh) 2018-09-27 2018-09-27 一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法

Country Status (1)

Country Link
CN (1) CN109000782A (zh)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111458395A (zh) * 2020-06-18 2020-07-28 北京英视睿达科技有限公司 一种改变q值的卡尔曼滤波方法、装置及可读存储介质
CN112305496A (zh) * 2020-10-26 2021-02-02 哈尔滨工程大学 一种被动测向通道相位校正方法
CN112737547A (zh) * 2020-12-18 2021-04-30 北京邮电大学 一种基于卡尔曼滤波的信号补偿方法及装置
CN112967256A (zh) * 2021-03-09 2021-06-15 扬州大学 一种基于空间分布的隧道椭圆化检测方法
CN113325452A (zh) * 2021-05-25 2021-08-31 哈尔滨工程大学 一种三星无源融合定位体制机动目标跟踪方法
CN114323092A (zh) * 2021-12-28 2022-04-12 中国人民解放军国防科技大学 一种计算与消除内调制pgc信号检测中伴生调幅的方法
CN114820721A (zh) * 2022-05-17 2022-07-29 苏州轻棹科技有限公司 一种卡尔曼滤波观测噪声的可视化调制方法和装置
CN115511979A (zh) * 2022-10-14 2022-12-23 南通伊诺精密塑胶导管有限公司 一种基于灰度变换的内窥镜主板控制系统
CN117201250A (zh) * 2023-11-07 2023-12-08 武汉能钠智能装备技术股份有限公司 一种相位生成载波调解方法、装置、电子设备及存储介质
CN117310455A (zh) * 2023-11-30 2023-12-29 上海采起电子技术服务有限公司 一种基于电信号的口腔扫描仪电路板故障检测方法
CN118362198A (zh) * 2024-06-20 2024-07-19 中山大学 一种分布式声波传感相位漂移消除方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110142282A1 (en) * 2009-12-14 2011-06-16 Indian Institute Of Technology Bombay Visual object tracking with scale and orientation adaptation
US20160097788A1 (en) * 2014-10-07 2016-04-07 Snappafras Corp. Pedestrian direction of motion determination system and method
CN106936513A (zh) * 2017-03-22 2017-07-07 北京邮电大学 一种基于卡尔曼滤波算法的载波相位恢复方法及装置
CN107707310A (zh) * 2017-09-20 2018-02-16 哈尔滨工业大学深圳研究生院 一种基于自适应卡尔曼的偏振解复用和载波相位恢复方法
CN107843189A (zh) * 2017-09-30 2018-03-27 浙江理工大学 正弦相位调制干涉仪pgc解调实时归一化修正装置及方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110142282A1 (en) * 2009-12-14 2011-06-16 Indian Institute Of Technology Bombay Visual object tracking with scale and orientation adaptation
US20160097788A1 (en) * 2014-10-07 2016-04-07 Snappafras Corp. Pedestrian direction of motion determination system and method
CN106936513A (zh) * 2017-03-22 2017-07-07 北京邮电大学 一种基于卡尔曼滤波算法的载波相位恢复方法及装置
CN107707310A (zh) * 2017-09-20 2018-02-16 哈尔滨工业大学深圳研究生院 一种基于自适应卡尔曼的偏振解复用和载波相位恢复方法
CN107843189A (zh) * 2017-09-30 2018-03-27 浙江理工大学 正弦相位调制干涉仪pgc解调实时归一化修正装置及方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
许洪飞: "多椭圆检测算法研究与实现", 《中国学位论文全文数据库》 *
陈明建: "卡尔曼滤波拟合椭圆在地铁隧道断面监测中的应用", 《市政技术》 *

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111458395B (zh) * 2020-06-18 2020-09-22 北京英视睿达科技有限公司 一种改变q值的卡尔曼滤波方法、装置及可读存储介质
CN111458395A (zh) * 2020-06-18 2020-07-28 北京英视睿达科技有限公司 一种改变q值的卡尔曼滤波方法、装置及可读存储介质
CN112305496A (zh) * 2020-10-26 2021-02-02 哈尔滨工程大学 一种被动测向通道相位校正方法
CN112305496B (zh) * 2020-10-26 2022-06-17 哈尔滨工程大学 一种被动测向通道相位校正方法
CN112737547A (zh) * 2020-12-18 2021-04-30 北京邮电大学 一种基于卡尔曼滤波的信号补偿方法及装置
CN112967256B (zh) * 2021-03-09 2023-11-24 扬州大学 一种基于空间分布的隧道椭圆化检测方法
CN112967256A (zh) * 2021-03-09 2021-06-15 扬州大学 一种基于空间分布的隧道椭圆化检测方法
CN113325452A (zh) * 2021-05-25 2021-08-31 哈尔滨工程大学 一种三星无源融合定位体制机动目标跟踪方法
CN114323092A (zh) * 2021-12-28 2022-04-12 中国人民解放军国防科技大学 一种计算与消除内调制pgc信号检测中伴生调幅的方法
CN114820721A (zh) * 2022-05-17 2022-07-29 苏州轻棹科技有限公司 一种卡尔曼滤波观测噪声的可视化调制方法和装置
CN114820721B (zh) * 2022-05-17 2024-03-26 苏州轻棹科技有限公司 一种卡尔曼滤波观测噪声的可视化调制方法和装置
CN115511979A (zh) * 2022-10-14 2022-12-23 南通伊诺精密塑胶导管有限公司 一种基于灰度变换的内窥镜主板控制系统
CN117201250A (zh) * 2023-11-07 2023-12-08 武汉能钠智能装备技术股份有限公司 一种相位生成载波调解方法、装置、电子设备及存储介质
CN117201250B (zh) * 2023-11-07 2024-01-23 武汉能钠智能装备技术股份有限公司 一种相位生成载波解调方法、装置、电子设备及存储介质
CN117310455A (zh) * 2023-11-30 2023-12-29 上海采起电子技术服务有限公司 一种基于电信号的口腔扫描仪电路板故障检测方法
CN117310455B (zh) * 2023-11-30 2024-02-09 上海采起电子技术服务有限公司 一种基于电信号的口腔扫描仪电路板故障检测方法
CN118362198A (zh) * 2024-06-20 2024-07-19 中山大学 一种分布式声波传感相位漂移消除方法

Similar Documents

Publication Publication Date Title
CN109000782A (zh) 一种基于卡尔曼滤波的椭圆拟合非线性误差校正方法
Carli et al. Quantized average consensus via dynamic coding/decoding schemes
Bolker Functions resembling quotients of measures
US8358719B2 (en) Method and device for the compensation of signal errors in IQ-modulators
CN110307780B (zh) 基于迭代计算的pgc相位解调误差实时补偿方法
US8761238B2 (en) Method and apparatus for correcting frequency offset
CN104076332B (zh) 一种雷达均匀线性阵列幅度和相位的估计方法
CN109100816B (zh) 一种重磁数据处理方法及系统
Paláncz et al. Application of Pareto optimality to linear models with errors-in-all-variables
Anttila et al. Blind signal estimation in widely-linear signal models with fourth-order circularity: Algorithms and application to receiver I/Q calibration
CN103259585B (zh) 基于收发机损耗的下行链路波束成形方法及其系统
CN104596651B (zh) 一种基于四象限二元相位调制的相位反演方法
CN108845975A (zh) 一种基于阵列天线的相位恢复方法
Xia et al. A regularised normalised augmented complex least mean square algorithm
CN109917658B (zh) 一种多径衰落信道下非线性网络化系统滤波器及设计方法
Stratonovich The quantum generalization of optimal statistical estimation and hypothesis testing
US10419118B2 (en) Optical transmission device and optical transmission method
Das et al. Distributed linear estimation of dynamic random fields
CN106209706B (zh) 一种混沌背景下弱谐波信号的检测方法
CN113055323B (zh) 一种通信系统的数字预失真处理的方法及系统
CN113054948A (zh) 并行卡尔曼滤波系统和方法
Saito et al. Discovery of relevant weights by minimizing cross-validation error
JP2022082435A (ja) 非線形通信システムの性能推定装置、方法及び電子機器
Cavalcante et al. Consensus acceleration in multiagent systems with the Chebyshev semi-iterative method
CN112787728B (zh) 模拟相干信号处理系统和方法

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20181214

RJ01 Rejection of invention patent application after publication