一种GPS/INS数据直接定向精度评估方法.pdf

本发明公开了一种GPS/INS数据直接定向精度评估方法，其以摄影测量的重点方程共线方程为基础，以光束法多片前方交会为理论依据，结合最小二乘原理平差计算得到地面控制点物方坐标，再经过与地面控制点的实测坐标比较，统计得到坐标的误差值来评估GPS/INS系统所得外方位元素直接定向的精度。本发明评估精度高、量测误差小，评估结果能直接应用于生产实践，能有效避免人工立体量测中增加的像点坐标量测误差，并且其算法不受像片数约束。

1.  一种GPS/INS数据直接定向精度评估方法，其特征在于该方法包括以下步骤：
步骤一、通过航空摄影方法按照航空摄影规范摄取待测区域的多张摄影像片并通过像片量测方法获取各张像片中地面控制点的像点坐标观测数据，同时通过GPS/INS系统获取各张像片摄影时的外方位元素，所述外方位元素包括外方位角元素ω和κ以及摄影中心点坐标(X_S，Y_S，Z_S)；并通过常规摄影测量方法获取所述各张像片的内方位元素x₀、y₀和f；所述多张摄影像片的数量为地面每个控制点在多张像片上成像的像片数，设为N张；
采用常规实地测量方法，测量出待测区域内所有控制点的实测坐标数据且将所有实测坐标数据传送至中央处理器；
步骤二、将所获取的N张摄影像片对应的像点坐标观测数据、内外方位元素、每个地面控制点对应的摄影像片数N均传送至中央处理器进行分析处理求取待测区域内所有地面点的最佳物方坐标及其坐标理论误差，对于待测区域内一地面控制点M1即物点而言，所述中央处理器的分析处理过程如下：
第一步、获取地面控制点M1所对应的摄影像片数N；
第二步、利用N张摄影像片对应的外方位元素ω和κ空间方位元素，计算各张像片的旋转矩阵，第i张像片的旋转矩阵 R i = a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 , ]]>其中，i＝1、2、3...N；
第三步、求出N张像片中N个同名点对应的像点坐标计算值和误差方程式的常数项：
首先，利用共线方程公式列出N张像片中N个同名点对应的N组共线方程和N组误差方程式，对于第i张像片而言，其同名点对应的共线方程和误差方程式分别为 x - x 0 = - f a 1 ( X - X S ) + b 1 ( Y - Y S ) + c 1 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) y - y 0 = - f a 2 ( X - X S ) + b 2 ( Y - Y S ) + c 2 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) - - - ( 1 ) ]]>
和 v s = ( x ) - x + ∂ x ∂ X ΔX + ∂ x ∂ Y ΔY + ∂ x ∂ Z ΔZ v y = ( y ) - y + ∂ y ∂ X ΔX + ∂ y ∂ Y ΔY + ∂ y ∂ Z ΔZ - - - ( 2 ) , ]]>式(1)和(2)中(X，Y，Z)为地面控制点M1对应的未知坐标数据，(x，y)在第i张像片中为地面控制点M1对应像点的像点坐标观测数据；(v_x，v_y)为像点坐标数据改正数；
假定(X，Y，Z)的初值为(0，0，0)，将(0，0，0)作为地面控制点M1坐标的近似值代入式(1)，求出第i张像片中同名点对应的像点坐标计算值(x)和(y)： ( x ) = - f a 1 ( X - X S ) + b 1 ( Y - Y S ) + c 1 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ( y ) = - f a 2 ( X - X S ) + b 2 ( Y - Y S ) + c 2 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) , ]]>并相应计算出第i张像片中同名点对应的误差方程式的常数项l_x和l_y： l x = x - ( x ) l y = y - ( y ) ; ]]>
第四步、逐点列出N张像片中N个同名点对应的N个误差方程式：
对于第i张像片而言，其误差方程式V_i＝A_iQ-L_i，其中V_i＝[v_x v_y]^T，L_i＝[l_x l_y]^T， A i = a 11 a 12 a 13 a 21 a 22 a 23 , ]]>Q＝[ΔX，ΔY，ΔZ]；
A_i中的各偏导数为：
a 11 = ∂ x ∂ X = - a 1 f + a 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 12 = ∂ x ∂ Y = - b 1 f + b 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 13 = ∂ x ∂ Z = - c 1 f + c 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 21 = ∂ y ∂ X = - a 2 f + a 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 22 = ∂ y ∂ Y = - b 2 f + b 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 23 = ∂ y ∂ Z = - c 2 f + c 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ; ]]>
第五步、将逐点列出的N个误差方程式合成为一总误差方程式：V＝AQ-L(3)，其中Q＝[ΔX，ΔY，ΔZ]，A＝[A₁A₂...A_N]^T，V＝A＝[V₁V₂...V_N]^T，L＝[L₁L₂...L_N]^T；
第六步、根据最小二乘间接平差原理，列出式(3)的法方程式：A^TPAQ＝A^TPL(4)，式(4)中，P为单位矩阵；求出式(4)解的表达式为：Q＝(A^TA)^-1A^TL(4)，从而可求出地面控制点M1坐标近似值的改正数ΔX₁，ΔY₁和ΔZ₁；
第七步、判断改正数ΔX₁，ΔY₁和ΔZ₁是否小于限差，若小于限差，计算此时地面控制点M1改正后的坐标(X₁，Y₁，Z₁)，其中X₁＝0+ΔX₁，Y₁＝0+ΔY₁和Z₁＝0+ΔZ₁，坐标(X₁，Y₁，Z₁)即为地面控制点M1的最佳物方坐标，地面控制点M1的分析处理过程完毕，转入第九步；
若改正数ΔX₁，ΔY₁和ΔZ₁不小于限差，则转入第八步；
第八步、以第七步改正后的坐标作为地面控制点M1坐标的近似值代入式(1)，重复第三步至第七步，经多次迭代运算后，直至第六步中所计算得出的改正数ΔX_k，ΔY_k和ΔZ_k小于限差，此时地面控制点M1改正后的坐标为(X_k，Y_k，Z_k)，其中X_k＝0+ΔX₁+…+ΔX_k，Y_k＝0+ΔY₁+…+ΔY_k，Z_k＝0+ΔZ₁+…+ΔZ_k，k为迭代的次数；坐标(X_k，Y_k，Z_k)即为地面控制点M1的最佳物方坐标，地面控制点M1的分析处理过程完毕，转入第九步；
第九步、通过第六步中所述法方程式中未知数矩阵Q的系数矩阵的逆矩阵(A^TA)^-1来求解地面控制点M1坐标的理论精度即理论误差，此时视地面控制点M1的像点坐标观测数据为等精度不相关观测值，由于矩阵(A^TA)^-1中第f个主对角线上的元素Q_ff就是法方程式中第f个未知数的权倒数，则第f个未知数的中误差 m f = Q ff m 0 , ]]>其中f＝1、2和3，m₀为单位权中误差且 m 0 = ± [ vv ] 2 N - 3 , ]]>其中 [ vv ] = Σ i = 1 N l ix 2 + l iy 2 , ]]>则地面控制点M1坐标的理论误差m_xM1＝m₁，m_yM1＝m₂和m_hM1＝m₃；
第十步、将计算得出的地面控制点M1的最佳物方坐标与地面控制点M1的实测坐标数据进行作差比较，得出相应差值ΔX_M1，ΔY_M1和ΔZ_M1；
步骤三、重复步骤二，计算下一地面控制点M2的最佳物方坐标及其坐标理论误差，直至求出待测区域内所有地面控制点的最佳物方坐标及其坐标理论误差m_xMj，m_yMj和m_hMj；并且将待测区域内所有地面控制点的最佳物方坐标与各自的实测坐标数据进行作差比较，并得出相应多组差值ΔX_Mj，ΔY_Mj和ΔZ_Mj，其中j＝1、2...t，t为待测区域内所有地面控制点的数量；
步骤四、计算直接定向精度结果，即所有地面控制点坐标的实际中误差mx、my、ms、mh、实际最大误差ΔX_max，ΔY_max、ΔS_max和ΔH_max和理论中误差mlx、mly、mls、mlh：
首先，计算多组差值ΔX_Mj，ΔY_Mj和ΔZ_Mj的坐标值中误差 ΔX M 1 2 + · · · + Δ X Mt 2 t = mx , ]]> ΔY M 1 2 + · · · + Δ Y Mt 2 t = my , ]]>和 ΔZ M 1 2 + · · · + Δ Z Mt 2 t = mz = mh ; ]]>同时，相应计算得出平面中误差 ms = mx 2 + my 2 ; ]]>
从ΔX_Mj，ΔY_Mj中计算所有地面控制点的平面误差值 Δ S Mj = Δ x Mj 2 + Δ Y Mj 2 , ]]>从ΔX_Mj，ΔY_Mj，ΔS_Mj和ΔZ_Mj筛选出最大的坐标值差值即坐标值最大误差：ΔX_max，ΔY_max，ΔS_max和ΔZ_max＝ΔH_max；
从所有地面控制点坐标的理论误差，计算得出总的理论中误差 m 1 x = m xM 1 2 + · · · + m xMt 2 t , ]]> m 1 y = m yM 1 2 + · · · + m yMt 2 t , ]]> m 1 h = m hM 1 2 + · · · + m hMt 2 t , ]]> m 1 s = mlx 2 + mly 2 . ]]>

2.  按照权利要求1所述的一种GPS/INS数据直接定向精度评估方法，其特征在于：步骤一中所述的常规实地测量方法为常规野外控制点测量方法。

一种GPS/INS数据直接定向精度评估方法
技术领域
本发明涉及摄影测量技术领域，尤其是涉及一种GPS/INS数据直接定向精度评估方法。
背景技术
获取影像的外方位元素是摄影测量的关键问题之一，很多摄影测量工作者一直为解决如何快速高效地获取影像的外方位元素的问题而努力着。长期以来，人们主要用空中三角测量技术，通过大量的外业控制点来间接解求外方位元素，上述方法既浪费经费，又延长了成图周期，严重影响国民经济建设所急需的资料。20世纪90年代以来，国际摄影测量界开始研究采用GPS/INS组合导航系统(即GPS/INS系统：GPS系统和INS惯性导航系统的组合系统)来获取像片摄影时的空间方位，以完全取代摄影测量加密，从而实现无地面控制的空中三角测量航空遥感已成为国内航空摄影测量的一种趋势。
而如何评估GPS/INS系统所得到的外方位元素精度成为目前迫切需要解决的问题，常用的评估方法有两种：一种是将GPS/INS系统所得的外方位元素数据即GPS/INS数据导入立体测图仪建立立体模型，人工立体观测获取地面检查点的物方坐标，再与实测检查点坐标对比来评估；另一种是利用高精度量测的检查点的像点坐标，经点投影系数法两片前交计算得到检查点的物方坐标，再与实测检查点坐标对比来评估。
第一种方法从常规立体测图入手，也是生产中常容易想到的评估外方位元素的方法。立体测图仪是成熟的用于数字摄影测量的系统，将外方位元素数据导入立体测图仪，即可在两张像片组成的立体模型中，人工观测地面控制点而获取地面控制点物方坐标。但是这种方法存在两个不足：1、由于GPS/INS系统获取的外方位元素存在误差，在立体测图仪中存在很大视差，很难进行立体测图，在视差的影响下，就不能正确量测到所需的地面控制点，导致像点坐标的量测存在较大误差，这样量测的地面控制点坐标的误差中有很大一部分来自像点坐标误差的影响，不能严密地评估外方位元素的数据质量。2、地面一个检查点会出现在多个立体像对中，需要人工逐个模型地观测，再将多个模型观测得到的坐标值平差得到最后的坐标值，同实测点坐标对比，才能比较客观地评估外方位元素质量，这是一个非常费工费时的工作，不仅浪费经费，而且也增加了质量评估的时间。正是因为上述问题的存在，导致本方法的弱点不可避免，很难很好地应用到实际的GIS/INS数据评估工作中。从实际用立体测图仪人工量测的GPS/INS数据直接定向的精度分析表可看出，由于人工成分太多，在存在视差的条件下，量测的像点坐标误差太多，无法正确评估直接定向精度，分析所得数据无法用于评估工作。
另一种方法是从摄影测量理论入手，进行摄影测量立体像对空间前方交会中的点投影系数法空间前方交会。用空三加密所得到的高精度的量测像点坐标，联合GPS/INS系统所得的外方位元素数据，代入点投影系数的空间前方交会公式，求出地面控制点坐标。空间前方交会公式如下：
X＝X_S1+NX₁＝X_S1+B_X+N′X₂
Y = Y S 1 + NY 1 = Y S 1 + B Y + N ′ Y 2 = Y S 1 + 1 2 ( NY 1 + N ′ Y 2 + B Y ) ]]>
Z＝Z_S1+NZ₁＝Z_S1+B_Z+N′Z₂                     式(1)
详细内容可参考摄影测量学教材中立体像对空间前方交会一节。这种方法有效避免了量测像点坐标的误差，理论上应该符合数学模型。但存在以下两种缺点：1、经过实验证明：这种方法不能符合实际摄影测量生产情况，只是数学理想化的双线前交，实际加密生产中的光束法平差是实现多个同名光束最佳地交会到一起，并不能真是交会到一起，只有由同一个地面控制点和左右片的外方位元素通过共线方程计算出来的理论的像点坐标数据，再和左右片的外方位元素数据计算前方交会才符合该方法的数学模型。即使使用高精度的空三加密成果数据，通过该方法前方交会计算，得到的地面控制点物方坐标的误差也很大，尤其是高程误差大，基本不符合实际情况，无法反映外方位元素数据的质量，所以不能正确评估GPS/INS的数据质量，算法不够严密。2、点投影系数法只限于两片前方交会，不能实现多光线的交会，受影像数约束。正是由于该方法的以上不足，所以只能作为理论的前方交会模型方法，不能真正用于实际的GPS/INS数据直接定向精度评估工作。据实际分析得出的常规空中三角测量所得的像片外方位元素定向的精度统计结果，其空三精度足以满足丘陵地区的规范要求；但用常规空中三角测量所得的外方位元素数据和像点坐标数据，通过点投影系数法两片前方交会所得的精度统计结果，其误差大大超过了规范要求，甚至不合常理，可见，用点投影系数法两片前方交会的方法不能正确评估外方位元素数据质量，不能用于GPS/INS所得外方位元素直接定向的精度评估工作。
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种GPS/INS数据直接定向精度评估方法，其评估精度高、量测误差小，评估结果能直接应用于生产实践，能有效避免人工立体量测中增加的像点坐标量测误差，并且其算法不受像片数约束。
为解决上述技术问题，本发明采用的技术方案是：一种GPS/INS数据直接定向精度评估方法，其特征在于该方法包括以下步骤：
步骤一、通过航空摄影方法按照航空摄影规范摄取待测区域的多张摄影像片并通过像片量测方法获取各张像片中地面控制点的像点坐标观测数据，同时通过GPS/INS系统获取各张像片摄影时的外方位元素，所述外方位元素包括外方位角元素、ω和κ以及摄影中心点坐标(X_S，Y_S，Z_S)；并通过常规摄影测量方法获取所述各张像片的内方位元素x₀、y₀和f；所述多张摄影像片的数量为地面每个控制点在多张像片上成像的像片数，设为N张；
采用常规实地测量方法，测量出待测区域内所有控制点的实测坐标数据且将所有实测坐标数据传送至中央处理器；
步骤二、将所获取的N张摄影像片对应的像点坐标观测数据、内外方位元素、每个地面控制点对应的摄影像片数N均传送至中央处理器进行分析处理求取待测区域内所有地面点的最佳物方坐标及其坐标理论误差，对于待测区域内一地面控制点M1即物点而言，所述中央处理器的分析处理过程如下：
第一步、获取地面控制点M1所对应的摄影像片数N；
第二步、利用N张摄影像片对应的外方位元素ω和κ空间方位元素，计算各张像片的旋转矩阵，第i张像片的旋转矩阵 R i = a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 , ]]>其中，i＝1、2、3...N；
第三步、求出N张像片中N个同名点对应的像点坐标计算值和误差方程式的常数项：
首先，利用共线方程公式列出N张像片中N个同名点对应的N组共线方程和N组误差方程式(因此，对于任何地面控制点而言，对应N(N≥2)个同名像点可列出2N(2n≥4)个共线方程和误差方程，但未知数只有三个，所以有2N-3个多余观测)，对于第i张像片而言，其同名点对应的共线方程和误差方程式分别为 x - x 0 = - f a 1 ( X - X S ) + b 1 ( Y - Y S ) + c 1 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) y - y 0 = - f a 2 ( X - X S ) + b 2 ( Y - Y S ) + c 2 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>(1)
和 v x = ( x ) - x + ∂ x ∂ X ΔX + ∂ x ∂ Y ΔY + ∂ x ∂ Z ΔZ v y = ( y ) - y + ∂ y ∂ X ΔX + ∂ y ∂ Y ΔY + ∂ y ∂ Z ΔZ - - - ( 2 ) , ]]>式(1)和(2)中(X，Y，Z)为地面控制点M1对应的未知坐标数据，(x，y)在第i张像片中为地面控制点M1对应像点的像点坐标观测数据；(v_x，v_y)为像点坐标数据改正数；
假定(X，Y，Z)的初值为(0，0，0)，将(0，0，0)作为地面控制点M1坐标的近似值代入式(1)，求出第i张像片中同名点对应的像点坐标计算值(x)和(y)： ( x ) = - f a 1 ( X - X S ) + b 1 ( Y - Y S ) + c 1 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ( y ) = - f a 2 ( X - X S ) + b 2 ( Y - Y S ) + c 2 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) , ]]>并相应计算出第i张像片中同名点对应的误差方程式的常数项l_x和l_y： l x = x - ( x ) l y = y - ( y ) ; ]]>
第四步、逐点列出N张像片中N个同名点对应的N个误差方程式：
对于第i张像片而言，其误差方程式V_i＝A_iQ-L_i，其中V_i＝[v_x v_y]^T，L_i＝[l_x l_y]^T， A i = a 11 a 12 a 13 a 21 a 22 a 23 , ]]>Q＝[ΔX，ΔY，Δz]；
A_i中的各偏导数为：
a 11 = ∂ x ∂ X = - a 1 f + a 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 12 = ∂ x ∂ Y = - b 1 f + b 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 13 = ∂ x ∂ Z = - c 1 f + c 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 21 = ∂ y ∂ X = - a 2 f + a 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 22 = ∂ y ∂ Y = - b 2 f + b 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 23 = ∂ y ∂ Z = - c 2 f + c 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ; ]]>
第五步、将逐点列出的N个误差方程式合成为一总误差方程式：V＝AQ-L  (3)，其中Q＝[ΔX，ΔY，ΔZ]，A＝[A₁A₂...A_N]^T，V＝A＝[V₁V₂...V_N]^T，L＝[L₁L₂...L_N]^T；
第六步、根据最小二乘间接平差原理，列出式(3)的法方程式：A^TPAQ＝A^TPL  (4)，式(4)中，P为单位矩阵；求出式(4)解的表达式为：Q＝(A^TA)^-1A^TL  (4)，从而可求出地面控制点M1坐标近似值的改正数ΔX₁，ΔY₁和ΔZ₁；
第七步、判断改正数ΔX₁，ΔY₁和ΔZ₁是否小于限差，若小于限差，计算此时地面控制点M1改正后的坐标(X₁，Y₁，Z₁)，其中X₁＝0+ΔX₁，Y₁＝0+ΔY₁和Z₁＝0+ΔZ₁，坐标(X₁，Y₁，Z₁)即为地面控制点M1的最佳物方坐标，地面控制点M1的分析处理过程完毕，转入第九步；
若改正数ΔX₁，ΔY₁和ΔZ₁不小于限差，则转入第八步；
第八步、以第七步改正后的坐标作为地面控制点M1坐标的近似值代入式(1)，重复第三步至第七步，经多次迭代运算后，直至第六步中所计算得出的改正数ΔX_k，ΔY_k和ΔZ_k小于限差，此时地面控制点M1改正后的坐标为(X_k，Y_k，Z_k)，其中X_k＝0+ΔX₁+...+ΔX_k，Y_k＝0+ΔY₁+...+ΔY_k，Z_k＝0+ΔZ₁+...+ΔZ_k，k为迭代的次数；坐标(X_k，Y_k，Z_k)即为地面控制点M1的最佳物方坐标，地面控制点M1的分析处理过程完毕，转入第九步；
第九步、通过第六步中所述法方程式中未知数矩阵Q的系数矩阵的逆矩阵(A^TA)^-1来求解地面控制点M1坐标的理论精度即理论误差，此时视地面控制点M1的像点坐标观测数据为等精度不相关观测值，由于矩阵(A^TA)^-1中第f个主对角线上的元素Q_ff就是法方程式中第f个未知数的权倒数，则第f个未知数的中误差 m f = Q ff m 0 , ]]>其中f＝1、2和3，m₀为单位权中误差且 m 0 = ± [ vv ] 2 N - 3 , ]]>其中 [ vv ] = Σ i = 1 N l ix 2 + l iy 2 , ]]>则地面控制点M1坐标的理论误差m_xM1＝m₁，m_yM1＝m₂和m_hM1＝m₃；
第十步、将计算得出的地面控制点M1的最佳物方坐标与地面控制点M1的实测坐标数据进行作差比较，得出相应差值ΔX_M1，ΔY_M1和ΔZ_M1；
步骤三、重复步骤二，计算下一地面控制点M2的最佳物方坐标及其坐标理论误差，直至求出待测区域内所有地面控制点的最佳物方坐标及其坐标理论误差m_xMj，m_yMj和m_hMj；并且将待测区域内所有地面控制点的最佳物方坐标与各自的实测坐标数据进行作差比较，并得出相应多组差值ΔX_Mj，ΔY_Mj和ΔZ_Mj，其中j＝1、2...t，t为待测区域内所有地面控制点的数量；
步骤四、计算直接定向精度结果，即所有地面控制点坐标的实际中误差mx、my、ms、mh、实际最大误差ΔX_max，ΔY_max、ΔS_max和ΔH_max和理论中误差mlx、mly、mls、mlh：
首先，计算多组差值ΔX_Mj，ΔY_Mj和ΔZ_Mj的坐标值中误差 ΔX M 1 2 + · · · + ΔX Mt 2 t = mx , ]]> ΔY M 1 2 + · · · + ΔY Mt 2 t = my ]]>和 ΔZ M 1 2 + · · · + ΔZ Mt 2 t = mz = mh ; ]]>同时，相应计算得出平面中误差 ms = mx 2 + my 2 ; ]]>
从ΔX_Mj，ΔY_Mj中计算所有地面控制点的平面误差值 ΔS Mj = ΔX Mj 2 + ΔY Mj 2 , ]]>从ΔX_Mj，ΔY_Mj，ΔS_Mj和ΔZ_Mj筛选出最大的坐标值差值即坐标值最大误差：ΔX_max，ΔY_max，ΔS_max和ΔZ_max＝ΔH_max；
从所有地面控制点坐标的理论误差，计算得出总的理论中误差 mlx = m xM 1 2 + · · · + m xMt 2 t , ]]> mly = m yM 1 2 + · · · + m yMt 2 t , ]]> mlh = m hM 1 2 + · · · + m hMt 2 t , ]]> mls = mlx 2 + mly 2 . ]]>
步骤一中所述的常规实地测量方法为常规野外控制点测量方法。
本发明与现有技术相比具有以下优点：1、设计新颖，本发明以摄影测量的重点方程共线方程为基础，以光束法多片前方交会为理论依据，结合最小二乘原理平差计算得到地面控制点物方坐标，再经过与地面控制点的实测坐标比较，统计得到坐标的误差值来评估GPS/INS系统所得外方位元素直接定向的精度。2、自动化程度高，不需要人工参与量测；而传统立体测图仪人工量测方法自动化程度低，人工参与成分大，手工工作量大，且手工工作难度较大，效率较低。3、本发明采用光束法原理进行计算处理，算法严密，能正确反映外方位元素质量；传统点投影系数法两片前交方法算法不严密，不能正确反映外方位元素质量。4、本发明采用多片最小二乘计算，即采用多像空间前方交会方法计算得出地面控制点的物方坐标，因而不受像片数约束；传统立体测图仪人工量测方法为两片立体模型量测，因而受像片数约束；传统点投影系数法两片前交方法为两片前方交会，因而也受像片数约束。5、本发明用加密的像点坐标，能有效避免像点坐标误差问题；传统立体测图仪人工量测方法存在视差，不能有效避免像点坐标误差。综上所述，本发明从摄影测量的核心方程共线方程入手，结合最小二乘原理，利用一点多度重叠的多余观测条件，实现地面控制点最佳值的计算，完成GPS/INS数据直接定向精度的评估，其不仅能够避免人工立体量测中增加的像点坐标量测误差，而且能充分利用一个地面控制点的多度重叠的条件和数据，增加平差条件，其算法不受像片数约束；同时，其也能结合实际生产现实，实现最小二乘光束法平差，算法理论严密、计算出地面控制点坐标的最佳值，能真正实现GPS/INS数据直接定向精度的评估。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明的流程框图。
具体实施方式
如图1所示，本发明所述的GPS/INS数据直接定向精度评估方法，包括以下步骤：
步骤一、通过航空摄影方法按照航空摄影规范摄取待测区域的多张摄影像片并通过像片量测方法获取各张像片中地面控制点的像点坐标观测数据，同时通过GPS/INS系统获取各张像片摄影时的外方位元素，所述外方位元素包括外方位角元素、ω和κ以及摄影中心点坐标(X_S，Y_S，Z_S)；并通过常规摄影测量方法获取所述各张像片的内方位元素x₀、y₀和f；所述多张摄影像片的数量为地面每个控制点在多张像片上成像的像片数，设为N张。
采用常规实地测量方法，测量出待测区域内所有控制点的实测坐标数据且将所有实测坐标数据传送至中央处理器。
本实施例中，所述常规实地测量方法为常规野外控制点测量方法。求取内方位元素x₀、y₀和f的方法为相机检校和内定向等方法。
步骤二、将所获取的N张摄影像片对应的像点坐标观测数据、内外方位元素、每个地面控制点对应的摄影像片数N均传送至中央处理器进行分析处理求取待测区域内所有地面点的最佳物方坐标及其坐标理论误差，对于待测区域内一地面控制点M1即物点而言，所述中央处理器的分析处理过程如下：
第一步、获取地面控制点M1所对应的摄影像片数N；
第二步、利用N张摄影像片对应的外方位元素ω和κ空间方位元素，计算各张像片的旋转矩阵，第i张像片的旋转矩阵 R i = a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 , ]]>其中，i＝1、2、3...N；
第三步、求出N张像片中N个同名点对应的像点坐标计算值和误差方程式的常数项：
首先，利用共线方程公式列出N张像片中N个同名点对应的N组共线方程和N组误差方程式，对于第i张像片而言，其同名点对应的共线方程和误差方程式分别为 x - x 0 = - f a 1 ( X - X S ) + b 1 ( Y - Y S ) + c 1 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) y - y 0 = - f a 2 ( X - X S ) + b 2 ( Y - Y S ) + c 2 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) - - - ( 1 ) ]]>和 v x = ( x ) - x + ∂ x ∂ X ΔX + ∂ x ∂ Y ΔY + ∂ x ∂ Z ΔZ v y = ( y ) - y + ∂ y ∂ X ΔX + ∂ y ∂ Y ΔY + ∂ y ∂ Z ΔZ - - - ( 2 ) , ]]>式(1)和(2)中(X，Y，Z)为地面控制点M1对应的未知坐标数据，(x，y)在第i张像片中为地面控制点M1对应像点的像点坐标观测数据；(v_x，v_y)为像点坐标数据改正数；
假定(X，Y，z)的初值为(0，0，0)，将(0，0，0)作为地面控制点M1坐标的近似值代入式(1)，求出第i张像片中同名点对应的像点坐标计算值(x)和(y)： ( x ) = - f a 1 ( X - X S ) + b 1 ( Y - Y S ) + c 1 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ( y ) = - f a 2 ( X - X S ) + b 2 ( Y - Y S ) + c 2 ( Z - Z S ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) , ]]>并相应计算出第i张像片中同名点对应的误差方程式的常数项l_x和l_y： l x = x - ( x ) l y = y - ( y ) ; ]]>
第四步、逐点列出N张像片中N个同名点对应的N个误差方程式：
对于第i张像片而言，其误差方程式V_i＝A_iQ-L_i，其中V_i＝[v_x v_y]^T，L_i＝[l_x l_y]^T， A i = a 11 a 12 a 13 a 21 a 22 a 23 , ]]>Q＝[ΔX，ΔY，Δz]；
A_i中的各偏导数为：
a 11 = ∂ x ∂ X = - a 1 f + a 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 12 = ∂ x ∂ Y = - b 1 f + b 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 13 = ∂ x ∂ Z = - c 1 f + c 3 ( x - x 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 21 = ∂ y ∂ X = - a 2 f + a 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 22 = ∂ y ∂ Y = - b 2 f + b 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ]]>
a 23 = ∂ y ∂ Z = - c 2 f + c 3 ( y - y 0 ) a 3 ( X - X S ) + b 3 ( Y - Y S ) + c 3 ( Z - Z S ) ; ]]>
第五步、将逐点列出的N个误差方程式合成为一总误差方程式：V＝AQ-L  (3)，其中Q＝[ΔX，ΔY，ΔZ]，A＝[A₁ A₂...A_N]^T，V＝A＝[V₁ V₂...V_N]^T，L＝[L₁ L₂...L_N]^T；
第六步、根据最小二乘间接平差原理，列出式(3)的法方程式：A^TPAQ＝A^TPL  (4)，式(4)中，P为单位矩阵；求出式(4)解的表达式为：Q＝(A^TA)^-1A^TL  (4)，从而可求出地面控制点M1坐标近似值的改正数ΔX₁，ΔY₁和ΔZ₁；
第七步、判断改正数ΔX₁，ΔY₁和ΔZ₁是否小于限差，若小于限差，计算此时地面控制点M1改正后的坐标(X₁，Y₁，Z₁)，其中X₁＝0+ΔX₁，Y₁＝0+Δy₁和Z₁＝0+ΔZ₁，坐标(X₁，Y₁，Z₁)即为地面控制点M1的最佳物方坐标，地面控制点M1的分析处理过程完毕，转入第九步；
若改正数ΔX₁，Δy₁和ΔZ₁不小于限差，则转入第八步；
第八步、以第七步改正后的坐标作为地面控制点M1坐标的近似值代入式(1)，重复第三步至第七步，经多次迭代运算后，直至第六步中所计算得出的改正数ΔX_k，ΔY_k和ΔZ_k小于限差，此时地面控制点M1改正后的坐标为(X_k，Y_k，Z_k)，其中X_k＝0+ΔX₁+...+ΔX_k，Y_k＝0+ΔY₁+...+ΔY_k，Z_k＝0+ΔZ₁+...+ΔZ_k，k为迭代的次数；坐标(X_k，Y_k，Z_k)即为地面控制点M1的最佳物方坐标，地面控制点M1的分析处理过程完毕，转入第九步；
第九步、通过第六步中所述法方程式中未知数矩阵Q的系数矩阵的逆矩阵(ATA)-¹来求解地面控制点M1坐标的理论精度即理论误差，此时视地面控制点M1的像点坐标观测数据为等精度不相关观测值，由于矩阵(A^TA)^-1中第f个主对角线上的元素Q_ff就是法方程式中第f个未知数的权倒数，则第f个未知数的中误差 m f = Q ff m 0 , ]]>其中f＝1、2和3，m₀为单位权中误差且 m 0 = ± [ vv ] 2 N - 3 , ]]>其中 [ vv ] = Σ i = 1 N l ix 2 + l iy 2 , ]]>则地面控制点M1坐标的理论误差m_xM1＝m₁，m_yM1＝m₂和m_hM1＝m₃；
同时，需注意的是：此处的l_ix和l_iy为第i张像片经k次迭代后的误差方程式的常数项，具体是第三步中计算得出的第i张像片中同名点对应的误差方程式的常数项l_x和l_y。
第十步、将计算得出的地面控制点M1的最佳物方坐标与地面控制点M1的实测坐标数据进行作差比较，得出相应差值Δ_XM1，Δ_YM1和Δ_ZM1。
步骤三、重复步骤二，计算下一地面控制点M2的最佳物方坐标及其坐标理论误差，直至求出待测区域内所有地面控制点的最佳物方坐标及其坐标理论误差m_xMj，m_yMj和m_hMj；并且将待测区域内所有地面控制点的最佳物方坐标与各自的实测坐标数据进行作差比较，并得出相应多组差值ΔX_Mj，ΔY_Mj和ΔZ_Mj，其中j＝1、2...t，t为待测区域内所有地面控制点的数量。
步骤四、计算直接定向精度结果，即所有地面控制点坐标的实际中误差mx、my、ms、mh、实际最大误差ΔX_max，ΔY_max、ΔS_max和ΔH_max和理论中误差mlx、mly、mls、mlh：
首先，计算多组差值ΔX_Mj，ΔY_Mj和ΔZ_Mj的坐标值中误差 ΔX M 1 2 + · · · + ΔX Mt 2 t = mx , ]]> ΔY M 1 2 + · · · + ΔY Mt 2 t = my ]]>和 ΔZ M 1 2 + · · · + ΔZ Mt 2 t = mz = mh ; ]]>同时，相应计算得出平面中误差 ms = mx 2 + my 2 ; ]]>
从ΔX_Mj，ΔY_Mj中计算所有地面控制点的平面误差值 ΔS Mj = ΔX Mj 2 + ΔY Mj 2 , ]]>从ΔX_Mj，ΔY_Mj，ΔS_Mj和ΔZ_Mj筛选出最大的坐标值差值即坐标值最大误差：ΔX_max，ΔY_max，ΔS_max和ΔZ_max＝ΔH_max；
从所有地面控制点坐标的理论误差，计算得出总的理论中误差 mlx = m xM 1 2 + · · · + m xMt 2 t , ]]> mly = m yM 1 2 + · · · + m yMt 2 t , ]]> mlh = m hM 1 2 + · · · + m hMt 2 t , ]]> mls = mlx 2 + mly 2 . ]]>
综上，要进行GPS/INS系统的数据直接定向的精度评估工作，就要将GPS/INS数据直接定向的结果和常规空中三角测量的定向结果进行比较，所以先进行常规的空中三角测量工作并且达到高精度要求，除了得到加密精度成果之外，从常规空中三角测量中，可以得到高精度的地面控制点的量测像点坐标，首先保证了前方交会中像点坐标的量测精度。
因而当将本发明用于GPS/INS应用的实验研究中时，先由GPS/INS组合导航系统来获取像片摄影时的空间方位，即像片的外方位元素X_S、Y_S、Z_S、ω、κ。因精度对比所需，在试验区进行了常规摄影测量即常规空中三角测量得到了一套满足精度要求的外方位元素数据和控制点像点坐标数据。从外方位元素数据和像点坐标数据这些已知数据出发，经上述步骤即步骤一到步骤四(在Microsoft Visual C++6.0平台上用C++语言开发具体实用程序)，可得到地面控制点的物方坐标值。
用本发明统计的常规空中三角测量所得外方位元素定向精度，如表1所示：
表1      本发明统计的常规空中三角测量成果直接定向精度

利用常规空中三角测量系统统计的常规空中三角测量所得的像片外方位元素定向的精度统计结果如表2所示：
表2

将表1和表2对比分析可知，表1与表2的统计结果基本一致，可见本发明能够正确客观地评估外方位元素的定向精度。
选用一个测区的574个控制点的地面控制点坐标，与实测地面控制点数据对比，即用本发明统计GPS/INS系统所得外方位元素直接定向精度结果如表3所示：
表3        本发明统计的GPS/INS数据直接定向精度

从表3可以看出，本发明能够完成GPS/INS数据的直接定向精度评估工作。
采用现有的立体测图仪人工量测的GPS/INS直接定向的精度分析结果，和通过点投影系数法两片前方交会统计的常规空中三角测量所得的像片外方位元素定向的精度结果分别如表4和表5所示：
表4：立体测图仪统计GPS/INS数据直接定向精度结果

表5：点投影系数法统计的常规空中三角测量成果直接定向精度结果

将表4和表5分别与表2对比分析可知，现有的立体测图仪人工量测和通过点投影系数法两片前方交会方法均不能正确客观地评估外方位元素的定向精度。
以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。