CN117606447B - 视场分割型航空面阵相机的试验场检校方法和装置 - Google Patents
视场分割型航空面阵相机的试验场检校方法和装置 Download PDFInfo
- Publication number
- CN117606447B CN117606447B CN202311242030.9A CN202311242030A CN117606447B CN 117606447 B CN117606447 B CN 117606447B CN 202311242030 A CN202311242030 A CN 202311242030A CN 117606447 B CN117606447 B CN 117606447B
- Authority
- CN
- China
- Prior art keywords
- camera
- test field
- image
- sub
- checked
- 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
- 238000012360 testing method Methods 0.000 title claims abstract description 126
- 238000000034 method Methods 0.000 title claims abstract description 50
- 230000011218 segmentation Effects 0.000 title abstract description 18
- 238000004364 calculation method Methods 0.000 claims abstract description 14
- 239000013598 vector Substances 0.000 claims description 78
- 238000012937 correction Methods 0.000 claims description 40
- 230000008859 change Effects 0.000 claims description 29
- 239000011159 matrix material Substances 0.000 claims description 27
- 230000014509 gene expression Effects 0.000 claims description 23
- 239000003550 marker Substances 0.000 claims description 16
- 238000005259 measurement Methods 0.000 claims description 10
- 238000012795 verification Methods 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 4
- 238000001514 detection method Methods 0.000 claims description 4
- 238000003860 storage Methods 0.000 claims description 3
- 230000001915 proofreading effect Effects 0.000 claims 1
- 230000009897 systematic effect Effects 0.000 abstract description 7
- 238000009826 distribution Methods 0.000 description 7
- 238000013461 design Methods 0.000 description 4
- 238000003384 imaging method Methods 0.000 description 4
- 238000013459 approach Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000005192 partition Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 208000011580 syndromic disease Diseases 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
- G01C11/02—Picture taking arrangements specially adapted for photogrammetry or photographic surveying, e.g. controlling overlapping of pictures
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Remote Sensing (AREA)
- Theoretical Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Algebra (AREA)
- Mathematical Optimization (AREA)
- Multimedia (AREA)
- Pure & Applied Mathematics (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Manufacturing & Machinery (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种视场分割型航空面阵相机的试验场检校方法和装置，所述方法包括：根据待校验相机信息和校验飞机飞行高度信息，确定试验场信息；利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息；利用所述试验场子影像集合信息和所述待校验相机信息，构建得到检校模型；对所述检校模型进行求解，得到所述待校验相机的校正参数。本发明方法能够克服附加参数自检校光束平差中附加参数自身的强的相关性，同时也能克服附加参数与内方位/外方位元素等其它未知数的强相关性所导致的计算不收敛难题，从而准确解算相机参数并完成***误差补偿，为生成无畸变全幅大面阵影像提供精准参数。
Description
技术领域
本发明涉及航空摄影测量技术领域，尤其涉及一种视场分割型航空面阵相机的试验场检校方法和装置。
背景技术
航空摄影测量***是获取大比例尺地理空间信息的主要手段，新一代航空测量相机技术发展要求同时满足大视场和高分辨率特性，视场分割集成面阵综合高分辨率成像体制是主要技术途径之一，此类相机的原始摄影数据为含有各类畸变的分块子影像，在利用分块子影像生成合成影像之前，需要获取相机的精准检校参数，以满足严苛的行业应用技术与规范要求。
视场分割型高分辨率航空面阵相机在完成出厂检校后，可交付用户使用。由于实际航飞作业的温湿压、震动及电磁等环境与厂商生产制造条件差异较大，相机参数往往会发生变化，这就需要用户对相机定期性的进行用户级检校。
试验场动态检校是用户采用的保证大幅面航空测量相机高精度量测性能的基本方法，对于新技术机制的相机***和实际地形条件需要设计不同的试验场检校策略。如何克服附加参数自检校光束平差中附加参数自身的强相关性，同时也克服附加参数与内方位/外方位元素等其它未知数的强相关性所导致的计算不收敛难题，从而准确解算相机参数并完成***误差补偿，为生成无畸变全幅大面阵影像提供精准参数，是当前需要解决的问题。
发明内容
针对如何克服附加参数自检校光束平差中附加参数自身的强的相关性，同时也克服附加参数与内方位/外方位元素等其它未知数的强相关性所导致的计算不收敛难题，从而准确解算相机参数并完成***误差补偿，为生成无畸变全幅大面阵影像提供精准参数，本发明公开了一种视场分割型航空面阵相机的试验场检校方法和装置。
本发明公开了一种视场分割型航空面阵相机的试验场检校方法，包括：
S1，根据待校验相机信息和校验飞机飞行高度信息，确定试验场信息；
S2，利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息；
S3，利用所述试验场子影像集合信息和所述待校验相机信息，构建得到检校模型；
S4，对所述检校模型进行求解，得到所述待校验相机的校正参数。
所述根据待校验相机信息和校验飞机飞行高度信息，确定试验场信息，包括：
S11，利用试验场信息生成模型，对所述待校验相机信息和校验飞机飞行高度信息进行处理，得到试验场场地信息；所述待校验相机信息，包括相机像幅宽度L、相机焦距f、相机像元尺寸Cs；所述试验场场地信息，包括试验场边长的下界、标志点尺寸值、标志点间隔值和标志点数目下限；
所述试验场信息生成模型，其表达式为：
W＝L×H/f，
GSD＝Cs×H/f，
其中，H为校验飞机飞行高度，W为试验场边长的下界；所述试验场为矩阵区域；GSD为试验场地面分辨率，R为试验场中的标志点尺寸值，S为试验场中的标志点间隔值，Tn为试验场中的标志点数目下限；所述标志点，为试验场中按照预设间隔S进行等间隔布设的标志物；
S12，根据所述试验场场地信息，在试验场中安装标志点；对每个安装后的标志点进行位置测量，得到标志点大地位置坐标；所述标志点大地位置坐标，为所述标志点在大地坐标系下的坐标；
S13，利用所述试验场场地信息和标志点大地位置坐标，构建得到试验场信息。
所述利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息，包括：
利用校验飞机搭载所述待校验相机，在校验飞机按照预设飞行高度进行飞行时，利用所述待校验相机对所述试验场进行拍摄，得到若干个子影像信息；每个子影像，对应于所述试验场的部分区域；所述子影像信息，包括子影像位置信息和子影像图像信息；
对每个子影像进行标志点提取，得到标志点图像位置坐标；所述标志点图像位置坐标，为所述标志点在像元坐标系下的位置坐标；
对所有子影像信息和其对应的标志点图像位置坐标进行合并处理，得到试验场子影像集合信息。
所述检校模型，包括i0×j0个校验表达式；每个校验表达式，与子影像的一个标志点相对应；对于第i个子影像的第j个标志点，其校验表达式为：
其中，i0为所述待校验相机采集得到的子影像数目，j0为一个子影像所包含的标志点数目，dx0、dy0、df分别为所述待校验相机主点在像元坐标系下的x轴位置坐标改正系数、y轴位置坐标改正系数、所述待校验相机的焦距改正系数，K1，K2，K3分别代表所述待校验相机的径向畸变改正系数，P1，P2分别代表所述待校验相机的切向畸变改正系数，b1，b2分别代表所述待校验相机的焦面探测器件的非正交变形改正系数；A0i，A1i，A2i，A3i，A4i，A5i，B0i，B1i，B2i，B3i，B4i，B5i为第i个子影像的像方二阶验后补偿系数；(xi，j，yi，j)为第i个子影像的第j个标志点的在像元坐标系下的位置坐标，x0、y0、f分别为所述待校验相机主点在像元坐标系下的x轴位置坐标、y轴位置坐标和焦距的初始值，(Xij，Yij，Zij)为所述第i个子影像的第j个标志点的在大地坐标系下的位置坐标，(XS，YS，ZS)为所述待校验相机对所述试验场进行拍摄时刻，所述待校验相机在大地坐标系下的位置坐标，(xi，yi)为第i个子影像的左上角点在所有子影像所构成的像元坐标系下的位置坐标，ri为第i个子影像的相对距离值，rij表示旋转矩阵的第i行、第j列的元素，i＝1，2，3，j＝1，2，3；所述旋转矩阵，根据待校验相机的摄影姿态角决定；在所述检校模型中，dx0、dy0、df、K1、K2、K3、P1、P2、b1、b2、(XS，YS，ZS)、r11、r12、r13、r21、r22、r23、r31、r32、r33、A0i、A1i、A2i、A3i、A4i、A5i、B0i、B1i、B2i、B3i、B4i、B5i为未知量，(xi，j，yi，j)、x0、y0、f、(Xij，Yij，Zij)、(xi，yi)、ri为已知量。
所述对所述检校方程进行求解，得到所述待校验相机的校正参数，包括：
利用所述检校模型的每个校验表达式，构建得到对应的校验方程；每个校验方程，与子影像的一个标志点相对应；
对于第i个子影像的第j个标志点，其对应的校验方程的表达式为：
其中，fij1和fij2为第i个子影像的第j个标志点的校验方程，v为未知量向量，v＝[A0i、A1i、A2i、A3i、A4i、A5i、B0i、B1i、B2i、B3i、B4i、B5i、r11、r12、r13、r21、r22、r23、r31、r32、r33、XS、YS、ZS、dx0、dy0、df、K1、K2、K3、P1、P2、b1、b2]；
对所有的校验方程进行合并处理，得到校验方程组；所述校验方程组，包括与每个标志点相对应的校验方程；
利用迭代求解模型，对所述校验方程组进行求解，得到所述待校验相机的校正参数。
所述利用迭代求解模型，对所述校验方程组进行求解，得到所述待校验相机的校正参数，包括：
S401，对未知量向量进行初始化，得到原始向量；预设变化向量；
S402，利用所述变化向量对所述原始向量进行叠加，得到变化叠加向量；
S403，利用校验方程组的每个校验方程，分别对所述变化叠加向量和原始向量进行计算处理，得到所述校验方程的第一输出值和第二输出值；将所述第一输出值减去第二输出值，得到所述校验方程的差值列向量；
S404，对校验方程组的每个校验方程，在所述原始向量处进行一阶导数计算，得到每个校验方程的一阶导数矩阵；
S405，对所有校验方程的差值列向量，按照列方向进行拼接，得到总差值列向量L；
对所有校验方程的一阶导数矩阵，按照列方向进行拼接，得到总导数矩阵A；
S406，利用总差值列向量L和总导数矩阵A，计算得到变化向量的更新值Δv1，其计算表达式为：
Δv1＝[ATA]-1[-ATL]，
判断Δv1的模值是否小于设定的判别阈值，若小于设定的判别阈值，利用变化向量的更新值与原始向量进行相加，得到相加结果向量，将所述相加结果向量作为所述待校验相机的校正参数；若不小于设定的判别阈值，利用变化向量的更新值，作为变化向量，利用变化叠加向量对原始向量进行更新，执行S402。
所述视场分割型航空面阵相机的试验场检校方法，在对所述检校方程进行求解，得到所述待校验相机的校正参数后，还包括：
利用影像位置坐标校正模型，对子影像的像素点的在像元坐标系下的位置坐标进行计算处理，得到子影像的像素点的校正坐标。
本发明公开了一种视场分割型航空面阵相机的试验场检校装置，所述装置包括：
存储有可执行程序代码的存储器；
与所述存储器耦合的处理器；
所述处理器调用所述存储器中存储的所述可执行程序代码，执行所述的视场分割型航空面阵相机的试验场检校方法。
本发明公开了一种计算机可存储介质，所述计算机存储介质存储有计算机指令，所述计算机指令被调用时，用于执行所述的视场分割型航空面阵相机的试验场检校方法。
本发明公开了一种信息数据处理终端，所述信息数据处理终端用于实现所述的视场分割型航空面阵相机的试验场检校方法。
本发明的有益效果为：
本发明公开了一种飞行条件下满足几何参数检校的试验场设计方案和视场分割型航空面阵相机的方法，能够克服附加参数自检校光束平差中附加参数自身的强的相关性，同时也能克服附加参数与内方位/外方位元素等其它未知数的强相关性所导致的计算不收敛难题，从而准确解算相机参数并完成***误差补偿，为生成无畸变全幅大面阵影像提供精准参数。
附图说明
图1为本发明方法的实施流程图；
图2为本发明的校验飞机搭载待校验相机进行拍摄的示意图；
图3为某检校场标志点分布图。
具体实施方式
为了更好的了解本发明内容，这里给出两个实施例。
图1为本发明方法的实施流程图；图2为本发明的校验飞机搭载待校验相机进行拍摄的示意图；图3为某检校场标志点分布图。
实施例一：
本发明公开了一种视场分割型航空面阵相机的试验场检校方法，包括：
S1，根据待校验相机信息和校验飞机飞行高度信息，确定试验场信息；
S2，利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息；
S3，利用所述试验场子影像集合信息和所述待校验相机信息，构建得到检校模型；
S4，对所述检校模型进行求解，得到所述待校验相机的校正参数。
所述根据待校验相机信息和校验飞机飞行高度信息，确定试验场信息，包括：
S11，利用试验场信息生成模型，对所述待校验相机信息和校验飞机飞行高度信息进行处理，得到试验场场地信息；所述待校验相机信息，包括相机像幅宽度L、相机焦距f、相机像元尺寸Cs；所述试验场场地信息，包括试验场边长的下界、标志点尺寸值、标志点间隔值和标志点数目下限；
所述试验场信息生成模型，其表达式为：
W＝L×H/f，
GSD＝Cs×H/f，
其中，H为校验飞机飞行高度，W为试验场边长的下界；所述试验场为矩阵区域；GSD为试验场地面分辨率，R为试验场中的标志点尺寸值，S为试验场中的标志点间隔值，Tn为试验场中的标志点数目下限；所述标志点，为试验场中按照预设间隔S进行等间隔布设的标志物，用于为视场分割型航空面阵相机的试验场检校提供位置已知的地面点。
S12，根据所述试验场场地信息，在试验场中安装标志点；对每个安装后的标志点进行位置测量，得到标志点大地位置坐标；所述标志点大地位置坐标，为所述标志点在大地坐标系下的坐标。
所述位置测量，可采用差分GPS定位方法；
S13，利用所述试验场场地信息和标志点大地位置坐标，构建得到试验场信息。
所述S2，包括：利用校验飞机搭载待校验相机，按照校验飞机预设航高进行飞行，飞行过程中，使用待校验相机拍摄采集得到试验场子影像集合；
所述利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息，包括：
利用校验飞机搭载所述待校验相机，在校验飞机按照预设飞行高度进行飞行时，利用所述待校验相机对所述试验场进行拍摄，得到若干个子影像信息；每个子影像，对应于所述试验场的部分区域；所述子影像信息，包括子影像位置信息和子影像图像信息；
对每个子影像进行标志点提取，得到标志点图像位置坐标；所述标志点图像位置坐标，为所述标志点在像元坐标系下的位置坐标；所述对每个子影像进行标志点提取，可采用灰度重心法实现；
对所有子影像信息和其对应的标志点图像位置坐标进行合并处理，得到试验场子影像集合信息。
所述对所有子影像信息和其对应的标志点图像位置坐标进行合并处理，得到试验场子影像集合信息，包括：
利用所有子影像信息和其对应的标志点图像位置坐标，构建得到试验场子影像集合信息。
所述检校模型，包括i0×j0个校验表达式；每个校验表达式，与子影像的一个标志点相对应；对于第i个子影像的第j个标志点，其校验表达式为：
其中，i0为所述待校验相机采集得到的子影像数目，j0为一个子影像所包含的标志点数目，dx0、dy0、df分别为所述待校验相机主点在像元坐标系下的x轴位置坐标改正系数、y轴位置坐标改正系数、所述待校验相机的焦距改正系数，K1，K2，K3分别代表所述待校验相机的径向畸变改正系数，P1，P2分别代表所述待校验相机的切向畸变改正系数，b1，b2分别代表所述待校验相机的焦面探测器件的非正交变形改正系数；A0i，A1i，A2i，A3i，A4i，A5i，B0i，B1i，B2i，B3i，B4i，B5i为第i个子影像的像方二阶验后补偿系数；(xi，j，yi，j)为第i个子影像的第j个标志点的在像元坐标系下的位置坐标，x0、y0、f分别为所述待校验相机主点在像元坐标系下的x轴位置坐标、y轴位置坐标和焦距的初始值，(Xij，Yij，Zij)为所述第i个子影像的第j个标志点的在大地坐标系下的位置坐标，(XS，YS，ZS)为所述待校验相机对所述试验场进行拍摄时刻，所述待校验相机在大地坐标系下的位置坐标，(xi，yi)为第i个子影像的左上角点在所有子影像所构成的像元坐标系下的位置坐标，ri为第i个子影像的相对距离值，rij表示旋转矩阵的第i行、第j列的元素，i＝1，2，3，j＝1，2，3；所述旋转矩阵，根据待校验相机的摄影姿态角决定；在所述检校模型中，dx0、dy0、df、K1、K2、K3、P1、P2、b1、b2、(XS，YS，ZS)、r11、r12、r13、r21、r22、r23、r31、r32、r33、A0i、A1i、A2i、A3i、A4i、A5I、B0i、B1i、B2i、B3i、B4i、B5i为未知量，(xI，j，yi，j)、x0、y0、f、(Xij，Yij，Zij)、(xi，yi)、ri为已知量。
所述第i个子影像的第j个标志点的在像元坐标系下的位置坐标，根据步骤S2得到的标志点图像位置坐标来确定；所述第i个子影像的第j个标志点的在大地坐标系下的位置坐标，根据步骤S1得到的标志点大地位置坐标确定。所述待校验相机对所述试验场进行拍摄时刻，所述待校验相机在大地坐标系下的位置坐标，利用所述待校验相机携带的定位装置得到，所述定位装置可以是GPS定位装置或北斗定位装置。所述摄影姿态角，利用所述待校验相机携带的姿态测量装置得到；所述姿态测量装置，可以是IMU测量装置。第i个子影像的左上角点在所有子影像所构成的像元坐标系下的位置坐标的获取，首先对所有子影像图像信息进行拼接，得到整体影像，根据每个子影像对应的成像子面阵在所述待校验相机的成像面阵上的位置信息，确定每个子影像在整体影像的位置信息；利用所有子影像所构成像元坐标系，根据每个子影像在整体影像的位置信息，确定所述子影像的左上角点在所有子影像所构成的像元坐标系下的位置坐标。
所述对所述检校方程进行求解，得到所述待校验相机的校正参数，包括：
利用所述检校模型的每个校验表达式，构建得到对应的校验方程；每个校验方程，与子影像的一个标志点相对应；
对于第i个子影像的第j个标志点，其对应的校验方程的表达式为：
其中，fij1和fij2为第i个子影像的第j个标志点的校验方程，v为未知量向量，v＝[A0i、A1i、A2i、A3i、A4i、A5i、B0i、B1i、B2i、B3i、B4i、B5i、r11、r12、r13、r21、r22、r23、r31、r32、r33、XS、YS、ZS、dx0、dy0、df、K1、K2、K3、P1、P2、b1、b2]；
对所有的校验方程进行合并处理，得到校验方程组；所述校验方程组，包括与每个标志点相对应的校验方程；
利用迭代求解模型，对所述校验方程组进行求解，得到所述待校验相机的校正参数。
所述利用迭代求解模型，对所述校验方程组进行求解，得到所述待校验相机的校正参数，包括：
S401，对未知量向量进行初始化，得到原始向量；预设变化向量；
S402，利用所述变化向量对所述原始向量进行叠加，得到变化叠加向量；
S403，利用校验方程组的每个校验方程，分别对所述变化叠加向量和原始向量进行计算处理，得到所述校验方程的第一输出值和第二输出值；将所述第一输出值减去第二输出值，得到所述校验方程的差值列向量；所述差值列向量的长度为2；
对第i个子影像的第j个标志点，其校验方程的差值列向量表示为[Δfij1，Δfij2]T；Δfij1＝fij1(v0+Δv)-fij1(v0)；Δfij2＝fij2(v0+Δv)-fij2(v0)；
其中，v0表示原始向量，Δv表示变化向量；
S404，对校验方程组的每个校验方程，在所述原始向量处进行一阶导数计算，得到每个校验方程的一阶导数矩阵；所述一阶导数矩阵，为2行28列的矩阵；
S405，对所有校验方程的差值列向量，按照列方向进行拼接，得到总差值列向量L；总差值列向量L的长度为2×i0×j0；
对所有校验方程的一阶导数矩阵，按照列方向进行拼接，得到总导数矩阵A；所述矩阵A的维度为2×i0×j0行28列；
S406，利用总差值列向量L和总导数矩阵A，计算得到变化向量的更新值Δv1，其计算表达式为：
Δv1＝[ATA]-1[-ATL]，
判断Δv1的模值是否小于设定的判别阈值，若小于设定的判别阈值，利用变化向量的更新值与原始向量进行相加，得到相加结果向量，将所述相加结果向量作为所述待校验相机的校正参数；若不小于设定的判别阈值，利用变化向量的更新值，作为变化向量，利用变化叠加向量对原始向量进行更新，执行S402；
所述利用变化叠加向量对原始向量进行更新，是将变化叠加向量的取值作为原始向量的值，继续执行S402；
对于第i个子影像的第j个标志点，其一阶导数矩阵A表示为：
A＝[A1，A2，A3]，
所述视场分割型航空面阵相机的试验场检校方法，在对所述检校方程进行求解，得到所述待校验相机的校正参数后，还包括：
利用影像位置坐标校正模型，对子影像的像素点的在像元坐标系下的位置坐标进行计算处理，得到子影像的像素点的校正坐标。
所述影像位置坐标校正模型，其计算表达式为：
x＝x0-Δx
y＝y0-Δy
其中：校正子影像i中的像素点像的位置坐标为(Ii，Ji)，x，y为子影像的像素点的校正坐标；x0，y0，f为相机的内方位元素(主点的图像坐标系的坐标，焦距)的初始值；a1，b1，c1，a2，b2，c2，a3，b3，c3为摄影姿态角决定的旋转矩阵的元素，a1，b1，c1为所述矩阵的第一行的元素，a2，b2，c2为所述矩阵的第二行的元素，a3，b3，c3为所述矩阵的第三行的元素；旋转矩阵R的表达式为：
其中，ω、k分别为待校验相机的俯仰角、侧滚角和旋偏角。X，Y，Z为第三特征点在大地坐标系下的位置坐标；Xs，Ys，Zs为摄影点在大地坐标系下的位置坐标；Δx，Δy为利用检校参数计算得到的像坐标改正量，主要包括两类：一类是子影像验后补偿***误差改正量Δx1，Δy1：
a0i，a1i，a2i，a3i，a4i，a5i，b0i，b1i，b2i，b3i，b4i为第i个子影像的像方二阶验后补偿系数，其通过相机出厂参数获得，或者，通过相机拍摄地面标校控制点，利用地面标校控制点的像方坐标和大地坐标，进行计算得到，i为子影像序号；a0i＝A0i，依次类推，b4i＝B4i；另一类是全视场附加参数***误差改正量Δx2，Δy2：
其中，Δx0，Δy0，Δf代表视场分割集成面阵相机***的内方位元素改正系数，通过相机出厂参数获得，或者，通过相机拍摄地面标校控制点，利用地面标校控制点的像方坐标和大地坐标，进行计算得到；K1，K2，K3代表视场分割集成面阵相机***的径向畸变改正系数(通过相机出厂参数获得，或者，通过相机拍摄地面标校控制点，利用地面标校控制点的像方坐标和大地坐标，进行计算得到)；P1，P2代表光学***切向畸变改正系数(通过相机出厂参数获得，或者，通过相机拍摄地面标校控制点，利用地面标校控制点的像方坐标和大地坐标，进行计算得到)；b10，b20代表光学传感器非正交变形改正系数(通过相机出厂参数获得，或者，通过相机拍摄地面标校控制点，利用地面标校控制点的像方坐标和大地坐标，进行计算得到)；x′，y′为最终无畸变影像的像坐标，即校正子影像的像素点在拼合后的全幅影像的位置坐标
为了便于第三方后续应用处理，需要消除与传感器相关因素形成的***误差(Δx，Δy)，即
x′-x0＝x-x0-Δx
y′-y0＝y-y0-Δy
Δx＝Δx1+Δx2
Δy＝Δy1+Δy2
对全幅影像位置坐标计算模型进行迭代运算，若前后两次迭代得到的像素点在拼合后的全幅影像的位置坐标的误差，小于设定阈值，则停止迭代，确定x′，y′为最终无畸变影像的像坐标，即校正子影像的像素点在拼合后的全幅影像的位置坐标。
本发明公开了一种视场分割型航空面阵相机的试验场检校装置，所述装置包括：
存储有可执行程序代码的存储器；
与所述存储器耦合的处理器；
所述处理器调用所述存储器中存储的所述可执行程序代码，执行所述的视场分割型航空面阵相机的试验场检校方法。
本发明公开了一种计算机可存储介质，所述计算机存储介质存储有计算机指令，所述计算机指令被调用时，用于执行所述的视场分割型航空面阵相机的试验场检校方法。
本发明公开了一种信息数据处理终端，所述信息数据处理终端用于实现所述的视场分割型航空面阵相机的试验场检校方法。
实施例二：
结合具体检校场景，对本方法实施流程进行了阐述。
对于试验场的确定，要着重考虑相机飞行高度H、像幅宽度L、影像分辨率R/地面采样间隔GSD、地面高程、检校参数、检定精度等因素，这里以国产视场分割型DMZ II相机选定的某通用航空机场区域的试验场为例进行说明：
1)首先，试验场范围必须大于像幅覆盖范围，即L×H/f；
2)考虑检定精度需要，要求航摄飞机工作相对高度1500米，地面采样间隔GSD为5厘米；
3)为了实现0.1像素水平级别检校精度，要求标志点大地坐标测量精度优于0.5厘米；
4)为了高精度量测影像上标志点，要求标志影像面积大小不少于25像素；
5)标志点间隔取决于像幅范围验后补偿分区大小和控制点要求个数。
最终设计结果为，试验场大小1200米×1200米，布设13×13共169个标志点，点间距约为100米，以满足DMZ II相机集成面阵分区补偿参数检校需求。同时由于GSD达到5厘米，因此对标志点位精度要求X、Y、Z三个方向分别为0.5、0.5、0.5厘米，以保证控制精度满足检校精度要求。最后某检校场标志点分布如图3所示。与设计分布的主要区别在于，根据实际点位布设区域地物分布情况，按需进行位置调整，调整幅度不影响检校要求。
试验场建设是在已设计好的地面范围内进行几何标志点的埋设固定、测量和数据处理。标志点测量，平面采用GNSS基站测量及差分动态基站数据采集，要求平面精度0.5cm；高程进行水准测量，要求达到三等高程水准测量0.5cm精度。标志点形状设计为白色圆形，标志直径0.45m，外接黑色正方形，以突出航空摄影影像上的标志点，边长0.6m，正方形、靶标中心标志和中心开孔为标志中心，开孔圆为直径10mm，便于固定并保持定位精度。
在确定试验场和标志点后，利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息；利用所述试验场子影像集合信息和所述待校验相机信息，构建得到检校模型；对所述检校模型进行求解，得到所述待校验相机的校正参数。
相对于仅带附加参数的自检校平差方法，本发明的顾及验后补偿的自检校法具有两方面优势，一是能够根据***误差特性灵活设置补偿***模型，可以从一阶线性补偿扩展到高阶拟合补偿，抑制不同特性***误差影响；二是能够较好规避自检校平差过度参数化问题，消除附加参数过多导致系数间强相关引起的不收敛或解算失败问题，对各类不同特性的***误差参数均能有效解算。
试验数据采用DMZ II相机航空摄影数据，摄影时间为2021年5月，摄影区域为试验场区域。相对航高1500米，地面分辨率为5厘米，其它摄影参数参见表1。分别采用交叉航线影像(片号1067#，1090#，1044#，1111#)进行检校，计算结果如下表所示。从表2中可见，摄影方位角90°的1067#影像，经过本方法的检校平差处理，量测误差从x和y两个方向0.943、1.14像素降低至0.415和0.472像素，采用三参数验后补偿达到0.157、0.174像素，采用四参数达到0.128和0.147像素，采用六参数达到0.107和0.123像素；摄影方位角180°的1090#影像，自检校后中误差两个方向从1.058、1.158像素降低至0.618和0.562像素，采用三参数验后补偿法量测误差达到0.125、0.165像素，采用四参数补偿达到0.119和0.147像素，采用六参数补偿达到0.100和0.118像素。270°和360°摄影方位角的1044#和1111#影像均能实现***误差有效检校，检校精度达到0.1像素水平，与标志点坐标量测精度一致，达到相机***检校理论精度水平。采用多片联合检校，1067#和1090#中误差分别为0.199和0.162像素，四片联合检校中误差为0.259和0.193像素，随着片数增加检校误差有所增大，主要因素是标志点量测误差累积所致，而检校结果基本一致、计算流程更为高效，具体比较可见后续分析。
表1试验场航摄设计参数表
/>
表2动态检校结果统计
根据1067#和1090#影像检校前后全视场范围标志点残差分布结果，未检校之前，各影像均存在明显***性成像几何误差，经过检校后，不同影像的控制点上残余误差随机分布，***性误差已经得到有效校正。
检校精度情况分析以1067#影像为例，整幅影像上138个标志点，自检校后标志点像坐标残差结果统计如表3所示，可见检校中误差达到1um精度水平。
表3动态检校精度统计(以1067#影像为例，单位：um)
从计算得到的误差分布来看，本发明方法能有效解算***误差，误差分布基本一致，反应了检校模型和检校技术途径的有效性。
本发明通过采用自检校结合验后补偿方法，有效实现了***误差模型参数和补偿参数的有效解算，检校精度达到0.1像素，检校误差达到1um水平。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。
Claims (7)
1.一种视场分割型航空面阵相机的试验场检校方法，其特征在于，包括：
S1，根据待校验相机信息和校验飞机飞行高度信息，确定试验场信息；
S2，利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息；
S3，利用所述试验场子影像集合信息和所述待校验相机信息，构建得到检校模型；
S4，对所述检校模型进行求解，得到所述待校验相机的校正参数；
所述根据待校验相机信息和校验飞机飞行高度信息，确定试验场信息，包括：
S11，利用试验场信息生成模型，对所述待校验相机信息和校验飞机飞行高度信息进行处理，得到试验场场地信息；所述待校验相机信息，包括相机像幅宽度L、相机焦距f、相机像元尺寸Cs；所述试验场场地信息，包括试验场边长的下界、标志点尺寸值、标志点间隔值和标志点数目下限；
所述试验场信息生成模型，其表达式为：
W＝L×H/f，
GSD＝Cs×H/f，
其中，H为校验飞机飞行高度，W为试验场边长的下界；所述试验场为矩阵区域；GSD为试验场地面分辨率，R为试验场中的标志点尺寸值，S为试验场中的标志点间隔值，Tn为试验场中的标志点数目下限；所述标志点，为试验场中按照预设间隔S进行等间隔布设的标志物；
S12，根据所述试验场场地信息，在试验场中安装标志点；对每个安装后的标志点进行位置测量，得到标志点大地位置坐标；所述标志点大地位置坐标，为所述标志点在大地坐标系下的坐标；
S13，利用所述试验场场地信息和标志点大地位置坐标，构建得到试验场信息；
所述检校模型，包括i0×j0个校验表达式；每个校验表达式，与子影像的一个标志点相对应；对于第i个子影像的第j个标志点，其校验表达式为：
其中，i0为所述待校验相机采集得到的子影像数目，j0为一个子影像所包含的标志点数目，dx0、dy0、df分别为所述待校验相机主点在像元坐标系下的x轴位置坐标改正系数、y轴位置坐标改正系数、所述待校验相机的焦距改正系数，K1，K2，K3分别代表所述待校验相机的径向畸变改正系数，P1，P2分别代表所述待校验相机的切向畸变改正系数，b1，b2分别代表所述待校验相机的焦面探测器件的非正交变形改正系数；A0i，A1i，A2i，A3i，A4i，A5i，B0i，B1i，B2i,B3i，B4i，B5i为第i个子影像的像方二阶验后补偿系数；(xi,j,yi,j)为第i个子影像的第j个标志点的在像元坐标系下的位置坐标，x0、y0、f分别为所述待校验相机主点在像元坐标系下的x轴位置坐标、y轴位置坐标和焦距的初始值，(Xij,Yij,Zij)为所述第i个子影像的第j个标志点的在大地坐标系下的位置坐标，(XS,YS,ZS)为所述待校验相机对所述试验场进行拍摄时刻，所述待校验相机在大地坐标系下的位置坐标，(xi,yi)为第i个子影像的左上角点在所有子影像所构成的像元坐标系下的位置坐标，ri为第i个子影像的相对距离值，rij表示旋转矩阵的第i行、第j列的元素，i＝1,2,3，j＝1,2,3；所述旋转矩阵，根据待校验相机的摄影姿态角决定；在所述检校模型中，dx0、dy0、df、K1、K2、K3、P1、P2、b1、b2、(XS,YS,ZS)、r11、r12、r13、r21、r22、r23、r31、r32、r33、A0i、A1i、A2i、A3i、A4i、A5i、B0i、B1i、B2i、B3i、B4i、B5i为未知量，(xi,j,yi,j)、x0、y0、f、(Xij,Yij,Zij)、(xi,yi)、ri为已知量。
2.如权利要求1所述的视场分割型航空面阵相机的试验场检校方法，其特征在于，所述利用校验飞机搭载所述待校验相机，对所述试验场进行影像采集和分析，得到试验场子影像集合信息，包括：
利用校验飞机搭载所述待校验相机，在校验飞机按照预设飞行高度进行飞行时，利用所述待校验相机对所述试验场进行拍摄，得到若干个子影像信息；每个子影像，对应于所述试验场的部分区域；所述子影像信息，包括子影像位置信息和子影像图像信息；
对每个子影像进行标志点提取，得到标志点图像位置坐标；所述标志点图像位置坐标，为所述标志点在像元坐标系下的位置坐标；
对所有子影像信息和其对应的标志点图像位置坐标进行合并处理，得到试验场子影像集合信息。
3.如权利要求1所述的视场分割型航空面阵相机的试验场检校方法，其特征在于，所述对所述检校方程进行求解，得到所述待校验相机的校正参数，包括：
利用所述检校模型的每个校验表达式，构建得到对应的校验方程；每个校验方程，与子影像的一个标志点相对应；
对于第i个子影像的第j个标志点，其对应的校验方程的表达式为：
其中，fij1和fij2为第i个子影像的第j个标志点的校验方程,v为未知量向量，v＝[A0i、A1i、A2i、A3i、A4i、A5i、B0i、B1i、B2i、B3i、B4i、B5i、r11、r12、r13、r21、r22、r23、r31、r32、r33、XS、YS、ZS、dx0、dy0、df、K1、K2、K3、P1、P2、b1、b2]；
对所有的校验方程进行合并处理，得到校验方程组；所述校验方程组，包括与每个标志点相对应的校验方程；
利用迭代求解模型，对所述校验方程组进行求解，得到所述待校验相机的校正参数。
4.如权利要求3所述的视场分割型航空面阵相机的试验场检校方法，其特征在于，所述利用迭代求解模型，对所述校验方程组进行求解，得到所述待校验相机的校正参数，包括：
S401，对未知量向量进行初始化，得到原始向量；预设变化向量；
S402，利用所述变化向量对所述原始向量进行叠加，得到变化叠加向量；
S403，利用校验方程组的每个校验方程，分别对所述变化叠加向量和原始向量进行计算处理，得到所述校验方程的第一输出值和第二输出值；将所述第一输出值减去第二输出值，得到所述校验方程的差值列向量；
S404，对校验方程组的每个校验方程，在所述原始向量处进行一阶导数计算，得到每个校验方程的一阶导数矩阵；
S405，对所有校验方程的差值列向量，按照列方向进行拼接，得到总差值列向量L；
对所有校验方程的一阶导数矩阵，按照列方向进行拼接，得到总导数矩阵A；
S406，利用总差值列向量L和总导数矩阵A，计算得到变化向量的更新值Δv1，其计算表达式为：
Δv1＝[ATA]-1[-ATL]，
判断Δv1的模值是否小于设定的判别阈值，若小于设定的判别阈值，利用变化向量的更新值与原始向量进行相加，得到相加结果向量，将所述相加结果向量作为所述待校验相机的校正参数；若不小于设定的判别阈值，利用变化向量的更新值，作为变化向量，利用变化叠加向量对原始向量进行更新，执行S402。
5.如权利要求3所述的视场分割型航空面阵相机的试验场检校方法，其特征在于，所述视场分割型航空面阵相机的试验场检校方法，在对所述检校方程进行求解，得到所述待校验相机的校正参数后，还包括：
利用影像位置坐标校正模型，对子影像的像素点的在像元坐标系下的位置坐标进行计算处理，得到子影像的像素点的校正坐标。
6.一种视场分割型航空面阵相机的试验场检校装置，其特征在于，所述装置包括：
存储有可执行程序代码的存储器；
与所述存储器耦合的处理器；
所述处理器调用所述存储器中存储的所述可执行程序代码，执行如权利要求1至5中任一项所述的视场分割型航空面阵相机的试验场检校方法。
7.一种计算机可存储介质，其特征在于，所述计算机存储介质存储有计算机指令，所述计算机指令被调用时，用于执行如权利要求1至5中任一项所述的视场分割型航空面阵相机的试验场检校方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311242030.9A CN117606447B (zh) | 2023-09-25 | 2023-09-25 | 视场分割型航空面阵相机的试验场检校方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311242030.9A CN117606447B (zh) | 2023-09-25 | 2023-09-25 | 视场分割型航空面阵相机的试验场检校方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117606447A CN117606447A (zh) | 2024-02-27 |
| CN117606447B true CN117606447B (zh) | 2024-06-18 |
Family
ID=89954057
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311242030.9A Active CN117606447B (zh) | 2023-09-25 | 2023-09-25 | 视场分割型航空面阵相机的试验场检校方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117606447B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103364012A (zh) * | 2013-07-18 | 2013-10-23 | 中测新图(北京)遥感技术有限责任公司 | 一种带约束条件的多面阵航摄仪平台检校方法 |
| CN110966991A (zh) * | 2019-11-21 | 2020-04-07 | 中国人民解放军国防科技大学 | 一种无控制点下的单幅无人机影像定位方法 |
Family Cites Families (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR100762891B1 (ko) * | 2006-03-23 | 2007-10-04 | 연세대학교 산학협력단 | Ｌｏｓ벡터 조정모델을 이용한 영상의 기하보정 방법 및 그장치 |
| CN101354250B (zh) * | 2008-09-08 | 2010-12-29 | 中国测绘科学研究院 | 具有自校验自稳定功能的组合宽角航空数码相机*** |
| CN101750619B (zh) * | 2010-01-18 | 2012-07-11 | 武汉大学 | 自检校pos直接对地目标定位方法 |
| CN102636159B (zh) * | 2012-04-18 | 2014-05-07 | 中国人民解放军总参谋部测绘研究所 | 多镜头航天线阵相机***在轨几何自检校方法 |
| CN103048872B (zh) * | 2013-01-19 | 2015-07-08 | 吴军 | 全自动数控相机检校*** |
| CN103363960B (zh) * | 2013-07-16 | 2016-03-16 | 中测新图(北京)遥感技术有限责任公司 | 一种基于dlt系数大幅面数字航摄仪影像拼接方法 |
| KR101346323B1 (ko) * | 2013-11-01 | 2014-01-02 | 주식회사 범아엔지니어링 | 지상기준점과 부가매개변수를 이용한 비측량용 디지털 카메라 촬영 시스템의 자체검정 방법 |
| CN106875449B (zh) * | 2017-02-17 | 2019-06-18 | 深圳飞马机器人科技有限公司 | 一种基于飞行数据的无人机非量测型相机检校方法 |
| CN109724625B (zh) * | 2019-01-22 | 2021-05-04 | 中国人民解放军61540部队 | 一种光学复合大面阵测绘相机的像差改正方法 |
| CN115393738A (zh) * | 2021-05-25 | 2022-11-25 | 北京天华航宇科技有限公司 | 一种基于无人机的papi飞行校验方法和*** |
| CN218824686U (zh) * | 2022-08-09 | 2023-04-07 | 中国人民解放军战略支援部队信息工程大学 | 激光雷达和光学相机联合检校靶标器 |
| CN116468621A (zh) * | 2023-03-10 | 2023-07-21 | 国网浙江省电力有限公司湖州供电公司 | 一种一键式数字化航空影像数据处理方法 |
-
2023
- 2023-09-25 CN CN202311242030.9A patent/CN117606447B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103364012A (zh) * | 2013-07-18 | 2013-10-23 | 中测新图(北京)遥感技术有限责任公司 | 一种带约束条件的多面阵航摄仪平台检校方法 |
| CN110966991A (zh) * | 2019-11-21 | 2020-04-07 | 中国人民解放军国防科技大学 | 一种无控制点下的单幅无人机影像定位方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117606447A (zh) | 2024-02-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| KR101965965B1 (ko) | 위성영상과 제공 rpc로부터 제작된 수치표고모델의 자동 오차보정 방법 | |
| US9194954B2 (en) | Method for geo-referencing an imaged area | |
| US8542947B2 (en) | Method for RPC refinement using ground control information | |
| KR20070096370A (ko) | Ｌｏｓ벡터 조정모델을 이용한 영상의 기하보정 방법 및 그장치 | |
| CN111473802A (zh) | 一种基于线阵推扫的光学传感器内方位元素定标方法 | |
| CN107967700B (zh) | 大视场宽工作距双目相机的在轨几何校正及精度验证方法 | |
| CN109724625B (zh) | 一种光学复合大面阵测绘相机的像差改正方法 | |
| Liu et al. | Accurate mapping method for UAV photogrammetry without ground control points in the map projection frame | |
| CN104807477A (zh) | 一种基于靶标控制点的卫星ccd阵列影像几何检校方法 | |
| Yuan | A novel method of systematic error compensation for a position and orientation system | |
| CN108447100A (zh) | 一种机载三线阵ccd相机的偏心矢量和视轴偏心角标定方法 | |
| CN112097794B (zh) | 遥感卫星载荷平台标定方法及*** | |
| Zhang et al. | Auto-calibration of GF-1 WFV images using flat terrain | |
| KR102167847B1 (ko) | 레이저 관측 장비를 이용한 ｍｍｓ의 캘리브레이션을 위한 장치 및 방법 | |
| CN117606447B (zh) | 视场分割型航空面阵相机的试验场检校方法和装置 | |
| CN105571598A (zh) | 一种卫星激光高度计足印相机姿态的测定方法 | |
| CN109696155B (zh) | 光线共面约束的弱交会光学卫星影像联合平差方法及*** | |
| CN111121729A (zh) | 基于平坦地形的光学遥感卫星内检校方法及*** | |
| Rokhmana et al. | Cadastral surveys with non-metric camera using UAV: a feasibility study | |
| CN117237192B (zh) | 视场分割集成面阵相机的全幅影像拼合方法和装置 | |
| CN116091546B (zh) | 光学相机推扫模式下的观测构建方法 | |
| CN116753918B (zh) | 基于空天线阵传感器的地面目标位置估计方法和装置 | |
| Seo et al. | Merge of sub-images from two PAN CCD lines of Kompsat-3 AEISS | |
| CN111623802B (zh) | 一种基于足印影像的星载激光测高仪在轨指向检校方法 | |
| Gao et al. | Application of Rodrigues Matrix in High Accuracy Geo-location for ZY-3 Panchromatic Imagery. |
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 |