一种卫星光学遥感相机内方元素在轨检校方法 技术领域 本发明属于遥感卫星在轨几何检校领域， 特别是涉及一种卫星光学遥感相机内方 元素在轨检校方法。
背景技术 随着国产遥感卫星空间分辨率的提高， 遥感影像的内部几何精度成为一项重要的 考察指标。线阵推扫式光学遥感相机是一种主要的光学遥感卫星载荷， 其内方位元素检校 能够提取光学畸变模型参数， 减小畸变， 提高影像内部几何精度。
目前国内卫星遥感领域， 在光学遥感相机的内方元素几何检校方面， 除了在三线 阵立体成像相机的在轨动态检校方面有一定的研究外， 其他还有对于面阵相机地面检校的 研究， 其中三线阵相机动态检校的重点是研究三线阵相机的主点和焦距及三个相机之间的 几何关系 ； 面阵相机地面检校则要借助地面测试设备如转台和光源发生器， 并校正的是面 阵模型， 但是对于单线阵推扫式相机在轨检校方面则缺少相应的研究。
发明内容 本发明的目的在于克服现有技术的上述不足， 提供一种卫星光学遥感相机内方元 素在轨检校方法， 该方法通过在轨检校光学遥感相机的几何参数， 提高了遥感卫星图像的 内部几何精度， 减小了光学畸变误差。
本发明的上述目的是通过如下技术方案予以实现的 ：
一种卫星光学遥感相机内方元素在轨检校方法， 包括如下步骤 ：
建立 OXYZ 三维坐标系， 其中相机在坐标系中的位置为 O1， 地球中心在坐标系中的 位置为 O2， 坐标系中 X 向为线阵推扫方向， Y 向为垂轨方向， Z 向为遥感相机视轴， 设地面某 控制点 P， 向量 OP 与地球相交于 P 点， OP 与焦平面相交于 P2 点， 由于光学部件设计与加工 缺陷的存在， 视向量 O1P 在焦平面的实际成像位置为 P1 点， 设 P2 点的 Y 坐标为 y2， P1 点的 Y 坐标为 y1， 则根据径向畸变模型公式得出 y2 与 y1 的关系式 ：
y2 ＝ c0+c1y1+c2y12+c3y13
计算模型参数 c0、 c1、 c2 和 c3， 将计算得到的模型参数 c0、 c1、 c2 和 c3 代入上式， 得 到 P2 点修正后的 Y 坐标 y2′， 进一步得到拟合残差 Δy2 ＝ y′ 2-y2， 调整相机焦距 f， 使得 拟合残差 Δy2 最小， 完成在轨检校， 其中计算模型参数 c0、 c1、 c2 和 c3 的步骤如下 ：
(1) 选择几何检校区， 卫星在几何检校区成像后得到 1 级图像， 同时选择与 1 级图 像具有相同分辨率或更高分辨率的正射影像和 DEM 高程图作为参考图， 在 1 级图像上选取 k 个控制点， 同时在所述参考图上选取 k 个同名点， 记录 1 级图像上 k 个控制点的坐标 (m， n)， m 表示列号， n 表示行号， 记录参考图上 k 个同名点的经纬度坐标和高度值 (Lon， Lat， h)， Lon 表示经度， Lat 表示纬度， h 表示高度， 其中 k 为正整数 ；
(2) 根据列号 m， 计算 OXYZ 三维坐标系下控制点的 y 坐标，
y ＝ c(m-M/2)， 其中 M 表示一行的像元个数， c 表示像元尺寸，
得到控制点线阵向量 w0(0， y， f)， 进一步归一化得到控制点线阵向量 wi ；
(3) 根据行号 n， 计算控制点所在行的成像时刻秒计数 t，
t ＝ n·d-t0， 其中 d 为积分时间， t0 为 1 级图像第一行秒计数，
根据秒计数 t， 利用卫星下传的 GPS 数据， 计算成像时刻卫星在地固系的位置坐标 向量 Pt， 并假设卫星和相机的坐标重合， 根据参考图上同名点的经纬度坐标和高度值 (Lon， Lat， h)， 计算得到地固系下控制点的坐标向量 Pj， 控制点坐标向量 Pj 减去卫星位置坐标向 量 Pt 得到控制点视向量 v0， 进一步归一化得到控制点视向量 vi ；
(4) 在 k 个控制点的线阵向量 wi 中选择最接近相机视轴的向量 wij， 计算向量 wij 与视轴的夹角 α ：
y0 表示离相机视轴最近的控制点的 y 坐标 ；分别以最接近相机视轴的向量 wij、 vij 为基准， 计算其它 k-1 个向量与所述向量 wij、 vij 的夹角， 得到 k-1 个控制点线阵向量间夹角 Awi 和 k-1 个控制点视向量间夹角 Avi， 每个夹角都加上一个夹角 α， 得到 Aw’ i 和 Av’ l ；
本步骤中以最接近相机视轴的向量 wij、 vij 为基准， 并假设该向量的畸变可以忽 vij 之间的夹角后， 两组夹角数据的差异反映了光学几何 略， 得到其它向量与基准向量 wij、 畸变的特征， 这是本方法中能够解耦卫星外方位元素影响的关键。
(5) 计算线阵坐标 Bwi 和 Bvi， 计算公式为 ：
Bwi ＝ f*tan(Aw′ i)
Bvi ＝ f*tan(Av′ i)， f 为相机焦距 ；
线阵坐标 Bwi 和 Bvi 分别对应公式 y2 ＝ c0+c1y1+c2y12+c3y13 中的 y1 和 y2， 将 Bwi 和 2 3 Bvi 代入公式 y2 ＝ c0+c1y1+c2y1 +c3y1 ， 用最小二乘法拟合得到模型参数 c0、 c1、 c2 和 c3。
步骤 (5) 中的公式是将两组夹角数据转换为线阵 Y 坐标方向的数值， 使得到的畸 变模型系数单位和 Y 坐标单位一致， 方便应用该模型 ； 对于直接利用角度进行几何校正的 程序， 则可以直接进行角度畸变计算 c0、 c1、 c2 和 c3， 而模型中的 y1 和 y2 也对应于角度值。
在上述卫星光学遥感相机内方元素在轨检校方法中， 步骤 (1) 中几何检校区选择 人工地物目标丰富， 道路交通发达的地区。
在上述卫星光学遥感相机内方元素在轨检校方法中， 步骤 (1) 中在 1 级图像上选 取沿 Y 方向均匀分布， 沿 X 方向较窄区域中的 k 个控制点， 并且 20 ＜ k ＜ 40。 本发明与现有技术相比具有如下优点 ：
(1) 本发明在轨检校方法根据通常的光学镜头畸变模型建立了 3 阶内方元素模 型， 并基于相对角度误差受外方元素误差影响小的原理， 从视向量之间的相对角度关系中 提取了内方元素的模型参数， 实现了内外方元素解耦 ；
(2) 本发明在轨检校方法通过对光学高阶畸变建模， 并通过最小二乘法拟合得到 模型参数， 从而得到修正后的线阵坐标， 并进一步得到拟合残差， 通过调整焦距， 使得拟合 残差最小， 完成几何检校， 本发明方法能够达到较高的检校精度 ；
(3) 本发明在轨检校方法采用 3 阶模型， 具有相当的柔韧性， 对于不同的焦距误差 均能予以吸收， 并通过模型参数平衡能够达到总体模型的高精度 ；
(4) 本发明在选择控制点时， 选取沿像元排列方向均匀分布的控制点， 并且控制点
的选择区域为沿 X 推扫方向的较窄区域， 并且越窄约好， 目的是使得时间引起的姿态误差 达到最小。 附图说明 图 1 为本发明单线阵推扫相机在轨成像模型 ；
图 2 为本发明内方元素几何检校流程图 ；
图 3 为本发明内方元素几何检校过程中控制点选取示意图 ；
图 4 为本发明实施例中相机畸变拟合前后散点图 ；
图 5 为本发明实施例中相机畸变拟合后残差图。
具体实时方式
下面结合附图通过具体实施例对本发明进行进一步详细的描述 ：
如图 1 所示为本发明单线阵推扫相机在轨成像模型， 建立 OXYZ 三维坐标系， 其中 航天相机在坐标系中的位置为 O1， 地球中心在坐标系中的位置为 O2， 图 1 中左图表示线阵推 扫的情况， 右图表示相机视向量和地球相交的情况。相机焦距为 f， X 为线阵推扫方向， Y为 垂轨方向， Z 为相机视轴。设地面某控制点 P， 向量 O1P 与地球相交于 P 点， 与焦平面 XY 相 交于 P2 点， 由于光学部件设计和加工的缺陷， 光学相机总是存在一定的几何畸变， 使得视向
量 O1P 在焦平面的实际成像位置为 P1 点。
通常光学畸变为径向畸变， 提出建立一个 3 阶多项式畸变模型来表达该径向畸 变。
r′＝ c0+c1r+c2r2+c3r3 (1)
其中， r 为畸变后实际成像距离即 OP1， r′为共线径向距离 OP2。
地球中心为 O2， 卫星和地球的距离为 O1O2。
向量 O1P2 在 OXYZ 相机坐标系下构成视向量 w(0， y2， f)， P2 点的 Y 坐标为 y2。向量 O1P1 在 OXYZ 相机坐标系下构成视向量 v(0， y1， f)， P1 点的 Y 坐标为 y1。
根据公式 (1) 的径向畸变模型， y1 和 y2 的关系为 ： 2 3
y2 ＝ c0+c1y1+c2y1 +c3y1 (2)
计算模型参数 c0、 c1、 c2 和 c3， 将计算得到的模型参数 c0、 c1、 c2 和 c3 代入公式 (2)， 得到 P2 点修正后的 Y 坐标 y2′， 进一步得到拟合残差 Δy2 ＝ y′ 2-y2， 调整相机焦距 f， 使 得拟合残差 Δy2 最小， 完成在轨检校。
其中计算模型参数 c0、 c1、 c2 和 c3 的步骤如下 ：
步骤一 ： 控制点数据采集。
选择一景人工地物目标丰富， 道路交通发达的地区作为几何检校区， 卫星在该区 域成像后得到 1 级图像， 同时选择与 1 级图像同分辨率或更高分辨率的正射影像和 DEM 高 程图作为参考图。在经过了辐射校正的 1 级图像上选取控制点， 同时在参考图上选取同名 点。记录 1 级图像上控制点坐标 (m， n)， m 表示列号， n 表示行号， 记录参考图上的经纬度坐 标和高度值 (Lon， Lat， h)， Lon 表示经度， Lat 表示纬度， h 表示高度。
如图 3 所示为本发明内方元素几何检校过程中控制点选取示意图， 每一个 “+” 表 示一个控制点， 选择清晰交叉点， 个数 k 个， 20 ＜ k ＜ 40。控制点沿像元排列方向选择， 均 匀分布 ； 控制点的选择区域为沿 X 推扫方向的较窄区域， 并且越窄约好， 使得时间引起的姿态误差达到最小。
步骤二 ： 数据处理
首先， 根据列号 m， 计算相机坐标系下控制点的 y 坐标 ：
y ＝ c(m-M/2) (3)
这里， M 表示一行的像元个数， c 表示像元尺寸。得到控制点线阵向量 w0(0， y， f)， 进一步归一化得到控制点线阵向量 wi。如图 2 所示为本发明内方元素几何检校流程图。
然后， 根据行号 n， 可以计算控制点所在行的成像时刻秒计数 ：
t ＝ n·d-t0 (4)
其中 d 为积分时间， t0 为该景第一行秒计数。
根据卫星秒计数可以推算出成像时刻的卫星位置在地固系下的向量 Pt， 并假设卫 星质心和相机位置 O1 重合。
根据地面控制点的经纬度和高度值 (Lon， Lat， h)， 计算得到地固系下控制点坐标 向量 Pj。
控制点坐标向量 Pj 减去卫星位置向量 Pt 得到控制点视向量 v0， 并对该向量归一化 得到控制点视向量 vi。 步骤三 ： 角度计算
上述步骤二分别得到了归一化线阵向量 wi 和归一化视向量 vi， 对于 k 个控制点， 有线阵向量 wi1 ～ k 和视向量 vi1 ～ k。
首先， 在一系列线阵向量 wi1 ～ k 中选择最接近相机视轴的向量 wij， 计算其与视轴的 夹角 a， 该夹角计算公式为 ：
其中 y0 表示离相机视轴最近的控制点的 y 坐标，
分别以最接近相机视轴的的向量 wij、 vij 为基准， 计算其他 k-1 个向量与所述向量 wij、 vij 的夹角， 得到 k-1 个控制点线阵向量间夹角 Awi 和 k-1 个控制点视向量间夹角 Avi， 每个夹角都加上一个常值角度 α， 得到新的 Aw’ i 和 Av’ l。
根据三角关系计算在 OXYZ 相机坐标系下的各个 y 坐标， 即线阵坐标， 得到 Bwi 和 Bvi， 计算公式为 ：
Bwi ＝ f*tan(Aw′ i)
Bvi ＝ f*tan(Av′ i) (6)
线阵坐标 Bwi 和 Bvi 分别对应公式 y2 ＝ c0+c1y1+c2y12+c3y13 中的 y1 和 y2， 将 Bwi 和 2 3 Bvi 代入 y2 ＝ c0+c1y1+c2y1 +c3y1 ， 用最小二乘法拟合得到模型参数 c0、 c1、 c2 和 c3。
基于模型参数 c0、 c1、 c2 和 c3， 可以得到修正后的 P2 点的 Y 坐标 y2′， 进一步得到 拟合残差 Δy2 ＝ y′ 2-y2。调整相机焦距 f， 使得拟合残差 Δy2 最小， 完成在轨检校。
下面列举一个具体的实施例
以我国环境减灾卫星 1B 的 CCD2 数据为参考， 抽取 2009 年 8 月 29 日位于安徽和江 苏北部地区的一景遥感影像， 进行试验计算。相机 CCD 像元尺寸 0.065mm， 对应地面 30 米。 影像参考图选择 landsat 的 ETM 正射影像， 高程数据为 STRM90。控制点共有 21 个， 平面精 度约为 50 米， 高程精度约为 5 米。
如图 4 所示为本发明实施例相机畸变拟合前后散点图， 由图 4 可知拟合前畸变误 差非常明显， 其中 ‘☆’ 表示控制点畸变误差坐标， ‘+’ 表示控制点误差拟合修正坐标， 本实 施例分别选择了 5 个焦距值， 得到 5 组结果， 见表 1 为畸变拟合结果图， 图 4 的相机畸变拟 合前后散点图为选取表 1 中序号 1 数据得到的结果图。
表1
序号 1. 2. 3. 4. 5.
焦距 /mm 140.8 140.9 141.0 141.1 141.2C0 0.00095899 0.0023983 0.0038375 0.0052768 0.006716C1 0.99776 0.99847 0.99918 0.99989 1.0006C2 1.7785E-5 1.7874E-5 1.7964E-5 1.8053E-5 1.8142E-5C3 2.0774E-5 2.0789E-5 2.0804E-5 2.0819E-5 2.0833E-5残差 /mm 0.015083 0.015094 0.015104 0.015115 0.015126从结果可以看出， 3 阶畸变模型得到的残差大约在 0.015mm， 等于 2.3 个像元， 残差 分布见图 5， 说明残差已经没有明显规律。考虑参考点精度也大约在 2 个像元左右， 因此说 明本发明在轨检校方法达到了预期精度， 同时， 本发明方法采用的 3 阶模型， 具有相当的柔 韧性， 对于不同的焦距差异， 均能予以吸收， 图 5 相机畸变拟合后各个控制点的残差图。
以上所述， 仅为本发明最佳的具体实施方式， 但本发明的保护范围并不局限于此， 任何熟悉本技术领域的技术人员在本发明揭露的技术范围内， 可轻易想到的变化或替换， 都应涵盖在本发明的保护范围之内。
本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。