CN100541232C - 无姿态信息条件下的航空多光谱扫描仪几何粗校正方法 - Google Patents

Abstract

本发明公开了一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法。包括如下步骤：1)对多光谱扫描仪数据进行坐标转换，将WGS-84大地坐标系下的测量观测值转换为高斯平面直角坐标下的值；2)根据理想飞行模型，利用1秒模式模拟滚动角；3)在竖直平面内根据高度数据拟合一条曲线，然后根据曲线的切向计算每一个更新数据点处的飞行俯仰角；4)利用中心投影构像方程求得机下点坐标；5)由扫描方式，根据飞行高度、扫描角、瞬时扫描角求取同一扫描行其它点的坐标；6)利用直接法生成粗校正后的图像。本发明提高了航空遥感的地理精度，是对无姿态信息条件下遥感技术的一项创新，提高时效性，使航空遥感技术能够更好的应用于国民经济生产生活中。

Description

无姿态信息条件下的航空多光谱扫描仪几何粗校正方法

技术领域

本发明涉及遥感技术领域，尤其涉及一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法。

背景技术

从60年代至今，机载多光谱扫描图像几何特性的研究一直为人们所关注。在无姿态信息条件下，传统的几何校正方法仅仅利用控制点和最小二乘法拟合求取空间转换多项式，经插值求得校正图像。该方法需耗费大量的人力于控制点的选取，精度受主、客观因素的限制。更主要的是此方法仅能在有限的子区内，在一定的程度上消除图像的整体畸变，对于由平台姿态、运动变化引起的高频畸变则无法消除。

1990年，B.J.Devcrqux等提出了匹配Delannay三角形的机载多光谱扫描图像几何畸变校正方法^[1]。其原理是在图像空间三角形ABC与校正空间三角行A‘B’C’对应，任意在三角形ABC中的象元X与顶点A可用矢量AX联系起来，表达为线性组合，待定系数α、β一旦确定，X′在校正空间的对应点存在A′X′＝αA′B′+βA′C′。与传统方法相比，该方法虽在精度上略有提高，但同样需选取大量的控制点以形成三角形联网。

另一类校正方法称非参数校正方法^[2]，将原始图像和校正图像空间坐标间的差异视为二维随机场。最典型的代表是算术平均法，其利用一组控制点，对每一点计算其原始空间坐标与校正空间坐标的差值，通过加权平均来完成几何校正，权因子与至每个GCP的距离成反比。

早在1962年Elms就提出了利用位置和方位参数对带式相机图像进行几何畸变校正，但当时获取参数的仪器缺乏或者过于昂贵以至于未能实现^[2]。随着惯导系统以及GPS定位技术的发展，从空间实现象元点几何定位已成为可能。美国、加拿大等国都根据自己所能获取的参数及手段建立了相应的校正模式。由于欧美国家具有先进的POS/DG(Position and Orientation Solutions for DirectGeoreferencing)等测量飞行姿态角的设备，对没有姿态角的情况研究较少。在国内，机载多光谱扫描图像几何特性的研究起步比较晚。其中遥感应用研究所的吴传庆等人，利用POS/DG获得的姿态角数据，完成了另一种机载遥感仪器——OMIS的几何校正^[3]。熊桢等人在不利用控制点和姿态值的情况下，对OMIS图像进行了航线校正，但在精度上仍需进一步提高^[4]。本专利申请人曾经发表文章《基于GPS的机载多通道扫描仪图像几何校正》^[5]，文中采用前人已有的方法模拟了飞行姿态进而几何校正，而本专利则根据理想飞行模型采用了新的1秒模式进行了模拟，具有更高的几何校正精度。

参考文献：

1.郭德方。遥感图象的计算机处理和模式识别。北京：电子工业出版杜，1984。

2.KindeLan M，Moreno V，Valverde A，Geometric Correction of AirborneMulti-Scanner Images.15th Intemational Symposium on Remote Sensing ofenvironment，1981：1539.

3.吴传庆。基于POS/DG的无稳定平台机载高光谱图像几何校正。中国科学院遥感应用研究所硕士论文。

4.熊桢，王向军，郑兰芬，童庆禧。基于GPS数据的OMIS图像航线校正研究。遥感技术与应用，2000，15(1)：1～5.

5.刘涛，毛志华，王迪峰，潘德炉。基于GPS的机载多通道扫描仪图像几何校正。仪器仪表学报，2006年S3期，2224-2226。

发明内容

本发明的目的是提供了一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法。

本发明包括如下步骤：

1)对多光谱扫描仪数据进行坐标转换，将WGS-84大地坐标系下的测量观测值转换为高斯平面直角坐标下的值，利于校正后遥感反演应用；

2)根据理想飞行模型，利用1秒模式模拟滚动角；

3)根据高度数据拟合飞机俯仰角，即在竖直平面内根据高度数据拟合一条曲线，然后根据曲线的切向计算每一个更新数据点处的飞行俯仰角；

4)利用中心投影构像方程求得机下点坐标；

5)由扫描方式，根据飞行高度、扫描角、瞬时扫描角求取同一扫描行其它点的坐标；

6)利用直接法生成粗校正后的图像。

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的根据理想飞行模型：理想飞行模型是把飞行近似看成半径不断变化的相对于水平面的圆周运动，由于飞机在空中方位角不停地摆动，导致飞行半径也不停地变化，所以时间间隔越短，飞机作同一个飞行半径的圆周运动的概率越大；因此可以假设飞机在一秒内作圆周运动。

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用1秒模式模拟飞行滚动角：假设P1，P2是飞行一秒前后空间点；其两点的地面高斯大地坐标为(x1，y1)，(x2，y2)；V1，V2为一秒前后的地面速率；θ表示V1与V2之间方向的夹角；记P1，P2之间的距离为Len，则有飞行轨迹的曲率半径：

$r=\frac{\mathrm{Len}}{2*\sin \left(\frac{\mathrm{\θ}}{2}\right)}$ 其中 $\mathrm{Len}=\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2}$

则飞机的滚动角即侧翻角ω：

其中：g为引力常数，m为飞机质量，v表示飞行速度，F为向心力。

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用中心投影构像方程求得机下点坐标：

在地面坐标系与传感器坐标系之间建立的转换关系称为通用构像方程。设地面点P在地面坐标系中的坐标为(X，Y，Z)_P，P在传感器坐标系中的坐标为(U，V，W)_P，传感器投影中心S在地面坐标系中的坐标为(X，Y，Z)_S，传感器的姿态角为

则通用构像方程为：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+A{\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P$

式中：A为传感器坐标系相对地面坐标系的旋转矩阵，是传感器姿态角的函数。根据中心投影特点，图象坐标(x，y，-f)和传感器系统坐标(U，V，W)_P之间有如下关系：

${\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P={\mathrm{\λ}}_p\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]$

λ_p为成像比例尺分母，f为扫描成像主距，中心投影像片坐标与地面点大地坐标的关系即中心投影构像方程为：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+{\mathrm{\λ}}_{p}A{\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]}_p$

其中：

$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$

具体表达式为：

a₂₁＝cosωsink

a₂₂＝cosωcosk

a₂₃＝-sinω

由像点坐标可以解算大地平面坐标，称为正算公式：

$X_P=X_S+(Z_P-Z_S)\frac{a_{11}x+a_{12}y-a_{13}f}{a_{31}x+a_{32}y-a_{33}f}$

$Y_P=Y_S+(Z_P-Z_S)\frac{a_{21}x+a_{22}y-a_{23}f}{a_{31}x+a_{32}y-a_{33}f}$

令x＝0，y＝0，Z_p＝0，得到扫描行中心点像元地面坐标公式如下：

$X_{\mathrm{mid}}=X_s-Z_s*\frac{a_{13}}{a_{33}}$

$Y_{\mathrm{mid}}=Y_s-Z_s*\frac{a_{23}}{a_{33}}$

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用直接法生成粗校正后的图像：

从原始图像阵列出发，按行列的顺序依次对每个原始像素点位求其在地面坐标系，即输出图像坐标系中的正确位置： $\left\{\begin{array}{c}X=F_x(x,y)\\ Y=F_y(x,y)\end{array}\right.$

式中F_x和F_y为直接校正变换函数，同时，把该像素的亮度值移置到由上式算得的输出图像中的相应点位上去，

这样生成的bmp图，会存在一些没有灰度值的点，

对那些有灰度值的点，把灰度值直接赋予相应的位置，对于没有灰度值的点，以5像元为半径向周围搜索有灰度值的点，把各点的灰度值以距离平方的倒数为权值叠加给所求的点，叠加公式为

$P=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i{z}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i}$

其中P为所求的灰度值， $W=\frac{1}{d_i^2},$ z_i为点i的灰度值，d_i为点i到所求点的距离。本发明的具有的有益效果：

航空遥感技术能够机动灵活的获取大范围地面信息，广泛地应用于农、林、地质、石油、测绘、城市规划、海洋执法等领域，对国民经济的发展有重要的意义。然而和卫星遥感器相比，机载成像光谱仪具有姿态稳定性差、飞行高度低、视场角度大等特点，这些因素的结合，使得所获图像的几何畸变更加复杂，所以航空遥感图像几何校正是一个十分重要的技术手段，是提高对地观测资料可信度的前提条件。本发明在水平面内，根据GPS数据拟合的曲线，计算每一个节点处的曲率半径r，然后根据理想飞行模型和动力学原理计算飞机的滚动角，进而模拟其它飞机姿态信息，实现了航空多光谱扫描仪的几何粗校正，提高了航空遥感的地理精度，为精校正的进行打下了良好基础，是对无姿态信息条件下遥感技术的一项创新，具有极大的实用价值。

本发明目的在于发明一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，提高航空多光谱扫描仪的几何粗校正精度，同时提高时效性，使航空遥感技术能够更好的应用于国民经济生产生活中。

附图说明

图1是无姿态信息条件下的航空多光谱扫描仪几何粗校正方法技术路线图；

图2是本发明的一秒模式航迹示意图；

图3是本发明的飞行半径求解示意图；

其中：P1，P2是飞行一秒前后空间点；其两点的地面高斯大地坐标为(x1，y1)，(x2，y2)；V1，V2为一秒前后的地面速率；θ表示V1与V2之间方向的夹角；记P1，P2之间的距离为Len；

图4是本发明的仰俯角计算示意图；

图5是构象方程中的坐标系示意图

图6是多光谱扫描仪扫描几何示意图；

图7是几何粗校正前的原始数据示例；

图8是完成粗校正后的最终图像。

具体实施方式

无姿态信息条件下的航空多光谱扫描仪几何粗校正方法包括如下步骤：

1)对多光谱扫描仪数据进行坐标转换，将WGS-84大地坐标系下的测量观测值转换为高斯平面直角坐标下的值，利于校正后遥感反演应用；

2)根据理想飞行模型，利用1秒模式模拟滚动角；

3)根据高度数据拟合飞机俯仰角，即在竖直平面内根据高度数据拟合一条曲线，然后根据曲线的切向计算每一个更新数据点处的飞行俯仰角；

4)利用中心投影构像方程求得机下点坐标；

5)由扫描方式，根据飞行高度、扫描角、瞬时扫描角求取同一扫描行其它点的坐标；

6)利用直接法生成粗校正后的图像。

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的根据理想飞行模型：理想飞行模型是把飞行近似看成半径不断变化的相对于水平面的圆周运动，由于飞机在空中方位角不停地摆动，导致飞行半径也不停地变化，所以时间间隔越短，飞机作同一个飞行半径的圆周运动的概率越大；因此可以假设飞机在一秒内作圆周运动。

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用1秒模式模拟飞行滚动角：假设P1，P2是飞行一秒前后空间点；其两点的地面高斯大地坐标为(x1，y1)，(x2，y2)；V1，V2为一秒前后的地面速率；θ表示V1与V2之间方向的夹角；记P1，P2之间的距离为Len，则有飞行轨迹的曲率半径：

$r=\frac{\mathrm{Len}}{2*\sin \left(\frac{\mathrm{\θ}}{2}\right)}$ 其中 $\mathrm{Len}=\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2}$

则飞机的滚动角即侧翻角ω：

其中：g为引力常数，m为飞机质量，v表示飞行速度，F为向心力。

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用中心投影构像方程求得机下点坐标：

在地面坐标系与传感器坐标系之间建立的转换关系称为通用构像方程。设地面点P在地面坐标系中的坐标为(X，Y，Z)_P，P在传感器坐标系中的坐标为(U，V，W)_P，传感器投影中心S在地面坐标系中的坐标为(X，Y，Z)_S，传感器的姿态角为

则通用构像方程为：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+A{\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P$

式中：A为传感器坐标系相对地面坐标系的旋转矩阵，是传感器姿态角的函数。根据中心投影特点，图象坐标(x，y，-f)和传感器系统坐标(U，V，W)_P之间有如下关系：

${\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P={\mathrm{\λ}}_p\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]$

λ_p为成像比例尺分母，f为扫描成像主距，中心投影像片坐标与地面点大地坐标的关系即中心投影构像方程为：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+{\mathrm{\λ}}_{p}A{\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]}_p$

其中：

$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$

具体表达式为：

a₂₁＝cosωsink

a₂₂＝cosωcosk

a₂₃＝-sinω

由像点坐标可以解算大地平面坐标，称为正算公式：

$X_P=X_S+(Z_P-Z_S)\frac{a_{11}x+a_{12}y-a_{13}f}{a_{31}x+a_{32}y-a_{33}f}$

$Y_P=Y_S+(Z_P-Z_S)\frac{a_{21}x+a_{22}y-a_{23}f}{a_{31}x+a_{32}y-a_{33}f}$

令x＝0，y＝0，Z_p＝0，得到扫描行中心点像元地面坐标公式如下：

$X_{\mathrm{mid}}=X_s-Z_s*\frac{a_{13}}{a_{33}}$

$Y_{\mathrm{mid}}=Y_s-Z_s*\frac{a_{23}}{a_{33}}$

所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用直接法生成粗校正后的图像：

从原始图像阵列出发，按行列的顺序依次对每个原始像素点位求其在地面坐标系，即输出图像坐标系中的正确位置： $\left\{\begin{array}{c}X=F_x(x,y)\\ Y=F_y(x,y)\end{array}\right.$

式中F_x和F_y为直接校正变换函数，同时，把该像素的亮度值移置到由上式算得的输出图像中的相应点位上去，

这样生成的bmp图，会存在一些没有灰度值的点，

对那些有灰度值的点，把灰度值直接赋予相应的位置，对于没有灰度值的点，以5像元为半径向周围搜索有灰度值的点，把各点的灰度值以距离平方的倒数为权值叠加给所求的点，叠加公式为

$P=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i{z}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i}$

其中P为所求的灰度值， $W=\frac{1}{d_i^2},$ z_i为点i的灰度值，d_i为点i到所求点的距离。

具体来说，几何校正所用数据来自于机载GPS，在GPS中，记录着飞机在飞行中每秒钟的基本信息，其中可以用来进行几何校正的有飞机所在位置大地坐标，飞行速度，飞行高度和飞行偏航角共4种信息。由以上4种信息，直接可以校正的仅有速高比变化引起的畸变和正切畸变。由飞机所在位置坐标、飞行速度和飞行方位角和合理的飞行模型可以反演出近似飞行滚动角和俯仰角，因而姿态角引起的畸变可以近似得以校正。总体上本发明的技术路线如图1所示。

其中包括如下关键步骤：

步骤1、对多光谱扫描仪数据进行坐标转换，将WGS-84大地坐标系下的测量观测值转换为高斯平面直角坐标下的值，利于校正后遥感反演应用；

GPS卫星全球定位系统的观测值是WGS-84大地坐标系下的测量观测值(纬度，经度，高程)，而在航空成像光谱遥感图像处理与分析中，所采用的坐标系是高斯平面直角坐标系，像元大小的单位是m。几何校正的方法也是基于高斯平面直角坐标系建立的.因此在对图像进行几何校正前，必须把GPS大地坐标观测值(纬度，经度，高度)转换为高斯平面直角坐标值(x，y，z)。

高斯投影坐标正算公式如下：

(1)高斯投影正算：已知椭球面上某点的大地坐标(L，B)，求该点在高斯投影平面上的直角坐标(x，y)，即 $(L,B)\⇒(x,y)$ 的坐标变换。

(2)投影变换必须满足的条件

●中央子午线投影后为直线；

●中央子午线投影后长度不变；

●投影具有正形性质，即正形投影条件；

(3)投影过程

在椭球面上有对称于中央子午线的两点P₁和P₂，它们的大地坐标分别为(L，B)及(l，B)，式中l为椭球面上P点的经度与中央子午线(L₀)的经度差：l＝L-L₀，P点在中央子午线之东，l为正，在西则为负，则投影后的平面坐标一定为P′₁(x，y)和P′₂(x，y)。

(4)计算公式

$x=l\left(B\right)+\frac{t}{2}N{\cos }^2{\mathrm{Bl}}^2+\frac{t}{24}N{\cos }^{4}B(5-t^2+9{\mathrm{\η}}^2)l^4$

$+\frac{t}{720}N{\cos }^{6}B(61-58t^2+t^4+270{\mathrm{\η}}^2-330t^2{\mathrm{\η}}^2)l^6$

$y=N\cos \mathrm{Bl}+\frac{1}{6}N{\cos }^{3}B(1-t^2+{\mathrm{\η}}^2)l^3$

$+\frac{1}{120}N{\cos }^5(5-18t^2+t^4+14{\mathrm{\η}}^2-58t^2{\mathrm{\η}}^2)l^5$

为了使得y永远大于零，令：y＝y+500000.

式中：l(B)为从赤道到投影点的子午线弧长； $N=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}$ 为卯酉圈半径；

t＝tanB；l＝L-L₀为经差；L₀为中央子午线经度； $e^{\′}=\sqrt{\frac{a^2-b^2}{a^2+b^2}};$ η＝e′cosB；a，b分别代表地球的长短轴半径。

l(B)＝α(B+βsin2B+γsin4B+δsin6B+εsin8B)

其中：

$\left\{\begin{array}{c}\mathrm{\α}=\frac{a+b}{2}(1+\frac{1}{4}{n}^2+\frac{1}{64}{n}^4)\\ \mathrm{\β}=-\frac{3}{2}n+\frac{9}{16}{n}^3-\frac{3}{32}{n}^5\\ \mathrm{\γ}=\frac{15}{16}{n}^2-\frac{15}{32}{n}^4\\ \mathrm{\δ}=-\frac{35}{48}{n}^3+\frac{105}{256}{n}^5\\ \mathrm{\ϵ}=\frac{315}{512}{n}^4\\ n=\frac{a-b}{a+b}\end{array}\right.$

另外，为了使坐标转换引起的几何变形最小，当坐标转换时，选取图像所在区域的平均子午线为投影的中央子午线。

步骤2、根据理想飞行模型，利用1秒模式模拟滚动角；

飞机飞行的精确模式很难得知，但在采集数据的飞行过程中，飞行员尽量保持其平稳飞行，因而可以把飞行近似看成半径不断变化的相对于水平面的圆周运动。三点可以确定一个圆。由于飞机在空中方位角不停地摆动，导致飞行半径也不停地变化，所以时间间隔越短，飞机作同一个飞行半径的圆周运动的概率越大；所以可以假设飞机在一秒内作圆周运动。示意图见图2。

其中飞行轨迹的曲率半径r计算原理请见图3。

则有以下推导：

$\mathrm{Len}=\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2})$

$r=\frac{\mathrm{Len}}{2*\sin \left(\frac{\mathrm{\θ}}{2}\right)}$

则飞机的滚动角即侧翻角ω：

其中：g为引力常数，m为飞机质量，v表示飞行速度，F为向心力。

步骤3、根据高度数据拟合飞机俯仰角，即在竖直平面内根据高度数据拟合一条曲线，然后根据曲线的切向计算每一个更新数据点处的飞行俯仰角；

在竖直平面内根据高度数据拟合一条曲线，然后根据相邻高度数据计算每一个更新数据点处的飞行俯仰角(见图4)。

x1与x3之间的曲线是拟合的曲线，求点x2的仰附角，可以由x1与x3之间的直线的斜率求得：

$\mathrm{\γ}=\arctan \left(\frac{H_3-H_1}{S}\right)$ 其中 $S={\left({(X_3-X_1)}^2{+(Y_3-Y_1)}^2\right)}^{\frac{1}{2}}$

式中，H₃，H₁分别是遥感平台在点x1，x3的高；X₁，X₃，Y₁，Y₃分别是遥感平台在点x1，x3的高斯X轴，Y轴坐标。

步骤4、利用中心投影构像方程求得机下点坐标；

前述步骤2和步骤3已获得飞机滚动角和俯仰角，而偏航角已经记录在数据的GPS信息中，这样我们就获得了飞机的所有三个姿态信息。然后可以进入图像数据点真实坐标的求解。

遥感图像的构像方程是指地物点在图像上的图像坐标(x，y)和其在地面对应点的大地坐标(X、Y、Z)之间的数学关系。根据摄影测量原理，这两个对应点和传感器成像中心成共线关系，可以用共线方程来表示。这个数学关系是对任何类型传感器成像进行几何纠正和对某些参量进行误差分析的基础。

传感器坐标系S-UVW，S为传感器投影中心，作为传感器坐标系的坐标原点，U轴的方向为遥感平台的飞行方向，V轴垂直于U，W轴则垂直于UV平面，该坐标系描述了像点在空间的位置。地面坐标系O-XYZ，主要采用地心坐标系统。当传感器对地成像时，Z轴与原点处的天顶方向一致，XY平面垂直于Z轴。图像(像点)坐标系o-xyf，(x，y)为像点在图像上的平面坐标，f为传感器成像时的等效焦距，其方向与S-UVW方向一致。参见图5。

上述坐标系统都是三维空间坐标系，而最基本的坐标系统是图像坐标系统o-xy和地图坐标系统O_m-X_mY_m，它们是二维的平面坐标系统，是遥感图像几何处理的出发点和归宿。

在地面坐标系与传感器坐标系之间建立的转换关系称为通用构像方程。设地面点P在地面坐标系中的坐标为(X，Y，Z)_P，P在传感器坐标系中的坐标为(U，V，W)_P，传感器投影中心S在地面坐标系中的坐标为(X，Y，Z)_S，传感器的姿态角为

则通用构像方程为：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+A{\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P$

式中：A为传感器坐标系相对地面坐标系的旋转矩阵，是传感器姿态角的函数。根据中心投影特点，图象坐标(x，y，-f)和传感器系统坐标(U，V，W)_P之间有如下关系：

${\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P={\mathrm{\λ}}_p\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]$

λ_p为成像比例尺分母，f为扫描成像主距，中心投影像片坐标与地面点大地坐标的关系即中心投影构像方程为：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+{\mathrm{\λ}}_{p}A{\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]}_p$

其中：

$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$

具体表达式为：

a₂₁＝cosωsink

a₂₂＝cosωcosk

a₂₃＝-sinω

由像点坐标可以解算大地平面坐标，称为正算公式：

$X_P=X_S+(Z_P-Z_S)\frac{a_{11}x+a_{12}y-a_{13}f}{a_{31}x+a_{32}y-a_{33}f}$

$Y_P=Y_S+(Z_P-Z_S)\frac{a_{21}x+a_{22}y-a_{23}f}{a_{31}x+a_{32}y-a_{33}f}$

令x＝0，y＝0，Z_p＝0，得到扫描行中心点像元地面坐标公式如下：

$X_{\mathrm{mid}}=X_s-Z_s*\frac{a_{13}}{a_{33}}$

$Y_{\mathrm{mid}}=Y_s-Z_s*\frac{a_{23}}{a_{33}}$

步骤5、由扫描方式，根据飞行高度、扫描角、瞬时扫描角求取同一扫描行其它点的坐标；

先利用中心投影构像方程就得机下点坐标，然后由扫描方式，求得同一扫描行其它点的坐标(扫描示意图请参见图6)。

其中：H海拔高度，γ仰俯角，h遥感器到扫描行中心点的距离，S_i第i个像元到扫描行中心点的距离，W_i第i个扫描线与中心点扫描线的夹角。

求取扫描行其它点大地坐标值：

Xi＝Xmid+(h*tan(Wi))*cos(fw)

Yi＝Ymid+(h*tan(Wi))*sin(fw)

其中：fw是多光谱扫描仪扫描方向角；

第i点的滚动角；Xmid，Ymid分别表示扫描行中心点的坐标；Xi，Yi分别为第i个扫描点的横纵坐标。

步骤6、利用直接法生成粗校正后的图像。

在输出图像边界及其坐标系统确立后，就可以按照选定的纠正变换函数把原始数字图像逐个像素变换到图像贮存空间中去。这里采用直接法，从原始图像阵列出发，按行列的顺序依次对每个原始像素点位求其在地面坐标系，即输出图像坐标系中的正确位置： $\left\{\begin{array}{c}X=F_x(x,y)\\ Y=F_y(x,y)\end{array}\right.$

式中F_x和F_y为直接校正变换函数，同时，把该像素的亮度值移置到由上式算得的输出图像中的相应点位上去，

这样生成的bmp图，会存在一些没有灰度值的点，

对那些有灰度值的点，把灰度值直接赋予相应的位置，对于没有灰度值的点，以5像元为半径向周围搜索有灰度值的点，把各点的灰度值以距离平方的倒数为权值叠加给所求的点，叠加公式为

$P=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i{z}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i}$

其中P为所求的灰度值， $W=\frac{1}{d_i^2},$ z_i为点i的灰度值，d_i为点i到所求点的距离。

Claims (5)

1.一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于包括如下步骤：

1)对多光谱扫描仪数据进行坐标转换，将WGS-84大地坐标系下的测量观测值转换为高斯平面直角坐标下的值，利于校正后遥感反演应用；

2)根据理想飞行模型，利用1秒模式模拟滚动角；

3)根据高度数据拟合飞机俯仰角，即在竖直平面内根据高度数据拟合一条曲线，然后根据曲线的切向计算每一个更新数据点处的飞行俯仰角；

4)利用中心投影构像方程求得机下点坐标；

5)由扫描方式，根据飞行高度、视场角、瞬时视场角求取同一扫描行其它点的坐标；

6)利用直接法生成粗校正后的图像。

2.根据权利要求1所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的理想飞行模型是把飞行近似看成半径不断变化的相对于水平面的圆周运动，由于飞机在空中方位角不停地摆动，导致飞行半径也不停地变化，所以时间间隔越短，飞机作同一个飞行半径的圆周运动的概率越大；因此可以假设飞机在一秒内作圆周运动。

3.根据权利要求1所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用1秒模式模拟飞行滚动角：假设P1，P2是飞行一秒前后空间点；其两点的地面高斯大地坐标为(x1，y1)，(x2，y2)；V1，V2为一秒前后的地面速率；θ表示V1与V2之间方向的夹角；记P1，P2之间的距离为Len，则有飞行轨迹的曲率半径：

$r=\frac{\mathrm{Len}}{2*\sin \left(\frac{\mathrm{\θ}}{2}\right)}$ 其中 $\mathrm{Len}=\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2}$

则飞机的滚动角即侧翻角ω：

$\mathrm{tg}\left(\mathrm{\ω}\right)=\frac{F}{G},$ 其中 $F=\frac{{\mathrm{mv}}^2}{r},$ G＝mg，因此可以得到 $\mathrm{\ω}=\mathrm{arctg}\left(\frac{v^2}{\mathrm{rg}}\right)$

其中：g为引力常数，m为飞机质量，v表示飞行速度，F为向心力。

4.根据权利要求1所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用中心投影构像方程求得机下点坐标：在地面坐标系与传感器坐标系之间建立的转换关系称为通用构像方程；设地面点P在地面坐标系中的坐标为(X，Y，Z)_P，P在传感器坐标系中的坐标为(U，V，W)_P，传感器投影中心S在地面坐标系中的坐标为(X，Y，Z)_s，传感器的姿态角为

则通用构像方程为：

${\left[\begin{array}{c}X\\ T\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+A{\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P$

式中：A为传感器坐标系相对地面坐标系的旋转矩阵，是传感器姿态角的函数；根据中心投影特点，图象坐标(x，y，-f)和传感器系统坐标(U，V，W)_P之间有如下关系：

${\left[\begin{array}{c}U\\ V\\ W\end{array}\right]}_P={\mathrm{\λ}}_p\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]$

λ_p为成像比例尺分母，f为扫描成像主距，中心投影像片坐标与地面点大地坐标的关系即中心投影构像方程为：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_P={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_S+{\mathrm{\λ}}_{p}A{\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]}_p$

其中：

$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$

具体表达式为：

a₂₁＝cosωsink

a₂₂＝cosωcosk

a₂₃＝-sinω

由像点坐标可以解算大地平面坐标，称为正算公式：

$X_P=X_S+(Z_P-Z_S)\frac{a_{11}x+a_{12}y-a_{13}f}{a_{31}x+a_{32}y-a_{33}f}$

$Y_P=Y_S+(Z_P-Z_S)\frac{a_{21}x+a_{22}y-a_{23}f}{a_{31}x+a_{32}y-a_{33}f}$

令x＝0，y＝0，Z_p＝0，得到扫描行中心点像元地面坐标公式如下：

$X_{\mathrm{mid}}=X_s-Z_s*\frac{a_{13}}{a_{33}}$

$Y_{\mathrm{mid}}=Y_s-Z_s*\frac{a_{23}}{a_{33}}$

5.根据权利要求1所述的一种无姿态信息条件下的航空多光谱扫描仪几何粗校正方法，其特征在于所述的利用直接法生成粗校正后的图像：

从原始图像阵列出发，按行列的顺序依次对每个原始像素点位求其在地面坐标系，即输出图像坐标系中的正确位置： $\left\{\begin{array}{c}X=F_x(x,y)\\ Y=F_y(x,y)\end{array}\right.$

式中F_x和F_y为直接校正变换函数，同时，把该像素的亮度值移置到由上式算得的输出图像中的相应点位上去，

这样生成的bmp图，会存在一些没有灰度值的点，

对那些有灰度值的点，把灰度值直接赋予相应的位置，对于没有灰度值的点，以5像元为半径向周围搜索有灰度值的点，把各点的灰度值以距离平方的倒数为权值叠加给所求的点，叠加公式为

$P=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i{z}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{W}_i}$

其中P为所求的灰度值， $W=\frac{1}{d_i^2},$ z_i为点i的灰度值，d_i为点i到所求点的距离。

Images