CN117480356A - 基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量结构位移的方法及其系统 - Google Patents
基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量结构位移的方法及其系统 Download PDFInfo
- Publication number
- CN117480356A CN117480356A CN202280041438.7A CN202280041438A CN117480356A CN 117480356 A CN117480356 A CN 117480356A CN 202280041438 A CN202280041438 A CN 202280041438A CN 117480356 A CN117480356 A CN 117480356A
- Authority
- CN
- China
- Prior art keywords
- feature point
- displacement
- image frame
- roi
- acceleration data
- 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
Links
- 238000006073 displacement reaction Methods 0.000 title claims abstract description 216
- 230000001133 acceleration Effects 0.000 title claims abstract description 117
- 238000000034 method Methods 0.000 title claims abstract description 86
- 238000005259 measurement Methods 0.000 title claims abstract description 81
- 230000004927 fusion Effects 0.000 title claims abstract description 21
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 69
- 238000005070 sampling Methods 0.000 claims abstract description 42
- 230000003044 adaptive effect Effects 0.000 claims abstract description 33
- 238000013519 translation Methods 0.000 claims description 50
- 239000011159 matrix material Substances 0.000 claims description 31
- 230000000007 visual effect Effects 0.000 claims description 21
- 230000008859 change Effects 0.000 claims description 20
- 238000004364 calculation method Methods 0.000 claims description 18
- 238000004590 computer program Methods 0.000 claims description 11
- 238000001914 filtration Methods 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 6
- 238000013459 approach Methods 0.000 claims description 3
- 230000001360 synchronised effect Effects 0.000 abstract description 3
- 230000014616 translation Effects 0.000 description 36
- 230000006870 function Effects 0.000 description 7
- 238000012545 processing Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 238000009434 installation Methods 0.000 description 3
- 238000012423 maintenance Methods 0.000 description 3
- 230000015654 memory Effects 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 239000000969 carrier Substances 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000000691 measurement method Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 241000282320 Panthera leo Species 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 238000010420 art technique Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 230000004044 response Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B21/00—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant
- G01B21/32—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant for measuring the deformation in a solid
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/16—Measuring arrangements characterised by the use of optical techniques for measuring the deformation in a solid, e.g. optical strain gauge
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
- G01D21/00—Measuring or testing not otherwise provided for
- G01D21/02—Measuring two or more variables by means not covered by a single other subclass
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
- G01P15/00—Measuring acceleration; Measuring deceleration; Measuring shock, i.e. sudden change of acceleration
- G01P15/02—Measuring acceleration; Measuring deceleration; Measuring shock, i.e. sudden change of acceleration by making use of inertia forces using solid seismic masses
- G01P15/08—Measuring acceleration; Measuring deceleration; Measuring shock, i.e. sudden change of acceleration by making use of inertia forces using solid seismic masses with conversion into electric or magnetic values
- G01P15/097—Measuring acceleration; Measuring deceleration; Measuring shock, i.e. sudden change of acceleration by making use of inertia forces using solid seismic masses with conversion into electric or magnetic values by vibratory elements
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0004—Industrial image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/50—Depth or shape recovery
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Quality & Reliability (AREA)
- Length Measuring Devices By Optical Means (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量所述结构的位移的方法和系统。比例因子(其是图像帧数据的像素分辨率与真实空间的距离分辨率之间的转换因子)是随时间使用由安装在所述结构上的相机以在所述结构外部的静止目标的预定第一采样频率获取的所述图像帧数据,以及由安装在所述结构上的加速度计以预定第二采样频率测量的加速度数据自动计算的。所述结构的位移值是使用改进的特征点匹配算法和所述比例因子进行估计的。所述结构的位移是通过使用自适应多速率卡尔曼滤波器来无缝融合异步的所述图像帧数据和所述加速度数据以使它们以减小的噪声彼此同步,同时使用所述估计的位移值来以所述加速度数据的采样频率进行估计。
Description
技术领域
本发明涉及位移测量技术的领域,并且更特别地涉及一种一起使用加速度数据和视觉测量数据来实时测量结构的位移的方法和系统。
背景技术
位移在民用基础设施的评估、监测和控制中发挥着重要作用,这是因为它提供了与结构完整性和健康相关的关键信息。有几种技术可用于直接位移测量。
大多数位移测量传感器需要在安装期间无法移动的刚性固定支撑点,这限制了它们在大型结构上的使用。线性可变差动变压器(LVDT)是一种活塞形传感器,其具有约50厘米长的细长圆柱形结构。它可以以这样一种方式安装,使得其一端附接到结构,并且另一端附接到固定的支撑点,即,其被放置在结构上的目标点和固定参考点之间以测量相对位移。LVDT可以以高准确性测量位移。然而,当安装在结构,诸如桥上时,往往很难在距测量点50cm的距离处固定支撑点,并且需要额外的安装装置,诸如临时结构,这对于实际应用来说并不方便。
使用实时运动全球导航卫星系统(RTK-GNSS)是测量位移的另一个选项。GPS-RTK是一种不需要安装固定支撑点并且利用两个GPS传感器的传感器。一组GNSS天线和接收器放置在待测结构上的目标点处,而另一组则放置在距结构预定距离内的固定参考点处。两个点处的天线同时接收来自多个卫星的载波,因此可以基于在两个测量点接收的载波之间的相位差计算距目标点的位移。然而,RTK-GNSS通常具有高达20Hz的低采样率,并且其准确性限于约7-10mm。当由传感器接收到的卫星数量少于四颗、天气条件恶劣或GNSS天线的视野被附近移动物体或碎屑阻挡时,测量准确性急剧下降。此外,其无法在GNSS拒止环境中工作。
视觉相机有时用于位移估计。视觉相机基于模板匹配技术、光流技术、特征匹配算法等来跟踪目标结构上的人工或自然目标的移动。在这种基于视觉的技术中,首先根据从视觉测量获得的图像帧以像素单位估计平移,并且然后使用比例因子将估计的平移转换为长度单位的位移。转换的位移对应于在真实空间中行进的长度。
然而,这些现有技术具有以下缺点或限制。可以通过预先识别目标的物理单位的大小来估计比例因子。然而,在现场手动测量目标的大小可能很困难。此外,该方法需要额外的维护,诸如上面安装有视觉相机的固定支撑点,并且可能需要安装单独的人工目标,这增加了成本。现有技术具有相对较高的计算成本,这阻碍了在高采样率下对位移的实时估计。另外,由于视觉相机的相对较低的采样频率,大部分的测量结果都受到噪声的污染,这可能会限制测量的准确性。
另一方面,也提出了基于数据融合的位移估计技术。这些技术结合了来自不同模态的多个测量,以产生更准确的位移估计。示例包括基于RTK-GNSS和加速度计融合的技术(参见下面的现有技术文献(1))以及基于应变计和加速度计融合的技术(参见下面的现有技术文献(2)和(3))。然而,前者的局限性在于它不能在GNSS拒止环境中使用,并且后者需要目标结构的模态形状的先验知识。
还尝试了融合视觉相机和加速度计(参见下面的现有技术文献(4)和(5))。通过多速率卡尔曼滤波器融合从视觉测量获得的低采样率位移和高采样率加速度测量,以估计高采样率位移。然而,如图1所示,这些技术涉及将人工目标放置在待测量的结构上,并且使用安装在远离结构的外部固定点处的视觉相机对人工目标进行成像。因此,该现有技术无法避免图1所示的基于视觉相机的位移测量方法的缺点或限制。也就是说,这些技术仍然需要目标物体的实际长度的先验知识,以便估计比例因子。其还受到限制,因为其只能在视觉测量和加速度测量之间存在同步的条件下使用。
与本发明相关的背景技术包括
(1)Kim,K.、Choi,J.、Koo,G.和Sohn,H.(2016),通过使用两级卡尔曼估计器来融合偏置高采样率加速度和低采样率位移测量进行动态位移估计,智能结构和系统,17(4),647-667。
(2)Park,J.W.、Sim,S.H.和Jung,H.J.(2013),使用多测量数据融合进行位移估计,IEEE/ASME机电技术论文集,18(6),1675-1682。
(3)Park,J.W.、Lee,K.C.、Sim,S.H.、Jung,H.J.和Spencer Jr,B.F.(2016),使用扩展的多传感器数据融合对铁路桥梁进行交通安全评估,计算机辅助土木和基础设施工程,31(10),749-760。
(4)Chang,C.C.和Xiao,X.H.(2010),一种用于结构位移和速度测量的集成视觉惯性技术,智能结构与系统,6(9),1025-1039。
(5)Xu,Y.、Brownjohn,J.M.W.和Huseynov,F.(2019),使用与相机和加速度计相结合的具有成本效益的感测系统对桥梁结构进行准确的变形监测:案例研究,桥梁工程杂志,24(1)。
(6)Smyth,A.和Wu,M.(2007),用于动态系统监测中位移和加速度响应测量的数据融合的多速率卡尔曼滤波,机械系统和信号处理,21(2),706-723。
发明内容
[技术任务]
本发明的目的是提供一种用于基于分别由安装在待测量的结构上的加速度计和视觉相机测量的结构的加速度测量数据与结构外部的自然目标的视频测量数据的融合的以低维护成本和高准确性测量结构的位移的方法和系统。
本发明的另一个目的是提供一种通过融合自然目标的异步视觉测量数据和结构的加速度测量数据来测量结构的位移的方法和系统,其中可以通过使用结构的加速度测量数据自动计算用于平移视觉测量数据的比例因子并且将改进的特征匹配算法应用于视觉测量数据来准确地测量结构的位移。
本发明的又一个目的是提供一种基于自然目标的异步视觉测量与结构的加速度数据的融合来测量结构的位移的方法和系统,其中异步低速采样视觉测量和高速采样加速度测量用自适应多速率卡尔曼滤波器融合,以便以高采样率和提高的准确性测量结构的位移。
本发明所要解决的问题并不限于上述的那些,并且可以在不脱离本发明的精神和范围的情况下以各种方式进行扩展。
[技术方案]
根据用于实现本发明的一个目的的实施例的一种基于自然目标的异步视觉测量数据和结构的加速度数据的融合的测量结构的位移的方法包括:使用随时间由安装在要测量的结构上的相机以在结构外部的静止目标的预定第一采样频率获取的图像帧数据,以及由安装在结构上的加速度计以预定第二采样频率测量的加速度数据,自动计算比例因子(α),比例因子(α)是图像帧数据的像素分辨率与真实空间的距离分辨率之间的转换因子;以及使用改进的特征点匹配算法和比例因子估计结构的位移值,并且通过使用自适应多速率卡尔曼滤波器来无缝融合异步的图像帧数据和加速度数据以使它们以减小的噪声彼此同步,同时使用估计的位移值来以加速度数据的采样频率计算结构的位移。
在示例性实施例中,“自动计算比例因子(α)”可以包括通过使用用于比较参考图像帧与比较图像帧之间的预定特征点的特征点匹配算法来匹配相同的特征点来估计特征点的逐像素平移;以及通过对加速度数据进行二重积分来估计真实空间中的纵长位移。
在示例性实施例中,“自动计算比例因子(α)”还可以包括对平移的获得的估计值和位移的估计值进行带通滤波处理以减少噪声。
在示例性实施例中,在带通滤波处理中,相机的采样频率可以比图像帧的频率分量的最高频率高至少10倍。
在示例性实施例中,“自动计算比例因子(α)”还可以包括对位移的带通滤波估计值进行下采样以使位移的估计值与平移的估计值同步。
在示例性实施例中,比例因子(α)可以是在通过比较平移的带通滤波估计值与位移的带通滤波估计值(uf(t))获得的两个值之间的斜率。
在示例性实施例中,比例因子(α)可以通过将最小二乘法或RANSAC(随机抽样一致)算法应用于带通滤波平移的估计值(df(t))和下采样位移的估计值(uf(t))获得。
在示例性实施例中,“估计结构的位移值”可以包括更新当前第i个图像帧的感兴趣区域(ROI);使用预定的特征点匹配算法将第1个图像帧的ROI(第1个ROI)中的N个特征点与第i个图像帧的感兴趣区域(第i个ROI)内的N个特征点进行匹配;使用失配拒绝算法拒绝特征点失配,以仅留下在N个特征点匹配中良好的特征点匹配,以用于可靠的位移测量;并且通过使用方程
来估计最终的平移,其中Ng表示良好的特征点匹配的数量,其中比例因子(α)被应用于良好特征点匹配的平移和ROI的移动的总和。
在示例性实施例中,“更新第i个图像帧的ROI”可以包括使用比例因子将位移的先验状态估计值从纵长的一个转换为逐像素的一个;通过使用转换的逐像素的先验状态估计值获得第i个图像帧的感兴趣目标(TOI)的逐像素位移(di TOI);使用轮函数获得第i个ROI的平移(di TOI);并且使用第i个ROI的获得的平移(di TOI)更新第i个ROI的位置。
在示例性实施例中,“拒绝特征点失配以仅留下良好的特征点匹配”可以包括使用第一特征点失配拒绝算法来拒绝特征点失配,其中第一特征点失配拒绝算法包括:从特征点匹配选择最佳的特征点匹配;以及基于最佳的特征点匹配通过交叉检查剩余的特征点来拒绝失配的特征点。
在示例性实施例中,“选择最佳的特征点匹配”可以包括:将第1个图像帧的ROI与第i个图像帧的ROI之间的N个特征点匹配两两配对为N/2或(N-1)/2个特征点匹配组;计算N/2或(N-1)/2个特征点匹配组中每一个的距离变化指数;在计算的N/2或(N-1)/2个距离变化指数中选择具有最小距离变化指数(μs)的匹配组作为最佳匹配组;并且选择属于最佳匹配组的两个特征点匹配中的一个作为最佳特征点匹配。
在示例性实施例中,“拒绝失配的特征点”可以包括:基于选择的最佳匹配来计算剩余特征点匹配中每一个的距离间隔,以确定计算的距离间隔是否满足失配条件;并且从特征点失配拒绝满足失配条件的特征点,从而仅留下良好的匹配。
在示例性实施例中,失配条件可以包括第一和第二距离间隔条件,第一间隔条件为“在第j个特征点与第1个图像帧的ROI内的最佳匹配的特征点之间的最大距离”大于或等于“在第j个特征点与第i个图像帧的ROI内的最佳匹配特征点之间的最小距离”,并且第二距离间隔条件为“在第j个特征点与第i个图像帧的ROI内的最佳匹配的特征点之间的最大距离”大于或等于“在第j个特征点与第1个图像帧的ROI内的最佳匹配特征点之间的最小距离”,其中j=1、2、......、N。
在示例性实施例中,距离变化指数可以通过方程计算,其中μs是第s个匹配组的距离变化指数,l1 s是在第1个图像帧的ROI内的两个特征点之间的距离,并且li s是在第i个图像帧的ROI内的两个特征点之间的距离。
在示例性实施例中,“拒绝特征点失配以仅留下良好的特征点匹配”可以包括使用第二特征点失配拒绝算法来拒绝特征点失配,并且其中第二特征点失配拒绝算法包括计算由第一特征点失配拒绝算法进行滤波的剩余特征点匹配中的每一个的平移并且将所计算的平移未落在由
方程表示的平移值的范围内的特征点匹配作为失配进行拒绝。
在示例性实施例中,使用自适应多速率卡尔曼滤波器进行的图像帧数据与加速度数据的融合可以通过以下实现:在类型I时间步骤中执行,其中在当前的加速度数据(ak)与先前的加速度数据(ak-1)之间的时间间隔[(k-2)Δta,(k-1)Δta]中,仅当前的加速度数据(ak)可用并且没有可用的图像帧数据,分别使用第一方程和第二方程获得位移的先验状态估计值/>以及先验状态估计值/>的误差协方差矩阵/>并且在时间(k-1)Δta,获得对应于先验状态估计值/>和先验状态估计值/>的误差协方差矩阵/>的后验状态估计值/>和后验状态估计值/>的后验误差协方差矩阵/>其中Δta表示在加速度数据之间的时间间隔,k表示加速度数据的第k个时间间隔,并且/>和ak-1分别表示在“前一时间步骤”和的后验状态估计值和加速度数据。
在示例性实施例中,使用自适应多速率卡尔曼滤波器进行的图像帧数据与加速度数据的融合可以通过以下实现:在类型II时间步骤中执行,其中在时间间隔[(k-2)Δta,(k-1)Δta]内仅图像帧数据(第i个帧)可用,并且加速度数据不可用,使用第一和第二方程估计位移的先验状态估计值和位移的先验状态估计值/>的先验误差协方差矩阵其中将时间间隔从Δta改变为(i-1)Δtd-(k-1)Δta;使用获得的先验状态估计值及其先验误差协方差矩阵/>更新第i个图像帧的ROI;通过基于第一和第二特征点失配拒绝算法将获得的先验状态估计值/>及其先验误差协方差矩阵/>应用于特征点失配拒绝来估计基于视觉的位移(ui);通过使用位移的获得的先验状态估计值/>其先验误差协方差矩阵/>以及基于视觉的位移(ui)估计值来估计基于视觉的位移(ui)的噪声方差(Ri);使用/>方程和基于视觉的位移(ui)的获得的噪声方差(Ri)来计算卡尔曼增益(K);以及使用方程/>和以及卡尔曼增益(K)来计算在时间t=(i-1)Δtd的位移的后验状态估计值/>及其误差协方差矩阵/>其中在上述方程中,x和P分别由y和G表示,H为向量[1,0],i是图像帧数据的第i个时间步骤,并且Δtd是图像帧数据之间的时间间隔。
在示例性实施例中,使用自适应多速率卡尔曼滤波器进行的图像帧数据与加速度数据的融合可以通过以下实现:在类型III时间步骤中执行,其中在当前加速度数据(ak)与后续加速度数据(ak+1)之间的时间间隔[t-Δta,t]内,加速度数据和图像帧数据两者是可用的,使用
方程获得位移的先验状态估计值在类型I时间步骤中的加速度数据(ak)、在类型II时间步骤中获得的位移的后验状态估计值/>以及后验状态估计值/>的误差协方差矩阵/>使用/>方程获得位移的先验状态估计值/>的误差协方差矩阵/>其中W1=A(kΔta(i-1)Δtd)(I-KH)A((i-1)Δtd-(k-1)Δta)、W2=A(kΔta-(i-1)Δtd)(I-KH)B((i-1)Δtd-(k-1)Δta)+B(kΔta-(i-1)Δtd)和W3=A(kΔta-(i-1)Δtd)K。
同时,根据用于实现本发明的一个目的的实施例的一种基于自然目标的异步视觉测量数据和结构的加速度数据的融合的测量结构的位移的系统包括安装在要测量的结构上的相机,相机被配置为以预定的第一采样频率对结构外部的静止目标进行成像;安装在结构上的相机附近的加速度计,加速度计被配置为以预定的第二采样频率测量加速度;以及位移估计单元,位移估计单元包括比例因子计算单元和自适应多速率卡尔曼滤波器单元,其中比例因子计算单元被配置为随时间分别使用由相机和加速度计测量的一系列图像帧数据和加速度数据来执行自动计算比例因子(α)的功能,比例因子(α)是在图像帧数据的像素分辨率与真实空间的距离分辨率之间的转换因子,并且其中自适应多速率卡尔曼滤波器单元被配置为使用改进的特征点匹配算法和比例因子来执行估计结构的位移值的功能;以及通过使用自适应多速率卡尔曼滤波器来无缝融合异步的加速度数据和图像帧数据,同时使用估计的位移值来以加速度数据的采样频率估计结构的位移。
在示例性实施例中,位移估计单元可以包括实施比例因子计算单元和自适应多速率卡尔曼滤波器单元的计算机程序,以及处理器单元,处理器单元用于执行计算机程序以执行自动计算比例因子(α)的任务和计算结构的位移估计值的任务。
在示例性实施例中,“自动计算比例因子(α)的功能”可以包括通过使用用于比较参考图像帧与比较图像帧之间的预定特征点的特征点匹配算法来匹配相同的特征点来估计特征点的逐像素平移;以及通过对加速度数据进行二重积分来估计真实空间中的纵长位移。
在示例性实施例中,“估计结构的位移值”的功能可以包括:更新当前第i个图像帧的感兴趣区域(ROI);使用预定的特征点匹配算法将第1个图像帧的ROI内的N个特征点与第i个图像帧的ROI内的N个特征点进行匹配;使用失配拒绝算法拒绝特征点失配,以仅留下在N个特征点匹配中良好的特征点匹配,以用于可靠的位移测量;并且通过使用
方程来估计最终的平移估计值,其中Ng表示良好的特征点匹配的数量,其中比例因子(α)被应用于良好特征点匹配的平移和ROI的移动的总和。
[发明效果]
根据本发明的示例性实施例,视觉相机和加速度计安装在待测量结构上,并且利用位于结构外部的自然目标作为视觉相机的目标。不需要用于安装视觉相机的固定支撑点,并且不需要预先测量关于目标的大小以及视觉相机与目标之间的距离的信息。与需要将视觉相机安装在待测量的结构外部的现有技术相比,本发明可以降低位移测量系统的维护成本并且使其更易于实施。
根据本发明的示例性实施例,可以自动计算将视觉测量数据中的位移转换为真实空间中的距离的缩放因子。此外,根据本发明的示例性实施例,可以融合异步加速度测量数据和视觉测量数据以通过使用自适应多速率卡尔曼滤波器以高采样率和准确性来估计结构的位移。
此外,根据本发明的示例性实施例,可以通过更新感兴趣的区域(ROI)来改进连续视觉图像中的自然目标的匹配,并且可以通过使用自动失配拒绝算法来提高特征匹配的准确性,以进一步提高结构的位移测量的准确性。
附图说明
图1示意性地示出了使用视觉相机的传统的基于特征点的位移估计方法;
图2示出了使用特征匹配算法的传统位移估计程序;
图3是示出典型的位移测量技术的一些限制的图;
图4示意性地示出了根据本发明的一个示例性实施例的用于执行结构位移估计方法的系统;
图5示意性地示出了根据本发明的示例性实施例的结构位移估计方法的概述;
图6是示出执行图5所示的比例因子计算步骤的详细过程的流程图;
图7a和图7b是示出根据本发明的一个示例性实施例的基于自适应多速率卡尔曼滤波器和改进的特征匹配算法的实时位移估计步骤S200的详细程序的流程图;
图8示出了根据本发明的一个示例性实施例的感兴趣区域(ROI)更新算法的工作原理;
图9示意性地示出了根据本发明的一个示例性实施例的用于从图像帧中的特征点匹配选择最佳匹配的方法;
图10示意性地示出了根据本发明的一个示例性实施例的用于对照最佳匹配交叉检查剩余匹配以消除失配的方法;
图11示意性地示出了根据本发明的一个示例性实施例的使用自适应多速率卡尔曼滤波器来融合异步加速度测量和视觉测量的概念的概述;
图12示出了根据本发明的一个示例性实施例的使用自适应多速率卡尔曼滤波器来融合异步加速度测量和视觉测量的示例性过程。
具体实施方式
下面将参考附图更详细地描述本发明的优选实施例。附图中相同的附图标记用于相同的部件,并且省略对相同部件的重复描述。
本发明中使用的术语仅用于描述具体实施例,并且不旨在限制本发明。单数表达包括复数,除非上下文另有明确指示。在本申请中,术语包括”或“具有”等旨在指示说明书中列举的特征、数字、步骤、动作、部件、零件或其组合的存在,并且不排除一个或多个其他特征、数字、步骤、动作、部件、零件或其组合的存在或添加的可能性。此外,诸如第一、第二等术语可以用于描述各种部件,但是部件不受此类术语的限制。术语仅用于区分一个部件与另一个。
首先,为了更好地理解本发明的概念,有必要理解使用视觉相机的基于特征点的位移估计方法的基本概念。图1示意性地示出了使用视觉相机的传统的基于特征匹配的位移估计方法。一般来说,视觉相机14通常安装在远离为目标结构的桥梁10的静止位置处,并且测量存在或安装在桥梁10上的自然或人工目标12。目标12可以在视觉相机14的视野(FOV)内。在视觉相机14拍摄的一系列图像帧中,可以利用第1个图像帧作为参考帧,并且可以将特征匹配算法应用于第1个和第i个帧,以估计在第i个帧处的目标相关于在第1个帧处的目标的相对位移。图2中示出了详细程序。
图2示出了使用特征匹配算法的传统位移估计程序。虽然图2中仅示出了竖直位移估计程序,但是所描述的位移估计技术可以很容易地扩展到水平位移估计。
参考图2,可以首先从待比较的两个图像帧201和20i选择感兴趣区域(ROI)22(参见图2的(a))。ROI 22的大小和位置可以包括第1个图像帧201的感兴趣目标(TOI)24,并且可以在后续图像帧内维持相同的ROI 22。
检测两个图像帧201、20i的感兴趣区域22中的特征点并且使其彼此匹配(参见图2的(b))。从第1个图像帧201的第1个ROI 22检测到的特征点M1 1、M1 2、......、M1 N与在第i个图像帧20i的第i个ROI 22内的对应特征点Mi 1、M1 2、......、Mi N进行匹配。在这里,两个匹配的特征点,诸如图2(b)中所描绘的M1 1和Mi 1被命名为特征匹配或简称为匹配。本领域已知有若干可以用于特征查找和匹配的算法(Lowe,2004;Bay等,2006;Rublee等,2011),其可以应用于本发明。特别地,加速鲁棒特征(SURF)(Bay、Tuytelaars和Van Gool,2006)技术由于其高准确性和计算速率可以用于本发明。
然后,可以计算第j个特征匹配在第1个图像帧201和第i个图像帧20i之间的相对移位/>N个特征点Mi 1、M1 2、......、Mi N的相对移位的平均相对平移di可以由以下方程进行计算(参见图2(c)),
在这里,N是匹配的特征点的数量。可以使用比例因子α将像素单位的平移转换为空间长度单位的位移ui(参见图2(d))。换句话说,比例因子α是用于将视频图像的像素分辨率转换为空间分辨率的转换因子。
ui=α×di......(2)
在这里,可以使用以下方程来估计比例因子α。
在这里,utarget表示人工或自然目标的实际物理尺寸,并且dtarget表示相同目标的像素单位的对应尺寸。
然而,这种传统的位移估计技术有若干的限制。图3是用于示出这一点的图。首先,应提前估计最大位移以确定ROI的大小。因为ROI的位置对于所有帧都是固定的,所以当TOI移动与ROI的尺寸相比太大时,TOI 32可能移动到ROI 30之外。这可以在图3(a)中看到。然而,如图3(b)所示,扩大ROI 34的尺寸可能需要更多的计算,并且增加了捕获ROI 34内的不需要的目标36的可能性。其次,传统的特征匹配算法经常产生失配,这损害了位移估计的准确性。当匹配两个不相关的特征(例如,图2(b)中的M1 1和Mi 2)以进行位移估计时,这种现象被称为失配。尽管已经提出了几种失配拒绝算法,但它们在实施中具有限制,诸如需要用户交互来建立特设阈值或是计算密集型的,从而使得实时估计变得困难。第三,必须提前知道目标的大小和长度以估计比例因子,这对于民用基础设施应用来说可能很麻烦。
另一方面,基于视觉的位移估计方法需要大量的计算,因此基于视觉的位移估计方法常常受限于低采样率。为了补偿这一点,已经尝试通过用多速率卡尔曼滤波器融合低采样率基于视觉的位移测量和高采样率加速度测量来以更高的采样率估计位移。多速率卡尔曼滤波器的工作原理被简单地总结如下。
令和xk分别表示在第k个时间步骤的真实速度和位移,用于加速度-位移关系的离散状态空间模型可以表示如下。
xk=A(Δta)xk-1+B(Δta)ak-1+B(Δta)wk-1......(4)
uk=Hxk+vk......(5),
在这里,xk是状态变量并且ak-1和uk分别是测量的加速度和基于视觉的位移。wk-1和vk分别是具有方差Q和R的对应噪声。Δta是加速度测量的时间间隔。H是向量[1,0]。A和B是加速度测量时间间隔Δta的函数。
基于该模型,制定了多速率卡尔曼滤波器,以用于使用以不同采样率但同步采样的加速度和基于视觉的位移进行位移估计。在这里,同步测量指出加速度测量的采样率是基于视觉的位移的采样率的整数倍。可以使用后验状态估计值和“前一时间步骤”处的加速度ak-1来获得先验状态估计值/>如下面的方程中所示。
先验状态估计值的误差协方差矩阵/>也通过以下获得:
在这里,表示/>的误差协方差矩阵。
当在第k个时间步骤基于视觉的位移uk和加速度都可用时,此时的卡尔曼增益(K)计算如下:
并且,后验状态估计值及其误差协方差矩阵/>按如下方式获得。
如果在第k个时间步骤仅加速度数据可用,但基于视觉的位移uk不可用,则位移的后验状态估计值等于先验状态估计值,如下面的两个方程中所示,并且这两个估计值的误差协方差矩阵也彼此相等。
该卡尔曼滤波器仅适用于当加速度和视觉测量数据同步并且其采样率之比为整数时。另外,卡尔曼滤波器的性能显著地依赖于Q和R的估计。虽然可以在实验室环境中估计用于加速度测量的Q值,但用于基于视觉的位移的R值的估计是很难的,这是因为该值深受光照条件、视觉相机与目标之间的距离以及目标的纹理的影响。需要一种用于噪声方差R的自适应估计技术。
为了解决上述限制,根据本发明的一个示例性实施例提供了一种新的位移估计方法。图4示意性地示出了根据本发明的一个示例性实施例的用于执行结构位移估计方法的系统。
参考图4,结构位移估计系统100可以包括视觉相机110、加速度计120和位移估计单元130。视觉相机110可以包括用于记录视频的相机和用于分析记录的视频的分析程序。视觉相机110和加速度计120可以设置在待测量位移的结构70上。即,具有相对低采样率的视觉相机110可以放置在结构70上的期望点处以跟踪结构70周围的自然目标80。加速度计120也可以放置在与视觉相机110相同的位置处,以按相对较高的采样率测量加速度。目标80可以选自存在于结构70外部并且不改变其位置的静止物体。
位移估计单元130可以被配置为在给定来自视觉相机110的目标80的视觉测量数据和来自加速度计120的结构70的加速度数据的情况下,计算结构70的位移。为此,位移估计单元130可以包括根据本文所述的算法实施的用于估计位移的计算机程序,以及能够执行计算机程序以执行计算处理以产生期望的位移估计值的硬件资源。用于位移估计单元130的硬件资源可以包括计算装置,其包括处理器132。除了处理器132之外,计算装置还可以包括存储器134、为非易失性存储装置的数据存储器136、输入/输出装置138等。例如,位移估计单元130的硬件可以包括通用或本发明专用的计算机装置、工作站装置等,其包括上述装置。
图5示意性地示出了根据本发明的一个示例性实施例的结构位移估计方法的一般概述。
参考图5,结构位移估计方法可以广泛地包括比例因子计算步骤S100,其中比例因子由比例因子计算单元150自动计算;以及位移估计步骤S200,其中基于自适应多速率卡尔曼滤波器单元160和改进的特征匹配算法实时地估计结构70的位移。
图6是示出执行图5所示的比例因子计算步骤S100的详细程序的流程图。
参考图6,由视觉相机110获取的视觉测量数据和由加速度计120测量的加速度数据用于计算比例因子。比例因子的计算可以由位移估计单元130的比例因子计算单元150执行。比例因子计算单元150可以被实施为要由处理器132执行的计算机程序。比例因子计算单元150可以使处理器132分别从视觉相机110和加速度计120接收由视觉相机110获取的视觉测量数据和由加速度计120测量的加速度数据(S102、S112)。视觉相机110的成像和加速度计120的测量可能以短时间间隔发生,并且视觉测量数据的采样率可能慢于加速度数据的采样率。
在视觉测量数据(即,随时间拍摄的一系列图像帧)中,可以通过应用特征匹配算法来估计特征点的逐像素平移,该特征匹配算法比较在参考图像帧(例如,第1个图像帧)与比较目标图像帧(例如,第i个图像帧)之间的预定特征点以匹配相同的特征点,如图2所述(S104、S106)。可以使用上述方程(1)等来估计特征点的平移。可以使用所测量的加速度数据随时间的二重积分来估计真实空间中的长度单位的位移(S114、S116)。
可以对先前获得的平移估计值和位移估计值进行带通滤波以去除噪声(S108、S118)。在这里,df(t)表示带通滤波后的平移并且uf(t)表示带通滤波后的位移(S110、S120)。在带通滤波中,带通滤波器的下截止频率flc可以设置得足够高,以消除加速度信号失真,即基于加速度的位移中的低频漂移,并且带通滤波器的上截止频率fuc可以考虑视觉相机110的采样率来进行设置。例如,上截止频率fuc可以设置为视觉相机110的采样率的约1/10,使得视觉相机110的采样率可以比信号的最高频率,即视觉测量数据的频率分量快至少10倍。换句话说,如果需要识别时域中的一个频段,则必须以快至少10倍的采样率获取数据。
滤波后的平移df(t)和滤波后的位移uf(t)与比例因子(α)之间存在以下关系。
uf(t)=α×df(t)......(14)
在一个示例性实施例中,可以使用最小二乘算法来估计比例因子α(S126)。作为另一个示例,可以使用随机抽样一致(RANSAC)算法来估计比例因子α。使用最小二乘算法或RANSAC算法,可以将平移估计值df(t)与位移估计值uf(t)进行比较,以获得两个值之间的斜率,其是期望的比例因子α。
为了匹配位移的采样率uf(t)与平移的采样率df(t),可以在应用最小二乘算法(S126)之前将带通滤波后的位移uf(t)下采样至视觉相机110的采样率。通过该处理,位移估计值可以与平移估计值同步。
所获得的比例因子α可以用于位移估计的步骤S200中(S128)。
图7a和图7b示出了根据本发明的一个示例性实施例的基于自适应多速率卡尔曼滤波器和改进的特征匹配算法的实时位移估计步骤S200的详细执行。
位移估计步骤S200可以在自适应多速率卡尔曼滤波器单元160中执行。自适应多速率卡尔曼滤波器单元160可以被实施为计算机程序并且由处理器132执行。
参考图7a和图7b,自适应多速率卡尔曼滤波器单元160可以使处理器132以逐帧增量读取由视觉相机110提供的并且存储在数据存储器136中的视觉测量数据,同时使帧索引i以1递增(S210、S212)。可以随着时间t以第一预定时间间隔Δtd拍摄该系列图像帧。
与传统的特征点匹配相比,本发明提出的方法的主要改进可以是增加了ROI更新和两种失配拒绝算法。由于例如结构70的移动,连续图像帧中的ROI的位置可以随时间改变。本发明使用基于相关性的模板匹配。除了被视为模板的ROI的初始选择之外,模板匹配几乎不依赖于用户干预。模板匹配方法对背景和照明变化敏感。因此,对于长期记录(例如,超过几个小时的)而言,可能需要定期更新ROI模板以减轻误差累积。改进的特征匹配算法可以使基于视觉的位移测量更加可靠。为了提高可靠性,可以更新当前第i个图像帧内的ROI(以下简称“第i个ROI')(S214和S216)。使用可以在图7b的步骤S234中获得的长度单位的位移的先验状态估计值可以更新第i个ROI。
图8示出了根据本发明的一个示例性实施例的ROI更新算法的工作原理。
参考图8,可以使用所估计的纵长位移的先验状态估计值(1)来获得感兴趣目标(TOI)的像素单位的平移di TOI。由于所估计的纵长位移的先验状态估计值/>(1)采用位移单位,其可以使用下面带比例因子的方程转换为像素单位。
长度单位的位移的先验状态估计值表示/>的第一项。第i个图像帧中的ROI(第i个ROI)使用采用TOI的像素单位的平移di TOI进行移位。第i个ROI的平移di TOI以逐像素为基础进行离散化。因此,可以使用如下的“轮”函数获得第i个ROI的di TOI。
第i个ROI的平移di TOI可以用于更新第i个图像帧的第i个ROI的位置。
在更新第i个图像帧的第i个ROI之后,可以使用预定的特征点匹配算法在第1个图像帧的ROI(第1个ROI)中的N个特征点与在更新的第i个图像帧的ROI(第i个ROI)中的N个特征点之间执行特征点匹配(S218)。可以使用传统的特征匹配算法。
然而,在N个特征点匹配中,可能存在一些正确匹配和一些失配,即不正确匹配的特征。对于图像帧中准确的平移估计而言,重要的是要确保相关联的特征点以高置信度彼此匹配,如图2b所示。然而,传统的特征匹配技术经常导致特征点失配。这降低了平移和比例因子的准确性。对于高置信度的位移测量而言,可以使用失配拒绝算法来从N个特征点匹配拒绝失配,从而仅留下良好的特征点匹配(S220、S222)。可以通过使用在步骤S234中产生的纵长位移的先验状态估计值及其误差协方差矩阵/>的两种拒绝算法来拒绝失配。最终,仅使用剩余的高置信度的特征点匹配来进行位移估计,从而提高位移测量的可靠性。
在一个示例性实施例中,可以对匹配的特征匹配相继执行第一和第二特征点失配拒绝算法。
首先,第一特征点失配拒绝算法可以包括两个步骤,其大致如下:(1)在特征匹配中选择最佳特征匹配(参见图8),以及(2)通过对照最佳特征匹配交叉检查剩余特征匹配来拒绝失配(参见图9)。在执行第一特征点失配拒绝算法之后,可以对剩余的特征匹配执行第二特征点失配拒绝算法。拒绝特征点失配的基本假设是在图像平面中的两个特征点之间的距离不会随着时间变化。
首先,下面将描述使用第一特征点失配拒绝算法的对特征点失配的第一阶段的拒绝。
图9示意性地示出了根据本发明的一个示例性实施例的用于从图像帧中的特征点匹配选择最佳匹配的方法。
参考图6和图9,选择最佳匹配可以是第一步骤。首先,可以将特征点匹配分成两组,如图9(A)所示。假设第1个图像帧的ROI(第1个ROI)和第i个图像帧的ROI(第i个ROI)各自包含N个特征点,通过应用特征匹配算法从第1个ROI和第i个ROI获得N个特征点匹配。在这里,每个特征点匹配可以由从第1个ROI选择的一个特征点和从第i个ROI选择的一个特征点组成。在N个特征点匹配中,可以随机选择两个特征点匹配来形成一组匹配。可以对所有N个特征点匹配重复该分组以创建N/2个匹配组。如果N为奇数,则分组时可以忽略剩余的特征点匹配,这会导致(N-1)/2个匹配组。每个匹配组包含第1个ROI中的两个特征点和第i个ROI中对应的两个特征点(参见图9(A))。
在如上构造N/2或(N-1)/2个匹配组后,就可以计算出每个匹配组的距离变化指数,如图9(B)所示。可以使用以下方程来计算距离变化指数。
在这里,μs表示第s个匹配组的距离变化指数(其中,s为小于N/2或(N-1)/2的自然数)。另外,l1 s为第1个图像帧的ROI内第s个匹配组的两个特征点之间的距离,并且为在第i个图像帧的ROI内第s个匹配组的两个特征点之间的距离。
对于具有两个良好匹配的组而言,距离变化指数的值应当接近于零。因此,如图9(c)所示,可以选择具有计算出的N/2或(N-1)/2个距离变化指数中的最小距离变化指数(μs)值的匹配组作为最佳匹配组。属于最佳匹配组的两个特征点匹配中的一个可以被定义为最佳匹配。在图9(c)中,例如,第r个特征匹配被例示为最佳匹配。
图10示意性地示出了根据本发明的一个示例性实施例的用于对照最佳匹配交叉检查剩余匹配以拒绝失配的方法。
图10示出了基于最佳匹配来计算剩余匹配的距离间隔。
参考图10,一旦选择了最佳匹配,例如,就可以基于所选择的最佳匹配来计算剩余特征匹配中的每一个的距离间隔,以确定并拒绝失配。具体地,在第1个图像帧中,在第1个ROI内的任意特征点,例如,/>与最佳匹配特征点M1 r之间的真实距离落在最大距离/>与最小距离/>之间,如图9所示。在这里,j=1,2,...,N且j≠r。
在第1个图像帧中,在第1个ROI内的任意特征点与最佳匹配特征点M1 r之间的最大距离/>和最小距离/>可以表示如下。
在这里,ε表示x轴和y轴方向上的最大像素离散化误差,并且在图9中,该值被示例性地设置为0.5像素。另外,和/>分别表示特征点/>和特征点M1 r的x和y坐标值之间的差。
类似地,在第i个图像帧的第i个ROI内的任意特征点与最佳匹配特征点/>之间的最大距离/>和最小距离/>还可以使用方程(18)计算。
最终,以这种方式,就可以计算出关于在第1个和第i个图像帧的每个ROI内的最佳匹配特征点至剩余特征点的最大和最小距离。
可以确定在第1个图像帧的ROI内的计算出的最大距离和最小距离/>以及在第i个图像帧的ROI内的计算出的最大距离/>和最小距离/>是否满足以下两个距离间隔条件,其中j=1,2,...,N且j≠r。
(i)距离间隔条件1:“在第1个图像帧的ROI内第j个特征点与最佳匹配特征点之间的最大距离”不小于“在第i个图像帧的ROI内的第j个特征点与最佳匹配特征点之间的最小距离”;以及
(ii)距离间隔条件2:“在第i个图像帧的ROI内第j个特征点与最佳匹配特征点之间的最大距离”不小于“在第1个图像帧的ROI内的第j个特征点与最佳匹配特征点之间的最小距离”。
下面的方程(19)示出了两个条件。作为判断的结果,不满足以下方程的特征点匹配被分类为失配并且可以被拒绝。
在通过应用第一特征点失配拒绝算法来拒绝特征点失配之后,可以应用第二特征点失配拒绝算法来拒绝特征点失配。改进的特征点匹配算法还可以提高基于视觉的位移测量的可靠性。
具体地,在由第一特征点失配拒绝算法进行滤波之后,可以为剩余特征点匹配中的每一个计算平移假设长度单位的位移的先验状态估计值/>中的误差具有带协方差的正态分布/>则可以估计平移/>的99.7%置信区间。平移/>的范围可以由以下方程表示。/>
在这里,协方差误差表示位移的先验状态估计值的协方差误差矩阵/>的第一项。可以在该置信区间的估计内考虑像素离散误差ε。不在该置信区间内的特征点匹配被视为失配并且被拒绝。
以这种方式,可以执行第一和第二特征点失配拒绝算法来拒绝被分类为失配的特征匹配,并且可以获得剩余的特征点匹配作为良好的匹配(S222)。
良好的匹配可以用于估计平移(S224)。最终平移被估计为所有剩余的良好匹配和ROI移动的平均平移的总和。
通过将在步骤S100中获得的缩放因子α应用于最终平移的总和,可以估计结构70的期望位移ui(S224)。在一个示例性实施例中,可以使用下面的方程来计算位移ui(S226)。
在这里,Ng是良好匹配,即,高置信度匹配的数量。
另一方面,如上面所提及的,假设在传统卡尔曼滤波器中用于数据融合的加速度数据和视觉测量数据是时间同步的。然而,如上面所提及的,视觉相机110和加速度计120的测量速率可以是异步的。根据示例性实施例的自适应多速率卡尔曼滤波器单元160可以无缝地融合异步的时间加速度数据和基于视觉的数据。另外,视觉相机110的每秒帧数(FPS)在真实拍摄中可能不会随时间精确地保持恒定(例如,30FPS),而是可能变化(例如,29.5FPS、30.3FPS等)。自适应多速率卡尔曼滤波器单元160允许相关于基于视觉的位移自适应地估计该噪声方差R。
图11示意性地示出了根据本发明的一个示例性实施例的使用自适应多速率卡尔曼滤波器单元来融合异步加速度测量数据和视觉测量数据的概念的概述。图12示出了根据本发明的一个示例性实施例的使用自适应多速率卡尔曼滤波器单元来融合异步加速度测量数据和视觉测量数据的示例性过程。
结合图7B参考图11和图12,由加速度计120随时间测量的加速度数据ak可以具有时间间隔Δta。在第k个时间步骤,加速度数据ak可以被输入到位移估计单元130的卡尔曼滤波器单元160(S230),其中k是加速度数据的索引。
在卡尔曼滤波器单元160中,可以根据加速度和视觉测量数据的可用性将所有时间步骤分为三种类型。可以保留加速步骤处的状态估计值,使得最终估计的位移具有与加速度测量相同的采样率,即,所得位移的采样频率可以等于相对较高的加速度数据的采样频率。
视觉相机110所捕获的图像帧数据(第i个帧)可以具有时间间隔Δtd。在时间尺度上,图像帧数据可能位于或可能不位于两个相邻加速度数据之间。根据一个示例性实施例的卡尔曼滤波器可以针对三种不同类型的时间步骤进行制定。
(i)类型I时间步骤(180-1):在该时间步骤中,在当前加速度数据ak和前一个加速度数据ak-1之间的时间间隔[t-Δta,t]内,仅当前的加速度数据ak可用并且没有可用的图像帧数据。
(ii)类型II时间步骤(180-2):在该时间步骤中,仅图像帧数据(第i个帧)可用,而加速度数据不可用。
(iii)类型III时间步骤(180-3):在该时间步骤中,在当前加速度数据(ak)和下一个加速度数据(ak+1)之间的时间间隔[t-Δta,t]内,加速度数据和图像帧数据均可用。
位移估计单元130的卡尔曼滤波器单元160可以确定在时间间隔[(k-2)Δta,(k-1)Δta]内是否存在图像帧数据(S232)。
如果在时间间隔[(k-2)Δta,(k-1)Δta]之间没有可用的图像帧数据,则其对应于类型I时间步骤180-1。在类型I时间步骤180-1中,由于仅加速度数据可用,可以分别使用前述方程(7)和(8)以及方程(12)和(13)来计算状态变量xk及其协方差误差矩阵这类似于传统的卡尔曼滤波器。具体地,在时间(k-1)Δta,可以使用方程(7)来计算位移的先验状态估计值/>并且可以使用方程(8)获得其先验误差协方差矩阵/>接下来,在时间(k-1)Δta,可以使用方程(10)来获得位移的后验状态估计值/>并且可以使用方程(11)获得其后验误差协方差矩阵/>如上面的方程(12)和(13)中所提及的,以这种方式获得的位移的先验和后验状态估计值/>和/>彼此相同,并且它们的误差协方差矩阵/>和/>也彼此相同。
如果在步骤S232中确定在时间间隔[(k-2)Δta,(k-1)Δta]之间存在图像帧数据,则其可以对应于类型II时间步骤180-2或类型III时间步骤180-3。
在类型II时间步骤180-2中,可以使用方程(7)和方程(8)来获得状态变量及其误差协方差矩阵。然而,在该时间,时间间隔只是从Δta变为((i-1)Δtd-(k-1)Δta),其中k是加速度数据的第k个时间步骤,并且i是视觉测量数据的第i个时间步骤。Δta是加速度数据之间的时间间隔,并且Δtd是视觉测量数据之间的时间间隔。在类型II时间步骤180-2中,可以获得在时间t=(i-1)Δtd处的位移的先验状态估计值及其误差协方差矩阵/>(S234)。
所获得的位移的先验状态估计值及其误差协方差矩阵/>可以用于更新第i个图像帧的ROI(S214、S216),如上面所提及的。此外,改进的第一和第二特征点失配拒绝算法可以应用于视觉测量数据以拒绝特征点失配(S218至S222)以获得基于视觉的位移ui的估计值(S224至S226)。
此外,可以使用获得的位移的先验状态估计值及其误差协方差矩阵/>以及估计的基于视觉的位移ui来估计基于视觉的位移ui的噪声方差Ri(S236至S238)。可以使用下面的方程自适应地估计基于视觉的位移ui的噪声方差Ri(S238)。该方差方程是通过已知的协方差匹配技术提出的(Mohamed和Schwarz,1999)。
在这里,Ri表示基于视觉的位移ui的噪声方差,并且E(·)表示期望操作。并且,创新度ηi,即测量值和预测值之间的差异可以定义如下。
尽管可以通过在移动时间窗口上对ηi 2进行平均来接近ηiηi T的期望值,但由于移动时间窗口的大小,这可能是计算密集型的。为了解决这个问题,可以使用遗忘因子β来自适应地估计基于视觉的位移的噪声方差Ri,如下面的方程中所示(Akhlaghi等,2017)。
在这里,遗忘因子β具有在范围0<β<1中的值。
所获得的基于视觉的位移的噪声方差Ri可以用于计算卡尔曼增益K(S240)。可以使用方程(9)来计算卡尔曼增益K。
在类型II时间步骤(180-2)中,即在时间t=(i-1)Δtd处,一旦获得卡尔曼增益K,其就可以用来获得位移的后验状态估计值及其误差协方差矩阵/>可以使用上述方程(10)和(11)来获得位移的后验状态估计值/>及其误差协方差矩阵/>(S242)。
接下来,可以获得类型III时间步骤(180-3)处的位移的先验状态估计值和后验状态估计值以及它们的误差协方差矩阵(S244、S246)。
在类型III时间步骤180-3中,其中在时间kΔta-(i-1)Δta期间在当前加速度数据ak与后续加速度数据ak+1之间的时间间隔[t-Δta,t]内加速度数据和图像帧数据均可用,可以使用类型I时间步骤180-1的加速度数据ak和在类型II时间步骤180-2中获得的位移的后验状态估计值及其误差协方差矩阵/>来获得位移的先验状态估计值。而且,在该时间步骤180-3中,位移的先验状态和后验状态估计值相等。可以如下估计这些。
根据在ui和ak之间的关系,方程(25)可以改写为:
在这里,W1、W2和W3分别表示如下。
W1=A(kΔta-(i-1)Δtd)(I-KH)A((i-1)Δtd-(k-1)Δta)......(27)
W2=A(kΔta-(i-1)Δtd)(I-KH)B((i-1)Δtd-(k-1)Δta)+B(kΔta-(i-1)Δtd)......(28)
W3=A(kΔta-(i–1)Δtd)K......(29)
由于ak和ui的值彼此独立,因此/>误差协方差矩阵计算如下。
通过如上所述用自适应多速率卡尔曼滤波器融合异步加速度数据和视觉数据,可以估计在时间t=(k-1)Δta处的位移的状态变量
如上所述,由于作为异构数据的图像帧和加速度以时间不同步的方式获取,因此可以通过应用卡尔曼滤波器来校正根据示例性实施例的结构位移估计方法。在这种情况下,自动计算和减少随分辨率和着色变化的基于视觉的位移噪声。
根据上述实施例的结构位移估计方法可以是以能够通过各种计算装置执行的程序指令的形式实施的软件。软件可以包括计算机程序、代码、指令或其一种或多种组合,并且可以配置处理装置以根据需要操作,或者可以独立地或共同地指示处理装置。程序指令可以记录在计算机可读介质上。计算机可读介质可以单独地或组合地包括程序指令、数据文件、数据结构等。记录在介质上的程序指令可以是为实施例专门设计和配置的,并且可以是计算机程序领域中的技术人员已知和可获得的。计算机可读记录介质的示例包括磁介质,诸如硬盘、软盘和磁带,光学介质,诸如CD-ROM和DVD,磁光介质,诸如光磁软盘,以及硬件装置,诸如ROM、RAM、闪存等,其被专门配置为存储和执行程序指令。
用于实施根据所描述的实施例的结构位移估计方法的计算机装置可以使用一个或多个通用或专用计算机来实现,诸如,例如,处理器、控制器、算术逻辑单元(ALU)、数字信号处理器、微型计算机、现场可编程阵列(FPA)、可编程逻辑单元(PLU)、微处理器或任何其他能够执行和响应指令的装置。
[工业实用性]
可以应用本发明以对安全性要求较高的结构,诸如桥梁、建筑物等的位移进行实时监测。
虽然已经借助于有限的附图示出了上述实施例,但是应当理解本领域的技术人员可以在不背离以下权利要求书中描述的本发明的构想和范围的情况下对本发明进行各种修改和改变。因此,其他实施方式、其他实施例以及与本专利的权利要求等同的事物也在以下专利权利要求的范围内。
Claims (24)
1.一种基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量所述结构的位移的方法,其包括:
使用随时间由安装在要测量的所述结构上的相机以在所述结构外部的静止目标的预定第一采样频率获取的图像帧数据,以及由安装在所述结构上的加速度计以预定第二采样频率测量的加速度数据,自动计算比例因子(α),所述比例因子(α)是所述图像帧数据的像素分辨率与真实空间的距离分辨率之间的转换因子;以及
使用改进的特征点匹配算法和所述比例因子估计所述结构的位移值,并且通过使用自适应多速率卡尔曼滤波器来无缝融合异步的所述图像帧数据和所述加速度数据以使它们以减小的噪声彼此同步,同时使用所述估计的位移值来以所述加速度数据的采样频率计算所述结构的位移。
2.根据权利要求1所述的方法,其中所述“自动计算比例因子(α)”包括通过使用用于比较参考图像帧与比较图像帧之间的预定特征点的特征点匹配算法来匹配所述相同的特征点来估计所述特征点的逐像素平移;以及通过对所述加速度数据进行二重积分来估计真实空间中的纵长位移。
3.根据权利要求2所述的方法,其中所述“自动计算比例因子(α)”还包括对所述平移的所述获得的估计值和所述位移的所述估计值进行带通滤波处理以减少噪声。
4.根据权利要求3所述的方法,其中在所述带通滤波处理中,所述相机的所述采样频率比所述图像帧的频率分量的最高频率高至少10倍。
5.根据权利要求3所述的方法,其中所述“自动计算比例因子(α)”还包括对所述位移的带通滤波估计值进行下采样以使所述位移的估计值与所述平移的所述估计值同步。
6.根据权利要求3所述的方法,其中所述比例因子(α)是在通过比较所述平移的所述带通滤波估计值与所述位移的所述带通滤波估计值(uf(t))获得的所述两个值之间的斜率。
7.根据权利要求5所述的方法,其中所述比例因子(α)是通过将最小二乘法或RANSAC(随机抽样一致)算法应用于所述带通滤波平移的估计值(df(t))和下采样位移的估计值(uf(t))获得的。
8.根据权利要求1所述的方法,其中所述“估计所述结构的位移值”包括更新当前第i个图像帧的感兴趣区域(ROI);使用预定的特征点匹配算法将第1个图像帧的ROI(第1个ROI)内的N个特征点与所述第i个图像帧的感兴趣区域(第i个ROI)内的N个特征点进行匹配;使用失配拒绝算法拒绝特征点失配,以仅留下在所述N个特征点匹配中良好的特征点匹配,以用于可靠的位移测量;并且通过使用方程
来估计最终的平移,其中Ng表示所述良好的特征点匹配的数量,其中所述比例因子(α)被应用于所述良好特征点匹配的所述平移和所述ROI的移动的总和。
9.根据权利要求8所述的方法,其中所述“更新第i个图像帧的ROI”包括使用所述比例因子将所述位移的先验状态估计值从纵长的一个转换为逐像素的一个;通过使用转换的逐像素的先验状态估计值获得所述第i个图像帧的感兴趣目标(TOI)的逐像素位移(di TOI);使用轮函数获得所述第i个ROI的平移(di TOI);并且使用所述第i个ROI的所述获得的平移(di TOI)更新所述第i个ROI的位置。
10.根据权利要求8所述的方法,其中所述“拒绝特征点失配以仅留下良好的特征点匹配”包括使用第一特征点失配拒绝算法来拒绝特征点失配,其中所述第一特征点失配拒绝算法包括:从所述特征点匹配选择最佳的特征点匹配;以及基于所述最佳的特征点匹配通过交叉检查剩余的特征点来拒绝失配的特征点。
11.根据权利要求10所述的方法,其中所述“选择最佳的特征点匹配”包括:将所述第1个图像帧的所述ROI与所述第i个图像帧的所述ROI之间的N个特征点匹配两两配对为N/2或(N-1)/2个特征点匹配组;计算所述N/2或(N-1)/2个特征点匹配组中每一个的距离变化指数;在所述计算的N/2或(N-1)/2个距离变化指数中选择具有最小距离变化指数(μs)的匹配组作为最佳匹配组;并且选择属于所述最佳匹配组的所述两个特征点匹配中的一个作为所述最佳特征点匹配。
12.根据权利要求11所述的方法,其中所述“拒绝失配的特征点”包括:基于所述选择的最佳匹配来计算所述剩余特征点匹配中每一个的距离间隔,以确定所述计算的距离间隔是否满足失配条件;并且从所述特征点失配拒绝满足所述失配条件的特征点,从而仅留下良好的匹配。
13.根据权利要求12所述的方法,其中所述失配条件包括第一和第二距离间隔条件,所述第一间隔条件为“在第j个特征点与所述第1个图像帧的所述ROI内的所述最佳匹配的特征点之间的最大距离”大于或等于“在所述第j个特征点与所述第i个图像帧的所述ROI内的所述最佳匹配特征点之间的最小距离”,并且所述第二距离间隔条件为“在第j个特征点与所述第i个图像帧的所述ROI内的所述最佳匹配的特征点之间的最大距离”大于或等于“在所述第j个特征点与所述第1个图像帧的所述ROI内的所述最佳匹配特征点之间的最小距离”,其中j=1、2、......、N。
14.根据权利要求11所述的方法,其中所述距离变化指数是通过方程计算的,其中μs是第s个匹配组的距离变化指数,l1 s是在所述第1个图像帧的所述ROI内的两个特征点之间的距离,并且li s是在所述第i个图像帧的所述ROI内的两个特征点之间的距离。
15.根据权利要求8所述的方法,其中所述“拒绝特征点失配以仅留下良好的特征点匹配”包括使用第二特征点失配拒绝算法来拒绝所述特征点失配,并且其中所述第二特征点失配拒绝算法包括计算由所述第一特征点失配拒绝算法进行滤波的剩余特征点匹配中的每一个的平移(di j);并且将所计算的平移未落在由
方程表示的平移值(di j)的范围内的特征点匹配作为失配进行拒绝。
16.根据权利要求1所述的方法,其中使用所述自适应多速率卡尔曼滤波器进行的所述图像帧数据与所述加速度数据的所述融合是通过以下实现的:在类型I时间步骤中执行,其中在所述当前的加速度数据(ak)与先前的加速度数据(ak-1)之间的时间间隔[(k-2)Δta,(k-1)Δta]中,仅当前的加速度数据(ak)可用并且没有可用的图像帧数据,分别使用第一方程和第二方程获得所述位移的先验状态估计值/>以及所述先验状态估计值/>的误差协方差矩阵/>并且在时间(k-1)Δta,获得对应于所述先验状态估计值/>和所述先验状态估计值/>的所述误差协方差矩阵/>的后验状态估计值/>和所述后验状态估计值/>的后验误差协方差矩阵/>其中Δta表示在所述加速度数据之间的时间间隔,k表示所述加速度数据的第k个时间间隔,并且/>和ak-1分别表示在
“前一时间步骤”和的所述后验状态估计值和加速度数据。
17.根据权利要求16所述的方法,其中使用所述自适应多速率卡尔曼滤波器进行的所述图像帧数据与所述加速度数据的所述融合是通过以下实现的:在类型II时间步骤中执行,其中在时间间隔[(k-2)Δta,(k-1)Δta]内仅所述图像帧数据(第i个帧)可用,并且所述加速度数据不可用,使用所述第一和第二方程估计所述位移的先验状态估计值和所述位移的所述先验状态估计值/>的先验误差协方差矩阵/>其中将所述时间间隔从Δta改变为(i-1)Δtd-(k-1)Δta;使用所述获得的先验状态估计值/>及其所述先验误差协方差矩阵/>更新所述第i个图像帧的所述ROI;通过基于所述第一和第二特征点失配拒绝算法将所述获得的先验状态估计值/>及其所述先验误差协方差矩阵/>应用于所述特征点失配拒绝来估计基于视觉的位移(ui);通过使用所述位移的所述获得的先验状态估计值/>其所述先验误差协方差矩阵/>以及所述基于视觉的位移(ui)估计值来估计所述基于视觉的位移(ui)的噪声方差(Ri);使用/>方程和所述基于视觉的位移(ui)的所述获得的噪声方差(Ri)来计算卡尔曼增益(K);以及使用方程和/>以及所述卡尔曼增益(K)来计算在时间t=(i-1)Δtd的所述位移的后验状态估计值/>及其误差协方差矩阵/>其中在所述上述方程中,x和P分别由y和G表示,H为向量[1,0],i是所述图像帧数据的第i个时间步骤,并且Δtd是所述图像帧数据之间的时间间隔。
18.根据权利要求17所述的方法,其中使用所述自适应多速率卡尔曼滤波器进行的所述图像帧数据与所述加速度数据的所述融合是通过以下实现的:在类型III时间步骤中执行,其中在所述当前加速度数据(ak)与后续加速度数据(ak+1)之间的时间间隔[t-Δta,t]内,所述加速度数据和所述图像帧数据两者是可用的,使用
方程获得所述位移的先验状态估计值在所述类型I时间步骤中的所述加速度数据(ak)、在所述类型II时间步骤中获得的所述位移的所述后验状态估计值/>以及所述后验状态估计值/>的所述误差协方差矩阵/>使用方程获得所述位移的所述先验状态估计值的误差协方差矩阵/>其中W1=A(kΔta-(i-1)Δtd)(I-KH)A((i-1)Δtd-(k-1)Δta)、W2=A(kΔta-(i-1)Δtd)(I-KH)B((i-1)Δtd-(k-1)Δta)+B(kΔta-(i-1)Δtd)和W3=A(kΔta-(i-1)Δtd)K。
19.一种计算机可读记录介质,其上记录有计算机程序,其用于执行根据权利要求1至18中任一项中所述的基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量所述结构的位移的方法。
20.一种存储在计算机可读记录介质上的计算机可执行程序,其用于执行根据权利要求1至18中任一项中所述的基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量所述结构的位移的方法。
21.一种基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量所述结构的位移的系统,其包括:
安装在要测量的所述结构上的相机,所述相机被配置为以预定的第一采样频率对所述结构外部的静止目标进行成像;
安装在所述结构上的所述相机附近的加速度计,所述加速度计被配置为以预定的第二采样频率测量加速度;以及
位移估计单元,所述位移估计单元包括比例因子计算单元和自适应多速率卡尔曼滤波器单元,其中所述比例因子计算单元被配置为随时间分别使用由所述相机和所述加速度计测量的一系列图像帧数据和加速度数据来执行自动计算比例因子(α)的功能,所述比例因子(α)是在所述图像帧数据的像素分辨率与真实空间的距离分辨率之间的转换因子,并且其中所述自适应多速率卡尔曼滤波器单元被配置为使用改进的特征点匹配算法和所述比例因子来执行估计所述结构的位移值的功能;以及通过使用自适应多速率卡尔曼滤波器来无缝融合异步的所述加速度数据和所述图像帧数据,同时使用所述估计的位移值来以所述加速度数据的采样频率估计所述结构的位移。
22.根据权利要求21所述的系统,其中所述位移估计单元包括实施所述比例因子计算单元和所述自适应多速率卡尔曼滤波器单元的计算机程序,以及处理器单元,所述处理器单元用于执行所述计算机程序以执行自动计算所述比例因子(α)的任务和计算所述结构的所述位移估计值的任务。
23.根据权利要求21所述的系统,其中所述“自动计算比例因子(α)的功能”包括通过使用用于比较参考图像帧与比较图像帧之间的预定特征点的特征点匹配算法来匹配所述相同的特征点来估计所述特征点的逐像素平移;以及通过对所述加速度数据进行二重积分来估计真实空间中的纵长位移。
24.根据权利要求21所述的系统,其中所述“估计所述结构的位移值”的功能包括:更新当前第i个图像帧的感兴趣区域(ROI);使用预定的特征点匹配算法将第1个图像帧的ROI内的N个特征点与所述第i个图像帧的所述ROI内的N个特征点进行匹配;使用失配拒绝算法拒绝特征点失配,以仅留下在所述N个特征点匹配中良好的特征点匹配,以用于可靠的位移测量;并且通过使用
方程来估计最终的平移估计值,其中Ng表示所述良好的特征点匹配的数量,其中所述比例因子(α)被应用于所述良好特征点匹配的所述平移和所述ROI的移动的总和。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
KR1020210047367A KR102565719B1 (ko) | 2021-04-12 | 2021-04-12 | 비동기 자연 표적 영상 계측데이터와 가속도 데이터의 융합에 기초한 구조물 변위 측정 방법 및 이를 위한 시스템 |
KR10-2021-0047367 | 2021-04-12 | ||
PCT/KR2022/003490 WO2022220414A1 (ko) | 2021-04-12 | 2022-03-11 | 비동기 자연 표적 영상 계측데이터와 가속도 데이터의 융합에 기초한 구조물 변위 측정 방법 및 이를 위한 시스템 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN117480356A true CN117480356A (zh) | 2024-01-30 |
Family
ID=83640395
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202280041438.7A Pending CN117480356A (zh) | 2021-04-12 | 2022-03-11 | 基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量结构位移的方法及其系统 |
Country Status (3)
Country | Link |
---|---|
KR (1) | KR102565719B1 (zh) |
CN (1) | CN117480356A (zh) |
WO (1) | WO2022220414A1 (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115790401B (zh) * | 2023-02-09 | 2023-06-16 | 西北工业大学 | 一种基于视觉测量的位移测量方法及相关设备 |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR101803503B1 (ko) * | 2017-02-06 | 2017-11-30 | 주식회사 풍산에프앤에스 | 구조물의 정밀 계측 시스템 및 그 방법 |
KR101996779B1 (ko) * | 2019-02-15 | 2019-07-04 | 화창건축사사무소 주식회사 | 초고층 건축요소의 치수 변형 예측 및 안정화 제어 시스템 |
KR102230397B1 (ko) * | 2020-06-15 | 2021-03-22 | 반석안전주식회사 | 변형률 및 가속도값에 기반한 구조물의 변위 추정방법 |
KR102200824B1 (ko) * | 2020-09-01 | 2021-01-12 | (주)영신디엔씨 | 카메라 및 가속도센서를 이용한 실시간 파일항타 관입량 및 리바운드량 자동측정시스템 |
-
2021
- 2021-04-12 KR KR1020210047367A patent/KR102565719B1/ko active IP Right Grant
-
2022
- 2022-03-11 CN CN202280041438.7A patent/CN117480356A/zh active Pending
- 2022-03-11 WO PCT/KR2022/003490 patent/WO2022220414A1/ko active Application Filing
Also Published As
Publication number | Publication date |
---|---|
KR102565719B1 (ko) | 2023-08-11 |
KR20220141160A (ko) | 2022-10-19 |
WO2022220414A1 (ko) | 2022-10-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Ma et al. | Real‐time structural displacement estimation by fusing asynchronous acceleration and computer vision measurements | |
JP6760114B2 (ja) | 情報処理装置、データ管理装置、データ管理システム、方法、及びプログラム | |
US10762643B2 (en) | Method for evaluating image data of a vehicle camera | |
JP6222353B2 (ja) | 物標検出装置及び物標検出方法 | |
EP2824425B1 (en) | Moving-object position/attitude estimation apparatus and method for estimating position/attitude of moving object | |
JP4874607B2 (ja) | 物体測位装置 | |
US9904841B2 (en) | Method and system for estimating finger movement | |
CN107941212B (zh) | 一种视觉与惯性联合定位方法 | |
CN109716062B (zh) | 姿势估计装置 | |
JP2009008662A (ja) | センサ及びビデオ三角測量を共同使用する物体の検出 | |
JP4887520B2 (ja) | 物体検出装置、物体検出方法及び物体検出プログラム | |
JP6751280B2 (ja) | 位置推定装置、位置検出方法及びプログラム | |
KR101737950B1 (ko) | 지형참조항법에서 영상 기반 항법해 추정 시스템 및 방법 | |
KR101885961B1 (ko) | 이미지를 기반으로 한 객체 위치 추정 방법 및 장치 | |
CN117480356A (zh) | 基于自然目标的异步视觉测量数据和结构的加速度数据的融合来测量结构位移的方法及其系统 | |
JP4361913B2 (ja) | 動き量計算装置 | |
CN116359905A (zh) | 基于4d毫米波雷达的位姿图slam计算方法及系统 | |
CN111723597B (zh) | 跟踪算法的精度检测方法、装置、计算机设备和存储介质 | |
JP2017130067A (ja) | 衛星映像の位置正確度改善のための自動映像処理システム及びその方法 | |
JP2019002790A (ja) | 形状推定装置、形状推定方法及びプログラム | |
CN112630729B (zh) | 一种基于热电堆传感器定位和跟踪室内人类目标的方法 | |
WO2021059765A1 (ja) | 撮像装置、画像処理システム、画像処理方法及びプログラム | |
JP5267100B2 (ja) | 運動推定装置及びプログラム | |
JP4935769B2 (ja) | 平面領域推定装置及びプログラム | |
WO2018139297A1 (ja) | カメラ装置 |
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 |