基于数字地图的数码像片方位元素获取方法

摘要

本发明提供了一种基于数字地图的数码像片方位元素获取方法，步骤1：拍摄数码像片和与该数码像片相对应的数字地图，该数码像片包括至少5个地物点，数字地图中包括与地物点对应的同名点，地物点及同名点统称为控制点，数码像片设有像片二维坐标系，数字地图设有数字地图坐标系；步骤2：获取控制点在数字地图坐标系中的坐标及在像片二维坐标系中的坐标和像片像主点在像片二维坐标系中的坐标，像片像主点为摄影物镜中心在像片上的垂足点；步骤3：建立摄影测量坐标系D-XDYDZD。本发明利用现有的数字地图为条件，计算普通数码图像的方位元素，为获取三维空间信息提供了一条简便、快捷、廉价的途径，外业工作量小，代价低。

权利要求

1.一种基于数字地图的数码像片方位元素获取方法，其特征在于：该方法包括如下步骤：步骤1：拍摄数码像片和与该数码像片相对应的数字地图，该数码像片包括至少5个地物点，数字地图中包括与地物点对应的同名点，地物点及同名点统称为控制点，数码像片设有像片二维坐标系，数字地图设有数字地图坐标系；

步骤2：获取控制点在数字地图坐标系中的坐标及在像片二维坐标系中的坐标和像片像主点在像片二维坐标系中的坐标，像片像主点为摄影物镜中心在像片上的垂足点；

步骤3：建立摄影测量坐标系D-XDYDZD，该摄影测量坐标系为以近似的拍摄点位置为原点、实地东方向为X轴、垂直向上方向为Y轴、实地南方向为Z轴的右手三维空间坐标系，将控制点在数字地图坐标系中的坐标转换为摄影测量坐标系中的坐标；

步骤4：利用控制点在摄影测量坐标系中的坐标和在像片二维坐标系中的坐标获取像片的方位元素，获取方法包括如下步骤：步骤41：利用控制点在摄影测量坐标系中的坐标和在像片二维坐标系中的坐标，根据二维直接线性变换原理计算像片二维直接线性变换参数，即：式1

其中， 其中，

(x、y)是控制点在像片二维坐标系中的坐标；(X、Y)是控制点在摄影测量坐标系中的坐标；k1代表物镜畸变参数；Δx、Δy是物镜畸变对像片二维坐标系的影响；l1、l2、l3、l4、l5、l6、l7、l8分别代表二维直接线性第一变换参数、第二变换参数、第三变换参数、第四变换参数、第五变换参数、第六变换参数、第七变换参数和第八变换参数；(x0、y0)为像片像主点在像片二维坐标系中的坐标，像片像主点为物镜中心在像片上的垂足点；

步骤42：引入旋转矩阵 设

为表述方便，引入过渡第一参数γ1、过渡第二参数γ2、过渡第三参数γ3，令可简化得：当x0＝0、y0＝0时得到

步骤43：由二维直接线性变换参数获取摄影物镜焦距，即其中，f为摄影物镜焦距；

进一步得到：

或

可以得到旋转矩阵中的六个参数：

同样可以得到：

步骤44：获取像片的方位元素，像片的方位元素包括摄影点在摄影测量坐标系中的坐标(XS、YS、ZS)和摄影的侧滚角 俯仰角ω和像片旋转角κ；

其中 ω、κ获取方法为：

摄影点在摄影测量坐标系中的坐标(XS、YS、ZS)获取方法为：f为摄影物镜焦距，x0、y0为像片像主点在像片二维坐标系中的坐标，像片像主点为物镜中心在像片上的垂足点，a1、b1、c1、a2、b2、c2、a3、b3、c3分别为旋转矩阵中不同的数值。

说明书

基于数字地图的数码像片方位元素获取方法

技术领域

[0001] 本发明涉及一种获取像片方位元素的方法，尤其涉及一种利用数字地图获取普通数码像片方位元素的方法。

背景技术

[0002] 计算像片的方位元素是摄影测量和计算机领域的重要内容，利用像片的方位元素可以进一步对摄影目标进行三维测量。目前，计算像片的方位元素的方法主要有两种方法：借助控制点计算和GPS/INS观测的方法。

[0003] 借助控制点计算的方法是利用一定数量、合理分布的三维控制点(计算机视觉中又称为黄金点)，建立与相应影像点的对应关系，计算像片的参数具有很广泛的研究和应用，该方法同样可以实现视频帧的内外方位元素的解算(张祖勋.数字摄影测量与计算机视觉[J].武汉大学学报·信息科学版，2004.12，12(29))。目前已广泛应用于摄影测量领域。

[0004] 利用GPS/INS观测的方法是近十多年来发展而来的新技术，采用GPS/INS组合系统(简称POS系统)来获取航空摄影时影像的空间方位(即利用GPS确定摄站的空间位置，利用IMU惯性测量装置获取影像的姿态角)，以直接用于航测内业的像片定向，目的是取代摄影测量加密工序(CannonM E，Sun H.Experimentalassessment of a non-dedicate GPS receiver system for airborne attitudedetermination[J].ISPRS Jour nal of Photogrammetry & Remote Sensing，1996，51(2)：99-108.)。同样，该技术也可以用于地面上，通过在机动车上装配GPS/CCD/INS/或航位推算系统等先进的传感器和设备，可以在车辆的高速行进之中快速采集道路及道路两旁地物的空间位置数据和属性数据，如：道路中心线或边线位置坐标，目标地物的位置坐标，路、车道、宽、桥、隧道高、交通标志、坡度、道路设施等数据同步存储在车载计算机系统中，经事后编辑处理，形成内容丰富的道路空间信息数据库(李德仁.移动测量技术及其应用[J].地理空间信息.2006年8月：1-5)。

[0005] 目前采用的技术方法尽管精度较高，但也有其不足。借助控制点计算的方法无法做到实时性的要求，同时工作量大。而GPS/INS观测高昂的价格一直制约着它的推广应用(何秀凤，徐勇，桑文刚.微小型IMU/GPS组合定位定向系统研究[J].武汉大学学报·信息科学版.第30卷第11期，2005年11月：991-994)。

发明内容

[0006] 技术问题：本发明要解决的技术问题是提供一利用数字地图计算普通数码像片方位元素的方法，使得获取三维空间信息更简便、快捷、廉价。

[0007] 技术方案：为解决上述技术问题，本发明采用的技术方案是：一种基于数字地图的数码像片方位元素获取方法，该方法包括如下步骤：

[0008] 步骤1：拍摄数码像片和与该数码像片相对应的数字地图，该数码像片包括至少5个地物点，数字地图中包括与地物点对应的同名点，地物点及同名点统称为控制点，数码像片设有像片二维坐标系，数字地图设有数字地图坐标系；

[0009] 步骤2：获取控制点在数字地图坐标系中的坐标及在像片二维坐标系中的坐标和像片像主点在像片二维坐标系中的坐标，像片像主点为摄影物镜中心在像片上的垂足点；

[0010] 步骤3：建立摄影测量坐标系D-XDYDZD，该摄影测量坐标系以近似的拍摄点位置为原点、实地东方向为X轴、垂直向上方向为Y轴、实地南方向为Z轴的右手三维空间坐标系，将控制点在数字地图坐标系中的坐标转换为摄影测量坐标系中的坐标；

[0011] 步骤4：利用控制点在摄影测量坐标系中的坐标和在像片二维坐标系中的坐标获取像片的方位元素，获取方法包括如下步骤：

[0012] 步骤41：利用控制点在摄影测量坐标系中的坐标和在像片二维坐标系中的坐标，根据二维直接线性变换原理计算像片二维直接线性变换第一参数、第二变换参数、第三变换参数、第四变换参数、第五变换参数、第六变换参数、第七变换参数和第八变换参数，即：

[0013] 其中，

[0014] (x、y)是控制点在像片二维坐标系中的坐标；(X、Y)是控制点在摄影测量坐标系中的坐标；k1代表物镜畸变参数；Δx、Δy是物镜畸变对像片二维坐标系的影响；l1、l2、l3、l4、l5、l6、l7、l8分别代表二维直接线性第一变换参数、第二变换参数、第三变换参数、第四变换参数、第五变换参数、第六变换参数、第七变换参数和第八变换参数；(x0、y0)为像片像主点的在像片二维坐标系中的坐标，像片像主点为物镜中心在像片上的垂足点；

[0015] 引入过渡第一参数γ1、过渡第二参数γ2、过渡第三参数γ3，令

[0016] 可以得到：

[0017]

[0018] 步骤42：由二维直接线性变换参数获取摄影物镜焦距，即

[0019] 其中，f为摄影物镜焦距；

[0020] 步骤43：引入旋转矩阵 设

[0021]

[0022] 可以得到：

[0023]

[0024] 步骤44：获取像片的方位元素，像片的方位元素包括摄影点在摄影测量坐标系中的坐标(XS、YS、ZS)和摄影的侧滚角、俯仰角和像片旋转角分别用 ω、κ表示；

[0025] 其中 κ获取方法为：

[0026]

[0027] 摄影点在摄影测量坐标系中的坐标(XS、YS、ZS)获取方法为：

[0028]

[0029]

[0030] f为摄影物镜焦距，x0、y0为像片像主点的图像坐标，像片像主点为物镜中心在像片上的垂足点。

[0031] 有益效果：本发明的有益效果在于：利用现有的数字地图为条件，计算普通数码图像的方位元素，为获取三维空间信息提供了一条简便、快捷、廉价的途径，外业工作量小，代价低。

附图说明

[0032] 图1为本发明提供的像片方位元素计算方法的流程图。

具体实施方式

[0033] 下面结合附图对本发明做进一步说明。

[0034] 本发明提供的计算方法能够实现对数码像片拍摄位置、摄影姿态、像片焦距和畸变的计算，可广泛应用于数字城市建设和三维建模，属于近景摄影测量和计算机视觉领域。本发明在利用二维直接线性变换公式，根据数码像片内方位元素中的x0、y0近似等于0的特点，提出了计算像片内方位元素的方法，参见图l，其步骤如下：

[0035] 步骤S11：获取数码像片和与该数码像片相对应的数字地图，该数码像片包括至少5个地物点，数字地图中包括与地物点对应的同名点，数码像片设有像片二维坐标系，数字地图有数字地图坐标系；

[0036] 步骤S12：利用Visual C#平台开发的软件获取同名点在数字地图坐标系中的坐标和地物点在像片二维坐标系中的坐标；

[0037] 步骤S13：建立摄影测量坐标系，该坐标系以近似的摄影点位置为原点、实地东方向为X轴、垂直向上方向为Y轴、实地南方向为Z轴的右手三维空间坐标系，将同名点在数字地图坐标系中的坐标转换为摄影测量坐标系中的坐标；

[0038] 步骤S14：利用同名点在摄影测量坐标系中的坐标和地物点在像片二维坐标系中的坐标计算像片的方位元素。

[0039] 步骤S14包括如下步骤：

[0040] 步骤S141：首先利用控制点在摄影测量坐标系中的坐标和在像片二维坐标系中的坐标计算像片二维直接线性变换参数和像片畸变参数；

[0041] 步骤S142：由二维直接线性变换参数计算摄影机的焦距；

[0042] 步骤S143：根据二维直接线性变换参数和摄影机拍摄时的焦距进一步计算像片的外方位元素(包括摄影位置参数XS、YS、ZS和摄影姿态参数 ω、κ)，然后将摄影位置参数XS、YS、ZS转换到地图坐标系中。

[0043] 也就是说，本发明提供的方法包括：

[0044] 第一步：获取研究范围的数字地图。

[0045] 第二步：用普通数码相机拍摄研究目标的数码像片。拍摄时要求位于同一水平面上的地物点的个数在5个以上，即利用数码相机进行拍摄时，要包含至少5个地物点，这些地物点在实地应处于同一水平面上，并且成像清晰，便于识别。同时这些地物点必须是数字地图存在的点。

[0046] 第三步：选择数字地图和数码像片上一一对应的同名点。通过编制相应的软件获取同名点在数字地图坐标系中的坐标和在像片坐标系中的二维坐标。

[0047] 第四步：建立摄影测量坐标系，具体方法见表1。然后将同名点在地图坐标系中的坐标转换为摄影测量坐标系中的坐标。

[0048] 表1 摄影测量坐标系与地图坐标系的关系

[0049]坐标系 X Y Z
地图坐标系 北 东 高程
摄影测量坐标系 东 高程 南
关系 X摄影测量＝Y地图 Y摄影测量＝Z地图 Z摄影测量＝-X地图

[0050] 第五步：利用5个以上同名点的摄影测量坐标和二维像片坐标计算像片的方位元素。

[0051] 计算像片的方位元素，步骤如下：

[0052] ①计算二维直接线性变换参数

[0053] 考虑影像畸变，即为二维直接线性变换公式：

[0054]

[0055] 其中 包含像片畸变参数k1。

[0056] 式中：x、y是控制点的图像坐标；

[0057] X、Y是控制点的摄影测量坐标；

[0058] Δx、Δy是物镜畸变对图像坐标的影响；

[0059] l1、l2、l3、l4、l5、l6、l7、l8是二维直接线性变换参数；

[0060] k1代表物镜畸变系数；

[0061] x0、y0 x0、y0为像片像主点的图像坐标，像片像主点为物镜中心在像片上的垂足点[0062] 根据最小二乘法原理，利用间接平差计算公式(1)中的二维直接线性变换参数l1……l8和畸变参数k1。

[0063] ②计算像片方位元素与二维直接线性变换之间的关系。

[0064] 根据平面点高程相同的条件，将摄影测量中的共线方程进行变换为：

[0065]

[0066] 式中：f是摄影物镜焦距；

[0067] al、bl、cl是旋转矩阵 中的各值。

[0068] 通分合并为：

[0069]

[0070] 为表述方便，引入过渡参数γ1、γ2、γ3，令

[0071]

[0072] 可简化得：

[0073]

[0074] 当x0＝0、y0＝0时得到

[0075] 式中：l1、l2、l3、l4、l5、l6、l7、l8是二维直接线性变换参数。

[0076] 对公式(3)的参数进行代数运算：

[0077]

[0078]

[0079] 由以上两式可以解算影像焦距：

[0080]

[0081] 进一步得到：

[0082] 或

[0083] 代人到(3)式，可以得到旋转矩阵中的六个参数：

[0084]

[0085] 同样可以得到：

[0086]

[0087]

[0088] 摄影姿态参数

[0089]

[0090] 式中： ω、κ是摄影瞬间的摄影姿态参数。

[0091] 其中 ω、κ分别表示摄影的侧滚角、俯仰角和像片旋转角，f为像片焦距，[0092] 摄影位置参数的解算：

[0093]

[0094]

[0095] ③利用②中的格式由二维直接线性变换参数计算像片方位元素。

[0096] 下面更进一步的说明本发明。

[0097] 选取所拍摄的像片中中位于同一地平面上的六个点，在像片上测量像素点坐标，在数字地图上测量二维地图坐标。各坐标系的算例观测数据见表2。

[0098] 表2 平面控制结构控制点坐标

[0099]

[0100] 计算步骤和结果如下：

[0101] (1)计算二维DLT参数和影像畸变

[0102] 计算结果为：

[0103] l1＝-0.0073896

[0104] l2＝-2.8547409

[0105] l2＝3159.7316

[0106] l4＝-0.3400813

[0107] l5＝0.0011163

[0108] l6＝353.6026

[0109] l7＝-0.000970

[0110] lg＝-0.0000002-10

[0111] k1＝-6.4×10

[0112] (2)计算焦距

[0113]

[0114] (3)计算γ1、γ2、γ3的值

[0115] 先计算

[0116]

[0117] 然后计算

[0118]

[0119] (4)计算旋转矩阵

[0120]

[0121] (5)计算外方位元素，然后转换回地图坐标系中，计算结果为：

[0122]

[0123] ω＝0.000209(弧度)

[0124] κ＝-0.000405(弧度)

[0125] XS＝1031.344(米)

[0126] YS＝1104.133(米)

[0127] ZS＝10.445(米)