基于直接线性变换算法的普通数码相机检校的应用研究.doc

基于直接线性变换算法的普通数码相机检校的应用研究孔 建 黄建魏 沈 周 (西南交通大学 四川成都 610031 中铁十局 山东济南 520000)摘要：本文采用直接线性变换（DLT）算法，完成了普通数码相机检校的应用研究。通过编程实验，解算普通数码相机在不同焦距情况下内方位元素（,）以及畸变参数（径向畸变系数，、偏心畸变系数，），同时对直接线性变换方法中初值的问题给出解决方案。提出了解决控制点布设在一个近似平面上解算系数初始值的方法，并且依据实验数据分析了在不同焦距下，相机内方位元素和光学畸变参数的变化情况。关键字：直接线性变换；相机检校；径向畸变；偏心畸变AbstractIn this paper, to complete a common application of digital camera calibration by using the direct linear transformation algorithm. This paper have solved different elements of interior orientation （,）and distortion parameters（Radinal Distortion ，,Decentering Distortion ，）of ordinary digital camera focal length by the programming experiments and meanwhile, put forward the solutions of the initial value problem in the direct linear transformation method. Proposed a solution in an approximate control points for solving plane initial value coefficient method, and analyzed the changes of the camera orientation elements and optical distortion parameters in the base of experimental data at different focal lengths.1 概述在数字摄影测量中，数字影像的获取，通常采用的是专业的摄影设备。这些专业设备的价格昂贵，对非专业部门是无法应用的。随着数码相机技术的发展与进步，普通数码相机在数字摄影测量领域中得到了广泛的应用，尤其是在近景数字摄影测量、无人机低空摄影测量的应用中，表现出了巨大的优势。普通数码相机不仅价格便宜，且操作方便，是专业摄影机不能比拟的。随着数码相机技术的不断进步，其像幅、分辨率不断提高，数码相机相对于专业摄影机的优势更加凸显。数码相机是非量测相机，应用到数字摄影测量中有其自身的缺陷，例如其像主点位置是未知的；相机镜头存在较大的光学畸变等。因此，为了保证影像的精度，获得可靠的测量结果，必须对数码相机的内方位元素及相关参数进行精确的检校确定。非量测数码相机的内方位元素及相关参数检校的算法有很多，其中一种常用的方法是直接线性变换（DLT）法。该方法解算过程中无需内、外方位元素的初始值，故特别适用于非量测相机的检校。其基本思想是在空间设置足够数量（大于等于6个）的控制点，用相机对控制点所在的目标进行摄影，最后将像点坐标和空间坐标代人直接线性变换的方程，引入光学畸变误差进行整体评差，即可解算出该相机的内外方位元素及相关参数。2 非量测相机检校的内容和原理2.1 相机检校内容非量测像机检校的主要内容有：（1） 像主点位置（）与主距的测定；（2） 光学畸变系数的测定，包括径向畸变系数（，）和偏心畸变系数（，）的测定；（3） 变焦后畸变差变化的测定；（4） 变焦后主距变化的测定；（5） 摄影机偏心常数的测定。本文对摄影机检校研究的内容主要是内方位元素的检校和光学畸变参数的解算。2.2非量测相机检校原理2.2.1 内方位元素的确定像片的内方位元素是恢复摄影时光束状态的要素；内方位元素确定了摄影中心S与所摄像片P相对位置关系，依据相对位置关系可以恢复摄影时光束的形状。借助内方位元素可以确定摄影中心与所摄像片间的位置关系，即恢复光束（光线，）在摄影时候的形状。由此可见，确定内方位元素对摄影测量来说是重要的一步。理论上，像主点在底片（模拟影像）或者CCD（数字影像）的中心，但由于相机设计制造时的偏差，因此无论是量测相机还是非量测相机，像主点通常不与框标的连线或者像片中心重合，因此需要解算出像主点的偏差值（）。如图21。若像主点在框标坐标系中的原点坐标为，量测出来的像点坐标化算到以像主点为原点的像平面坐标系的坐标为（）。图21 内方位元素2.2.2 非量测相机的光学畸变差（1） 径向畸变差（Radial Distortion）径向畸变差使得构像点沿向径方向偏离其准确理想位置。设在主距为的标准像片上，物方点A的标准位置为，实际构像于点，而且像方构像角与物方角不等，如图23。由图可知，径向畸变差与物方点的入射角有关，影像上不同点位的畸变差，因其角度的不一样而不同。图22 径向畸变像点的径向畸变可以用如下奇次多项式来表达： （21）其中，和是以像主点为原点的像片面坐标。 本文进行径向畸变改正时候只考虑了，因为与焦距具有相同的影响，因此不能与焦距校准同时发生1。，的精度取决于控制点的精度。控制点像素坐标的量测精度应该高于0.2个像素才能用于相机校准1。（2） 偏心畸变差（Decentering Distortion）物镜系统的各单元透镜，因装配和震动偏离了轴线或者歪斜，从而引起像点偏离其理想位置。如图24，当轴线偏离其设计轴线或者旋转了一个角度，影响将会产生变形。图23 引起偏心误差的原因在近景数字摄影测量应用中，偏心畸变会随着焦距D的变化而变化；而且不在焦距D上的物点也存在偏心畸变差。一般情况下，偏心畸变比径向畸变小，但是这样的畸变仍然不能忽视。偏心畸变差的表达式如下： （22）式中：焦距为D时的偏心畸变差分量； 焦距为D时的摄影时的主距； 偏心畸变系数； r 像点向径 r=； 是像主点坐标；图24即为径向畸变差和偏心畸变差对像点的共同影响。图24 径向畸变和偏心畸变12.2.3 不垂直性误差d和比例尺不一误差ds当对数码相机进行量测时，像点坐标的改正数（）中还包含了坐标轴不垂直性误差d和比例尺不一误差ds。如图25，像素坐标系是非直角坐标系，其两坐标轴之间的不垂直度为d。以像主点为原点有两个坐标系，分别是直角坐标系和非直角坐标系。像主点的坐标为（）。像点在像素坐标系中的坐标是（），在非直角坐标系中的坐标为（）。由于受和的影响，此坐标包含线性误差。点p是点的理想位置，它在中的坐标是（）；其中，。假设向无比例误差，而方向比例尺规划系数为1+。此时向的主距为，方向的主距为： （23）图25 坐标不垂直性误差和比例误差 比例尺误差可以认为是所用的像素坐标系轴和轴的单位长度不一致及摄影材料的不均匀变形引起的。不正交性误差认为像素坐标轴和轴的不垂直性引起的。因此和的改正应为： （24）3 直接线性变换解法 非量测相机检校，常用的方法有光学实验室检校法，试验场检校法和在任检校法。其中后面两种方法皆是依据物方空间分布一群合理的高质量控制点，并依据单片空间后方交会解法或者多片空间后方交会解法，解求像片内外方位元素及各类光学畸变系数。本文所采用的方法是试验场检校法。单片后方交会的解法主要有：基于共线条件方程式的单像空间后方交会、基于共面条件方程的空间后方交会解法、直接线性变化解法和基于角锥原理的空间后方交会解法。基于共线条件方程式的单像空间后方交会解法缺点为解算前必须知道外方位元素的初始值，因而较适用于量测相机；基于共面条件方程的后方交会解法利用了像片间内在的几何关系，因而对控制点的要求比较少，但是这种解法精度稍低2。基于角锥原理的空间后方交会，能够很好的避免控制点的构想范围较小时候，外方位元素的不稳定，但是这种解法比较复杂。综合上述考虑，本文采用直接线性变（DLT）的解法，该方法建立了像平面坐标和空间坐标之间的关系，计算中无需内外方位元素的初始值，因此特别适用于非量测相机的检校。直接线性变换的模型如下： （31）引入径向畸变差、偏心畸变差、坐标轴不垂直误差和比例尺误差并线性化后，其误差方程为： （32）根据间接平差的原理，此误差方程式及其相应的法方程式可以写成： （33）这个运算过程是一个迭代的过程，以相邻两次迭代运算的差值是否小于0.01mm作为判断。3.1 存在的问题及解决方案3.1.1 存在的问题 若初始解算的控制点分布在一个平面上，则式(33)中的矩阵不是满秩矩阵2。直接线性变换解法不能用于控制点分布或者近似在一个平面的场合。因此当用非量测相机对近似在一个平面上的地物（建筑物）进行摄影时，使用DLT解法就出现了问题。在使用十一个方程直接解算系数时，的初值无法求解。在平差理论中，如果有多余观测值，在解求观测值或者待定值的平差值时，通常会在等式两边同时乘以系数矩阵的转秩矩阵。如果在解求的初值时，也在方程两边同时乘以系数矩阵，是否能够解决初值的问题呢？3.1.2 解决方案首先，用n（n=6）个控制点列出2n个的方程式，即： （34）所以， （35）式（31）中有11个系数，解算系数的初值至少需要十一个方程，即需要6个控制点。如果直接使用十一个方程解算，由于矩阵（1,1）、（3,3）、（5,5）、（7,7）都等于0；因此使用一般的矩阵求逆方法不能求解的逆矩阵。使用（35）式，将矩阵与相乘，能够避免矩阵M中的上述四个元素为0，因而一般矩阵求逆方法也能使用。由于直接线性变化是一个迭代的过程，若所选择的6个控制点在一个近似的平面上，计算的的初值不稳定，但是将所有控制点都代入计算，在迭代的过程中对的值进行修正，因而最后解算的结果比较稳定。3.2 数据分析控制点的像点坐标和三维坐标如下： 控制点的像素和三维坐标 表31点号x（像素）y（像素）X(m)Y(m)Z(m)20001514.13040.4497.5045446.2537295.687320002496.42413.5497.4540446.2532299.894820003487.21758.5497.4357446.2199304.299420004505.21136.5497.5883446.2048308.5020200051287.82947.5504.8410454.1712297.0586200061191.21576.2504.4446454.3183304.8016200081867.92694.2508.8083457.6446298.6852200091836.71883.5508.7136457.6433302.8787200101845.81077.9508.8224457.6420307.0718200111622.1906.5506.0524449.9470309.2177200122718.23052.2512.9311454.3017296.39912001327032322.5512.9222454.3034300.4483200142677.51561.7512.8625454.3071304.6817200163926.42322.1520.3013446.2188300.1135200173779.11682.9519.4427446.2059304.2976200205116.32943527.9277446.2150295.9279200215124.62297.3528.0793446.2067300.1111200225046.31653527.6832446.2027304.3088200235038.01011.0527.7622446.1945308.5126 本文实验摄站的三维坐标的初始值为（500,500,300）；比较表31中的Y坐标值，可知最大值与最小值之差=457.6446-446.2027=11.4491（m）；摄站的Y坐标与控制点的平均Y坐标值的约为500-450=50（m）。因此控制点相对于摄站并不能看作在一个平面上或者近似的一个平面上。控制点的分布及定义的独立坐标系O-XYZ如图： 图31 控制点的分布及定义的独立坐标系O-XYZ下面分四种种情况进行比较：1）由于受目标点的限制，在布设控制点时候，只有四个平面，如表31中的Y坐标所示。因此近似选择六个Y坐标相差较大的点（2001、2002、2005、2010、2011、2020）用来解算的初始值，最后将所有的控制点代入方程，解算值和， ，；2）用控制点（2003、2006、2010、2011、2020、2023）列十一个方程，解算的初值，最后将所有的控制点代入方程，解算值和， ，；3）任意代入六个控制点（如2001、2002、2003、2004、2005、2006）解算的初始值，最后将所有的控制点代入方程，解算的值和， ，；4）所选择控制点（2001、2002、2003、2020、2021、2022、2023）近似在一个平面上，解算的初始值，最后将所有的控制点代入方程，解算值和， ，。上面四种情况所解算的外方位元素、内方位元素及值如表32： 各个待求参数的最终解 表32控制点不在一个平面上十一个方程解算任意六个控制点 控制点在一 个 近似的平面最大值与最小值之差单位516.2523516.2510516.2559516.25280.0013m498.9267 498.9266 498.9262 498.9180 0.0082m301.5005 301.5005 301.4973 301.5000 0.0027m1.087465 1.088100 1.095114 1.085520 0.001946rad1.504859 1.505283 1.502772 1.504463 0.001691rad1.107982 1.108615 1.115540 1.106022 0.009518rad16.88 19.58 0.67 14.71 18.91像素12.95 14.79 9.79 10.73 5.00像素8051.31 8051.27 8049.73 8048.50 2.81像素8050.47 8050.15 8050.44 8047.79 2.69像素8.187E-056.843E-051.015E-049.938E-053.305E-05m1.043E-041.381E-04-8.725E-058.837E-052.253E-04m1.253E-091.243E-091.185E-091.135E-091.177E-10m3.135E-173.300E-173.692E-174.601E-171.466E-17m-2.519E-09-1.415E-086.490E-081.113E-087.905E-08m3.109E-082.189E-086.503E-082.927E-084.315E-08m分析表32可知： 1）比较表32的前面六行的各列数据，发现外方位元素解算值非常的接近，其中摄站坐标中相差最大的是为0.0082m；外方位角元素相差最大的是为0.009518（rad）=03243，不到一个厘米，这是空间后方交会自身的缺点导致的；内方位元素中、最大值与最小值的差值最大的是=18.919（像素）=18.9196.41=116.2（m），后面的六个参数的值相差更小。 2）表32种各个参数变换比较大的主要是第四列（任意六个控制点），其原因主要是2001到2006这六个点集中分布在像片的左半边（如图31），没有均匀分布在像片的构想范围内。 3）如果采用式（35），即使用多于方程解算系数的初值，虽然解算出来系数的初值相差比较大，不够稳定，但是如果有足够多的控制点，即使起算控制点都集中分布在一个近似的平面上面，对最后内方位元素、外方位元素、解算的结果影响不是很大，但是应当注意初始起算数据的控制点最好均匀在整幅像片中。4 数据分析以往对于非量测相机精度的探讨和分析都是基于同一个焦距下利用非量测相机获取的像片进行研究。如果将非量测相机直接用在各项工程中，由于野外环境比较恶劣，如摄影时候光线较差，摄站点与目标间的距离受限，在对目标摄影过程中，目标在影像上很小等，如果还是使用同一个主距，在应用中势必受到限制。因此有必要研究非量测相机在不同的焦距下，其内方位元素和物镜构像畸变参数是否会发生改变。4.1试验方案为了探讨数码相机检校的精度，现在利用Canon EOS-1Ds Mark III对西南交大测量试验场进行摄影，采用佳能24-75毫米的变焦镜头，摄影主距分别为24mm，50mm，和60mm。每个主距摄影三张以上的像片。为了使试验数据稳定，采用6个均匀分布在像片上的控制点（2001、2006、2010、2016、2020、2023）作为起算数据，像点坐标的量测采用的是更加精确的LISA进行量测，其标称精度为0.2个像素。在表411中像片1808与1803、1804大致在同一个位置拍摄，但在拍摄1808时，首先将主距调到50mm拍摄另外一组像片后在调回24mm，这样能够检验当相机主距重新调整后在调回去，内方位元素是否会发生改变。像片1830是同1826、1827、1828大致在同一个位置拍摄，但1830使用的是60mm主距。 解算参数值（主距24mm） 表411180818031804单位XS508.1348508.3799508.3757mYS499.3719498.8445498.9346mZS301.5944301.5934301.599mphi1.2995-0.8602-0.741radomega1.56281.52811.5279radkappa-1.28840.87330.7538radx0-22.33-44.11-34.43像素y029.9927.5519.89像素fx3816.543797.673799.67像素fy3814.983796.443799.02像素ds-1.14E-043.48E-044.13E-04Dbeta4.09E-043.24E-041.71E-04k11.12E-088.44E-098.57E-09mk2-5.64E-163.34E-161.55E-16mP11.99E-077.60E-077.04E-07mP2-3.54E-07-5.24E-07-4.34E-07m 解算参数值（主距50mm） 表412182618271828单位XS503.4188503.4285503.4041mYS501.5148501.4438501.5247mZS301.6393301.6225301.6176mphi-1.3776-1.4594-1.5185radomega1.43221.43721.4449radkappa1.37711.4621.522radx0-10.23-12.21-8.93像素y013.3733.7246.19像素fx7614.637606.87624.53像素fy7613.977608.427624.67像素ds4.34E-043.36E-042.42E-04Dbeta8.62E-05-2.14E-04-1.81E-05k1-1.55E-09-1.98E-09-1.57E-09mk2-1.14E-16-6.00E-17-1.04E-16mP12.67E-072.53E-071.44E-07mP2-3.83E-07-4.85E-08-2.06E-07m 解算参数值（主距60mm） 表413183018321833单位XS503.4332520.0471520.0455mYS501.424501.8119501.7926mZS301.624301.5953301.593mphi-1.41941.50091.5499radomega1.4181.44371.4523radkappa1.4262-1.5046-1.5461radx0-18.95-26.67-15.86像素y055.7869.5881.8像素fx9464.029490.799480.35像素fy9461.729493.989480.7像素ds1.84E-042.57E-043.61E-04Dbeta2.44E-04-3.35E-04-3.66E-05k1-2.61E-09-2.66E-09-2.55E-09mk2-3.41E-188.45E-18-7.33E-18mP12.32E-072.62E-071.90E-07mP2-2.67E-072.75E-07-7.79E-07m4.2 解算精度分析4.2.1 摄站坐标精度分析 由于摄影时候是手持摄影，对每一张像片摄影的过程中，不可能保持同一个位置。而1826、1827、1828、1830是在间隔更短时间内拍摄的，拍摄瞬间位置基本上没变，因此这几张影像的摄站坐标（XS,YS,ZS）能够更好的反应解算的精度。其中XS、YS、ZS的最大值与最小值的偏差分别为3cm、10cm、2cm。其中YS的偏差最大是由于后方交会解法本身的缺陷所致，其解算的精度受到交会角大小的限制。4.2.2 内方位元素精度分析 内方位元素的均值和方差（单位：像素） 表4224mm50mm60mm均值中误差均值中误差均值中误差-33.62 8.91 -12.36 1.35 -20.49 4.55 25.81 4.30 15.06 13.53 69.05 10.63 3804.63 8.46 7615.32 7.25 9478.39 11.02 3803.48 8.20 7615.69 6.74 9478.80 13.24 在表42中内方位元素（，）最大中误差=13.53像素=0.0867mm。比较不同主距的，可知，在不同的主距下，相机的内方位元素是变化的，因此使用对同一个相机，使用不同的主距，其具有不同的内方位元素。比较表411中1808与1803、1804内方位元素（，），可知重新调焦后对相机的内方位元素（，）的影响比较大。4.2.3 相机参数变化 各项畸变参数均值和方差 表4324mm50mm60mm均值中误差均值中误差均值中误差ds2.16E-042.35E-042.25E-047.84E-052.67E-047.27E-05Dbeta3.01E-049.86E-052.00E-041.24E-04-4.25E-052.36E-04k19.40E-091.27E-095.34E-091.98E-10-2.61E-094.58E-11k2-2.51E-173.88E-161.81E-162.35E-17-7.62E-196.71E-18P15.54E-072.52E-074.03E-075.50E-082.28E-072.96E-08P2-4.37E-076.94E-08-1.84E-071.37E-07-2.57E-074.30E-07（a） ds、Dbeta随主距变化图（b） k1随主距变化图（c） P1、P2随焦距变化图图41 各项畸变参数随焦距变化图图41的（a）、（b）、（c）三幅图可以看出，随着焦距的变化，的值也是变化的。，随着焦距的增加而变大；，随着主距的变化而减小；由于参数的值太小，在主距发生改变时候，可以看作是不变的。从上面