具体实施方式

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本申请的一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

本申请中的术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”、“第三”的特征可以明示或者隐含地包括至少一个该特征。本申请的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。本申请实施例中所有方向性指示（诸如上、下、左、右、前、后……）仅用于解释在某一特定姿态（如附图所示）下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、系统、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤或单元。

在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。

图1是本发明实施例的机载面阵相机时空标定方法的流程示意图。需注意的是，若有实质上相同的结果，本发明的方法并不以图1所示的流程顺序为限。如图1所示，该机载面阵相机时空标定方法包括步骤：

步骤S101：在预先构建的控制场中进行覆盖成像，并基于所成影像中提取到的第一SIFT特征进行三维重构，得到SIFT特征控制点，控制场的尺度预先设定且重力加速度方向预先获知。

需要说明的是，控制场是实现相机定位与定向的保障系统。为实现该目的，其应具备以下特征：为保证成像清晰度，成像距离（即相机对控制场成像位置到控制场平均距离）一般应大于30米；控制场景深与成像距离之比大于1/30；控制场范围与主流相机成像范围之比大于2；控制场需布设有控制点标志，该标志既可以被全站仪精确测量位置，也可以在成像影像上准确获取位置；保证主流相机成像时至少有8个控制标志可被成像，并较均匀分布于影像上。在完成控制场布设后，在适当位置架设高精度全站仪，调平后对控制场中所有控制点坐标进行精确测量，显然所获取坐标尺度为1且重力方向已知。

本实施例通过三维重构技术（Structure from motion,简称SFM）实现SIFT特征控制场构建。具体地，该步骤S101具体包括以下步骤：

1、利用高分辨率定焦相机对控制场进行覆盖成像，得到控制场影像。

具体地，在进行覆盖成像时，首先，选择3个成像站点，站点间距大于5米，可对控制场清晰成像；然后，在其中一个成像站点，手持相机利用仿航线法对控制场成像，保证航带内、航带间影像重叠度均大于50%，所有影像可覆盖整个检校场，重复执行该过程，以对三个成像站点均完成成像。

2、将预先获取的控制场中的控制点坐标与控制场影像上控制点对应成像处的成像点坐标进行关联。

3、基于控制点坐标和成像点坐标的关联关系，利用相机初始内参得到相机的初始外方位元素。

具体地，利用控制点坐标与成像点坐标间的关联关系，基于相机初始内参进行直接线性变换（如kneip-P3P法），以获取相机的初始外方位元素（即相机在检校场中的位置及姿态）。

4、分别从每张控制场影像中提取第一SIFT特征。

具体地，从控制场影像中提取第一SIFT特征以及其对应的描述子。

5、将不同控制场影像间的第一SIFT特征进行匹配，且保留匹配后关联特征点数量满足预设数量阈值的控制场影像对。

具体地，利用描述子实现不同控制场影像的第一SIFT特征匹配，并利用初始外方位元素构建的基本矩阵对错误匹配进行过滤，保留匹配后关联特征点数量满足预设数量阈值的控制场影像对，该预设数量阈值优选设置为8个。

6、基于初始外方位元素、相机初始内参和控制场影像对中成功匹配的关联特征点进行三维重构，得到SIFT特征控制点。

具体地，首先进行三角测量，即基于初始外方位元素、相机初始内参对匹配成功的关联特征点估计初始三维坐标。然后，分别进行无相机内参光束法平差和带有相机内参的光束法平差；无相机内参光束法平差，即通过优化方法、利用控制点约束及SIFT特征关联约束对初始外方位元素及初始三维坐标进行优化；带有相机内参的光束法平差，即通过优化方法、利用控制点约束及SIFT特征关联约束对相机内参、初始外方位元素及初始三维坐标进行优化。其次，进行空三质量评估，即检查平差结果中是否存在异常的关联特征点（如三角测量角过小点、投影误差较大点）、刺点误差较大点；若不存在则认为数据合格，否则需调整后重新进行无相机内参光束法平差和带有相机内参的光束法平差。最后，输出SIFT特征控制点，若空三质量评估中无异常情况则导出SIFT三维重构点坐标及其描述子，作为SIFT特征控制点，典型重构结果如图2所示。

需要说明的是，从理论上而言，检校场中一点（控制点或SIFT特征重构点）在像平面上的像素坐标可表达为：

其中，为控制场中点在像空间坐标系下坐标，、为相机坐标系到检校场坐标系的旋转及平移；为相机焦距，为相机主点。由公式（1）~（2）可构建检校场中一点与成像像素点间残差式，由此实现以上步骤中的光束法平差，并获取各SIFT特征控制点在控制场坐标系中坐标值精确值。

步骤S102：基于控制场，结合预设方式进行运动成像，并采集标定数据，标定数据包括多帧面阵影像和IMU数据。

具体地，在预先构建好大尺度的控制场后，基于该控制场，结合预设方式进行运动成像，并采集标定数据，标定数据包括多帧面阵影像和IMU数据。

需要说明的是，现有地面检校方案主要用于视觉—惯性里程计（Visual-inertialodometry，VIO）外方位元素及IMU初始化标定，适配于高频触发相机（10Hz以上），故存在检校期间不同时刻IMU角速度偏差、加速度偏差不变假设，这与常用航空相机（一般触发频率低于1Hz）性能并不相符，也无法使用相关假设。因此，本实施例通过“静—动—静”完成相机与IMU同步标定数据的采集，不仅保证初始时刻及结束时刻速度为零，并使初始时刻及结束时刻IMU角速度偏差亦可通过IMU数据直接读取，从而利用这些约束的增加估计结果的可靠性。具体地，该步骤S102具体包括以下步骤：

1、在控制载荷设备航向朝向天空方向、成像方向朝向控制场的情况下，将载荷设备置于初始静置状态，采集至少第一预设数量的面阵影像且同时采集IMU数据。

具体地，本实施例可根据设备重量等手持设备或置于运动平台上进行数据采集，采集数据时载荷航向朝向天顶方向、成像方向朝向控制场，然后，将载荷设备置于初始静置状态，采集至少第一预设数量的面阵影像且同时采集IMU数据，该第一预设数量优选设置为至少2组。

2、控制载荷设备分别绕自身X、Y、Z轴进行顺时针和逆时针旋转，且在绕每个轴旋转时采集第二预设数量的面阵影像且同时采集IMU数据。

具体地，控制载荷设备分别绕自身X、Y、Z轴进行顺时针和逆时针旋转，旋转角度大于20度，且每个轴采集至少第二预设数量的面阵影像和IMU数据，该第二预设数量优选设置为10组。

3、将载荷设备置于结束静置状态，采集至少第一预设数量的面阵影像且同时采集IMU数据。

本实施例通过“静—动—静”数据采集模式构建前置约束条件（初始时刻及结束时刻IMU角速度偏差已知，初始时刻及结束时刻速度为0）完成时空标定，且视重力加速度方向已知、大小未知，这些更加准确的假设将提高检校精度，使得最终的标定结果更为准确。

步骤S103：分别从每帧面阵影像中提取第二SIFT特征并与SIFT特征控制点进行匹配、且将不同帧面阵影像间的第二SIFT特征进行匹配，并根据匹配结果得到各帧面阵影像成像时相机的位姿信息。

具体地，在得到面阵影像后，对所有的面阵影像进行第二SIFT特征提取，再将该第二SIFT特征与SIFT特征控制点进行匹配，以及将不同帧的面阵影像的第二SIFT特征进行匹配，最后结合匹配结果得到各帧面阵影像成像时相机的位姿信息。

进一步，步骤S103具体包括以下步骤：

1、分别从每帧面阵影像中提取第二SIFT特征。

2、将第二SIFT特征与SIFT特征控制点进行匹配，确认影像外方位元素。

具体地，将第二SIFT特征与SIFT特征控制点进行匹配，通过kneip-P3P+RANSAC（Random Sample Consensus，随机一致性采样）确认影像外方位元素。

3、将不同帧面阵影像间的第二SIFT特征进行匹配，并利用初始外方位元素对匹配的特征点进行过滤，得到匹配特征点。

4、根据面阵影像间匹配特征点的像素位移确认处于初始静置状态及结束静置状态的静置面阵影像。

具体地，根据面阵影像间匹配特征点的像素位移（即各匹配特征点视差）判断处于初始静置状态及结束静置状态的面阵影像，剔除部分静置面阵影像，使初始阶段及结束阶段各保留两张静态面阵影像。

5、基于光束法平差对非静置的面阵影像进行位姿估计，最终得到各帧面阵影像成像时相机的位姿信息。

具体地，影像位姿数据作为时空标定的输入数据，其中第帧影像的位姿为：，其中，、为第帧影像成像时相机在控制场坐标系下的旋转参数及位置参数。

步骤S104：对各帧面阵影像间的IMU数据进行预积分处理，得到IMU的旋转、速度和位移的状态描述。

需要说明的是，原则上，给定初始位姿及速度，IMU位姿可由IMU测量获取角速度及加速度值积分获取，但其测量值不仅包含高斯白噪声和游走偏差，且加速度测量值应减去重力加速度的影响。由此，IMU测量模型可描述为：

其中，及为IMU在全局帧下的角速度及线加速度值，为IMU在世界坐标系下的旋转矩阵，和分别为角速度的高斯白噪声和游走偏差，和分别为线加速度值的高斯白噪声和游走偏差。本实施例采用预积分对相机帧间的IMU数据进行预处理，设某连续相机相邻帧间对应的IMU数据为第至第时刻，则IMU预积分的旋转、速度及位移分别可描述为：

其中，表示第至第时刻的旋转差，表示第至第时刻的速度差，表示第至第时刻的位移差，表示第至第时刻之间的时刻，表示至第时刻的时间差。

若考虑到预积分量中角速度偏差及加速度偏差，则IMU积分获取状态为：

其中：、为第时刻IMU在世界坐标系下的旋转矩阵及平移向量；、为第时刻IMU在世界坐标系下的旋转矩阵及平移向量；及为第及时刻IMU在世界坐标系下的速度；为至第时刻的时间差；为旋转预积分量对的导数；、为速度预积量对、的导数；、为位移预积分量对、的导数，以上导数可由预积分过程获取。

进一步的，如图3所示，由于时钟不同步、传输延迟、传感器响应和操作系统开销，这些延时的存在致使传感器采样时刻与测量时刻间存在未对齐情况。若在时刻同时对相机和IMU数据进行采样，则相机和IMU测量值实际时间戳分别为：

其中，、分别为相机和IMU的实际时间戳，、分别为相机和IMU的延时。由此定义相机及IMU间时间偏移：

由此可获取相机及IMU间位姿存在如下关系：

其中，为采样时刻相机位姿，为采样时刻IMU的位姿，为相机在IMU坐标系下的位姿。显然，所谓时空标定即对时间偏移量、相机与IMU间相对变换矩阵（包含角元素、线元素）进行估计的过程。因此，步骤S104具体包括以下步骤：

1、构建时间偏移构建相机和IMU之间的时空同步状态。

具体地，本发明将时刻IMU在世界坐标系下旋转矩阵及平移向量分别进行线性内插为：

其中，为时刻IMU在体坐标系下角速度，为时刻IMU在世界坐标系下线速度。由此可获取相机与IMU时空同步状态可描述为：

2、基于时空同步状态对各帧面阵影像间的IMU数据进行预积分处理，得到面阵影像的旋转、速度和位移的状态描述。

具体地，面阵影像旋转、速度、位移状态描述分别为：

；（13）

；（14）

；（15）

其中，表示第帧影像成像时相机在控制场坐标系下的旋转预积分量，表示第帧影像成像时相机在控制场坐标系下的旋转预积分量，表示角元素，表示指数函数，表示时刻IMU在体坐标系下角速度，表示时刻IMU在体坐标系下角速度，表示相机与IMU间的时间偏移，表示第至第时刻的旋转差，表示旋转预积分量对的导数，表示角速度偏差，表示第时刻IMU在世界坐标系下速度预积分量，表示第时刻IMU在世界坐标系下速度预积分量，表示重力加速度，表示第至第时刻的时间差，表示第至第时刻的速度差，表示速度预积分量对的导数，表示速度预积分量对的导数，表示加速度偏差，表示第帧影像成像时相机在控制场坐标系下的位移预积分量，表示第帧影像成像时相机在控制场坐标系下的位移预积分量，表示第至第时刻的位移差，表示位移预积分量对的导数，表示位移预积分量对的导数，表示线元素。

步骤S105：基于位姿信息和旋转的状态描述对相机和IMU间的角元素和时间偏差进行标定，且基于位姿信息、速度和位移的状态描述对相机和IMU间的线元素进行标定。

具体地，在得到位姿信息后，利用位姿信息和旋转预计分量即可对相机和IMU间的角元素和时间偏差进行标定，利用位姿信息、速度预积分量和位移预积分量即可对相机和IMU间的线元素进行标定。

进一步的，基于位姿信息和旋转的状态描述对相机和IMU间的角元素和时间偏差进行标定的步骤，具体包括：

1、基于旋转的状态描述构建残差式。

具体地，由于本实施例中待处理连续影像间时间差较大，故假设不同时刻IMU角速度偏差不同。在本实施例中，通过公式（13）构建残差式实现标定。公式（13）中旋转矩阵基于四元数可表达为：

其中，表示相机在IMU坐标系下角元素的四元数表示法；为的共轭四元数，表示的逆变换。

四元数可定义为：

其中，分别为四元数实部、方向虚部、方向虚部、方向虚部，。

则四元数左乘运算和右乘运算可被定义为：

其中为取向量反对称阵，表示3×3大小的单位对角阵。

2、初始化不同时刻的角速度偏差。

具体地，本实施例所采集数据中第帧面阵影像间、第帧面阵影像均处于静态，故该期间IMU角速度真值为0，根据IMU测量模型（见公式（3））可知其测量值为角速度游走偏差与测量噪声之和，取第帧面阵影像间邻近第2帧面阵影像的5帧IMU角速度数据、帧面阵影像间邻近第帧面阵影像的5帧IMU角速度数据，分别取平均值、，作为第2、时刻IMU的角速度偏差，通过线性内插获取第时刻IMU的角速度偏差初值，利用该角速度偏差进行影像帧间IMU数据预积分。需要理解的是，在现有方案中，采用的相机频率高， 因此不同时刻相机的偏差可以认定为相同，但是，本发明中，所采用相机不再局限于高频率相机，还可以采用低频率相机，而低频率相机是不能简单的认定不同时刻的偏差相同的，因此，本发明利用初始及结束时刻的静态数据来对不同时刻的重力加速度偏差进行初始化。

3、基于不同时刻的角速度偏差，利用残差式分别进行角元素标定、时间偏差标定和不同时刻角速度偏差标定。

具体地，针对角元素进行标定，通过公式（16）移项可获取第面阵影像间数据存在等式：

在条件下，利用最小二乘法对的估计：

其中，，为相机在IMU坐标系下角元素的共轭四元数表示法，将其转换为旋转矩阵即完成标定。

针对时间偏差进行标定，通过公式（16）移项可获取第面阵影像间数据存在等式：

若令，则整理公式（21）中四元数虚部可得：

利用最小二乘法对时间偏差进行标定：

；（23）

针对不同时刻角速度偏差进行标定：

通过公式（16）移项可获取第面阵影像间数据存在等式：

若令：

；

则：

利用最小二乘可得第时刻IMU角速度偏差修正量为：

由此修正第时刻角速度偏差，从而对不同时刻角速度偏差进行标定。

3、判断时间偏差是否超过预设时间偏差阈值。

具体地，该预设时间偏差阈值预先设置。

4、若未超过，则确认角元素、时间偏差和不同时刻角速度偏差标定完成。

5、若超过，则循环利用标定得到的时间偏差对IMU数据再次进行预积分处理后进行标定，并将新的标定结果中的时间偏差再次与预设时间偏差阈值进行比较，直至时间偏差低于预设时间偏差阈值时为止。进一步的，基于位姿信息、速度和位移的状态描述对相机和IMU间的线元素进行标定的步骤，包括：

1、基于角元素、时间偏差和不同时刻角速度偏差的完成标定后的标定数据对IMU数据重新进行预积分，得到新的速度和新的位移的状态描述。

具体地，在对线元素进行标定前，利用已标定的角元素、时间偏差和不同时刻角速度偏差重新进行预积分，公式（14）和（15）可变形为：

本实施例基于该式进行线元素标定，除该元素外本实施例中还存在未知数重力加速度、加速度偏差，而在本实施例中控制场重力方向已知，且可描述为，故仅对重力加速度范数进行估计，由此可得第时刻公式（27）~（28）表达为：

其中：

。

2、初始化不同时刻的重力加速度偏差。

具体地，本实施例所采集数据中第帧面阵影像间、第帧面阵影像均处于静态，故该期间IMU加速度真值为0，根据IMU测量模型（见公式（3））可知其测量值为加速度游走偏差、测量噪声及重力加速度之和，取第帧面阵影像间邻近第2帧面阵影像的5帧IMU加速度数据平均值、帧面阵影像间邻近第帧面阵影像的5帧IMU加速度数据平均值后，可计算得第2、时刻IMU的重力加速度偏差：

其中：及为第2、时刻相机在世界坐标系中旋转矩阵的转置，是对角元素标定得到的。通过线性内插获取第时刻IMU的重力加速度偏差初值，利用该重力加速度偏差进行面阵影像帧间IMU数据预积分。需要理解的是，在现有方案中，采用的相机频率高，因此不同时刻相机的偏差可以认定为相同，但是，本发明中，所采用相机不再局限于高频率相机，还可以采用低频率相机，而低频率相机是不能简单的认定不同时刻的偏差相同的，因此，本发明利用初始及结束时刻的静态数据来对不同时刻的重力加速度偏差进行初始化。

3、在控制场重力加速度方向已知的情况，结合重力加速度方向、不同时刻的重力加速度偏差、新的速度和新的位移的状态描述分别对线元素、重力加速度、速度和加速度偏差进行标定。

具体地，针对线元素、重力加速度及速度进行标定：

将中对应列置0后利用最小二乘法对未知参数进行估计：

其中，表示重力加速度的共轭四元数，表示速度的共轭四元数，由于本实施例数据采集的特殊性，其初始速度、结束速度为0，在以上优化中设置相应列为0即可。

针对加速度偏差进行标定：

将中、、对应列置0后，利用最小二乘法对未知参数进行估计：

其中，表示。

3、判断重力加速度是否超过预设加速度阈值。

具体地，该预设加速度阈值预先设置。

4、若未超过，则确认线元素、重力加速度、速度和加速度偏差标定完成。

5、若超过，则循环利用标定得到的加速度偏差对IMU数据再次进行预积分处理后进行标定，并将新的标定结果中的重力加速度再次与预设加速度阈值进行比较，直至重力加速度低于预设加速度阈值时为止。

可以理解的是，在现有方案中，因采用高频相机进行数据采集，因此，对不同时刻IMU的角速度偏差和加速度偏差均可认定为相同。而本发明采用低频相机采集数据，因此，不同时刻IMU的角速度偏差和加速度偏差不能直接认定为相同，导致未知数增加，而为了减少标定过程中的未知数，本发明在进行数据采集时，对数据的采集过程进行约束，使得SIFT特征控制场的尺度、重力方向已知，再分步骤对不同待估计值求解的方法，实现对所有值得标定。

本实施例的机载面阵相机时空标定方法通过在预先设定尺度且重力加速度方向已知的控制场中进行覆盖成像，以得到SIFT特征控制点，再结合控制场采集标定数据，标定数据包括多帧面阵影像和IMU数据，再分别从每帧面阵影像中提取各帧面阵影像成像时相机的位姿信息，再对各帧面阵影像间的IMU数据进行预积分处理，得到面阵影像的旋转、速度和位移的状态描述，最后基于位姿信息和旋转、速度和位移的状态描述对相机和IMU间的角元素、时间偏差、线元素进行标定，基于预先构建的大尺度SIFT特征控制场，其具有重力加速度方向已知、尺度已知的特定，从而可保证航空影像清晰成像质量，在标定过程中不再需要对重力加速度方向进行估计，从而减少标定过程中需要估计的参数量，从而提高标定的精度。

图4是本发明实施例的机载面阵相机时空标定装置的功能模块示意图。如图4所示，该机载面阵相机时空标定装置20包括三维重构模块21、采集模块22、位姿估计模块23、预积分处理模块24和标定模块25。

三维重构模块21，用于在预先构建的控制场中进行覆盖成像，并基于所成影像中提取到的第一SIFT特征进行三维重构，得到SIFT特征控制点，控制场的尺度预先设定且重力加速度方向预先获知；

采集模块22，用于基于控制场，结合预设方式进行运动成像，并采集标定数据，标定数据包括多帧面阵影像和IMU数据；

位姿估计模块23，用于分别从每帧面阵影像中提取第二SIFT特征并与SIFT特征控制点进行匹配、且将不同帧面阵影像间的第二SIFT特征进行匹配，并根据匹配结果得到各帧面阵影像成像时相机的位姿信息；

预积分处理模块24，用于对各帧面阵影像间的IMU数据进行预积分处理，得到面阵影像的旋转、速度和位移的状态描述；

标定模块25，用于基于位姿信息和旋转的状态描述对相机和IMU间的角元素和时间偏差进行标定，且基于位姿信息、速度和位移的状态描述对相机和IMU间的线元素进行标定。

可选地，三维重构模块21执行在预先构建的控制场中进行覆盖成像，并基于所成影像中提取到的第一SIFT特征进行三维重构，得到SIFT特征控制点的操作，具体包括：

利用高分辨率定焦相机对控制场进行覆盖成像，得到控制场影像；

将预先获取的控制场中的控制点坐标与控制场影像上控制点对应成像处的成像点坐标进行关联；

基于控制点坐标和成像点坐标的关联关系，利用相机初始内参得到相机的初始外方位元素；

分别从每张控制场影像中提取第一SIFT特征；

将不同控制场影像间的第一SIFT特征进行匹配，且保留匹配后关联特征点数量满足预设数量阈值的控制场影像对；

基于初始外方位元素、相机初始内参和控制场影像对中成功匹配的关联特征点进行三维重构，得到SIFT特征控制点。

可选地，采集模块22执行基于控制场，结合预设方式进行运动成像，并采集标定数据的操作，具体包括：

在控制载荷设备航向朝向天空方向、成像方向朝向控制场的情况下，将载荷设备置于初始静置状态，采集至少第一预设数量的面阵影像且同时采集IMU数据；

控制载荷设备分别绕自身X、Y、Z轴进行顺时针和逆时针旋转，且在绕每个轴旋转时采集第二预设数量的面阵影像且同时采集IMU数据；

将载荷设备置于结束静置状态，采集至少第一预设数量的面阵影像且同时采集IMU数据。

可选地，位姿估计模块23执行分别从每帧面阵影像中提取第二SIFT特征并与SIFT特征控制点进行匹配、且将不同帧面阵影像间的第二SIFT特征进行匹配，并根据匹配结果得到各帧面阵影像成像时相机的位姿信息的操作，具体包括：

分别从每帧面阵影像中提取第二SIFT特征；

将第二SIFT特征与SIFT特征控制点进行匹配，确认影像外方位元素；

将不同帧面阵影像间的第二SIFT特征进行匹配，并利用初始外方位元素对匹配的特征点进行过滤，得到匹配特征点；

根据面阵影像间匹配特征点的像素位移确认处于初始静置状态及结束静置状态的静置面阵影像；

基于光束法平差对非静置的面阵影像进行位姿估计，最终得到各帧面阵影像成像时相机的位姿信息。

可选地，预积分处理模块24执行对各帧面阵影像间的IMU数据进行预积分处理，得到面阵影像的旋转、速度和位移的状态描述的操作，具体包括；

构建时间偏移构建相机和IMU之间的时空同步状态；

基于时空同步状态对各帧面阵影像间的IMU数据进行预积分处理，得到面阵影像的旋转、速度和位移的状态描述。

可选地，标定模块25执行基于位姿信息和旋转的状态描述对相机和IMU间的角元素和时间偏差进行标定的操作，具体包括：

基于旋转的状态描述构建残差式；

初始化不同时刻的角速度偏差；

基于不同时刻的角速度偏差，利用残差式分别进行角元素标定、时间偏差标定和不同时刻角速度偏差标定；

判断时间偏差是否超过预设时间偏差阈值；

若未超过，则确认角元素、时间偏差和不同时刻角速度偏差标定完成；

若超过，则循环利用标定得到的时间偏差对IMU数据再次进行预积分处理后进行标定，并将新的标定结果中的时间偏差再次与预设时间偏差阈值进行比较，直至时间偏差低于预设时间偏差阈值时为止。

可选地，标定模块25执行基于位姿信息、速度和位移的状态描述对相机和IMU间的线元素进行标定的操作，具体包括：

基于角元素、时间偏差和不同时刻角速度偏差的完成标定后的标定数据对IMU数据重新进行预积分，得到新的速度和新的位移的状态描述；

初始化不同时刻的重力加速度偏差；

在控制场重力加速度方向已知的情况，结合重力加速度方向、不同时刻的重力加速度偏差、新的速度和新的位移的状态描述分别对线元素、重力加速度、速度和加速度偏差进行标定；

判断重力加速度是否超过预设加速度阈值；

若未超过，则确认线元素、重力加速度、速度和加速度偏差标定完成；

若超过，则循环利用标定得到的加速度偏差对IMU数据再次进行预积分处理后进行标定，并将新的标定结果中的重力加速度再次与预设加速度阈值进行比较，直至重力加速度低于预设加速度阈值时为止。

关于上述实施例机载面阵相机时空标定装置中各模块实现技术方案的其他细节，可参见上述实施例中的机载面阵相机时空标定方法中的描述，此处不再赘述。

需要说明的是，本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。对于装置类实施例而言，由于其与方法实施例基本相似，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

请参阅图4，图4为本发明实施例的计算机设备的结构示意图。如图4所示，该计算机设备30包括处理器31及和处理器31耦接的存储器32，存储器32中存储有程序指令，程序指令被处理器31执行时，使得处理器31执行上述任一实施例所述的机载面阵相机时空标定方法步骤。

其中，处理器31还可以称为资源（Central Processing Unit，中央处理单元）。处理器31可能是一种集成电路芯片，具有信号的处理能力。处理器31还可以是通用处理器、数字信号处理器（DSP）、专用集成电路（ASIC）、现场可编程门阵列（FPGA）或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。

参阅图5，图5为本发明实施例的存储介质的结构示意图。本发明实施例的存储介质存储有能够实现上述机载面阵相机时空标定方法的程序指令41，其中，该程序指令41可以以软件产品的形式存储在上述存储介质中，包括若干指令用以使得一台计算机设备（可以是个人计算机，服务器，或者网络设备等）或处理器（processor）执行本申请各个实施方式所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器（ROM，Read-Only Memory）、随机存取存储器（RAM，Random Access Memory）、磁碟或者光盘等各种可以存储程序代码的介质，或者是计算机、服务器、手机、平板等计算机设备设备。

在本申请所提供的几个实施例中，应该理解到，所揭露的计算机设备，装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。以上仅为本申请的实施方式，并非因此限制本申请的专利范围，凡是利用本申请说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本申请的专利保护范围内。

The specific embodiment

The technical solution in the embodiment of the present application will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present application, and it is obvious that the described embodiments are only part of the embodiments of the present application, not all embodiments. Based on the embodiments in the present application, all other embodiments obtained by a person skilled in the art without creative work fall within the scope of protection of the present application.

The terms "first", "second", "third" in the present application are for descriptive purposes only and should not be construed as indicating or implying relative importance or implying the number of technical features indicated. Thus, defining the "first", "second", and "third" features may explicitly or implicitly include at least one of these features. In the description of this application, "plurality" means at least two, such as two, three, etc., unless otherwise expressly and specifically qualified. All directional indications (such as up, down, left, right, front, rear...... in the embodiment of the present application ) is only used to explain the relative position relationship between the parts under a specific attitude (as shown in the attached drawing), the movement, etc., and if the specific attitude changes, the directional indication will change accordingly. In addition, the terms "including" and "having" and any variation of them, are intended to cover non-exclusive inclusions. For example, a process, method, system, product, or apparatus that contains a series of steps or elements is not limited to a listed procedure or unit, but may optionally include a step or unit that is not listed, or may optionally include other steps or elements that are inherent to those processes, methods, products, or equipment.

A reference to "embodiment" in this document means that a particular feature, structure or characteristic described in conjunction with an embodiment may be included in at least one embodiment of the present application. The occurrence of the phrase in various places in the description does not necessarily refer to the same embodiment, nor is it an independent or alternative embodiment that is mutually exclusive with other embodiments. It is understood explicitly and implicitly by those skilled in the art that the embodiments described herein may be combined with other embodiments.

Fig. 1 is a schematic diagram of the flow of the airborne area array camera spatiotemporal calibration method of the embodiment of the present invention. It should be noted that if there is substantially the same result, the method of the present invention is not limited to the sequence of processes shown in Figure 1. As shown in Figure 1, the spatiotemporal calibration method of the airborne area scan camera comprises steps:

Step S101: carrying out coverage imaging in a pre-constructed control field, and performing three-dimensional reconstruction based on the first SIFT feature extracted from the finished image, obtaining a SIFT feature control point, the scale of the control field is preset and the direction of gravitational acceleration is known in advance.

It should be noted that the control field is a guarantee system for camera positioning and orientation. In order to achieve this purpose, it should have the following characteristics: in order to ensure the imaging clarity, the imaging distance (i.e., the average distance from the camera's imaging position to the control field) should generally be greater than 30 meters, and the ratio of the depth of the control scene to the imaging distance should be greater than 1/ 30. The ratio of the control field range to the imaging range of the mainstream camera is greater than 2; the control field needs to be equipped with a control point mark, which can be accurately measured by the total station, and can also accurately obtain the position on the imaging image; ensure that at least 8 control signs can be imaged when the mainstream camera imaging, and be evenly distributed on the image. After the control field is laid out, a high-precision total station is set up at the appropriate position, and the coordinates of all control points in the control field are accurately measured after leveling, and the obtained coordinate scale is obviously 1 and the gravity direction is known.

In this embodiment, the SIFT feature control field construction is realized through a three-dimensional reconstruction technique (Structure from motion, SFM). Specifically, the step S101 specifically comprises the following steps:

1. The high-resolution fixed-focus camera is used to cover the control field to obtain the control field image.

Specifically, when carrying out coverage imaging, firstly, three imaging stations are selected, and the distance between the stations is greater than 5 meters, so that the control field can be clearly imaged; then, at one of the imaging stations, the handheld camera uses the imitation route method to image the control field, and the overlap degree of images in and between the flight bands is guaranteed to be greater than 50%, and all images can cover the entire inspection field, and the process is repeated to complete the imaging of the three imaging stations.

2. The coordinates of the control point in the control field obtained in advance are associated with the coordinates of the imaging point corresponding to the control point on the control field image.

3. Based on the correlation between the coordinates of the control point and the coordinates of the imaging point, the initial external orientation element of the camera is obtained by using the initial internal parameters of the camera.

Specifically, the correlation between the coordinates of the control point and the coordinates of the imaging point is used to perform a direct linear transformation (such as the kneip-P3P method) based on the initial internal parameters of the camera to obtain the initial external orientation element of the camera (i.e., the position and attitude of the camera in the calibration field).

4. Extract the first SIFT feature from each control field image separately.

Specifically, the first SIFT feature and its corresponding descriptor are extracted from the control field image.

5. Match the first SIFT features between different control field images, and retain the control field image pairs whose number of associated feature points meets the preset number threshold after matching.

Specifically, the descriptor is used to realize the first SIFT feature matching of different control field images, and the basic matrix constructed by the initial external orientation element is used to filter the error matching, and the number of associated feature points after matching meets the preset number threshold of the control field image pairs is retained, and the preset number threshold is preferably set to 8.

6. The SIFT feature control points were obtained by 3D reconstruction based on the associated feature points that were successfully matched in the initial external azimuth element, the initial internal parameters of the camera and the control field image pair.

Specifically, triangulation is carried out firstly, that is, the initial three-dimensional coordinates are estimated based on the associated feature points that are successfully matched by the initial external azimuth element and the initial internal reference pair of the camera. Then, the beam adjustment without camera internal control and the beam adjustment with camera internal parameters were carried out, respectively, the initial external azimuth element and the initial 3D coordinate were optimized by the optimization method, the control point constraint and the SIFT feature association constraint, and the beam adjustment with camera internal control was optimized by the optimization method, the control point constraint and the SIFT feature association constraint. Secondly, the aerotriangulation quality assessment is carried out, that is, whether there are abnormal associated feature points (such as the point where the triangulation angle is too small, the point with large projection error) and the point error with a large puncture point in the adjustment results; if not, the data is considered to be qualified, otherwise it is necessary to adjust and re-perform the beam adjustment without camera internal control and the beam adjustment with camera internal control. Finally, the SIFT feature control points are output, and if there is no abnormality in the aerotriangulation quality assessment, the SIFT 3D reconstruction point coordinates and their descriptors are exported as SIFT feature control points, and the typical reconstruction results are shown in Figure 2.

It should be noted that, theoretically, the pixel coordinates of a point (control point or SIFT feature reconstruction point) in the calibration field on the image plane can be expressed as:

Among them, it is the coordinate of the control field midpoint in the image space coordinate system, the rotation and translation from the camera coordinate system to the inspection field coordinate system, and the camera focal length and the camera principal point. From the equation (1)~(2), the residual formula between a point in the calibration field and the imaging pixel can be constructed, so as to realize the beam adjustment in the above steps, and obtain the exact value of the coordinate value of each SIFT feature control point in the control field coordinate system.

Step S102: based on the control field, combined with the preset mode for motion imaging, and collecting calibration data, the calibration data comprises multi-frame area array images and IMU data.

Specifically, after the large-scale control field is constructed in advance, motion imaging is carried out based on the control field and combined with the preset mode, and the calibration data includes multi-frame area array images and IMU data.

It should be noted that the existing ground calibration scheme is mainly used for the external azimuth element of Visual-inertialodometry (VIO) and IMU initialization calibration, and is suitable for high-frequency triggering cameras (above 10Hz), so there is the assumption that the IMU angular velocity deviation and acceleration deviation remain unchanged at different times during the calibration period, which is not consistent with the performance of commonly used aerial cameras (generally the trigger frequency is less than 1Hz), and the relevant assumptions cannot be used. Therefore, the present embodiment completes the acquisition of synchronous calibration data between the camera and the IMU through "static-dynamic-static", which not only ensures that the velocity at the initial moment and the end moment is zero, but also makes the IMU angular velocity deviation at the initial moment and the end moment can be directly read through the IMU data, so as to increase the reliability of the estimation results by using these constraints. Specifically, the step S102 specifically comprises the following steps:

1. Under the condition that the heading of the payload equipment is facing the sky direction and the imaging direction is facing the control field, the payload equipment is placed in the initial static state, and at least the first preset number of area array images are collected and the IMU data is collected at the same time.

Specifically, the present embodiment may carry out data acquisition according to a handheld device such as the weight of the equipment or placed on a motion platform, and the load heading direction is toward the zenith direction and the imaging direction is towards the control field when the data is collected, and then the load equipment is placed in an initial static state, and at least the first preset number of area array images are collected and IMU data is collected simultaneously, and the first preset number is preferably set to at least 2 groups.

2. The control load equipment rotates clockwise and counterclockwise around its own X, Y, and Z axes respectively, and collects a second preset number of area array images and collects IMU data at the same time when rotating around each axis.

Specifically, the control load equipment rotates clockwise and counterclockwise around its own X, Y, and Z axes respectively, with a rotation angle greater than 20 degrees, and each axis collects at least a second preset number of area array images and IMU data, and the second preset number is preferably set to 10 groups.

3. Put the load equipment in the end static state, collect at least the first preset number of area array images, and collect IMU data at the same time.

In this embodiment, the spatiotemporal calibration is completed by constructing the pre-constraints (the angular velocity deviation of the IMU at the initial moment and the end time is known, and the velocity at the initial time and the end time is 0) through the "static-dynamic-static" data acquisition mode, and the direction and magnitude of the apparent gravitational acceleration are known, and these more accurate assumptions will improve the calibration accuracy and make the final calibration results more accurate.

step S103: extracting a second SIFT feature from each frame area array image and matching it with the SIFT feature control point, matching the second SIFT feature between different frame area array images, and obtaining the pose information of the camera when each frame area array image is imaged according to the matching result.

Specifically, after obtaining the area array images, the second SIFT feature is extracted from all the area array images, and then the second SIFT feature is matched with the SIFT feature control point, and the second SIFT feature of the area array images of different frames is matched, and finally the pose information of the camera during the imaging of the area array image of each frame is obtained by combining the matching results.

Further, step S103 specifically comprises the following steps:

1. Extract the second SIFT feature from each frame of area array image.

2. Match the second SIFT feature with the SIFT feature control point to confirm the orientation element outside the image.

Specifically, the second SIFT feature is matched with the SIFT feature control point, and the external orientation elements of the image are confirmed by kneip-P3P+RANSAC (Random Sample Consensus).

3. The second SIFT features between different frame area array images are matched, and the matching feature points are filtered by using the initial external orientation element to obtain the matching feature points.

4. According to the pixel displacement of the matching feature points between the area array images, confirm the static area array images in the initial static state and the end static state.

Specifically, according to the pixel displacement of the matching feature points between the area array images (i.e., the parallax of each matching feature point), the area array images in the initial static state and the end static state are judged, and part of the static area array images are eliminated, so that two static area array images are retained in the initial stage and the end stage.

5. Based on the beam adjustment, the pose estimation of the non-static area array image is carried out, and finally the pose information of the camera during the imaging of each frame area array image is obtained.

Specifically, the image pose data is used as the input data for spatiotemporal calibration, wherein the pose of the first frame image is: , where , is the rotation parameter and position parameter of the camera in the control field coordinate system when the first frame image is imaging.

Step S104: perform pre-integration processing on the IMU data between the area array images of each frame to obtain the state description of the rotation, velocity and displacement of the IMU.

It should be noted that, in principle, given the initial pose and velocity, the IMU pose can be obtained by the integral of the angular velocity and acceleration value obtained by the IMU measurement, but the measured value not only includes Gaussian white noise and wandering deviation, and the acceleration measurement value should subtract the influence of gravitational acceleration. From this, the IMU measurement model can be described as:

where and are the angular velocity and linear acceleration values of the IMU in the global frame, the rotation matrix of the IMU in the world coordinate system, and the Gaussian white noise and wander deviation of the angular velocity, and the Gaussian white noise and wander deviation of the linear acceleration values, respectively. In this embodiment, the IMU data between camera frames is preprocessed by pre-integral, and the corresponding IMU data between adjacent frames of a continuous camera is the first to the first moment, then the rotation, velocity and displacement of the IMU pre-integration can be described as:

Wherein, the rotation difference from the first to the first time, the velocity difference from the first to the first time, the displacement difference from the first to the first time, the time between the first and the first time, and the time difference between the first and first moments.

If the angular velocity deviation and acceleration deviation in the pre-integration are taken into account, the IMU integration acquisition status is:

Wherein: is the rotation matrix and translation vector of the IMU in the world coordinate system at the first time, the rotation matrix and translation vector of the IMU in the world coordinate system at the first time, and the velocity of the IMU in the world coordinate system at the first time, and the time difference at the first time;is the derivative of the rotation preintegration pair, is the velocity preintegration pair, is the derivative of the displacement preintegration pair, and the above derivatives can be obtained from the preintegration process.

Further, as shown in Figure 3, there is a misalignment between the sensor sampling time and the measurement time due to clock out-of-sync, propagation delays, sensor responses, and operating system overhead. If the camera and IMU data are sampled at the same time, the actual timestamps of the camera and IMU measurements are:

where , , is the actual timestamp of the camera and IMU, and is the time-lapse of the camera and IMU, respectively. This defines the time offset between the camera and the IMU:

From this, it can be obtained that there is the following relationship between the pose between the camera and the IMU:

wherein, is the camera pose at the sampling time, the IMU at the sampling time, and the camera pose in the IMU coordinate system. Obviously, the so-called spatiotemporal calibration is the process of estimating the time offset and the relative transformation matrix between the camera and the IMU (including angular elements and line elements). Therefore, step S104 specifically comprises the following steps:

1. Build the time offset to build the spatio-temporal synchronization state between the camera and the IMU.

Specifically, the present invention linearly interpolates the rotation matrix and the translation vector of the IMU at the moment in the world coordinate system as follows:

wherein, is the angular velocity of the IMU in the body coordinate system at the time, and the linear velocity of the IMU in the world coordinate system at the time. From this, the spatiotemporal synchronization status of the camera and the IMU can be described as:

2. Based on the spatio-temporal synchronization state, the IMU data between each frame of area array images are pre-integrated to obtain the state description of the rotation, velocity and displacement of area array images.

Specifically, the rotation, velocity, and displacement states of the area array image are described as follows:

；（13）

；（14）

；（15）

Wherein, it represents the rotation preintegration of the camera in the control field coordinate system when the first frame image is imaging, represents the rotation preintegration of the camera in the control field coordinate system when the first frame image is imaging, represents the angular element, represents the exponential function, and represents the angular velocity of the IMU in the body coordinate system at the momentThe angular velocity of the IMU in the body coordinate system at the moment represents the time offset between the camera and the IMU, the rotation difference from the first to the first moment, the derivative of the rotation preintegration pair, the angular velocity deviation, and the velocity preintegration of the IMU in the world coordinate system at the first timeThe velocity preintegration of the IMU at the moment in the world coordinate system represents the acceleration of gravity, the time difference between the first and first moments, the velocity difference between the first and first moments, the derivative of the velocity preintegration pair, and the derivative of the velocity preintegration pairRepresents the acceleration deviation, represents the displacement preintegration of the camera in the control field coordinate system when the first frame image is imaging, represents the displacement preintegration of the camera in the control field coordinate system during the first frame image imaging, represents the displacement difference from the first to the first time, represents the derivative of the displacement preintegration pair, and represents the displacement preintegration pair, which represents the line element.

Step S105: calibrating the angular elements and time deviations between the camera and the IMU based on the pose information and the state description of rotation, and calibrating the line elements between the camera and the IMU based on the state descriptions of the pose information, velocity and displacement.

Specifically, after obtaining the pose information, the angular element and time deviation between the camera and the IMU can be calibrated by using the pose information and the rotation prediction component, and the line element between the camera and the IMU can be calibrated by using the pose information, velocity preintegration and displacement preintegration.

Further, the steps to calibrate the angular element and time deviation between the camera and the IMU based on the pose information and the state description of the rotation include:

1. Construct the residual formula based on the state description of rotation.

Specifically, because the time difference between the continuous images to be processed in this embodiment is large, it is assumed that the angular velocity deviation of the IMU at different times is different. In the present embodiment, the residuals are constructed through equation (13) to achieve calibration. In Eq. (13), the rotation matrix can be expressed based on quaternions as:

where is the quaternion representation of the lower corner element of the camera in the IMU coordinate system, and the conjugate quaternion represents the inverse transformation.

Quaternions can be defined as:

Among them, they are the real part of the quaternion, the imaginary part of the direction, the imaginary part of the direction, and the imaginary part of the direction.

Then the quaternion left-multiplication and right-multiplication operations can be defined as:

Among them is the vector-oriented antisymmetry matrix, which represents the diagonal matrix of units of 3×3 size.

2. Initialize the angular velocity deviation at different moments.

Specifically, in the data collected in this embodiment, the first frame area array image and the first frame area array image are static, so the true value of the IMU angular velocity during this period is 0, and the measured value can be known to be the sum of the angular velocity wandering deviation and the measured noise according to the IMU measurement model (see equation (3)), and the 5 frames of IMU angular velocity data adjacent to the second frame area array image between the first frame area array images are takenThe angular velocity data of the 5 frames of IMU adjacent to the first frame of the area array images between the frame area array images are taken as the average value and the angular velocity deviation of the IMU at the second and second moments, and the first is obtained by linear interpolationThe initial value of the angular velocity deviation of the IMU at the moment is used to pre-integrate the IMU data between image frames. It needs to be understood that, in the existing scheme, the camera frequency adopted is high, so the deviation of the camera at different moments can be identified as the same, but, in the present invention, the camera adopted is no longer limited to the high-frequency camera, and the low-frequency camera can also be adopted, and the low-frequency camera cannot simply identify the deviation of different moments is the same, therefore, the present invention utilizes the static data of the initial and end moments to initialize the gravitational acceleration deviation at different moments.

3. Based on the angular velocity deviation at different moments, the residual formula is used to calibrate the angular element, the time deviation and the angular velocity deviation at different times.

Specifically, for the calibration of the angular elements, the equation of data existence between the first surface array images can be obtained by shifting the term in equation (16):

Under the conditions, the estimation using the least squares pairs:

, is the conjugate quaternion representation of the lower corner elements of the camera in the IMU coordinate system, which is converted into a rotation matrix to complete the calibration.

For the calibration of the time deviation, the equation between the data between the first area array images can be obtained by shifting the term of equation (16):

If it is ordered, then the imaginary part of the quaternion in equation (21) can be arranged:

The time deviation is calibrated using the least squares method:

；（23）

Calibration: Calibration for angular velocity deviation at different moments:

The equation of data existence between the first area array images can be obtained by shifting the term in equation (16):

If ordered:

；

Rule:

Using the least squares, the correction of the angular velocity deviation of the IMU at the first moment is as follows:

In this way, the angular velocity deviation at the first moment is corrected, so that the angular velocity deviation at different moments can be calibrated.

3. Judge whether the time deviation exceeds the preset time deviation threshold.

Specifically, this preset time deviation threshold is preset.

4. If it is not exceeded, confirm that the angular element, time deviation and angular velocity deviation at different times are calibrated.

5. If it is exceeded, the time deviation obtained by the calibration is recycled to pre-integrate the IMU data again and then calibrated again, and the time deviation in the new calibration result is compared with the preset time deviation threshold again until the time deviation is lower than the preset time deviation threshold. Further, the steps to calibrate the line elements between the camera and the IMU based on the posture information, velocity, and displacement state descriptions include:

1. The IMU data is pre-integrated based on the calibration data after the completion of the calibration of angular elements, time deviations and angular velocity deviations at different moments to obtain new velocity and new displacement state descriptions.

Specifically, before calibrating the line elements, the pre-integration is re-integrated by using the calibrated angular elements, time deviations, and angular velocity deviations at different times, and equations (14) and (15) can be deformed as:

The present embodiment carries out the calibration of the line element based on this formula, in addition to the element, there are unknown gravitational acceleration, acceleration deviation in the present embodiment, and in the present embodiment, the gravitational direction of the control field is known, and can be described as, so only the gravitational acceleration norm is estimated, and the formula (27)~(28) of the first moment can be expressed as:

Thereinto:

。

2. Initialize the gravitational acceleration deviation at different moments.

Specifically, in the data collected in this embodiment, the first frame area array image and the first frame area array image are static, so the true value of IMU acceleration during this period is 0, and the measured value is the sum of acceleration drift deviation, measurement noise and gravitational acceleration according to the IMU measurement model (see equation (3)), and the average value of the five IMU acceleration data adjacent to the second frame area array image between the first area array images is takenAfter the average value of the 5 frames of IMU acceleration data adjacent to the first frame area array image between the frame area array images, the gravitational acceleration deviation of the IMU at the second and time can be calculated:

Wherein: and is the transposition of the rotation matrix of the camera in the world coordinate system at the second time, which is obtained by calibrating the diagonal element. The initial value of the gravitational acceleration deviation of the IMU at the first moment is obtained by linear interpolation, and the gravitational acceleration deviation is used to preintegrate the IMU data between the area array image frames. It needs to be understood that, in the existing scheme, the frequency of the camera adopted is high, so the deviation of the camera at different moments can be identified as the same, but, in the present invention, the camera adopted is no longer limited to the high-frequency camera, and the low-frequency camera can also be used, and the low-frequency camera cannot simply identify the deviation of different moments is the same, therefore, the present invention utilizes the static data of the initial and end moments to initialize the gravitational acceleration deviation at different moments.

3. When the direction of gravitational acceleration in the control field is known, the line element, gravitational acceleration, velocity and acceleration deviation are calibrated respectively in combination with the direction of gravitational acceleration, the gravitational acceleration deviation at different moments, the new velocity and the state description of the new displacement.

Specifically, calibration is performed for line elements, gravitational acceleration, and velocity:

After setting the corresponding column to 0, the least squares method is used to estimate the unknown parameters:

Wherein, the conjugate quaternion representing gravitational acceleration and the conjugate quaternion representing velocity are 0 due to the particularity of the data acquisition in this embodiment, and the corresponding column can be set to 0 in the above optimization.

Calibration for acceleration deviation:

After setting the middle, and corresponding columns to 0, the least squares method is used to estimate the unknown parameters:

Among them, it is expressed.

3. Determine whether the gravitational acceleration exceeds the preset acceleration threshold.

Specifically, this preset acceleration threshold is preset.

4. If it is not exceeded, confirm that the calibration of line elements, gravitational acceleration, velocity and acceleration deviation is completed.

5. If it is exceeded, the acceleration deviation obtained by the calibration is recycled to pre-integrate the IMU data again and then calibrate, and the gravitational acceleration in the new calibration result is compared with the preset acceleration threshold again until the gravitational acceleration is lower than the preset acceleration threshold.

It is understandable that in the existing scheme, the angular velocity deviation and acceleration deviation of the IMU at different times can be considered the same because the high-frequency camera is used for data acquisition. And the present invention adopts a low-frequency camera to collect data, therefore, the angular velocity deviation and acceleration deviation of IMU at different moments can not be directly identified as the same, resulting in the increase of unknowns, and in order to reduce the unknowns in the calibration process, the present invention restricts the data acquisition process when carrying out data acquisition, so that the scale and gravity direction of the SIFT characteristic control field are known, and then the method of solving different values to be estimated in steps realizes the calibration of all worthiness.

The spatiotemporal calibration method of the airborne area array camera in the present embodiment obtains the SIFT feature control point by carrying out coverage imaging in the control field with a preset scale and the direction of gravitational acceleration is known, and then the calibration data is collected in combination with the control field, and the calibration data comprises multiple frames of area array images and IMU data, and then the pose information of the camera when each frame area array image is imaged is extracted from each frame area array image respectively, and then the IMU data between each frame area array image are pre-integrated to obtain the rotation, The state description of velocity and displacement, and finally the angular elements, time deviation and line elements between the camera and the IMU are calibrated based on the posture information and the state description of rotation, velocity and displacement, and based on the pre-constructed large-scale SIFT feature control field, which has the specificity of knowing the direction of gravitational acceleration and the known scale, so as to ensure the clear imaging quality of aerial images, and no longer need to estimate the direction of gravitational acceleration in the calibration process, so as to reduce the number of parameters that need to be estimated in the calibration process, so as to improve the accuracy of calibration.

Fig. 4 is a schematic diagram of the functional module of the airborne area array camera spatiotemporal calibration device of the embodiment of the present invention. As shown in Fig. 4, the airborne area array camera spatiotemporal calibration device 20 comprises a three-dimensional reconstruction module 21, an acquisition module 22, a pose estimation module 23, a preintegration processing module 24 and a calibration module 25.

the three-dimensional reconstruction module 21 is used for carrying out coverage imaging in the preconstructed control field, and carrying out three-dimensional reconstruction based on the first SIFT feature extracted from the finished image, obtaining the SIFT feature control point, the scale of the control field is preset and the gravitational acceleration direction is known in advance;

the acquisition module 22 is used for carrying out motion imaging based on the control field, combining with a preset mode, and collecting calibration data, and the calibration data comprises a multi-frame area array image and IMU data;

the pose estimation module 23 is used for extracting a second SIFT feature from each frame area array image and matching it with the SIFT feature control point, and matching the second SIFT feature between different frame area array images, and obtaining the pose information of the camera when each frame area array image is imaged according to the matching result;

the pre-integration processing module 24 is used for pre-integrating the IMU data between each frame area array image to obtain the state description of rotation, velocity and displacement of area array image;

The calibration module 25 is used for calibrating the angular element and the time deviation between the camera and the IMU based on the posture information and the state description of the rotation, and calibrating the line element between the camera and the IMU based on the state description of the pose information, velocity and displacement.

Optionally, the 3D reconstruction module 21 performs an operation of carrying out coverage imaging in a preconstructed control field and performing 3D reconstruction based on the first SIFT feature extracted from the finished image to obtain a SIFT feature control point, comprising:

The high-resolution fixed-focus camera was used to cover the control field to obtain the control field image.

The coordinates of the control points in the pre-obtained control field are associated with the coordinates of the imaging points corresponding to the control points on the control field image.

Based on the correlation between the coordinates of the control point and the coordinates of the imaging point, the initial external azimuth element of the camera is obtained by using the initial internal parameters of the camera.

The first SIFT feature was extracted from each control field image respectively.

The first SIFT features between different control field images are matched, and the number of associated feature points after matching is retained.

The SIFT feature control points were obtained by 3D reconstruction based on the initial external orientation element, the initial internal control of the camera and the successfully matched associated feature points in the control field image pair.

Optionally, the acquisition module 22 performs operations based on the control field, combined with a preset mode for motion imaging, and collects calibration data, comprising:

Under the condition that the heading direction of the payload device is facing the sky direction and the imaging direction is facing the control field, the payload device is placed in the initial static state, and at least the first preset number of area array images are collected and the IMU data is collected at the same time;

the control load equipment rotates clockwise and counterclockwise around its own X, Y, and Z axes respectively, and collects a second preset number of area array images and collects IMU data simultaneously when rotating around each axis;

The payload device is placed in the end static state, and at least the first preset number of area array images are acquired, and the IMU data is collected at the same time.

Optionally, the pose estimation module 23 performs operations that extract the second SIFT feature from each frame area array image and match it with the SIFT feature control point, match the second SIFT feature between different frame area array images, and obtain the pose information of the camera when each frame area array image is imaged according to the matching result, specifically comprising:

The second SIFT feature was extracted from each frame of area array image.

The second SIFT feature is matched with the SIFT feature control point to confirm the orientation element outside the image.

The second SIFT features between different frame area array images are matched, and the matching feature points are filtered by using the initial external orientation element to obtain the matching feature points.

According to the pixel displacement of the matching feature points between the area array images, the static area array images in the initial static state and the end static state are confirmed.

Based on the beam adjustment, the pose estimation of the non-static area array image is carried out, and finally the pose information of the camera during the imaging of each frame of area array image is obtained.

Optionally, the pre-integration processing module 24 performs pre-integration processing on the IMU data between each frame of area array images to obtain the state description of rotation, velocity and displacement of area array images, specifically comprising;

Construct time offset: construct spatiotemporal synchronization state between camera and IMU;

Based on the spatiotemporal synchronization state, the IMU data between the area array images of each frame are pre-integrated to obtain the state description of the rotation, velocity and displacement of the area array images.

Optionally, the calibration module 25 performs operations for calibrating angular elements and temporal deviations between the camera and the IMU based on the pose information and the state description of the rotation, comprising:

The residual formula is constructed based on the state description of rotation.

Initialize angular velocity deviations at different moments;

Based on the angular velocity deviation at different moments, the residual formula is used to calibrate the angular element, the time deviation and the angular velocity deviation at different times.

Determine whether the time deviation exceeds the preset time deviation threshold;

If it is not exceeded, it is confirmed that the angular element, time deviation and angular velocity deviation at different times are calibrated;

If it is exceeded, the IMU data is pre-integrated again by recycling the time deviation obtained from the calibration and then calibrated again, and the time deviation in the new calibration result is compared with the preset time deviation threshold again until the time deviation is lower than the preset time deviation threshold.

Optionally, the calibration module 25 performs the operation of calibrating the line element between the camera and the IMU based on the state description of the pose information, velocity and displacement, comprising:

Based on the calibration data of angular elements, time deviations and angular velocity deviations at different times, the IMU data are pre-integrated to obtain new velocity and new displacement state descriptions.

Initialize the gravitational acceleration deviation at different moments;

When the direction of gravitational acceleration in the control field is known, the line element, gravitational acceleration, velocity and acceleration deviation are calibrated by combining the direction of gravitational acceleration, the gravitational acceleration deviation at different moments, the new velocity and the state description of the new displacement, respectively.

Determine whether the gravitational acceleration exceeds the preset acceleration threshold;

If it is not exceeded, it is confirmed that the calibration of line element, gravitational acceleration, velocity and acceleration deviation is completed;

If it is exceeded, the acceleration deviation obtained by the calibration is recycled to pre-integrate the IMU data and calibrated again, and the gravitational acceleration in the new calibration result is compared with the preset acceleration threshold again until the gravitational acceleration is lower than the preset acceleration threshold.

For other details of the technical scheme for the implementation of each module in the airborne area scan camera spatio-temporal calibration device of the above embodiment, please refer to the description in the spatiotemporal calibration method of the airborne area scan camera in the above embodiment, and will not be repeated here.

It should be noted that each embodiment in this specification is described in a progressive manner, and each embodiment focuses on the differences from other embodiments, and the same and similar parts between each embodiment can refer to each other. For device embodiments, because they are basically similar to the method embodiments, the description is relatively simple, and the relevant places can be described in part of the method embodiments.

Refer to FIG. 4, FIG. 4 is a schematic diagram of the structure of the computer equipment of the embodiment of the present invention. As shown in Figure 4, the computer device 30 comprises a processor 31 and a memory 32 coupled with the processor 31, a program instruction is stored in the memory 32, and when the program instruction is executed by the processor 31, the processor 31 is made to perform the spatiotemporal calibration method step of the airborne area array camera described in any embodiment.

Wherein, the processor 31 can also be called a resource (Central Processing Unit, central processing unit). The processor 31 may be an integrated circuit chip with the ability to process signals. The processor 31 may also be a general-purpose processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA) or other programmable logic device, a discrete gate or transistor logic device, or a discrete hardware component. A general-purpose processor can be a microprocessor, or the processor can be any regular processor, etc.

Refer to FIG. 5, FIG. 5 is a schematic diagram of the structure of the storage medium of the embodiment of the present invention. The storage medium of the embodiment of the present invention stores a program instruction 41 that can realize the spatiotemporal calibration method of the airborne area array camera, wherein the program instruction 41 can be stored in the storage medium in the form of a software product, and comprises a plurality of instructions to enable a computer device (which may be a personal computer, a server, or a network device, etc.) or a processor (processor) to perform all or part of the steps of the method described in each embodiment of the present application. The aforementioned storage media include: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM), disk disks or optical discs and other media that can store program code, or computers, servers, mobile phones, tablets and other computer equipment and equipment.

In several embodiments provided in the present application, it should be understood that the disclosed computer equipment, apparatus and methods may be realized by other means. For example, the device embodiments described above are only schematic, for example, the division of units is only a logical function division, and there can be other ways of division in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. On the other hand, the coupling or direct coupling or communication connection between each other shown or discussed may be indirect coupling or communication connection through some interface, device or unit, which may be in electrical, mechanical or other form.

In addition, each functional unit in each embodiment of the present invention may be integrated in a processing unit, or each unit may exist separately physically, or two or more units may be integrated in a single unit. The above-mentioned integrated units can be implemented either in the form of hardware or in the form of software function units. The above is only an embodiment of the present application, and does not limit the scope of the patent of the present application, and all equivalent structural or equivalent process transformations made by using the contents of the description and drawings of the present application, or directly or indirectly applied in other related technical fields, are similarly included in the scope of patent protection of the present application.