CN105910584B - 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 - Google Patents
大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 Download PDFInfo
- Publication number
- CN105910584B CN105910584B CN201610297902.5A CN201610297902A CN105910584B CN 105910584 B CN105910584 B CN 105910584B CN 201610297902 A CN201610297902 A CN 201610297902A CN 105910584 B CN105910584 B CN 105910584B
- Authority
- CN
- China
- Prior art keywords
- mtd
- mtr
- msub
- mtable
- orientation
- 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.)
- Active
Links
- 238000011156 evaluation Methods 0.000 title claims abstract description 30
- 238000005259 measurement Methods 0.000 claims abstract description 66
- 238000000034 method Methods 0.000 claims description 69
- 239000011159 matrix material Substances 0.000 claims description 43
- 239000000243 solution Substances 0.000 claims description 38
- 238000003384 imaging method Methods 0.000 claims description 26
- 230000036544 posture Effects 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 9
- 230000008859 change Effects 0.000 claims description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 230000003287 optical effect Effects 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 239000003637 basic solution Substances 0.000 claims description 2
- 238000000354 decomposition reaction Methods 0.000 claims description 2
- 230000000875 corresponding effect Effects 0.000 description 16
- 230000008569 process Effects 0.000 description 12
- 238000002474 experimental method Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 2
- 230000001276 controlling effect Effects 0.000 description 2
- 238000013519 translation Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 229920000049 Carbon (fiber) Polymers 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 239000004917 carbon fiber Substances 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000003550 marker Substances 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000004643 material aging Methods 0.000 description 1
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C15/00—Surveying instruments or accessories not provided for in groups G01C1/00 - G01C13/00
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Multimedia (AREA)
- Length Measuring Devices By Optical Means (AREA)
- Image Processing (AREA)
Abstract
本发明提供了一种大尺寸动态摄影测量系统的高精度定向及定向精度评价方法,包括步骤:a)选定基准尺,在所述基准尺两端布置编码点对其进行长度测量;b)在测量空间,以第一方式移动所述基准尺并调整所述基准尺姿态,以第二方式对其拍摄图片;c)对所述测量系统进行初步相对定位定向;d)使用多长度约束限制自标定光束平差,平差参数包括主点、主距、径向畸变、偏心畸变、面内畸变、外方位参数及空间点坐标;e)对所述测量系统的定向精度进行可溯源评价。
Description
技术领域
本发明涉及高精度定向及定向精度评价方法,尤其涉及大尺寸动态摄影测量系统的高精度定向及定向精度评价方法。
背景技术
动态摄影测量系统由两台或多台可以高速采集图像的相机组成,相机同步采集被测目标点的图像。通过标定各相机的投影成像模型,可以由像面点恢复空间光线,进而交会计算被测点空间坐标。相机的成像模型包括内方位和外方位参数两类,内方位参数包括主点、主距及各种类型、不同阶次的畸变参数,外方位参数包括相机位于空间坐标系中的位置和指向角。对于外方位参数的求解就是测量系统的定位定向。
目前公认自标定光束平差方法可以得到最高的定位定向精度,因为这种技术在平差过程中考虑内方位参数的变化对于测量结果的影响,得到的内外方位参数结果是同测量环境及测量网络相适应的最优结果。在立体视觉测量领域应用最为广泛的自标定方法就是Tsai的两步标定法和张正友的平面标定法,常用的标定物有立体标定物、平面棋盘格标定板、平面圆点标定板等。这一类标定方法本质是通过标定物和相机之间的相对移动,产生具有若干控制点的测量网络,实现单相机的自标定光束平差。然后在同一坐标系下求解两相机之间的方位关系。
这种方法,不仅机器视觉在使用,结构光测量在使用,所谓的高精度大尺寸立体视觉测量系统也在使用。这一类标定物自标定方法适合于较小的测量范围(通常小于1m×1m),因为标定物的精度同其尺寸是反比关系,当用于小空间的结构光测量或显微测量时,可以得到满意的标定精度。
这种标定物自标定方法自身的缺陷在于:平差过程的重点是内方位参数,对于相机之间的方位关系没有直接的优化和评价;标定板上作为控制点的目标点,其坐标精度会带来系统误差;标定板目标点同测量点在材质和成像质量的差异会带来系统误差;除了像面误差或空间坐标误差估计外,没有客观可靠的空间评价指标。非常不适合应用于具有高精度要求及精度可溯源的大尺寸测量场合。
发明内容
为了解决上述技术问题,本发明提供一种大尺寸动态摄影测量系统的高精度定向及定向精度评价方法,包括步骤:a)选定基准尺,在所述基准尺两端布置编码点对其进行长度计量;b)在测量空间,以第一方式移动所述基准尺并调整所述基准尺姿态,以第二方式对其拍摄图片;c)对所述测量系统进行初步相对定位定向;d)使用多长度约束限制自标定光束平差,平差参数包括主点、主距、径向畸变、偏心畸变、面内畸变、外方位参数及空间点坐标;e)对所述测量系统的定向精度进行可溯源评价。
优选地,所述步骤b的第一方式为:所述编码点法向方向与相机光轴交会角平分线平行,并且所述基准尺位置覆盖被测物体积空间且有深度方向的变化。
优选地,所述步骤b的第二方式为:
b1)均匀划分被测空间,至少具有15个用于划分空间的节点;
b2)移动基准尺至某个节点位置;
b3)在每个节点位置,旋转基准尺,当基准尺与水平方向具有0°,45°,90°和135°时,短暂保持当前姿态;
b4)控制两台相机同时拍摄基准尺图像;
b5)重复步骤b3)和步骤b4),遍历所有空间节点位置,获取所述基准尺不同位置、不同姿态的全局平差测量网络图片;
优选地,所述步骤b5)中重复步骤b3)和步骤b4)至少60次,以遍历所有空间节点位置。
优选地,所述步骤c包括步骤:
c1)利用五点法确定基础矩阵的解空间。自动选择所述基准尺位置中的5个成像点,所述成像点位置为所述拍摄图片的中心和四个顶点;所述自动选择的算法为:
abs(xli)+abs(yli)+abs(xri)+abs(yri)→min
(xli)+(yli)+(xri)+(yri)→max
(-xli)+(yli)+(-xri)+(yri)→max
(-xli)+(-yli)+(-xri)+(-yri)→max
(xli)+(-yli)+(xri)+(-yri)→max
其中,[xli yli]表示左相机采集到的第i个点的像面坐标,[xri yri]表示右相机采集到的第i个点的像点坐标。
c2)利用消元法获取本质矩阵的解。
优选地,步骤c2)中获取本质矩阵的解包含以下三个关键步骤:
c21)建立关于未知数w的10次多项式;
c22)求解步骤c21中所述10次多项式的实数根;
c23)判断本质矩阵的解。
优选地,所述步骤c21建立关于w的10次多项式的方法中,使用Toeplitz矩阵:
c211)将所有多项式描述为数组形式;
amxm+am-1xm-1+…+a2x2+a1x+a0 A=[am am-1 … a2 a1 a0]T
bnxn+bn-1xn-1+…+b2x2+b1x+b0 B=[bn bn-1 … b2 b1 b0]T
c212)利用Toeplitz矩阵进行编程计算多项式乘法,公式为:
其中,T是对应多项式A的Toeplitz矩阵,T的列数等于B中元素个数,其行数等于(m+1)+(n+1)-1。
优选地,所述步骤c212还可以利用Toeplitz矩阵进行编程计算多项式除法,方法步骤为:
c213)利用d的Toeplitz矩阵描述n为:n=T(d)q+r
其中,设分子多项式为n,分母多项式为d,商为q,余项为r;
c214)求得多项式除法的商的最优解:q=(TTT)-1TTn。
优选地,所述步骤c22求解10次多项式的实数根的方法为:
c221)建立十次多项式的Sturm序列,公式为:
其中,rem表示对多项式求余数,P(x)为已知十次多项式。
c222)通过一种递归的二分法,结合Sturm定理,搜索多项式所有的单根区间。
c223)在得到单根区间后,在各单根区间上使用二分法求解|C(w)|=0的实数根的数值解。
优选地,所述步骤d包括步骤为:
d1)获取相机成像模型;
d2)对所述摄影测量系统建立误差方程;
d3)对误差方程项进行自适应比例调整;
d4)对法方程进行分块运算;
d5)获取带长度约束的自标定光束平差;
d6)对所述测量系统的精度进行评估。
优选地,所述步骤e对所述测量系统的定向精度进行可溯源评价,采用直接线性变换解法,步骤为:
e1)获取单个成像点的误差方程:
e2)当测量系统由两台相机组成时,获取对应一空间点[X Y Z]的误差方程组:
e3)获取相应的法方程及坐标改正量,并通过多次迭代优化空间点坐标:
e4)获取所有方位基准尺长度的测量值:
e5)得到所有方位基准尺长度测量结果的平均值:
e6)以平均长度与计量长度之间的比值为定向结果进行比例缩放;对长度测量结果进行不确定度分析,获得定向精度评价:
e7)获取测量仪器的品质参数:
e8)计算各长度的相对误差:
e9)计算长度测量相对误差的不确定度。
总结上述描述,本发明的大尺寸多相机摄影测量的高精度定位定向方法,测量系统保持稳定,定位定向在测量现场完成,标定物加工和计量简便,标定过程简单。定向过程完成时就可以给出测量系统长度测量误差的可溯源评价,获得了以下技术效果:
1、采用一根带有两个回光反射目标点(或球形回光反射标志点)的基准尺进行多相机摄影测量系统的定位定向,基准尺在测量范围内的多个位置、以不同姿态成像,形成覆盖整个测量范围的长度定向场,不需要知道定向场中的空间点坐标;
2、利用五点法进行相对定向的数值计算方法及利用空间正交长度约束进行合理解判断的技术;
3、多相机自标定光束平差是一种具有多个限制条件的间接平差。多个位置和姿态的基准长度作为三维约束条件控制整个平差过程,可以解耦内外方位参数之间的相关性,实现全参数的自标定光束平差(包括主点、主距、径向畸变、偏心畸变、面内畸变、外方位参数以及空间点坐标);
4、定向过程完成时,通过多个长度测量误差的统计结果对整个定位定向过程进行可溯源的精度评价。
附图说明
为了更清楚的说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单的介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他附图。
图1为根据本发明的大尺寸多相机摄影测量的高精度定位定向方法的实验总成;
图2为所述系统的实验示意图;
图3为所述系统中某个基准尺在两台相机中的像;
图4为位于空间不同位置及姿态的基准尺与两台相机之间的相对位置图。
具体实施方式
通过参考示范性实施例,本发明的目的和功能以及用于实现这些目的和功能的方法将得以阐明。然而,本发明并不受限于以下所公开的示范性实施例,可以通过不同形式来对其加以实现。说明书的实质仅仅是帮助相关领域技术人员综合理解本发明的具体细节。
本发明提供了一种大尺寸动态摄影测量系统,所述系统包括相机101、电源102、网线103、信号线104以及上位机105:所述电源102为所述相机1、相机2以及上位机105供电,网线103a连接于相机1与上位机105之间、网线103b连接于相机2与上位机105之间,用于相机与上位机间的通讯;所述信号线104a同样连接相机1与上位机105、信号线104b连接相机2与上位机105,用于传输上位机105与相机101间的信号;所述上位机105还配置有USB信号卡,上位机105通过USB信号卡106,发出同步信号,通过信号线104,控制两台相机101同步采集图像)
所述系统的实验示意图如图2所示:测量空间201主要由两台相机202和基准尺203组成,该测量空间201提供了该系统的高精度定向及定向精度评价方法的实验验证环境。
本发明提供的一种大尺寸动态摄影测量系统的高精度定向及定向精度评价方法,包括以下步骤:
a)选定基准尺203,在所述基准尺203两端布置编码点对其进行长度测量;
b)在测量空间201,以第一方式移动所述基准尺203并调整所述基准尺姿态,用相机202以第二方式对其拍摄图片;
c)对所述测量系统进行初步相对定位定向;
d)使用多长度约束限制自标定光束平差,平差参数包括主点、主距、径向畸变、偏心畸变、面内畸变、外方位参数及空间点坐标;
e)对所述测量系统的定向精度进行可溯源评价。
该方法的具体步骤描述如下:
步骤b:
根据本发明的一个实施例,所述步骤b的第一方式为:所述编码点法向方向与相机光轴交会角平分线平行,并且所述基准尺位置覆盖被测物体积空间且有深度方向的变化;
根据本发明的一个实施例,所述步骤b的第二方式为:b1)划分被测空间区域,形成至少15个划分节点;b2)移动基准尺至某个节点位置;b3)将所述基准尺更换姿态,在0°,45°,90°,135°角度短暂保持稳定;b4)在每个稳定位置,控制两台相机同时采集基准尺图像;b5)重复步骤b3)和步骤b4),遍历所有空间节点位置,获取所述基准尺不同位置、不同姿态的全局平差测量网络图片;
其中,步骤b5)中重复步骤b3)和步骤b4)至少60次,以遍历所有空间节点位置。
基于步骤b的方法,本发明的基准尺设计及标定操作如下:
基准尺主体为碳纤维杆,杆两端固定两个回光反射编码点作为计量长度的端点。标定过程中,某个位置的基准尺在两台相机中的成像,如图3所示,其中白色点区域为贴于所述基准尺的两端的回光反射点301;
手持并移动基准尺至测量空间的不同位置,调整基准尺姿态,使得编码点法线方向与两相机光轴交会角平分线大致平行,位于公共视场,拍摄基准尺图片。基准尺位置应覆盖被测空间,并有深度方向的变化,如图4所示,三维空间中包含两台相机的坐标402和多个不同位置的基准尺401,回光反射点403贴于基准尺的两端用于标定基准尺长度。拍摄图片数目应不少于60幅。
步骤c:
根据本发明的一个实施例,步骤c包括以下步骤:
c1)利用五点法确定基础矩阵的解空间:自动选择所述基准尺位置中的5个成像点,所述成像点位置为所述拍摄图片的中心和四个顶点;所述自动选择的算法为:
其中,[xli yli]表示左相机采集到的第i个点的像面坐标,[xri yri]表示右相机采集到的第i个点的像点坐标。
满足上面各式描述的点就是与预期位置最为接近的点。
由于对应像点满足:
X′TEX=0 (3.2)
其中,E是本质矩阵,X和X′分别是左右图像上的对应点由各自内参数归一化之后的坐标。
将选择的五点带入上式,得到如下的方程组:
通过SVD分解获得上式的基础解系,从而扩展出本质矩阵的解空间:
E=wEw+xEx+yEy+Ez (3.4)
c2)利用消元法获取本质矩阵的解;
本质矩阵具有如下两个性质:
利用上面的式子可以得到关于未知数w,x,y的10个方程,将其中的x,y,x2,xy,y2,x3,x2y,xy2,y3作为未知量,w作为已知量,改写这10个方程,消去了x和y,得到如下的方程组:
其中,C(w)是关于w的多项式矩阵。
(3.6)式描述的齐次线性方程组有非零解,因此C(w)的行列式为零。理论上,描述出C(w)的行列式表达的多项式,并求解多项式的根,就得到w以及C(w)的数值解,进一步求解(3.6)式的非零解就可以得到x,y的值。其中必然涉及到10阶符号矩阵行列式的多项式描述,以及如何求解高次多项式的实根。
根据本发明的一个实施例,步骤c2)中获取本质矩阵的解包含以下关键步骤;
c21)建立关于未知数w的10次多项式;
行列式的求解通过下面的方法编程实现,对于符号多项式或者数值运算都适合:
上面的算法涉及了多项式的乘法和除法;
根据本发明的一个实施例,步骤c21建立关于w的10次多项式的方法中,使用Toeplitz矩阵进行编程计算,方法为:
c211)将所有多项式描述为数组形式;
c212)利用Toeplitz矩阵进行编程计算多项式乘法,公式为:
其中,T是对应多项式A的Toeplitz矩阵,T的列数等于B中元素个数,其行数等于(m+1)+(n+1)-1。
根据本发明的一个实施例,步骤c212还可以利用Toeplitz矩阵进行编程计算多项式除法,方法步骤为:
c213)利用d的Toeplitz矩阵描述n为:n=T(d)q+r
其中,设分子多项式为n,分母多项式为d,商为q,余项为r;
c214)求得多项式除法的商的最优解:q=(TTT)-1TTn。
利用Toeplitz矩阵求解多项式除法的方法不存在长除法的误差累计现象,而且可以避免大数除小数的病态问题,得到的结果是一种最小二乘意义下的最优解。
c22)求解步骤c21中所述10次多项式的实数根;
高次多项式实数根的求解只能依靠数值解法,例如牛顿法、两分法等等,但是关键问题在于确定单根区间。
根据本发明的一个实施例,本发明中利用Sturm定理进行单根区间的搜索,已知一个多项式P(x),可以确定它的Sturm序列为:
其中,rem表示对多项式求余数。
Sturm定理描述了在区间(a,b]内,多项式P(x)实根的个数。设在端点a处Sturm序列函数值符号变化次数为c(a)次,在端点b处Sturm序列函数值符号变化次数为c(b)次,那么多项式P(x)在此半开区间中实根的个数为c(a)-c(b)。本发明提出了一种递归的二分技术,在某个很宽的数值范围(a,b]内搜索所有的单根区间,相应函数的伪代码如下所示:
确定多项式P(x)的Sturm序列sturmSequence;
function[rootIntervals]=rootIntervalsDetect(a,b,sturmSequence)
通过Sturm定理确定(a,b]内实根的个数numRoot;
在得到单根区间后,在各单根区间上使用二分法求解|C(w)|=0的实数根。将求解得到的w带入公式(3.6)中C(w)的表达式,通过SVD分解,求解关于x,y的未知项,并进一步得到x和y的值。将解得的w,x,y带入公式(3.4),得到本质矩阵E,其解的数目同w实根的数目相同。
c23)判断本质矩阵的解。
每个w都存在一组相机位置关系以及对应的重建空间点坐标,在特殊情况下,例如对平面场景成像,传统的像面误差最小的标准不能用来判别合理解,因为存在像面误差近似的两个解。为了从根本上解决根的判别问题,本发明利用空间信息作为约束,以两个正交位置的基准尺长度比例作为判别标准。重建的空间长度比最接近1的本质矩阵是最终的合理解。
进一步通过空间长度与重建长度之间的比例对平移量和空间坐标进行尺度缩放。得到的外方位参数、空间坐标可以作为下一步自标定光束平差的参数初值。
步骤d:
根据本发明的一个实施例,步骤d包括步骤为:
d1)获取相机成像模型;
通过[X0 Y0 Z0 Az El Ro]T描述相机的方位,空间点坐标为[X Y Z]T在相机坐标系下的坐标[X Y Z]T描述为:
其中,
此空间点的线性投影坐标为:
它相对于成像系统的物距为:
研究表明,一物距为s′的空间点在成像系统像面上的径向畸变参数为:
其中,s是成像系统对焦距离,以及是两个距离上径向畸变参数的标定结果,和分别是两个距离对应高斯成像公式的像距。
令相机的主点偏移为xp和yp,对应此空间点的偏心畸变及面内畸变参数为P1,P2,B1,B2。位于(xl,yl)处像点畸变量为:
于是空间点在此站位相机像面上最终的像点坐标为:
外方位参数与空间坐标通过畸变参数通过s′与径向畸变参数产生相关性。
d2)对所述摄影测量系统建立误差方程;
对于一双站位摄影测量系统,用(xij,yij)表示第i张相片对于空间第j点的成像坐标,那么误差方程为:
其中,Xij 0是所有与成像点(xij,yij)相关的参数,(vxij,vyij)是残余误差,Aij和Bij是成像模型(4.7)对于第i图片的外方位参数、第j空间点坐标以及相机成像参数的雅可比矩阵。
Aij中各项通过下面的过程求解:
由公式(4.4)描述的关系,可知:
通过公式(4.3)的描述,可以求解xl和yl对于X0的偏导数为:
于是,结合公式(4.7),观测值对于外方位参数中X0项的偏导数为:
同理,可以得到观测值对于其他外方位参数的偏导数:
线性项xl和yl对于各角度的偏导数描述比较复杂,可以参照相关文献。这里仅就Δx和Δy对于各角度的偏导数进行分析。
本发明使用单边自标定光束平差,参与平差的内参数是xp,yp,f,P1,P2,B1,B2以及S2单边径向畸变参数
Bij中各项通过下面的过程求解:
d3)对法方程进行分块运算;
设一个双相机测量系统对于n个基准尺位置进行摄影测量,那么误差方程组为:
式中,符号下标i,jk表示第i台相机,第j个基准长度上的第k个端点(i=1,2,j=1,2,...n,k=1,2)。
令vTv→min的待求参数增量δ通过下式求解:
鉴于A的稀疏性,ATA及ATl都可以通过对于图片和目标点的索引进行有规则的分块描述:
d4)获取带长度约束的自标定光束平差;
限制条件通过下面的公式描述:
式中,L是长度计量值,Lj是测量系统重建出的第j个位置的基准尺长度,长度对于相应端点空间坐标偏导数的求解方法为:
那么,公式(4.31)在公式(4.32)描述的约束条件下的最小二乘解及相应的联系数向量为:
d5)对所述测量系统的精度进行评估;
测量系统的单位权中误差为:
为了保证定向精度,推荐在60个以上位置摆放基准尺。
两相机的内外方位参数、基准尺特征点空间坐标的误差协方差矩阵为:
d6)对误差方程项进行自适应比例调整;
待求未知量在数量级上相差巨大,尤其是内参数,例如K1、P1、P2、B1和B2项一般具有10-5,K2项具有10-7,而K3具有10-11,这种数量级上的巨大差异带来了误差方程(4.31)中矩阵A的病态问题,由于数值计算中普遍存在的机器精度、大数吃小数以及大数除小数等问题,会严重影响法方程的计算结果,往往会出现矩阵缺秩等问题。为了统一误差方程中各项的数量级,这里设计了自适应的比例调整。
在每次平差迭代过程开始,首先统计误差方程中各列的数量级Mj及调整系数kj,计算方法如下:
然后将每一列的计算结果乘以相应列的比例系数,得到调整后的误差方程:
使用公式(4.31)计算得到的未知量与需要求解的之间的关系为:
仿真和实验都证明,这种自适应比例调整消除了法方程在数值计算方面的病态,带来了更好的平差稳定性及平差精度。
步骤e
定向过程完成,便固定两台相机的内外方位参数,进而通过最小二乘求解待测点的空间坐标。空间坐标求解过程中,由于被测点之间不再有相关性,因此可以逐点求解,避免大型矩阵运算带来的时间消耗,提高动态测量的性能,本质上是一种直接线性变换解法。
根据本发明的一个实施例,对所述测量系统的定向精度进行可溯源评价的具体方法为:
e1)描述单个成像点的误差方程:
稀疏矩阵Bij是像面坐标对于空间坐标的偏导数矩阵,通过公式(4.30)和(4.31)求解。
e2)当测量系统由两台相机组成时,对应一空间点[X Y Z]的误差方程组为:
e3)相应的法方程及坐标改正量为:
e4)定向精度评价是一种长度测量的评价,使用的是定向过程中的数据。利用公式(5.3)重建2n个基准尺端点坐标,进而求解对应的基准尺长度的测量值:
L1 L2 ... Ln
e5)测量结果的平均值为:
e6)以平均长度与计量长度之间的比值为定向结果(外方位参数的平移量)进行比例缩放:
L1′=sL1 L2′=sL2 ... Ln′=sLn (5.6)
对长度测量结果进行不确定度分析,获得定向精度评价:
e7)测量仪器的品质参数由长度不确定度的整数倍描述:
B=ku(Li′) (5.8)
e8)各长度的相对误差为:
e9)长度测量相对误差等于未进行缩放时长度测量的相对误差,因此长度测量相对误差的不确定度也就等于原始测量长度相对误差的不确定度,为:
根据本发明的大尺寸动态摄影测量系统的高精度定向及定向精度评价方法,进行了以下验证性实验:
用于实验的测量系统由两台工业相机组成,如图3所示,用于定向的是一根长度约为950mm的基准尺,基准长度通过尺两端的回光反射编码标志点指示。测量环境的体积约为3.6m×2.4m×2.4m,在测量环境中移动基准尺并成像,利用本发明中的方法完成测量系统的自标定定位定向,并利用权利要求11中所述的评价方法,评估系统的定向精度。三次定向实验的结果如表1和2所示,表1展示了单纯以像面误差最小作为自标定光束平差目标的定向结果,表2展示了本发明介绍的多长度约束的自标定光束平差的定向结果。
表1传统自标定光束平差的定向结果
表2多长度约束的自标定光束平差的定向结果
对比表1和表2三次定向实验的像面误差标准差可以看出,传统的自标定方法将像面x方向的误差减小至可忽略的状态,小于另一个方向的像面误差一个数量级。这明显是不正确的,因为成像在两个方向上应具有近似的误差等级。这种平差方法导致了自标定结果及定位定向结果的错误,相应的长度误差很大。而本发明中介绍的多长度约束的自标定平差技术,由于受到了空间长度的约束,像面误差并非决定平差结果的唯一因素,像面误差更加合理,相应的长度误差水平明显降低。
实验结果表明,以长度测量相对误差作为评价的指标,可以忽略基准长度计量误差带来的分析复杂性。如果长度计量具有可溯源性,那么测量系统也可以通过绝对测量误差进行精度评价,其评价结果也将具有可溯源性。
综上,本发明利用单根基准尺完成双相机大尺寸动态摄影测量系统的现场定位定向。将五点法同空间长度比结合,进行初步定向,能够克服特殊情况下的丢解问题,同时从多解中准确确定合理解。同时,能够利用多个空间长度约束自标定光束平差,克服参数之间的相关性,提高了测量系统整体的定向精度。通过长度测量误差的统计结果提供可溯源的客观评价指标,实验结果标定本发明的长度测量相对误差达到了1/15000。
以上只是本发明较佳的实例,并非来限制本发明实施范围,故凡依本发明申请专利范围所述的构造、特征及原理所做的等效变化或修饰,均应包括于本发明专利申请范围内。
Claims (9)
1.一种大尺寸动态摄影测量系统的高精度定向及定向精度评价方法,包括步骤:
a)选定基准尺,在所述基准尺两端布置编码点对其进行长度测量;
b)在测量空间,以第一方式移动所述基准尺并调整所述基准尺姿态,以第二方式对其拍摄图片,其中所述第一方式为:所述编码点法向方向与相机光轴交会角平分线平行,并且所述基准尺位置覆盖被测物体积空间且有深度方向的变化;
c)对所述测量系统进行初步相对定位定向;
d)使用多长度约束限制自标定光束平差,平差参数包括主点、主距、径向畸变、偏心畸变、面内畸变、外方位参数及空间点坐标;
e)对所述测量系统的定向精度进行可溯源评价。
2.根据权利要求1所述的方法,其特征在于:所述步骤b的第二方式为:
b1)均匀划分被测空间,至少具有15个用于划分空间的节点;
b2)移动基准尺至某个节点位置;
b3)在每个节点位置,旋转基准尺,当基准尺与水平方向具有0°,45°,90°和135°时,短暂保持当前姿态;
b4)控制两台相机同时拍摄基准尺图像;
b5)重复步骤b3)和步骤b4),遍历所有空间节点位置,获取所述基准尺不同位置、不同姿态的全局平差测量网络图片。
3.根据权利要求1所述的方法,其特征在于:所述步骤c包括以下步骤:
c1)利用五点法确定基础矩阵的解空间;自动选择所述基准尺位置中的5个成像点,所述成像点位置为所述拍摄图片的中心和四个顶点;所述自动选择的算法为:
abs(xli)+abs(yli)+abs(xri)+abs(yri)→min
(xli)+(yli)+(xri)+(yri)→max
(-xli)+(yli)+(-xri)+(yri)→max
(-xli)+(-yli)+(-xri)+(-yri)→max
(xli)+(-yli)+(xri)+(-yri)→max
其中,[xli yli]表示左相机采集到的第i个点的像面坐标,[xri yri]表示右相机采集到的第i个点的像点坐标;
c2)利用消元法获取本质矩阵的解。
4.根据权利要求3所述的方法,其特征在于:所述步骤c2)中获取本质矩阵的解包含以下关键步骤;
c21)建立关于未知数w的10次多项式,其中w为c1)算法中通过SVD分解获得的基础解系,扩展出本质矩阵的解空间:E=wEw+xEx+yEy+Ez,中的未知数;
c22)求解步骤c21中所述10次多项式的实数根;
c23)判断本质矩阵的解。
5.根据权利要求4所述的方法,其特征在于:所述步骤c21建立关于w的10次多项式的方法中,使用Toeplitz矩阵:
c211)将所有多项式描述为数组形式;
amxm+am-1xm-1+…+a2x2+a1x+a0A=[am am-1 … a2 a1 a0]T
bnxn+bn-1xn-1+…+b2x2+b1x+b0B=[bn bn-1 … b2 b1 b0]T
c212)利用Toeplitz矩阵进行编程计算多项式乘法,公式为:
<mrow>
<mi>C</mi>
<mo>=</mo>
<mi>T</mi>
<mo>*</mo>
<mi>B</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>a</mi>
<mi>m</mi>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mo>...</mo>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>a</mi>
<mrow>
<mi>m</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
<mtd>
<msub>
<mi>a</mi>
<mi>m</mi>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mo>...</mo>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>a</mi>
<mrow>
<mi>m</mi>
<mo>-</mo>
<mn>2</mn>
</mrow>
</msub>
</mtd>
<mtd>
<msub>
<mi>a</mi>
<mrow>
<mi>m</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
<mtd>
<msub>
<mi>a</mi>
<mi>m</mi>
</msub>
</mtd>
<mtd>
<mo>...</mo>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mi>a</mi>
<mn>0</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<mo>...</mo>
</mtd>
<mtd>
<msub>
<mi>a</mi>
<mi>m</mi>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
<mtd>
<mtable>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
</mtable>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mi>a</mi>
<mn>0</mn>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>b</mi>
<mi>n</mi>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>b</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>b</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>b</mi>
<mn>0</mn>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中,T是对应多项式A的Toeplitz矩阵,T的列数等于B中元素个数,其行数等于(m+1)+(n+1)‐1。
6.根据权利要求5所述的方法,其特征在于:所述步骤c212还可以利用Toeplitz矩阵进行编程计算多项式除法,方法步骤为:
c213)利用d的Toeplitz矩阵描述n为:n=T(d)q+r
其中,设分子多项式为n,分母多项式为d,商为q,余项为r;
c214)求得多项式除法的商的最优解:q=(TTT)-1TTn。
7.根据权利要求4所述的方法,其特征在于:所述步骤c22中求解10次多项式的实数根的方法为:
c221)建立十次多项式的Sturm序列,公式为:
p0(x)=p(x)
p1(x)=p′(x)
p2(x)=-rem(p0,p1)
p3(x)=-rem(p1,p2)
0=-rem(pm-1,pm)
其中,rem表示对多项式求余数,P(x)为已知十次多项式;
c222)通过一种递归的二分法,结合Sturm定理,搜索多项式所有的单根区间;
c223)在得到单根区间后,在各单根区间上使用二分法求解|C(w)|=0的实数根的数值解。
8.根据权利要求1所述的方法,其特征在于:所述步骤d包括步骤为:
d1)获取相机成像模型;
d2)对所述摄影测量系统建立误差方程;
d3)对误差方程项进行自适应比例调整;
d4)对法方程进行分块运算;
d5)获取带长度约束的自标定光束平差;
d6)对所述测量系统的精度进行评估。
9.根据权利要求1所述的方法,其特征是:所述步骤e对所述测量系统的定向精度进行可溯源评价,采用直接线性变换解法,步骤为:
e1)获取单个成像点的误差方程:
e2)当测量系统由两台相机组成时,获取对应一空间点[X Y Z]的误差方程组:
e3)获取相应的法方程及坐标改正量,并通过多次迭代优化空间点坐标:
e4)获取所有方位基准尺长度的测量值:
e5)得到所有方位基准尺长度测量结果的平均值:
e6)以平均长度与计量长度之间的比值为定向结果进行比例缩放;对长度测量结果进行不确定度分析,获得定向精度评价:
e7)获取测量仪器的品质参数:
e8)计算各长度的相对误差:
e9)计算长度测量相对误差的不确定度。
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610297902.5A CN105910584B (zh) | 2016-05-06 | 2016-05-06 | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 |
CN201810319004.4A CN108801218B (zh) | 2016-05-06 | 2016-05-06 | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610297902.5A CN105910584B (zh) | 2016-05-06 | 2016-05-06 | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 |
Related Child Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810319004.4A Division CN108801218B (zh) | 2016-05-06 | 2016-05-06 | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105910584A CN105910584A (zh) | 2016-08-31 |
CN105910584B true CN105910584B (zh) | 2018-05-04 |
Family
ID=56748538
Family Applications (2)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610297902.5A Active CN105910584B (zh) | 2016-05-06 | 2016-05-06 | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 |
CN201810319004.4A Active CN108801218B (zh) | 2016-05-06 | 2016-05-06 | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 |
Family Applications After (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810319004.4A Active CN108801218B (zh) | 2016-05-06 | 2016-05-06 | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 |
Country Status (1)
Country | Link |
---|---|
CN (2) | CN105910584B (zh) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108182709A (zh) * | 2017-12-28 | 2018-06-19 | 北京信息科技大学 | 一种相机标定和相对定向的方法 |
CN109559355B (zh) * | 2018-12-04 | 2021-08-10 | 北京航空航天大学 | 一种基于相机组的无公共视场的多相机全局标定装置及方法 |
CN110285827B (zh) * | 2019-04-28 | 2023-04-07 | 武汉大学 | 一种距离约束的摄影测量高精度目标定位方法 |
CN110260822B (zh) * | 2019-06-18 | 2021-01-19 | 西安交通大学 | 一种多目结构光系统高精度标定方法 |
CN112985455B (zh) * | 2019-12-16 | 2022-04-26 | 武汉四维图新科技有限公司 | 定位定姿系统的精度评定方法、装置及存储介质 |
CN114459345B (zh) * | 2021-12-22 | 2024-05-24 | 上海智能制造功能平台有限公司 | 基于视觉空间定位的飞机机身位置姿态检测系统及方法 |
Family Cites Families (15)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7609846B2 (en) * | 2004-07-13 | 2009-10-27 | Eastman Kodak Company | Matching of digital images to acquisition devices |
US20060125920A1 (en) * | 2004-12-10 | 2006-06-15 | Microsoft Corporation | Matching un-synchronized image portions |
CN101644563B (zh) * | 2009-08-18 | 2010-12-08 | 北京信息科技大学 | 基于距离约束拟合点的视觉测量系统不确定度评价方法 |
CN101694370B (zh) * | 2009-09-15 | 2011-09-21 | 北京信息科技大学 | 大尺寸工业摄影测量系统的空间误差场获取方法及基准装置 |
CN101727670B (zh) * | 2009-11-10 | 2012-01-04 | 西安交通大学 | 一种可变幅面多相机系统柔性标定方法及装置 |
CN102175261B (zh) * | 2011-01-10 | 2013-03-20 | 深圳大学 | 一种基于自适应标靶的视觉测量系统及其标定方法 |
JP5991821B2 (ja) * | 2012-02-03 | 2016-09-14 | 三菱電機株式会社 | 写真測量装置 |
CN102721376B (zh) * | 2012-06-20 | 2014-12-31 | 北京航空航天大学 | 一种大视场三维视觉传感器的标定方法 |
JP5715735B2 (ja) * | 2012-06-29 | 2015-05-13 | 富士フイルム株式会社 | 3次元測定方法、装置、及びシステム、並びに画像処理装置 |
CN102889882B (zh) * | 2012-09-03 | 2014-11-12 | 北京信息科技大学 | 一种基于光束平差的三维重建方法 |
CN103389072B (zh) * | 2013-07-22 | 2015-04-22 | 北京信息科技大学 | 一种基于直线拟合的像点定位精度评估方法 |
CN103544699B (zh) * | 2013-10-11 | 2017-02-01 | 国家电网公司 | 一种基于单图片三圆模版的摄像头标定方法 |
CN104036542B (zh) * | 2014-05-21 | 2017-01-25 | 北京信息科技大学 | 一种基于空间光线聚集性的像面特征点匹配方法 |
CN104019799B (zh) * | 2014-05-23 | 2016-01-13 | 北京信息科技大学 | 一种利用局部参数优化计算基础矩阵的相对定向方法 |
CN104091345B (zh) * | 2014-07-24 | 2017-01-25 | 中国空气动力研究与发展中心高速空气动力研究所 | 基于前方交会约束的五点相对定向方法 |
-
2016
- 2016-05-06 CN CN201610297902.5A patent/CN105910584B/zh active Active
- 2016-05-06 CN CN201810319004.4A patent/CN108801218B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN108801218B (zh) | 2021-07-02 |
CN108801218A (zh) | 2018-11-13 |
CN105910584A (zh) | 2016-08-31 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105910584B (zh) | 大尺寸动态摄影测量系统的高精度定向及定向精度评价方法 | |
US9857172B1 (en) | Method for implementing high-precision orientation and evaluating orientation precision of large-scale dynamic photogrammetry system | |
CN107883870A (zh) | 基于双目视觉系统和激光跟踪仪测量系统的全局标定方法 | |
CN108844459A (zh) | 一种叶片数字化样板检测系统的标定方法及装置 | |
CN108648242B (zh) | 基于激光测距仪辅助无公共视场的两相机标定方法及装置 | |
CN108198219B (zh) | 用于摄影测量的相机标定参数的误差补偿方法 | |
CN104677277B (zh) | 一种测量物体几何属性或距离的方法及系统 | |
CN108198224A (zh) | 一种用于立体视觉测量的线阵相机标定装置及标定方法 | |
CN110398208A (zh) | 基于摄影全站仪系统的大数据变形监测方法 | |
CN105066962A (zh) | 一种多分辨率大视场角高精度摄影测量装置 | |
CN112902930B (zh) | 一种自动计算光束法区域网平差初值的方法 | |
CN112229321B (zh) | 基于lasso算法求解三坐标测量机21项几何误差的方法 | |
CN106671081B (zh) | 一种基于单目视觉的少自由度机器人运动学标定方法 | |
CN114170321A (zh) | 一种基于测距的相机自标定方法及系统 | |
CN103258327A (zh) | 一种基于二自由度摄像机的单点标定方法 | |
CN105528788B (zh) | 相对位姿参数的标定方法、装置和确定三维形状的装置 | |
CN112116665A (zh) | 一种结构光传感器标定方法 | |
CN106023146B (zh) | 用于摄影测量中的场相关单边自标定光束平差方法 | |
CN108871204B (zh) | 摄影测量中长度测量相对误差的绝对评价方法 | |
CN208061260U (zh) | 一种用于立体视觉测量的线阵相机标定装置 | |
Ma et al. | Research on a precision and accuracy estimation method for close-range photogrammetry | |
KR100823070B1 (ko) | 3차원 계측을 위한 다수 교정 영상을 이용한 카메라 교정방법 | |
Tate et al. | Mastcam-z geometric calibration: an alternative approach based on photogrammetric and affine solutions specific to filter, focus, and zoom | |
Borghese et al. | Statistical comparison of DLT versus ILSSC in the calibration of a photogrammetric stereo-system | |
Sun et al. | 3D Estimation of Single Image based on Homography Transformation |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |