CN103310448A - 用于das的摄像装置姿态角估计及实时生成合成图的方法 - Google Patents

Abstract

本发明公开了一种用于DAS的摄像装置姿态角估计及实时生成合成图的方法，属于图像处理技术领域。本发明的姿态角估计方法为：每相邻的两个摄像装置拍摄一对图像对{I_i，I_j}，得到M对参考图像对；构建姿态角列向量Ω，其元素依次对应M个摄像装置的姿态角；基于M对参考图像对、M个摄像装置的内参数矩阵，通过数值优化计算参量e(Ω)取得最小时所对应的列向量Ω的取值，得到各摄像装置的工作姿态角；参量e(Ω)为：累加M个参考图像对的各图像对{I_i，I_j}的图像I_j的各像素点与在图像I_i上的映射点的像素值之差的平方和。基于所述姿态角估计方法，本发明还提出了一种实时合成图生成方法。本发明易于在通用产品上实时生成任意指定视场范围的合成图像。

Description

用于DAS的摄像装置姿态角估计及实时生成合成图的方法

技术领域

本发明属于图像处理技术领域，具体涉及一种对分布式孔径系统的相机姿态角的估计方法、实时生成合成图的方法。

背景技术

可见光相机是当前获取图像信息所使用的最常见摄像设备。实时获取周围环境的场景图像用于场景监视安保等领域具有重要的应用需求和实用价值，例如，执行警备任务的车辆需要监控车辆周围场景状况，直升机驾驶员需要了解以飞机为中心的周围空域中任意视角的情况等等。然而由于单个相机视场角有限，所以单个相机只能提供局部视场范围内的场景图像。为了解360度全视场范围内的场景情况，需要采用分布式孔径系统(Distributed ApertureSystem-DAS)，即通过在不同视角安装的多个相机来分别提供图像(称为孔径图像)，相邻相机之间的视场范围互相重叠，从而形成一个360度全视场，进而采用图像拼接方式生成任意视场范围内的合成图。

在DAS中，对于一个相机，采用内参数矩阵K和外参数矩阵R对其进行描述。内参数矩阵的形式为： $K=\left[\begin{array}{ccc}{f}_x& 0& p_x\\ 0& f_y& p_y\\ 0& 0& 1\end{array}\right],$ 其中f_x和f_y分别为该相机在成像平面上水平(x)方向和垂直(y)方向的焦距，[p_x,p_y]^T为成像平面上主点的坐标，通过相机定标的方法可以计算出内参数矩阵K中的参数，外参数矩阵R由相机的姿态角(方位角、俯仰角和绕轴线的旋转角)可通过Rodriguez（罗德里格斯）旋转公式确定，Rodriguez旋转公式具体可参考文献：R.M.Murray,Z.Li,S.S.Sastry著，”A Mathematical Introduction to Robotic Manipulation”,CRC Press出版，1994年。基于DAS中各相机的内参数矩阵K和外参数矩阵R，可以生成任意指定视角的场景的合成图，即等效为任意指定一个内参数矩阵为K_v、外参数矩阵为R_v的虚拟相机，该虚拟相机输出即为合成图。若生成合成图的像素点对应多个真实相机的成像平面中的像素，则可以对这些像素点通过羽化（feathering）的方法，对像素值进行融合，使得生成的合成图在视觉上可以接受。

在DAS的推广应用中，需要一些特点的限制（硬件的支持，或者是复杂的坐标变换计算）方能实现合成图的生成，使得其通用性和普适性差。如公布号为US2004169663(A1)的美国专利申请所公开的用于飞行员头盔的机载DAS系统，需要建立一种几何映射曲面用于投影多成像传感器传出的图像。然而，该几何映射曲面的半径需要飞机上装载的导航系统提供的数据来计算得到，不适合其他应用场合，例如安保监控，该方案的通用性差。

公布号为US2006066730(A1)的美国专利申请所公开的DAS系统中，对生成的虚拟图像的二维像素坐标首先需要转换为方位角、俯仰角坐标，再转换为三维FRD（Forward-Right-Down）坐标；FRD坐标通过两次坐标变换对应到孔径图像的相应像素的FRD坐标，然后再转换为孔径图像的实际二维像素坐标。这一过程涉及较多的变换计算，需要较强的硬件支持，通用性和普适性差。

在相机姿态估计方面，现有方法采用从多个相机图像中分别提取特征，再对特征点进行匹配的方法。如文献“M.Brown,D.G.Lowe,“Automatic panoramic image stitching usinginvariant features,”International Journal of Computer Vision,74(1):59-73,2007.”采用SIFT特征提取和描述方法估计对应每幅图像的相机姿态。然而因为这类方法采用的信息是稀疏的图像特征点，所以对于少纹理区域提取的特征点数量和质量都较差，会导致匹配失败。

发明内容

本发明的发明目的在于：提供一种易于实现、在图像拼接时精确估计相机的姿态的方法。

本发明的用于DAS的摄像装置姿态角估计方法，包括下列步骤：

每相邻的两个摄像装置拍摄一对图像对{I_i，I_j}，得到M对参考图像对，其中i，j为M个摄像装置的标识；

构建姿态角列向量Ω，所述列向量Ω的元素依次对应M个摄像装置的姿态角；

获取M个摄像装置的内参数矩阵；

基于M对参考图像对、M个摄像装置的内参数矩阵，通过数值优化计算参量e(Ω)取得最小时所对应的列向量Ω的取值，得到各摄像装置的工作姿态角；

所述参量e(Ω)为：累加M个图像对中，各图像对{I_i，I_j}的图像I_j的各像素点与在图像I_i上的映射点的像素值之差的平方和；其中图像I_i的像素点q_i与在图像I_j的像素点q_j的映射关系为：K_i和K_j分别为摄像装置i、j的内参数矩阵，R_i和R_j为摄像装置i、j的外参数矩阵，基于列向量Ω中对应摄像装置i、j的姿态角确定。

与现有方法相比，本发明的摄像装置的工作姿态角采用的信息是相邻摄像装置所拍摄图像重合区内所有像素的灰度值，故对于少纹理区域仍然有效。

为了简化处理，降低运算复杂度，本发明还包括对构建的列向量Ω进行初始化：

指定任一摄像装置为参考系，并设置参考系对应的摄像装置的初始姿态角为0；基于所述参考系，根据非参考系的摄像装置所处的位置对其姿态角进行预估，得到各摄像装置的初始姿态角；基于各摄像装置的初始姿态角对列向量Ω进行初始化。

基于本发明的摄像装置的工作姿态角的估计方法，本发明还公开了一种易于在通用产品上实时生成任意指定视场范围的合成图像的方法。本发明的用于DAS的实时生成合成图的方法，包括下列步骤：

获取摄像装置的基本参数：包括内参数矩阵K_r和镜头畸变参数，以及外参数矩阵R_r；

获取用户指定的预合成图参数：包括虚拟内参数矩阵K_v、虚拟外参数矩阵R_v、虚拟成像平面尺寸；

基于所述摄像装置的基本参数和预合成图参数生成合成图：

根据反投影公式对虚拟成像平面的各像素点坐标x_v进行反投影，得到被拍摄场景点坐标X；

根据投影公式x_r=K_rR_rX，将各被拍摄场景点坐标X，投影到各摄像装置的成像平面上，获取与虚拟成像平面像素点坐标x_v对应的像素值，生成合成图，其中真实像素坐标x_r为摄像装置的成像平面的像素点坐标。

综上所述，由于采用了上述技术方案，本发明的有益效果是：用于图像拼接时的摄像装置姿态角估计精确、易于在通用产品上实时生成任意指定视场范围的合成图像。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是本发明具体实施方式中，DAS系统示意图，其中实线表示相机，虚线表示对应的视场角；

图2是本发明具体实施方式中，矩阵J和向量r的示意图，以包含6个摄像装置为例，J中的方块为3*3大小子矩阵，r中的方块为3*1大小的子矩阵；

图3是本发明实施例中，合成图的示意图，其中I_r1、I_r2、I_r3为参与生成合图的来自不同相机的图像，O₁₂和O₁₃为相邻相机图像在合成图中重叠区域。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书（包括任何附加权利要求、摘要和附图）中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

参见图1（图中是以M的取值为6所给出的示例），用于实现本发明的DAS的摄像装置采用M（M的个数要确保所述DAS摄像装置经过摆放得以覆盖360度视场范围）个敏通彩色可见光相机，配备广角镜头，系统包括两个单元：成像单元和数据处理单元。其中成像单元负责实际图像的采集，包括可见光相机、图像采集卡、以及控制相机进行同步图像采集的电路；数据处理单元负责对各个相机送出的图像数据进行处理并生成实时的合成图。

在DAS中，为了得到任意指定视角的合成图，需要计算虚拟相机（对应于所指定的视角）的成像平面（虚拟的）上每一个像素点的像素值。首先对该虚拟成像平面上的每一个像素点，需要计算与之对应的位于真实相机成像平面的像素点位置（二维坐标）。如果多个真实相机成像平面的像素点对应虚拟相机成像平面上的同一个像素点，则需要对这些像素点的像素值进行羽化融合。设置虚拟相机成像平面上某个像素点坐标为x_v，真实相机r(其内、外参数矩阵分别为K_r和R_r)的成像平面上与x_v对应的像素点坐标为x_r，则x_v和x_r之间的关系可以用公式(1)表示。

$x_v=K_v{R}_v{R}_r^{-1}{K}_r^{-1}{x}_r---\left(1\right)$

在公式(1)中，为了实现计算，x_v和x_r为对应的二维坐标的三维齐次坐标，

和

分别代表R_r和K_r的逆矩阵。

根据公式(1)可知，为了计算与点x_v对应的像素点坐标为x_r，需要计算R_r和K_r，K_r可以通过对真实相机r定标的方法得到，而R_r则需要通过估算得到。由于R_r由相机r的姿态角(方位角、俯仰角和绕轴线的旋转角)通过Rodriguez公式确定，所以首先需要计算相机r的姿态角。对于M个相机而言，假设每个真实相机r的外参数矩阵为R_ri(i=1,2,…，M)，对应的姿态角分别为ω_i=[ω_i1,ω_i2,ω_i3]^T（上标T表示矩阵的转置，下同），把所有ω_i=[ω_i1,ω_i2,ω_i3]^T串联后放入一个3M维的列向量Ω中，可基于数值优化的方法计算列向量Ω。

为此，我们首先利用M个真实相机r拍摄M对图像，其中每两个相邻相机拍摄一对图像，用{I_i，I_j}表示，其中i，j为M个摄像装置的标识编号，i,j∈{1,2,...,M},且i≠j，如用1号和2号相机拍摄一对图像，记为{I₁，I₂}，用2号和3号相机拍摄一对图像，记为{I₂，I₃}，…，用M号和1号相机拍摄一对图像，记为{I_M，I₁}。

基于所拍摄的M个图像对{I_i，I_j}，通过数值优化的方法计算列向量Ω的目的是通过迭代的方式使得下面公式中的e(Ω)最小。

$e\left(\mathrm{\Ω}\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{{\tilde{q}}_i\∈I_i,{\tilde{q}}_j\∈I_j}{\underset{({\tilde{q}}_i=W({\tilde{q}}_j;\mathrm{\Ω}),{\tilde{q}}_j),}{\mathrm{\Σ}}}{(I_j\left(W({\tilde{q}}_j;\mathrm{\Ω})\right)-I_j\left({\tilde{q}}_j\right))}^2;$

其中， $J=\left\{\begin{array}{ccc}(i+1),& \mathrm{if}& i\≤M-1\\ 1,& \mathrm{if}& i=M\end{array}\right.---\left(2\right)$

在公式(2)中，

和

是分别位于图像I_i和I_j中的对应像素点二维坐标（以下用

和q分别表示像素点的二维坐标和对应的齐次坐标，下同），是满足

的映射关系，即其中K_i和K_j分别为摄像装置i、j的内参数矩阵，可以采用当前公开的相机标定方法计算得到；R_i和R_j为摄像装置i、j的外参数矩阵，由列向量Ω中对应的ω_i=[ω_i1,ω_i2,ω_i3]^T和ω_j=[ω_j1,ω_j2,ω_j3]^T通过Rodriguez公式确定。

和

分别是图像I_i和I_j中对应像素点

和的图像灰度值，和的确定方法如下：

1）对于图像I_j中的整数位置像素点

（假设其坐标为

首先将其扩充为齐次坐标形式 $q_j={({\tilde{q}}_{\mathrm{jx}},{\tilde{q}}_{\mathrm{jy}},1)}^T;$

2）确定映射关系

时，通过计算齐次坐标

确定。

对于齐次坐标q_i，假设其坐标为(q_ix,q_iy,q_iz)^T，通过

得到对应的二维像素坐标如果

的坐标是非整数，则采用双线性插值方法得到对应的图像灰度值

如果对于图像I_j中的某个整数位置像素点

所映射对应的坐标

不在图像I_i的范围之内，则

和

不参与公式(2)的计算。

在本发明中，基于M个图像对{I_i，I_j}，可以通过下述迭代方式，得到参量e(Ω)取得最小时所对应的Ω的取值，也可以采用其他数值优化方式来获取，本发明对具体的数值优化方式并不限制。

通过迭代方式计算参量e(Ω)取得最小时所对应的Ω的取值的处理过程如下：

对于给定Ω的当前值，我们求解Ω的增量ΔΩ使得下一个值为Ω←Ω+ΔΩ。计算过程如下：

输入：列向量Ω的初始值（基于对各相机的姿态角的预估得到列向量Ω的初始值），M个相机的内参数矩阵{K_i}_i=1～M，M个图像对{I_i，I_j}，i,j∈{1,2,...,M},且i≠j,迭代次数D；

输出：经过优化计算后的向量Ω；

其中，D次迭代过程具体如下：

当前迭代次数d初始为0，每次迭代后自行加1，若d小于迭代次数D，则对每个图像对{I_i，I_j}进行以下操作：

1）计算图像I_i的梯度图像▽I_i={(▽I_i)_x,(▽I_i)_y}，其中(▽I_i)_x和(▽I_i)_y分别为I_i在水平方向和垂直方向的梯度图像，在计算▽I_i时，采用数字图像处理中常用的Sobel算子作用于图像I_i，得到对应的梯度图像(▽I_i)_x和(▽I_i)_y；

对图像I_j中的每个整数位置像素点

执行步骤2）～6）

2）计算

其中 ${\tilde{q}}_i=W({\tilde{q}}_j;\mathrm{\Ω});$

3）计算 $I_j\left({\tilde{q}}_i\right)-I_j\left({\tilde{q}}_j\right);$

4）基于图像(▽I_i)_x和(▽I_i)_y，采用双线性插值方法得到在位置

的元素值和

5）基于公式（3）计算下面由偏导数构成的二维列向量：

$\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i1}},\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i2}},\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i3}},\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j1}},\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j2}}$ 和

6）构造三维行向量

$V_i\left({\tilde{q}}_j\right)=\left(\begin{array}{cc}{\left({\▿I}_i\right)}_x\left({\tilde{q}}_i\right)& {\left({\▿I}_i\right)}_y\left({\tilde{q}}_i\right)\end{array}\right)\left(\begin{array}{ccc}\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i1}}& \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i2}}& \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i3}}\end{array}\right)$

$V_j\left({\tilde{q}}_j\right)=\left(\begin{array}{cc}{\left({\▿I}_i\right)}_x\left({\tilde{q}}_i\right)& {\left({\▿I}_i\right)}_y\left({\tilde{q}}_i\right)\end{array}\right)\left(\begin{array}{ccc}\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j1}}& \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j2}}& \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j3}}\end{array}\right)$

7）基于公式（4）计算3M×3M大小的矩阵J和3M维列向量r；

基于公式（5）计算ΔΩ；

更新列向量Ω：Ω←Ω+ΔΩ；

更新当前迭代次数d：d←d+1。

以下是完成上述次迭代过程所涉及的公式：

$\left\{\begin{array}{c}\frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i1}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}\frac{\∂q_i}{\∂{\mathrm{\ω}}_{i1}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}{K}_i\frac{\∂\left(R_i\right)}{\∂{\mathrm{\ω}}_{i1}}{R}_j^{-1}{K}_j^{-1}{q}_j\\ \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i2}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}\frac{\∂q_i}{\∂{\mathrm{\ω}}_{i2}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}{K}_i\frac{\∂\left(R_i\right)}{\∂{\mathrm{\ω}}_{i2}}{R}_j^{-1}{K}_j^{-1}{q}_j\\ \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{i3}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}\frac{\∂q_i}{\∂{\mathrm{\ω}}_{i3}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}{K}_i\frac{\∂\left(R_i\right)}{\∂{\mathrm{\ω}}_{i3}}{R}_j^{-1}{K}_j^{-1}{q}_j\\ \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j1}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}\frac{\∂q_i}{\∂{\mathrm{\ω}}_{j1}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}{K}_i{R}_i\frac{\∂\left(R_j^{-1}\right)}{\∂{\mathrm{\ω}}_{j1}}{K}_j^{-1}{q}_j\\ \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j2}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}\frac{\∂q_i}{\∂{\mathrm{\ω}}_{j2}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}{K}_i{R}_i\frac{\∂\left(R_j^{-1}\right)}{\∂{\mathrm{\ω}}_{j2}}{K}_j^{-1}{q}_j\\ \frac{\∂{\tilde{q}}_i}{\∂{\mathrm{\ω}}_{j3}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}\frac{\∂q_i}{\∂{\mathrm{\ω}}_{j3}}=\frac{\∂{\tilde{q}}_i}{\∂q_i}{K}_i{R}_i\frac{\∂\left(R_j^{-1}\right)}{\∂{\mathrm{\ω}}_{j3}}{K}_j^{-1}{q}_i\end{array}\right.---\left(3\right)$

在公式（3）中，用q_ix、q_iy和q_iz表示齐次坐标q_i的三个分量，则 $\frac{\∂{\tilde{q}}_i}{\∂q_i}=\left[\begin{array}{ccc}1/q_{\mathrm{iz}}& 0& -q_{\mathrm{ix}}/{\left(q_{\mathrm{iz}}\right)}^2\\ 0& 1/q_{\mathrm{iz}}& -q_{\mathrm{iy}}/{\left(q_{\mathrm{iz}}\right)}^2\end{array}\right];$ 再根据Rodriguez公式，对姿态角列向量Ω中对应相机i和j的分量ω_i=[ω_i1,ω_i2,ω_i3]^T和ω_j=[ω_j1,ω_j2,ω_j3]^T，记

记向量

的三个分量分别为

和

构造3×3大小的矩阵 ${\left[{{\hat{n}}_{}}_i\right]}_{\×}=\left[\begin{array}{ccc}0& -{\hat{n}}_{i,z}& {\hat{n}}_{i,y}\\ {\hat{n}}_{i,z}& 0& -{\hat{n}}_{i,x}\\ -{\hat{n}}_{i,y}& {\hat{n}}_{i,x}& 0\end{array}\right],$ 则 $\frac{\∂\left(R_i\right)}{\∂{\mathrm{\ω}}_{i1}}=e^{{\left[{\hat{n}}_i\right]}_{\×}}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -1\\ 0& 1& 0\end{array}\right],$ $\frac{\∂\left(R_i\right)}{\∂{\mathrm{\ω}}_{i2}}=e^{{\left[{\hat{n}}_i\right]}_{\×}}\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ -1& 0& 0\end{array}\right],$ $\frac{\∂\left(R_i\right)}{\∂{\mathrm{\ω}}_{i3}}=e^{{\left[{\hat{n}}_i\right]}_{\×}}\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right],$ 其中

是一个3×3的矩阵，可以通过MATLAB函数expm(X)，将矩阵

作为expm(X)的变量计算得到；同理，可计算得到

和

基于三维行向量

计算3M×3M大小的矩阵J和3M维列向量r的公式为：

$\left\{\begin{array}{c}{J}_{\mathrm{ii}}=\underset{\tilde{x}}{\mathrm{\Σ}}{\left(V_i\left({\tilde{q}}_j\right)\right)}^T{V}_i\left({\tilde{q}}_j\right)\\ J_{\mathrm{jj}}=\underset{\tilde{x}}{\mathrm{\Σ}}{\left(V_j\left({\tilde{q}}_j\right)\right)}^T{V}_j\left({\tilde{q}}_j\right)\\ J_{\mathrm{ij}}=\underset{\tilde{x}}{\mathrm{\Σ}}{\left(V_i\left({\tilde{q}}_j\right)\right)}^T{V}_j\left({\tilde{q}}_j\right)\\ r_i=\underset{\tilde{x}}{\mathrm{\Σ}}{\left(V_i\left({\tilde{q}}_j\right)\right)}^T(I_j\left({\tilde{q}}_j\right)-I_i\left({\tilde{q}}_i\right))\\ r_j=\underset{\tilde{x}}{\mathrm{\Σ}}{\left(V_j\left({\tilde{q}}_j\right)\right)}^T(I_j\left({\tilde{q}}_j\right)-I_j\left({\tilde{q}}_i\right))\end{array},\right.i,j=1~M---\left(4\right)$

将大小为3M×3M的矩阵J看成是由3×3大小的分块构成的M×M分块矩阵，则公式(4)中，J_ii、J_jj和J_ij分别是矩阵J中的第(i,i)、(j,j)和(i,j)分块。同理将大小为3M维的矩阵r看成是由3×1大小的分块构成的M×1分块矩阵，则r_i和r_j分别是r中的第i和第j个分块，如图2所示（图中是以M的取值为6所给出的示例）。

因此，列向量Ω的增量ΔΩ是由(J+λdiag(J))ΔΩ=r构成的线性方程组求解得到，如公式（5）所示：

ΔΩ=(J+λdiag(J))^-1r   (5)

在公式(5)中，diag(J)是由J的主对角线元素构成的对角矩阵，选择标量λ使得e(Ω+ΔΩ)<e(Ω)满足，得到当前迭代次数所对应的增量ΔΩ的取值。具体计算时，可将标量λ的初始值选为一个较小的数值，如0.1，如果e(Ω+ΔΩ)<e(Ω)不成立，则增加λ的取值直到e(Ω+ΔΩ)<e(Ω)满足为止。

本发明中，迭代次数D的取值根据具体应用情况设定，原则是使得e(Ω+ΔΩ)与e(Ω)的差值小于满足具体应用需求的范围为止。若对要求e(Ω+ΔΩ)与e(Ω)的差值小，则D可以设得较大，例如10及以上；反之，可以将其设为小于10的正整数。

为了进一步简化处理复杂度，可指定其中的一个相机为参考系，例如1号相机，这时ω₁=[0,0,0]^T。此外，将上述计算过程结合图像金字塔方式进行：对每对图像{I_i，I_j}构建分辨率由高到低的金字塔式的数据结构，金字塔底部图像分辨率最高，顶部图像分辨率最低，同级的图像分辨率相同。上述计算过程从金字塔顶部图像开始，迭代D次以后转入金字塔的下一级重新计算。这一过程持续进行直到金字塔底部图像完成计算为止，得到Ω的最终结果。

本发明中，在生成合成图时，既可以是基于本发明的姿态角评估方法来获取各真实相机r的外参数矩阵R_r，也可以是基于其他惯用的姿态角评估方法来获取。

在得到了M个相机的姿态角后，可以生成任意指定视角的场景的图像，用户指定的预合成图参数（包括虚拟相机内参数矩阵K_v和外参数矩阵R_v，以及成像平面的尺寸E_v×F_v）后，实现以下两个操作：

(i)对于虚拟相机成像平面上每一个像素点，快速计算与之对应的位于真实摄像装置成像平面的像素点位置；

(ii)如果多个真实摄像装置成像平面的像素点对应虚拟相机成像平面上的同一个像素点，则对这些像素点的像素值进行羽化融合。

对于(i)可以通过以下两个检验来实现，以确保获取到与合成图相关的像素值：

检验（1）：由于虚拟相机的视线方向可以任意指定，所以可以通过虚拟相机的视线（光轴）方向与M个真实摄像装置的光轴方向（将各摄像装置、虚拟相机的外参数矩阵分别乘以[0,0,1]^T得到对应的光轴方向）的夹角来确定参与合成图生成的T个摄像装置。

在实施过程中，分别计算M个真实摄像装置的光轴方向与指定的虚拟相机的视线方向之间的点积，对这些计算得到的点积按照从大到小的顺序排序后，挑出前T个（T的取值为T个相邻摄像装置所覆盖的视场角达到180度）最大点积对应的真实摄像装置。然后对这T个真实摄像装置分别应用检验（2）；

检验（2）：用x_v表示虚拟相机成像平面上的任一像素点坐标（表示形式为三维齐次坐标），对某个真实摄像装置r而言(其内、外参数矩阵分别为K_r和R_r)，目的是确定虚拟成像平面上与x_v对应的像素点坐标x_r。根据相机成像模型，首先反投影x_v以计算被拍摄场景的场景点X，

$X=R_v^{-1}{K}_v^{-1}{x}_v---\left(6\right)$

然后，将X投影到摄像装置r的成像平面上得到x_r，

x_r=K_rR_rX   (7)

根据摄像装置r的成像平面的尺寸检验x_r是否位于摄像装置r的视场范围之内：对坐标为(x_r,x,x_r,y,x_r,z)^T的齐次坐标x_r，基于公式

可得到与x_r对应的二维像素坐标

即可判别是否位于摄像装置r的视场范围之内。如果x_r被判断为位于摄像装置r的视场范围之内，则进一步反投影x_r，计算被拍摄场景点如果X′与X的之间的标量积X·X′大于0，则可判断x_v对应x_r。

在得到了虚拟相机成像平面上每个像素点对应的真实摄像装置中的对应像素点之后，对于(ii)，我们可以采用传统的羽化方式实现这些像素点的像素值融合；也可以采用本专利提出的增益补充后再羽化的方式实现这些像素点的像素值融合，使得在图像的重叠区过渡更自然和平滑。增益补偿的具体操作如下：

设参与合成图生成的摄像装置的图像数目为S,其中S≤T，对应的摄像装置的图像为I_rx，其中ξ∈{1,2,...,S}，设对应的增益补偿系数为a_ξ，增益补偿系数a_ξ的选取为使得e2取得最小值时所对应的取值：

$e2=\underset{\tilde{x},\tilde{y}\&Element;O_{\mathrm{\&mu;\&nu;}}}{\mathrm{\&Sigma;}}\{c_1{(a_{\mathrm{\&mu;}}{I}_{\mathrm{r\&mu;}}\left(\tilde{y}\right)-a_{\mathrm{\&nu;}}{I}_{\mathrm{r\&nu;}}\left(\tilde{x}\right))}^2+c_2[{(1-a_{\mathrm{\&mu;}})}^2+{(1-a_{\mathrm{\&nu;}})}^2\left]\right\}---\left(8\right)$

增益补偿系数a_ξ的值可以通过令偏导数为零计算得到。其中，O_μν为相邻摄像装置的图像I_rμ和I_rν在合成图像中的重叠区域，μ,ν∈{1,2,...,S}，

和

分别是I_rμ和I_rν在O_μν中对应像素点

和

的像素值；c₁和c₂是两个正实数，起到和[(1-a_μ)²+(1-a_ν)²]二者之间平衡的作用，由于[(1-a_μ)²+(1-a_ν)²]远小于

所以需要设定c₂/c₁＞＞1，通常设置c₁小于c₂至少2个数量级，c₁和c₂的具体数值可由使用者自行调整，直到取得满意的图像融合结果为止。经过充分实验发现c₁=0.01，c₂=100时效果较好。

在实现时，

和

数目的取值范围是从1到重合区O_μν中所有像素的数目之间，具体数值可由使用者根据设备的计算性能自行调整：如果计算资源丰富，可以取为重合区O_μν中所有像素的数目；反之，如果计算资源相对有限，但为了保证合成图生成的实时性，可以从重合区O_μν中所有像素随机抽取一定数目的像素。

在基于公式（8）得到增益补偿系数后，将图像I_rξ分别乘以增益补偿系数a_ξ，以得到增益补偿后的参与合成图生成的图像，再基于此提取对应的像素值（当然，也可以先提取对应的像素值，仅对获取的各像素值进行增益补充，即基于所获取的各像素值所对应的图像I_rξ，将像素值乘以对应的增益补偿系数a_ξ后），生成合成图。在生成合成图像时，因位于相邻摄像装置图像在合成景象图像中的重叠区域（例如图3所示的O₁₂和O₁₃）中的像素点x_v对应多个摄像装置成像平面中的像素，为此需要对这些像素点的像素值进行羽化处理，以生成一幅视觉上可以接受的合成图。

实施例1

参见图1，采用6个敏通彩色可见光相机（配备广角镜头）构成本实施例的DAS，6个相机可以固定在一个安装支架上，呈圆周均布，相邻相机之间按照60度夹角放置，由于安装支架加工时存在误差，以及考虑到相机在加工时光轴与理想位置之间存在偏差，所以60度夹角只是一个粗略的估计值，需要基于本发明的姿态角的估计来进一步细化。

6路相机的输出为PAL制式的模拟视频，送入两块微视图像采集卡中，从图像采集卡输出6路25帧/秒数字图像序列。6个相机采用可编程CPLD进行同步控制，由计算机控制采集卡发出TTL(Transistor-Transistor Logic)脉冲信号触发控制相机之间的同步。

系统运行时，首先获取摄像装置的基本参数：包括内参数矩阵、镜头畸变参数，以及外参数矩阵；用户指定的预合成图参数：包括虚拟内参数矩阵K_v、虚拟外参数矩阵R_v、虚拟成像平面尺寸。基于摄像装置的基本参数和预合成图参数生成合成图：

S1：根据每个摄像装置的内参数矩阵和镜头的畸变参数对每个摄像装置输出的每帧视频图像校正操作；

S2：根据虚拟相机的虚拟光轴方向，分别计算各摄像装置的光轴方向与虚拟光轴方向的点积，仅对前T个最大的点积所对应的相机输出的视频图像，执行步骤S3；

S3：基于公式（6），对x_v进行反投影得到X，基于公式（7）对X正投影得到x_r，

然后根据真实摄像装置输出图像的尺寸，判断得到的x_r是否位于真实摄像装置视场范围内，若否，则跳过对不在真实摄像装置视场范围内的像素点坐标的处理；若是，则进入步骤S4；

S4：根据

反投影计算拍摄场景点X′，若坐标X与X′的标量积大于0，则记录该像素点坐标；否则跳过该像素点坐标，处理下一个像素点坐标；

S5：基于摄像装置输出的经校正的图像，确定参与当前合成图生成的图像为I_rξ（ξ=1,2,...,S），基于公式（8）求解增益补偿系数a_ξ，将图像I_ri乘以对应的增益补偿系数a_ξ后，获取所记录的各像素点坐标所对应的像素值；

S6：当一个x_v对应多个像素值时，则通过羽化操作进行融合，使得虚拟成像平面的各像素点位置具备唯一像素值，生成合成图像并输出。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims (10)

1.用于DAS的摄像装置姿态角估计方法，其特征在于，包括下列步骤：

每相邻的两个摄像装置拍摄一对图像对{I_i，I_j}，得到M对参考图像对，其中i，j为M个摄像装置的标识；

构建姿态角列向量Ω，所述列向量Ω的元素依次对应M个摄像装置的姿态角；

获取M个摄像装置的内参数矩阵；

基于M对参考图像对、M个摄像装置的内参数矩阵，通过数值优化计算参量e(Ω)取得最小时所对应的列向量Ω的取值，得到各摄像装置的工作姿态角；

所述参量e(Ω)为：累加M个图像对中，各图像对{I_i，I_j}的图像I_j的各像素点与在图像I_i上的映射点的像素值之差的平方和；其中图像I_i的像素点q_i与在图像I_j的像素点q_j的映射关系为：

K_i和K_j分别为摄像装置i、j的内参数矩阵，R_i和R_j为摄像装置i、j的外参数矩阵，基于列向量Ω中对应摄像装置i、j的姿态角确定。

2.如权利要求1所述的姿态角估计方法，其特征在于，还包括对构建的列向量Ω进行初始化：

指定任一摄像装置为参考系，并设置参考系对应的摄像装置的初始姿态角为零；

基于所述参考系，根据非参考系的摄像装置所处的位置对其姿态角进行预估，得到各摄像装置的初始姿态角；

基于各摄像装置的初始姿态角对列向量Ω进行初始化。

3.如权利要求1或2所述的姿态角估计方法，其特征在于，对所述图像对{I_i，I_j}，构建分辨率由高到低的金字塔式的数据结构，金字塔底部对应的图像分辨率最高，顶部对应的图像分辨率最低，同级的图像分辨率相同；

则在通过数值优化计算参量e(Ω)取得最小时所对应的列向量Ω的操作时，从所述金字塔顶部往底部的方向进行优化。

4.用于DAS的实时生成合成图的方法，其特征在于，包括下列步骤：

获取摄像装置的基本参数：包括内参数矩阵K_r和镜头畸变参数，以及外参数矩阵R_r；

获取用户指定的预合成图参数：包括虚拟内参数矩阵K_v、虚拟外参数矩阵R_v、虚拟成像平面尺寸；

基于所述摄像装置的基本参数和预合成图参数生成合成图：

根据反投影公式

对虚拟成像平面的各像素点坐标x_v进行反投影，得到被拍摄场景点坐标X；

根据投影公式x_r=K_rR_rX，将各被拍摄场景点坐标X，投影到各摄像装置的成像平面上，获取与虚拟成像平面像素点坐标x_v对应的像素值，生成合成图，其中真实像素坐标x_r为摄像装置的成像平面的像素点坐标。

5.如权利要求4所述的方法，其特征在于，基于权利要求1所述的姿态角估计方法同时确定每个摄像装置的外参数矩阵R_r。

6.如权利要求4或5所述的方法，其特征在于，在对虚拟成像平面的各像素点坐标x_v进行反投影之前还包括：

将各摄像装置、虚拟相机的外参数矩阵分别乘以[0,0,1]^T得到对应的光轴方向；

分别计算各摄像装置的光轴方向与虚拟相机的光轴方向的点积，仅对前T个最大点积所对应的摄像装置执行步骤4；

所述T的取值为T个相邻摄像装置所覆盖的视场角达到180度。

7.如权利要4、5或6所述的方法，其特征在于，获取与虚拟成像平面像素点坐标x_v对应的像素值，生成合成图具体为：

（a）根据摄像装置的成像平面的尺寸，检验得到的各真实像素坐标x_r是否位于摄像装置的视场范围内，若否，则跳过该真实像素坐标；

若是，则再根据反投影公式

反投影计算被拍摄场景点坐标X′，若X与X′的标量积大于0，则记录该真实像素坐标；否则跳过该真实像素坐标；

（b）基于摄像装置输出的经校正的图像，获取所记录的各真实像素坐标所对应的像素值；

（c）根据投影公式

基于所获得的真实像素坐标x_r的像素值，得到虚拟成像平面的各像素点的像素值，若一个x_v对应多个像素值，则通过羽化进行融合，使得虚拟成像平面的各像素点坐标具备唯一像素值，生成合成图像。

8.如权利要求7所述的方法，其特征在于，对所述校正的图像进行增益补偿后，再获取所记录的各真实像素坐标所对应的像素值，所述增益补偿为：

设参与合成图生成的摄像装置的图像数目为S,各图像表示为I_rξ，其中ξ∈{1,2,...,S}，S≤T；

设与图像I_rξ对应的增益补偿系数为a_ξ；

将所述图像I_rξ分别与对应的增益补偿系数a_ξ相乘；

其中，所述增益补偿系数a_ξ为参量e2取得最小值时的对应值：

$e2=\underset{\tilde{x},\tilde{y}\&Element;O_{\mathrm{\&mu;\&nu;}}}{\mathrm{\&Sigma;}}\{c_1{(a_{\mathrm{\&mu;}}{I}_{\mathrm{r\&mu;}}\left(\tilde{y}\right)-a_{\mathrm{\&nu;}}{I}_{\mathrm{r\&nu;}}\left(\tilde{x}\right))}^2+c_2[{(1-a_{\mathrm{\&mu;}})}^2+{(1-a_{\mathrm{\&nu;}})}^2\left]\right\},$ 其中O_μν为相邻摄像装置的图像I_rμ和I_rν在合成图像中的重叠区域，μ,ν∈{1,2,...,S}，

和分别是图像I_rμ和I_rν在O_μν中对应像素点

和

的像素值；正实数c₁、c₂用于平衡

和[(1-a_μ)²+(1-a_ν)²]，所述c₁小于c₂，且c₁与c₂至少相差2个数量级。

9.如权利要求8所述的方法，其特征在于，所述c₁=0.01，c₂=100，N_x=N_y=100。

10.如权利要求7所述的方法，其特征在于，所述步骤（b）中，对所获取的各像素值进行增益补偿后，再执行步骤（c），所述增益补偿为：

设参与合成图生成的摄像装置的图像数目为S,各图像表示为I_rξ，其中ξ∈{1,2,...,S}，S≤T；

设与图像I_rξ对应的增益补偿系数为a_ξ；

基于所获取的各像素值所对应的图像I_rξ，将像素值乘以增益补偿系数a_ξ；

其中，所述增益补偿系数a_ξ为参量e2取得最小值时的对应值：

$e2=\underset{\tilde{x},\tilde{y}\&Element;O_{\mathrm{\&mu;\&nu;}}}{\mathrm{\&Sigma;}}\{c_1{(a_{\mathrm{\&mu;}}{I}_{\mathrm{r\&mu;}}\left(\tilde{y}\right)-a_{\mathrm{\&nu;}}{I}_{\mathrm{r\&nu;}}\left(\tilde{x}\right))}^2+c_2[{(1-a_{\mathrm{\&mu;}})}^2+{(1-a_{\mathrm{\&nu;}})}^2\left]\right\},$ 其中O_μν为相邻摄像装置的图像I_rμ和I_rν在合成图像中的重叠区域，μ,ν∈{1,2,...,S}，

和分别是图像I_rμ和I_rν在O_μν中对应像素点

和

的像素值；正实数c₁、c₂用于平衡

和[(1-a_μ)²+(1-a_ν)²]，所述c₁小于c₂，且c₁与c₂至少相差2个数量级。

Images