CN103345737A - 一种基于误差补偿的uav高分辨率影像几何校正方法 - Google Patents

Abstract

本发明公开了一种基于误差补偿的UAV高分辨率影像几何校正方法，包括以下步骤：第一步，确定机载相机属性参数，计算内参数矩阵A_in；第二步，获取遥感影像外方位参数，计算外参数矩阵A_out；第三步，基于地面控制点计算误差补偿矩阵A_{error_cmp}；第四步，整合校正矩阵；第五步，基于全局变换矩阵A，对遥感影像进行几何校正；第六步，遥感影像边界地理坐标定位；本发明综合考虑了引起无人机遥感影像的几何畸变因素，并通过误差补偿机制有效提高了本发明方法精度；本发明方法基本实现全程自动化，某些参数值一旦计算出来就可重复利用，大大提高了工作效率；本发明方法通过引入高度标定参数，可根据需要获取不同分辨率的几何校正图像，提高了本发明方法的灵活性。

Description

一种基于误差补偿的UAV高分辨率影像几何校正方法

技术领域

本发明属于数字图像处理技术领域，具体涉及一种基于误差补偿的无人机（UAV）高分辨率遥感影像几何校正方法。

背景技术

无人机因其低空遥感数据具有高分辨率、高灵活性、高效率和低成本的优势而被广泛应用于自然灾害区域评估、战场侦察、环境监测等领域。然而，相对于无人机的快速发展，无人机遥感影像处理技术相对落后；与传统航天影像和航空影像相比，无人机遥感影像数据倾角大且无规律，并且无人机应用于资源勘察、灾害监测时，难以建立野外实测地面控制点（GCP）。这些特点都给无人机遥感影像，尤其是高分辨率遥感影像的精确几何校正带来了困难，进而影响到后期诸如目标定位、图像镶嵌等信息的有效提取。

影响无人机遥感影像几何畸变的因素很多，可简单分为两大类：

（1）外方位因素

遥感相机被置于无人机机载平台上，从不同的高度对地进行观测。这时由于大气环境、地球起伏、平台状况等多种因素的影响，使得获取的遥感图像存在种种扭曲变形。产生无人机遥感影像几何形变的外方位因素主要有3个：

1、地形起伏的影响：

无人机遥感多属中低空遥感，飞行高度与地形起伏有一定的关系，而且视场角较大，所以获取的图像不可避免的存在因地形起伏所造成的投影差。

2、大气折射的影响：

电磁波在大气层中传播会受到不同高度的大气密度不同所引起的折射影响，使其传播路径发生偏折，在成像时造成像点位移，在一定程度上引起遥感影像的几何失真。

3：遥感器（平台）姿态的影响

由于无人机通常在中低空飞行，易受大气环境影响，机体稳定性差，即使应用稳定传感器平台也不能保证遥感器的绝对平稳。分辨率越高，遥感影像畸变效果越明显。飞机姿态的影响体现在：侧滚角和俯仰角导致成像中心偏移和地物扭曲变形；偏航角导致遥感影像整体旋转；航高变化产生图像缩放。

（2）内方位因素

相机属性决定了遥感器成像模型，不同的成像模型成像原理有着本质区别。另外，当遥感器自身性能、技术指标偏离标称值时，就会造成遥感影像固有的几何变形，这些固有的畸变因素统称为内方位因素。

另外，由于测量误差、安装定位不准确等各方面因素影响，都会对在不同程度上产生几何畸变。传统的无人机遥感影像几何校正多采用基于地面控制点的方法，该方法准确度高，但需要人工采集地面控制点，工作量大，费时费力，而且在有些地方人工无法采集到控制点，如海平面、地震后的危险地带等。另外，如果有地图作为基准匹配，也可迅速建立基于控制点的遥感影像几何校正。然而这些方法都需要过多的先验知识及人工干预，实时性低，操作不方便，难以满足无人机地面控制站的实时性要求。

发明内容

本发明的目的是为了解决上述传统无人机高分辨率遥感影像几何校正方法的不足，提出了一种基于误差补偿的高分辨率无人机遥感影像的几何校正方法。实验结果表明，本发明的方法在实时性和精度方面都得到了明显提高。

本发明的一种基于误差补偿的UAV高分辨率影像几何校正方法，包括以下几个步骤：

第一步，确定相机属性，计算内参数矩阵A_in；

相机属性参数包括：相机焦距f、像元尺寸u、图像尺寸（Weight×Height）。

相机内参数矩阵A_in是由相机定标得来，对每一个相机而言，在焦距保持恒定的情况下内参数矩阵是固定的，不依赖场景的视图，一旦计算出来，就可以重复使用。

第二步，获取遥感影像外方位参数，计算外参数矩阵A_out；

每一幅遥感影像都有唯一的外方位参数，通常由专门的惯性导航系统获取得到，且与遥感影像获取时间同步。本发明需要的遥感影像外方位参数包括：飞行高度H、俯仰角

滚转角ω、偏航角κ、成像经纬度（Latitude，Longitude）。综合考虑外方位因素对遥感影像几何失真的影像，计算外参数矩阵A_out，为后续遥感影像几何校正做准备。

第三步，基于地面控制点计算误差补偿矩阵A_{error_cmp}。

基于A_in和A_out校正后的遥感影像与真实地图匹配时依然有较大差距，为提高几何校正精度，通过数据拟合计算误差补偿矩阵A_{error_cmp}，消除各测量误差对遥感影像几何失真的影响。

第四步，综合考虑各因素影响，整合校正矩阵。

综合考虑引入遥感影像引入几何失真的先后顺序，整合A_in、A_out和A_{error_cmp}得到全局变换矩阵A，为后续一步几何校正算法做准备。

第五步，基于全局变换矩阵A，对遥感影像进行几何校正。

第六步，遥感影像边界地理坐标定位；

基于相机属性、外方位因素和误差补偿矩阵，完成成像经纬度（Latitude，Longitude）在图像坐标中的定位问题。在中心定位和几何校正后图像的基础上，计算遥感影像边界的地理坐标，即东经、西经、北纬和南纬。

本发明具有以下优点：

1)综合考虑了引起无人机遥感影像的几何畸变因素，并通过误差补偿机制有效提高了本发明方法精度；

2)本发明方法基本实现全程自动化，某些参数值一旦计算出来就可重复利用，大大提高了工作效率；

3)本发明方法通过引入高度标定参数，可根据需要获取不同分辨率的几何校正图像，提高了本发明方法的灵活性。

4)本发明方法综合考虑了引入误差的先后顺序，几何校正矩阵一步完成，大大提高了运算效率高，适用于对实时性要求高的场合。

附图说明

图1为本发明方法总体流程图；

图2为无人机原始遥感影像；

图3为中心投影构像原理示意图；

图4为正射投影示意图；

图5为地形起伏对遥感影像的影响；

图6为大气折射对遥感影像的影响；

图7未经误差补偿的几何校正图像；

图8经误差补偿后的几何校正图像；

图9为图像边界地理坐标。

具体实施方式

下面结合附图对本发明的具体实施方法进行详细说明。

本发明以针孔摄像机模型作为研究对象，综合考虑引起无人机遥感影像成像畸变因素，确定了以全局透视变换为主的遥感图像几何校正解决方案，并解决了畸变图像几何校正后的地理坐标定位问题。

本发明是一种针对高分辨率无人机遥感影像的几何校正方法，在综合考虑内方位因素和外方位因素的基础上，通过引入误差补偿机制进行系统级几何校正，大大提高了算法校正精度。另外，本发明通过矩阵整合，实现了几何校正一步到位，有效缩短了运算时间。最后实现了校正后图像的地理边界定位问题。

本发明的一种基于误差补偿的UAV高分辨率影像几何校正方法，总体流程图如图1所示，图2为无人机原始遥感影像，图像尺寸为4096*3072，成像高度为1939.91米，成像位置：纬度：37.4396212、经度：118.5698322，成像姿态（单位：弧度）：俯仰角：-0.0196686、滚转角：-0.0128441、偏航角：-0.4469624，具体实施方法包括以下步骤：

第一步，确定相机属性，计算内参数矩阵A_in；

机载相机属性参数很多，本发明仅需要如下3个参数：

相机焦距：通常以毫米mm作为单位，记为f；

像元尺寸：通常情况下横向像元尺寸和径向像元尺寸相等，单位为微米um，记为u；

图像尺寸：图像尺寸为图像的像素尺寸，即图像的宽Width和高Height，无单位。

相机内参数矩阵A_in不依赖场景的视图，在焦距固定的情况下一旦被计算出来，就可以重复利用。获取内参数矩阵有专门的摄像机定标算法，不属于本发明研究内容。内参数矩阵如下式（1）所示：

$A_{\mathrm{in}}=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]---\left(1\right)$

其中，x、y代表像素坐标系方向，(c_x,c_y)是基准点（通常在图像中心）；f_x，f_y是以像素为单位的焦距。如果因为某些因素对来自于摄像机的一幅图像进行升采样或将采样，所有这些参数（c_x,c_y,f_x,f_y）都将被缩放（乘或除）同样尺度。

第二步，获取遥感影像外方位参数，计算外参数矩阵A_out；

每一幅遥感影像都有唯一的外方位参数，通常由专门的惯性导航系统获取得到，且与遥感影像获取时间同步。本发明需要的遥感影像外方位参数包括：飞行高度H、俯仰角

、滚转角ω、偏航角κ、成像经纬度（Latitude，Longitude）。当标定成像高度H′时，可直接计算得到A_out。下面就方法原理予以说明：

当机载相机满足中心投影构像原理时，即成像时各光线汇聚于物镜中心S(即摄影中心)，遥感影像属于中心投影瞬间一次成像，同一幅影像上的所有像点共用一个摄影中心和一个像片平面，共用一组外方位元素，如图3所示，其中S为物镜中心，f为相机焦距，H为成像高度，a、b为地物点A、B在正负成像平面上的投影。

一幅影像的外方位因素包括六个参数：三个线参数和三个角参数。当三个角参数呈（分别为俯仰角，滚转角和偏航角）转角系统时，可求得三个角元素

、ω、κ的方向余弦a_i、b_i、c_i（i＝1,2,3）为：

b₁＝cosωsinκ

b₂＝cosωcosκ    （2）

b₃＝-sinω

设投影中心S与地面一物点A在地面坐标系中的坐标分别为(X_S,Y_S,Z_S)和(X_A,Y_A,Z_A)，物点A的像点a在像平面坐标系中的坐标为(x,y)，由中心投影的构像几何关系可得：

$x-x_0=-f\frac{a_1(X_A-X_S)+b_1(Y_A-Y_S)+c_1(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)}$ $\left(3\right)$

$y-y_0=-f\frac{a_2(X_A-X_S)+b_2(Y_A-Y_S)+c_2(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)}$

其中，f为相机焦距，(x₀,y₀)为像主点偏移误差，属内方位因素。

由此可以看出，无人机遥感影像发生几何形变源于：1）：传感器自身内方位因素影响，存在系统误差；2）：传感器外方位变化，二者在不同程度上引起影像的几何形变。中心投影的共线方程为航空遥感影像的几何畸变分析和校正提供了理论基础。

传感器的外方位元素是用于描述传感器成像瞬时的位置和姿态，根据摄影测量学，在严格正射投影情况下，如图4所示，其中S(X_S,Y_S,Z_S)、A(X,Y,Z)分别为物镜中心与物点A在地面直角坐标系(X,Y,Z)内的坐标，(x,y)为物点A在像平面内的投影坐标，由相似三角形关系得：

$\frac{X-X_S}{x}=\frac{{Y-Y}_S}{y}=\frac{{Z-Z}_S}{-f}---\left(4\right)$

当外方位元素偏离理想位置而出现变动时，就会使所获取的图像产生畸变，体现在影像上地物相对位置的坐标关系发生改变。这种类型的图像畸变可通过对传感器的构象方程得以解析。由中心投影定理可知，在成像瞬间传感器的空间位置和姿态决定了遥感图像的几何畸变特性。另外，地表起伏和大气状况也会在一定程度上引起遥感影像几何畸变。

地球表面的地形起伏将使影像中的像点产生位移。在高差为正值的情况下，地形起伏在中心投影影像上造成的像点位移是远离原点向外移动的。如图5所示，摄影中心为S，地物点A相对周边地物的高度差为h，其对应像点为a，距像主点o的距离为r，地物点A′为A在地面上的投影，其对应像点为a′，则

为地形起伏引起的像点位移。则由中心投影构像几何关系可得，

$\mathrm{\Δr}=\stackrel{\‾}{a^{\′}a}=r\frac{h}{H}---\left(5\right)$

通常情况下，成像高度Hh>>h时,像点位移Δr≈0。另外，地球曲率可理解为某种程度上的地形起伏。无人机遥感影像区别于卫星图像，属中低空遥感，可忽略地球曲率对影像几何失真的影响。

大气是引起图像几何畸变的另一个因素，由于大气层是一个非均匀的介质层，大气密度随高度增加而减小，导致电磁波在大气中传播的折射率也随高度增加而改变，从而使电磁波传播的路径变成了一条曲线，引起了像点的位移。如图6所示，摄影中心为S，当无大气折射影响时，地物点A通过直线光线

成像与a点；当有大气折射影响时，A点通过曲线光线AS成像于点a′，因此，在遥感影像上引起像点位移

。无人机航空影像质量与拍摄时刻局部天气情况有很大的关系，由于介质层密度不均匀，大气折射引起的像点位移具有很大随机性，不易得到校正，在条件允许的情况下，应尽可能选择天气晴好的中午时分进行拍摄，以减小由于大气折射引起的图像几何畸变。

通常情况下，遥感影像成像姿态角满足ω<<κ，基于共线方程(3)，在满足转角系统

时，可得关系式（6）。

$A_{\mathrm{out}}\&ap;\left[\begin{array}{ccc}\frac{H}{H^{\&prime;}}& 0& 0\\ 0& \frac{H}{H^{\&prime;}}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& 0\\ -\frac{a_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{a_1{c}_2}{a_1{b}_2-a_2{b}_1}& 0\\ 0& 0& 1\end{array}\right]---\left(6\right)$

其中，a_i,b_i,c_i(i＝1,2,3)为姿态角的余弦值，具体定义见式（2）。H为实际成像高度，H′为遥感影像校正标定高度，为变量。取定不同的H′，基于式（6）对图像进行整体坐标变换即可得到不同分辨率的校正后图像。

第三步，基于地面控制点计算误差补偿矩阵A_{error_cmp}。

通常情况下，基于A_in和A_out校正后的遥感影像与基于地面控制点校正后的遥感影像相比仍然有较大差距，这主要是由于相机安装定位不准确、外方位测量出现偏差等因素造成，这些误差因素通常在一定的时期内维持稳定，即同一批无人机遥感图像有相同的误差值。

基于上述分析，可以通过将经过A_in和A_out校正后的遥感影像与基于地面控制点校正后的遥感影像相比较，获取误差补偿，进一步校正几何畸变图像，使校正后的遥感影像逼近真实图像。通常情况下，遥感影像几何畸变满足透视变换，即几何畸变误差补偿矩阵A_{error_cmp}满足如下形式：

$A_{\mathrm{error}\_\mathrm{cmp}}=\left[\begin{array}{ccc}{\mathrm{\&zeta;}}_1& {\mathrm{\&zeta;}}_2& {\mathrm{\&zeta;}}_3\\ {\mathrm{\&zeta;}}_4& {\mathrm{\&zeta;}}_5& {\mathrm{\&zeta;}}_6\\ {\mathrm{\&zeta;}}_7& {\mathrm{\&zeta;}}_8& 1\end{array}\right]---\left(8\right)$

其中ζ_i(i＝1,2,3,...8)为未知数，8个未知数需要最少需要8个方程求解。

设原始遥感影像四个角点(0,0)、(Width，0)、(Width，Height)和(0,Height)经内方位矩阵A_in和外方位矩阵A_out校正后的坐标为(x₀,y₀)、(x₁,y₁)、(x₂,y₂)和(x₃,y₃)，经地面控制点校正后的坐标为(x₀′,y₀′)、(x₁′,y₁′)、(x₂′,y₂′)和(x₃′,y₃′)。设图像畸变满足透视变换关系，则满足关系：

$\left\{\begin{array}{c}{x_i}^{\&prime;}=\frac{{\mathrm{\&zeta;}}_1\&CenterDot;x_i+{\mathrm{\&zeta;}}_2\&CenterDot;y_i+{\mathrm{\&zeta;}}_3}{{\mathrm{\&zeta;}}_7\&CenterDot;x_i+{\mathrm{\&zeta;}}_8\&CenterDot;y_i+1}\\ {y_i}^{\&prime;}=\frac{{\mathrm{\&zeta;}}_4\&CenterDot;x_i+{\mathrm{\&zeta;}}_5\&CenterDot;y_i+{\mathrm{\&zeta;}}_6}{{\mathrm{\&zeta;}}_7\&CenterDot;x_i+{\mathrm{\&zeta;}}_8\&CenterDot;y_i+1}\end{array}\right.(i=\mathrm{0,1,2,3})---\left(9\right)$

通过解方程（9），即可求得误差补偿矩阵A_{error_cmp}，一旦计算得到A_{error_cmp}就可对本批次遥感影像重复使用，无需重复计算。经误差补偿矩阵校正后的遥感影像精度会大大提高。

第四步，综合考虑各因素影响，整合校正矩阵。

如果按内方位因素、外方位因素、误差补偿矩阵依次对遥感影像进行几何校正，会造成计算效率低下，且计算过程中又会积累新的计算误差。考虑到A_in、A_out和A_{error_cmp}形式上的一致性，可对其进行整合，合成一个变换矩阵A，这样不仅在形式上简化了运算，也可是整个遥感影像几何坐标变换一次完成，大大提高了运算效率，且避免了逐步运算过程中的误差积累，提高了运算精度。

按内方位校正、外方位校正和误差补偿依次进行，可得总体变换矩阵A如下：

                   A＝A_{error_cmp}·A_out·A_in                  （10）

即：

$A=\left[\begin{array}{ccc}{\mathrm{\&zeta;}}_1& {\mathrm{\&zeta;}}_2& {\mathrm{\&zeta;}}_3\\ {\mathrm{\&zeta;}}_4& {\mathrm{\&zeta;}}_5& {\mathrm{\&zeta;}}_6\\ {\mathrm{\&zeta;}}_7& {\mathrm{\&zeta;}}_8& 1\end{array}\right]\left[\begin{array}{ccc}\frac{H}{H^{\&prime;}}& 0& 0\\ 0& \frac{H}{H^{\&prime;}}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& 0\\ -\frac{a_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{a_1{c}_2}{a_1{b}_2-a_2{b}_1}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]---\left(11\right)$

化简后如下：

$A=\left[\begin{array}{ccc}\frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_1{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_2{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_3\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_4{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_5{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_6\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_7{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_8{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& 1\end{array}\right]---\left(12\right)$

其中，a_i,b_i,c_i(i＝1,2,3)为姿态角的余弦值，具体定义见式（2）。H为实际成像高度，H′为遥感影像校正标定高度。(c_x,c_y)是基准点（通常在图像中心）；f_x，f_y是以像素为单位的焦距。ζ_i(i＝1,2,3...8)为误差补偿矩阵参数。

第五步，基于全局变换矩阵A，对遥感影像进行几何校正。

获得全局变换矩阵A后，可直接矩阵A对遥感影像进行全局透视变换，消除遥感影像固有的几何畸变。校正公式如下：

$\left[\begin{array}{c}{x}_{\mathrm{ref}}\\ y_{\mathrm{ref}}\\ 1\end{array}\right]=\mathrm{\&lambda;}\&times;A\&times;\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\mathrm{\&lambda;}\left[\begin{array}{ccc}\frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_1{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_2{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_3\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_4{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_5{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_6\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_7{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_8{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]---\left(13\right)$

其中(x,y)为原始遥感影像坐标，(x_ref,_yref)为校正后的图像坐标。

第六步，遥感影像边界地理坐标定位；

理想情况（正射投影）下，遥感影像记录的地理经纬度坐标(Lat,Lng)对应于图像中心坐标，即

当相机存在固有几何畸变及俯仰角

和滚转角ω不为零时在非理想情况下，经纬度(Lat,Lng)对应的图像坐标(x_C,y_C)会发生偏移；

基于全局透视变换矩阵和成像模型得坐标(x_C,y_C)计算公式如下：

$\left[\begin{array}{c}{x}_C\\ y_C\\ 1\end{array}\right]=\mathrm{\&lambda;}\left[\begin{array}{ccc}\frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_1{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_2{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_3\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_4{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_5{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_6\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_7{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_8{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& 1\end{array}\right]\left[\begin{array}{c}\frac{\mathrm{Width}}{2}\\ \frac{\mathrm{Height}}{2}\\ 1\end{array}\right]---\left(14\right)$

其中，(Width,Height)原始遥感影像图像尺寸。在得到(x_C,y_C)的基础上，为得到校正后图像边界的地理坐标，需要计算校正后图像的横向分辨率R_H和径向分辨率R_V（通常情况下R_H＝R_V）。基于相机模型和成像模型可得：

$R_H=R_V=\frac{u}{f}{H}^{\&prime;}---\left(15\right)$

其中，u为像元尺寸，f为相机焦距，H′为遥感影像标定高度。通常情况下，遥感区域远远小于地球表面，在遥感区域内可近似认为经度分辨率r_Latitude和纬度分辨率r_Longtitude为定值，计算公式如下：

$r_{\mathrm{Latitude}}=\frac{{\mathrm{\&Delta;}}_{\mathrm{Latitude}}}{f(\mathrm{Lat}-\frac{{\mathrm{\&Delta;}}_{\mathrm{Latitude}}}{2},\mathrm{Lng},\mathrm{Lat}+\frac{{\mathrm{\&Delta;}}_{\mathrm{Latitude}}}{2},\mathrm{Lng})}$ $\left(16\right)$

$r_{\mathrm{Longtitude}}=\frac{{\mathrm{\&Delta;}}_{\mathrm{Longtitude}}}{f(\mathrm{Lat},\mathrm{Lng}-\frac{{\mathrm{\&Delta;}}_{\mathrm{Longtitude}}}{2},\mathrm{Lat},\mathrm{Lng}+\frac{{\mathrm{\&Delta;}}_{\mathrm{Longtitude}}}{2})}$

单位：度/米。其中，Δ_Latitude、Δ_Longtitude为常数，默认值为0.02，f为距离计算公式，定义如下：设地面任意两点p₁(Lat₁,Lng₁)、p₂(Lat₂,Lng₂)，则两点之间距离公式：

$f({\mathrm{Lat}}_1,{\mathrm{Lng}}_1,{{\mathrm{Lat}}_2,\mathrm{Lng}}_2)=p_1{p}_2=2*\mathrm{EarthRadius}*\arctan \left(\frac{\sqrt{a}}{\sqrt{1-a}}\right)---\left(17\right)$

其中：

$a={\sin }^2(\mathrm{Lat}2-\mathrm{Lat}1)+\cos \left(\mathrm{Latl}\right)\cos \left(\mathrm{Lat}2\right){\sin }^2\left(\frac{\mathrm{Lng}2-\mathrm{Lng}1}{2}\right)$

EarthRadius＝6378137

设校正后图像尺寸为(W,H)，则校正后图像四个边界的地理坐标具体计算方法如下：

北纬： $\mathrm{Lat}+y_C\&times;\frac{r_{\mathrm{Latitude}}}{R_H}$

南纬： $\mathrm{Lat}-(H-y_C)\&times;\frac{r_{\mathrm{Latitude}}}{R_H}$

西经： $\mathrm{Lng}-x_C\&times;\frac{r_{\mathrm{Longtitude}}}{R_V}$

东经： $\mathrm{Lng}+(W-x_C)\&times;\frac{r_{\mathrm{Longtitude}}}{R_V}$

通过上述步骤，从原始图像得到几何校正后的图像，并得到校正后图像的地理坐标信息。

依次完成上述各步骤，即可从原始图像得到几何校正后的图像，并得到校正后图像的地理坐标信息。附图7给出了原始无人机遥感图像（图2）未经误差补偿后的几何校正影像，图8给出了经本发明算法校正后的图像，图9为图8边界经纬度地理坐标，经实际测量验证，校正图像与实测图像误差在5米以内。

Claims (1)

1.一种基于误差补偿的UAV高分辨率影像几何校正方法，包括以下步骤：

第一步，确定机载相机属性参数，计算内参数矩阵A_in；

机载相机属性参数为：相机焦距f、像元尺寸u和图像尺寸，图像尺寸为图像的宽Width和高Height；

内参数矩阵A_in为：

$A_{\mathrm{in}}=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]---\left(1\right)$

其中，x、y代表像素坐标系方向，(c_x,c_y)是基准点；f_x，f_y是以像素为单位的焦距；

第二步，获取遥感影像外方位参数，计算外参数矩阵A_out；

遥感影像外方位参数为：飞行高度H、俯仰角

、滚转角ω、偏航角κ、成像经纬度（Latitude，Longitude）；

外参数矩阵A_out为：

$A_{\mathrm{out}}\&ap;\left[\begin{array}{ccc}\frac{H}{H^{\&prime;}}& 0& 0\\ 0& \frac{H}{H^{\&prime;}}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& 0\\ -\frac{a_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{a_1{c}_2}{a_1{b}_2-a_2{b}_1}& 0\\ 0& 0& 1\end{array}\right]---\left(2\right)$

其中：H为实际成像高度，H′为遥感影像校正标定高度，a_i,b_i,c_i为姿态角的余弦值，i＝1,2,3，具体为：

b₁＝cosωsinκ

b₂＝cosωcosκ    （3）

b₃＝-sinω

第三步，基于地面控制点计算误差补偿矩阵A_{error_cmp}；

误差补偿矩阵A_{error_cmp}为：

$A_{\mathrm{error}\_\mathrm{cmp}}=\left[\begin{array}{ccc}{\mathrm{\&zeta;}}_1& {\mathrm{\&zeta;}}_2& {\mathrm{\&zeta;}}_3\\ {\mathrm{\&zeta;}}_4& {\mathrm{\&zeta;}}_5& {\mathrm{\&zeta;}}_6\\ {\mathrm{\&zeta;}}_7& {\mathrm{\&zeta;}}_8& 1\end{array}\right]---\left(4\right)$

其中：ζ_i为误差补偿矩阵参数，为未知数，i＝1,2,3,...8；

设原始遥感影像四个角点(0,0)、(Width，0)、(Width，Height)和(0,Height)经内方位矩阵A_in和外方位矩阵A_out校正后的坐标为(x₀,y₀)、(x₁,y₁)、(x₂,y₂)和(x₃,y₃)，经地面控制点校正后的坐标为(x₀′,y₀′)、(x₁′,y₁′)、(x₂′,y₂′)和(x₃′,y₃′)；设图像畸变满足透视变换关系，则满足关系：

$\left\{\begin{array}{c}{x_i}^{\&prime;}=\frac{{\mathrm{\&zeta;}}_1\&CenterDot;x_i+{\mathrm{\&zeta;}}_2\&CenterDot;y_i+{\mathrm{\&zeta;}}_3}{{\mathrm{\&zeta;}}_7\&CenterDot;x_i+{\mathrm{\&zeta;}}_8\&CenterDot;y_i+1}\\ {y_i}^{\&prime;}=\frac{{\mathrm{\&zeta;}}_4\&CenterDot;x_i+{\mathrm{\&zeta;}}_5\&CenterDot;y_i+{\mathrm{\&zeta;}}_6}{{\mathrm{\&zeta;}}_7\&CenterDot;x_i+{\mathrm{\&zeta;}}_8\&CenterDot;y_i+1}\end{array}\right.(i=\mathrm{0,1,2,3})---\left(5\right)$

通过解方程（5），得到误差补偿矩阵A_{error_cmp}；

第四步，获取全局变换矩阵A；

全局变换矩阵A为：

               A＝A_{error_cmp}·A_out·A_in                        （6）

即：

$A=\left[\begin{array}{ccc}{\mathrm{\&zeta;}}_1& {\mathrm{\&zeta;}}_2& {\mathrm{\&zeta;}}_3\\ {\mathrm{\&zeta;}}_4& {\mathrm{\&zeta;}}_5& {\mathrm{\&zeta;}}_6\\ {\mathrm{\&zeta;}}_7& {\mathrm{\&zeta;}}_8& 1\end{array}\right]\left[\begin{array}{ccc}\frac{H}{H^{\&prime;}}& 0& 0\\ 0& \frac{H}{H^{\&prime;}}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{b_2{c}_3}{a_1{b}_2-a_2{b}_1}& 0\\ -\frac{a_2{c}_3}{a_1{b}_2-a_2{b}_1}& \frac{a_1{c}_2}{a_1{b}_2-a_2{b}_1}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]---\left(7\right)$

化简后如下：

$A=\left[\begin{array}{ccc}\frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_1{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_2{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_3\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_4{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_5{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_6\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_7{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_8{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& 1\end{array}\right]---\left(8\right)$

第五步，基于全局变换矩阵A，对遥感影像进行几何校正；

基于全局变换矩阵A，对原始遥感影像进行全局透视变换，校正公式如下：

$\left[\begin{array}{c}{x}_{\mathrm{ref}}\\ y_{\mathrm{ref}}\\ 1\end{array}\right]=\mathrm{\&lambda;}\&times;A\&times;\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\mathrm{\&lambda;}\left[\begin{array}{ccc}\frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_1{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_2{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_3\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_4{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_5{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_6\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_7{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_8{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]---\left(9\right)$

其中：(x,y)为原始遥感影像坐标，(x_ref,y_ref)为校正后的图像坐标；

第六步，遥感影像边界地理坐标定位；

正射投影下，遥感影像记录的地理经纬度坐标(Lat,Lng)对应于图像中心坐标，即当相机存在固有几何畸变及俯仰角和滚转角ω不为零时在非理想情况下，经纬度(Lat,Lng)对应的图像坐标(x_C,y_C)会发生偏移；

成像模型坐标(x_C,y_C)为：

$\left[\begin{array}{c}{x}_C\\ y_C\\ 1\end{array}\right]=\mathrm{\&lambda;}\left[\begin{array}{ccc}\frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_1{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_2{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_3\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_4{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_5{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& {\mathrm{\&zeta;}}_6\\ \frac{H{b}_2{c}_3{\mathrm{\&zeta;}}_7{f}_x}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& \frac{H{a}_1{c}_2{\mathrm{\&zeta;}}_8{f}_y}{H^{\&prime;}(a_1{b}_2-a_2{b}_1)}& 1\end{array}\right]\left[\begin{array}{c}\frac{\mathrm{Width}}{2}\\ \frac{\mathrm{Height}}{2}\\ 1\end{array}\right]---\left(10\right)$

校正后图像的横向分辨率R_H和径向分辨率R_V为：

$R_H=R_V=\frac{u}{f}{H}^{\&prime;}---\left(11\right)$

在遥感区域内近似认为经度分辨率r_Latitude和纬度分辨率r_Longtitude为定值，计算公式如下：

$r_{\mathrm{Latitude}}=\frac{{\mathrm{\&Delta;}}_{\mathrm{Latitude}}}{f(\mathrm{Lat}-\frac{{\mathrm{\&Delta;}}_{\mathrm{Latitude}}}{2},\mathrm{Lng},\mathrm{Lat}+\frac{{\mathrm{\&Delta;}}_{\mathrm{Latitude}}}{2},\mathrm{Lng})}$ $\left(12\right)$

$r_{\mathrm{Longtitude}}=\frac{{\mathrm{\&Delta;}}_{\mathrm{Longtitude}}}{f(\mathrm{Lat},\mathrm{Lng}-\frac{{\mathrm{\&Delta;}}_{\mathrm{Longtitude}}}{2},\mathrm{Lat},\mathrm{Lng}+\frac{{\mathrm{\&Delta;}}_{\mathrm{Longtitude}}}{2})}$

单位：度/米，其中，Δ_Latitude、Δ_Longtitude为常数，f为距离计算公式，定义如下：设地面任意两点p₁(Lat₁,Lng₁)、p₂(Lat₂,Lng₂)，则两点之间距离公式：

$f({\mathrm{Lat}}_1,{\mathrm{Lng}}_1,{{\mathrm{Lat}}_2,\mathrm{Lng}}_2)=p_1{p}_2=2*\mathrm{EarthRadius}*\arctan \left(\frac{\sqrt{a}}{\sqrt{1-a}}\right)---\left(13\right)$

其中：

$a={\sin }^2(\mathrm{Lat}2-\mathrm{Lat}1)+\cos \left(\mathrm{Latl}\right)\cos \left(\mathrm{Lat}2\right){\sin }^2\left(\frac{\mathrm{Lng}2-\mathrm{Lng}1}{2}\right)$

EarthRadius＝6378137

设校正后图像尺寸为(W,H)，则校正后图像四个边界的地理坐标具体计算方法如下：

北纬： $\mathrm{Lat}+y_C\&times;\frac{r_{\mathrm{Latitude}}}{R_H}$

南纬： $\mathrm{Lat}-(H-y_C)\&times;\frac{r_{\mathrm{Latitude}}}{R_H}$

西经： $\mathrm{Lng}-x_C\&times;\frac{r_{\mathrm{Longtitude}}}{R_V}$

东经： $\mathrm{Lng}+(W-x_C)\&times;\frac{r_{\mathrm{Longtitude}}}{R_V}$

通过上述步骤，从原始图像得到几何校正后的图像，并得到校正后图像的地理坐标信息。

Images