CN110111274A - 一种星载推扫式光学传感器外方位元素定标方法 - Google Patents
一种星载推扫式光学传感器外方位元素定标方法 Download PDFInfo
- Publication number
- CN110111274A CN110111274A CN201910351621.7A CN201910351621A CN110111274A CN 110111274 A CN110111274 A CN 110111274A CN 201910351621 A CN201910351621 A CN 201910351621A CN 110111274 A CN110111274 A CN 110111274A
- Authority
- CN
- China
- Prior art keywords
- control point
- image
- error
- point
- model
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 107
- 230000003287 optical effect Effects 0.000 title claims abstract description 23
- 239000011159 matrix material Substances 0.000 claims abstract description 75
- 230000009194 climbing Effects 0.000 claims description 32
- 238000003384 imaging method Methods 0.000 claims description 27
- 230000008569 process Effects 0.000 claims description 19
- 238000009434 installation Methods 0.000 claims description 15
- 238000013507 mapping Methods 0.000 claims description 14
- 230000008859 change Effects 0.000 claims description 12
- 238000013519 translation Methods 0.000 claims description 11
- 238000012937 correction Methods 0.000 claims description 10
- VAYOSLLFUXYJDT-RDTXWAMCSA-N Lysergic acid diethylamide Chemical compound C1=CC(C=2[C@H](N(C)C[C@@H](C=2)C(=O)N(CC)CC)C2)=C3C2=CNC3=C1 VAYOSLLFUXYJDT-RDTXWAMCSA-N 0.000 claims description 8
- 238000005259 measurement Methods 0.000 claims description 6
- 230000007704 transition Effects 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 4
- 238000001514 detection method Methods 0.000 claims description 4
- 230000006641 stabilisation Effects 0.000 claims description 4
- 238000011105 stabilization Methods 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 3
- 230000007547 defect Effects 0.000 claims description 3
- 230000002950 deficient Effects 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 238000005516 engineering process Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 230000036244 malformation Effects 0.000 claims description 3
- VMXUWOKSQNHOCA-UKTHLTGXSA-N ranitidine Chemical compound [O-][N+](=O)\C=C(/NC)NCCSCC1=CC=C(CN(C)C)O1 VMXUWOKSQNHOCA-UKTHLTGXSA-N 0.000 claims description 3
- 238000011524 similarity measure Methods 0.000 claims description 3
- 238000012706 support-vector machine Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- -1 OXY, oxy Chemical group 0.000 claims description 2
- 238000013499 data model Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 4
- 230000009897 systematic effect Effects 0.000 abstract description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000007796 conventional method Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000012876 topography Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 238000011545 laboratory measurement Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/80—Geometric correction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明涉及一种相机的几何定标方法，特别涉及一种星载推扫式光学传感器外方位元素定标方法，属于航空航天领域。一种星载推扫式光学传感器外方位元素定标方法，采用顾及误差时间特性的偏置矩阵模型，利用控制点基于最小二乘原理解求未知数，实现面阵相机的高精度外方位元素定标。本发明基于高精度配准控制点实现卫星面阵相机几何定标，避免系统误差对视频卫星几何定位精度的影响，提升卫星几何定位精度，对保障卫星在动态观测领域的应用效果具有十分重要的意义。
Description
技术领域
本发明涉及一种相机的几何定标方法，特别涉及一种星载推扫式光学传感器外方位元素定标方法，属于航空航天领域。
背景技术
在建立成像几何模型的过程中，影像外方位元素的准确性将直接影响模型的精度。外定向过程是由相机坐标系下探元指向恢复其在地心地固坐标系下的指向。高精度几何检校是保障卫星几何质量的关键技术，该技术利用地面控制数据消除星上成像系统误差，提升影像几何定位精度。高精度几何检校是实现卫星高精度定位的关键环节，也是卫星探测和消除高频误差的基础。目前，国外研究的都是利用几何定标场控制数据实现高精度几何定标。2011年法国发射Pleiades卫星后，CNES的研究人员利用该卫星的高敏捷成像能力，提出了不依赖于几何定标场的系列自检校方法，在不依赖外部高精度控制数据的条件下实现了系统误差的消除补偿，并得到了不逊于常规定标场方法的定标精度。然而，该方法的实现，依赖Pleiades卫星的高敏捷、高稳定性。国内对线阵推扫光学卫星的在轨几何检校研究起步较晚，如，徐建艳等人提出了利用偏移矩阵补偿CBERS卫星系统误差；张过等人分析了不同姿态角对几何定位的影响，提出了利用不同成像条件影像分别求解偏置矩阵角度的方法；刘楚斌等人研究了ALOS PRISM影像的姿态系统误差补偿，利用1个控制点改正姿态角常差等。以上方法实现了姿态系统误差补偿，对提升卫星影像的系统定位精度有一定效果，但是精度不高，且无法考虑由于平台稳定性不够以及陀螺漂移对姿态确定精度的影响。
因此研究星载推扫式光学传感器外方位元素定标方法，提出常量偏置矩阵模型和顾及误差时间特性的偏置矩阵模型，以解决系统误差对卫星几何定位精度的影响问题，提升卫星几何定位精度，对保障卫星在动态观测领域的应用效果具有十分重要的意义。
发明内容
本发的目的是提供一种星载推扫式光学传感器外方位元素定标方法，该方法基于高精度配准控制点实现卫星面阵相机几何定标，避免系统误差对视频卫星几何定位精度的影响，提升卫星几何定位精度，对保障卫星在动态观测领域的应用效果具有十分重要的意义。
为实现上述目的，本发明采取以下技术方案。
一种星载推扫式光学传感器外方位元素定标方法，采用顾及误差时间特性的偏置矩阵模型，利用控制点基于最小二乘原理解求未知数，实现面阵相机的高精度外方位元素定标。
一种星载推扫式光学传感器外方位元素定标方法，包括如下步骤：
步骤1，控制点量测。
遥感卫星应用需要控制点。采用控制点影像库的卫星遥感影像控制点的半自动选取技术，准确而又快速地量测出控制点，且能够保证控制点影像量测的精度；半自动选取技术主要包括控制点影像库的建立与应用。
步骤1.1，建立控制点数据库。
建立地面控制点影像库首先需要在地面应用系统的数据库中创建有关地面控制点列表(TABLE)，然后通过提取地面控制点，将有关控制点的信息存入数据库中，地面控制点影像库表结构如表1。
表1地面控制点影像库表结构
提取地面控制点，建立控制点数据库的流程为：(1)在经扫描纠正的地形图上利用选取控制点的工具选取地面控制点；(2)在原始的卫星遥感影像上找到同名点，并选择同名点周围的邻域，即窗口图像，将所述窗口图像的元数据和影像存入控制点数据库中。重复此操作，控制点数据库在几何纠正处理过程中逐渐建立和丰富，最终形成完整的、覆盖全球的控制点数据库。
若发生变化的地面特征为控制点数据库中控制点所对应的地面区域，则需要及时更新控制点库；
对建立后的控制点数据库进行更新的方法为：
(1)当地面特征发生明显变化时，及时更新控制点库；
(2)当地面特征发生非明显变化时，通过预先设置阈值，然后将观察日志文件中控制点的匹配相关系数与阈值进行对比，当达到阈值时，及时更新控制点库；
步骤1.2，调用控制点数据库。
由于控制点数据库的不断更新，可能一个地物在数据库中对应着多个同名控制点的信息。当进行卫星遥感影像几何纠正时，系统为了快速从控制点数据库中检索到最合适的控制点，遵循下述原则:选择与待纠正影像采集时间最接近的控制点进行自动匹配。具体流程如下：从控制点数据库中检索位于待纠正景或启算景影像范围内的控制点经纬度以及控制带影像拍摄条件。若检索的结果中有相同经纬度的控制点，则读取控制点的拍摄时间与待纠正景或启算景影像时间最接近的控制点数据；若检索的结果中无相同经纬度的控制点，则只读取当前控制点数据。根据以上的读取结果判断能否获取读取范围内所有控制点。若可以，则用控制点影像与原始卫星遥感影像匹配；若不可以，则重复本步骤。
步骤1.3，相似性测度。
采用相关系数作为步骤1.2影响匹配的匹配测度，其定义为：
式中，
其中，m、n为模板大小，且N＝m*n，相关系数的取值在[-1，1]，越接近1则影像间相似程度越高；(i，j)为目标区中的像元行列号，(c，r)为搜索区中心的坐标；gi，j为待纠正影响控制点(i，j)的灰度值，g′i，j为影像数据库中与待纠正影像采集时间最接近的控制点的灰度值；gi+r，j+c为待纠正影响控制点(i+r，j+c)的灰度值，g′i+r，j+c为影像数据库中与待纠正影像采集时间最接近的控制点的灰度值。
步骤1.4，消除影像库中影像的旋转和比例畸变。
影像的几何变形主要表现为旋转关系和平移关系的组合，比例关系也涵盖在旋转关系中，而平移关系能够由影像匹配来解决，因此通过旋转、缩放控制点影像块的方式来解决旋转和比例畸变带来的问题。首先根据严密成像模型或RPC模型计算控制点影像的北方向与原始影像沿着轨道方向的差异，估计卫星的旋转角度和比例系数，方法如下：
分别建立原始影像和控制点影像的影像坐标OXY、oxy，在控制点影像的影像坐标上选取两个点P1(x0，y0)、P2(x1，y1)，并使P1和P2的连线与oy轴平行。通过反转换，得出P1点在原始影像的影像坐标上对应的点为P1′(X0，Y0)，P2点在原始影像的影像坐标上对应的点为P2′(X1,Y1′)。过P1′作OY轴的平行线与过P2′作OX轴的平行线交与一点P3′(X2,Y2)。此时线段P1′P2′与线段P2′P3′的夹角称为旋转角，令旋转角为：
比例系数m为：
式中，P'2P'3＝X0-X1；P'2P'3＝Y1-Y0；P1(x0,y0)为控制点影像中心，P2(x1,y1)为控制点影像块下边缘中点。P'1(X0,Y0)和P'2(X1,Y'1)是控制点影像上的对应点利用变换模型投影到原始影像上的投影点，D(P1,P2)和D(P′1,P′2)分别为两点之间的像素距离。
然后将控制点影像模板绕卫星本体坐标的质心进行旋转，在旋转过程中，加入比例系数改正，在旋转后的模板影像当中取一个内接的矩形区域作为新的控制点影像模板，此时，需保证控制点数据库中的控制点仍然在新的模板中心。控制点数据库中的控制点根据新的控制点影像模板生成新的控制点影像。即完成了旋转和比例畸变的消除，生成新的控制点影像。
步骤1.5，采用改进的爬山法搜索策略实现快速控制点匹配，即对步骤1.4生成的新的控制点影像的控制点匹配
步骤1.5.1：利用严密成像模型或RPC模型，根据匹配点预测的最大视差确定搜索窗口大小。
减小控制点匹配的搜索范围，不仅能够减小计算量，而且能够减少系统在读取影像过程中的开销。根据匹配点预测的最大视差来决定搜索窗口的大小。根据严密成像模型或RPC模型能够显著提高控制点的预测精度，减小搜索窗口的大小。
步骤1.5.2：利用严密成像模型或RPC模型的预测点位精度确定金字塔匹配的层数。
首先，由于生成金字塔时采用低通滤波，使得高层上的影像保留原始影像本身大的整体结构特征，滤掉局部小的非平稳信息和噪声，从而在高层上减少匹配的不确定性；其次，高层匹配结果能够作为初值向下传递，减少下一层匹配的搜索范围，降低匹配过程中的计算复杂度。一幅影像既含有高频信息又含有低频信息。高频信息在影像中占的比例小，低频信息占的比例大。利用金字塔匹配模式，在较粗的分辨率下获得一个比较接近控制点实际位置的初始值，然后再在原始分辨率下计算影像相关系数。
金字塔匹配中层数确定是关键问题，采用根据严密成像模型或RPC模型的预测点位精度来确定金字塔层数。
步骤1.5.3：改进的爬山法搜索策略。
在每一层金字塔的搜索过程中，采用改进的爬山法搜索策略，以减少计算相关系数的次数，并且保证搜索的极值点为整体搜索区域内的极值点。
所述改进的爬山法搜索策略为：在任意一个空间直角坐标系中，预测点P0位于非主峰的区域内，根据爬山法的原理，找到一个峰顶S1，在S1位置沿着X方向就是P0A方向做一次整体搜索，找到X方向的整体峰值点P1，然后沿着X方向就是P1B方向做整体搜索，搜索到Y方向的整体峰值点P2，然后将峰顶S1与整体峰值点P1、整体峰值点P2对比得出最高峰值；
传统的爬山法搜索属于启发式搜索策略，采用启发信息引导，使搜索朝目标方向进行，避免走弯路，从而提高求解的效率。
但是有其内在的局限性，从而使爬山法搜索的区域容易陷入局部极值点。由于有启发信息的引导，爬山法使问题能够尽快得到解。但爬山法也有局限性，即当山只有单峰时，爬山法能够保证找到峰顶。如果山有多峰，且爬山的起始点处于非主峰的区域时，爬山法就难以找到真正的顶峰。
在任意一个空间直角坐标系中，预测点P0位于非主峰的区域内，根据爬山法的原理，找到一个峰顶S1，并认为这就是最高峰，但实际上最高峰是S2，而传统爬山法无法区别。即使知道S1并不是最高峰，还要再搜索，但由于当前所处位置S1是一个局部最高点，再往东、南、西、北各个方向走一步都会导致高度下降，从而无法确定下一步该往哪个方向走。因此在这种情况下爬山法是无法找到最高峰的。由此可见，即使问题有解，爬山法也难以找到解，这主要是因为爬山法的起始搜索点在非主峰区域内，在多峰的情况下，只能爬到局部极值点。
因此，改进起始点的位置使起始点在主峰的区域内是爬山法能迅速达到最佳极值点的有效保证。改进的爬山法搜索策略是在起始点位置沿着X方向就是P0A方向做一次整体搜索，找到X方向的整体峰值点P1，然后沿着P1B方向做整体搜索，搜索到Y方向的整体峰值点P2，该点在主峰区域内的概率比较高。因为该重新搜索的初始位置只要尽可能保证在主峰的区域内即可，因此在X方向搜索和Y方向搜索过程中，按照一定的步距进行，没有必要每个点都计算相关系数。然后在重新搜索的初始位置开始，用传统的爬山法做启发式搜索，能够快速准确搜索到主峰S2。
步骤1.6，计算控制点匹配精度。
采用相邻域点相关系数拟合椭球抛物面来获得该亚像元的二维匹配精度。
二维匹配精度的方法为：以步骤1.4得到的新的控制点影像为目标窗口，在原始影像的搜索窗口中取和目标窗口相同大小的多个影像窗口进行运算，求得代表各个影像窗口中心象素的相关关系的相关系数。根据改进的爬山法搜索策略，对每个影像窗口进行搜索，得到各自的最大相关系数点；最大相关系数点就是所求的待匹配影像控制点的同名点。
当每个影像窗口的同名象素的中心点是真正的同名点时，通过计算同名点即可求得控制点匹配精度；
由于控制点影像模板与原始影像拍摄时间的差异，每个影像窗口的同名象素的中心点不是真正的同名点时。根据偏离情况，按照如下方法计算控制点匹配精度；
情况一、真正的同名点在两个方向上偏离象素中心半个象素之内，则产生了误差，所述误差的大小为象素级。显然，该项误差在两个方向分别均服从内的均匀分布，因此单方向的相关精度为其中
因此，
同理σy＝0.29Δ
因此即二维整象素匹配的理论精度为0.4个象素。
结合理论精度，将最大相关系数位置的邻域点中的相关系数与平差函数联系起来，从而将平差函数的极值点处作为寻求的同名点；
情况二、其他任何偏离情况均采用下述方法计算；
椭圆抛物面方程拟合式：
式中，(x,y)为5×5邻域内每点相对于中心的位置；z为控制点影像块的相关系数；p、q为椭圆抛物面的参数，pq＞0；p、q、x0和y0用间接观测平差求得，当p≠q时，X和Y两个方向的匹配精度不一样，当p＝q时，椭圆抛物面退化为旋转抛物面，两个方向的匹配精度一致，此时椭圆抛物面顶点(x0,y0)作为寻求的亚像元匹配的同名点。由相关系数椭圆抛物面拟合在信噪比较高时使相关精度达到0.1个象素的精度。
由式(4)求得的椭圆抛物面的顶点即为亚像元精度的同名点。原因如下：同名点与周围5×5邻域符合椭圆抛物面，因此根据同名点以及周围5×5的范围内相关系数，结合式(4)求得同名点。
步骤2，计算最小二乘偏置矩阵。
步骤2.1，建立常量偏置矩阵模型。
设备安装误差与姿轨测量误差等效，因此仅需根据姿轨测量误差特性构建补偿模型。对于视场角较小的高分光学卫星而言，轨道位置误差引起的几何定位误差为平移误差，与俯仰角误差、滚动角误差具有等效性；对沿卫星飞行方向和垂直卫星飞行方向对卫星的轨道位置误差进行分析可知，轨道位置误差与姿态误差的等效性。卫星真实位置与卫星带误差位置之间的距离能够转换成俯仰角误差
将平移误差等效为姿态角误差，采用偏置矩阵Ru补偿姿态误差，修正真实光线指向与带误差光线指向间的偏差，则由式(5)得常量偏置矩阵模型式(6)：
式中，为t时刻GPS相位中心在WGS84坐标系下的位置矢量；为t时刻J2000坐标系相对于WGS84坐标系的转换矩阵；为t时刻body坐标系相对于J2000坐标系的转换矩阵；为camera坐标系相对body坐标系的转换矩阵；(Dx Dy Dz)T为GPS相位中心在本体坐标系的坐标；为相机坐标系原点与本体坐标系原点偏移；代表相机主距；ψ为CCD阵列沿轨向偏航角；为影像坐标；为主视轴即过主点垂直于CCS线阵的垂点对应的位置；为探元大小；m为比例系数。
Ru为正交旋转矩阵，Ru的作用是使得与之相乘的坐标系绕坐标系本身的y轴、x轴、z轴旋转角度ωu、κu。即：
则式(6)写成：
令
则：
显然，(Xb Yb Zb)T是由地面点坐标确定的光线本体系下的指向；而(xb xb xb)T是由像方坐标确定的光线本体系下的指向；Ru用于修正两者的偏差从而实现姿态误差补偿。
展开式(9)，变化为：
其中：
对式(11)进行线性化并构建误差方程：
v误差＝A误差x误差-L误差；p误差 (12)
其中，x误差为L误差为根据初值计算的p误差为观测权值，A误差为系数矩阵具体为：
则：
偏置矩阵中待求未知数为三个偏置角，而一个平高控制点只能列两个方程。因此，已知两个控制点即可解求偏置矩阵。
步骤2.2，建立顾及误差时间特性的最小二乘偏置矩阵模型。
针对姿态系统误差引起的平移误差和旋转误差，常量偏置矩阵能很好的补偿。但是，由于平台稳定性不够以及陀螺漂移对姿态确定精度的影响，几何定位模型中常常存在着漂移误差。
常量偏置矩阵无法补偿姿态漂移误差，为此，在常量偏置矩阵中引入偏置角的时间变化率：
将式(14)带入公式(6)以替换Ru，得到最小二乘偏置矩阵模型：
步骤3，评估偏置矩阵的精度，以及引入偏置角的时间变化率偏置矩阵的精度。
步骤3.1，偏置补偿矩阵RU。
考虑到偏置补偿矩阵，式(5)变化为：
RU表示待消除的外方位元素误差；所述外方位元素误差包括轨道误差、姿态误差和安装误差。RU定义成式(17)：
其中，ω、κ分别是相对于X、Y、Z轴旋转的欧拉角。偏置补偿矩阵补偿的是一景短周期范围内的外方位元素误差。
步骤3.2，估计长周期变化模型。
上述偏置补偿矩阵仅仅能对一景时间范围内的外方位元素误差进行建模。实际上，由于卫星姿态精度受到低频漂移误差和结构变形的而影响，存在一个长周期的变化，其变化模型目前来说仍是未知的。目前提高无控定位精度的方法为大区域网平差方法，但由于长周期变化模型的存在，会破坏区域网平差模型中的正态分布，所以不适用于长周期变化。
因此需要对外方位元素的长周期变化模型估计，长周期变化模型表示如下：
其中，θ1,…,θn,t为备选的与长周期变化相关的物理和几何参数，所述物理参数包括：太阳角度、季节以及时间；所述几何参数包括：地理位置；
对于上述长周期变化模型，传统的策略是首先采用地面控制点求出某些时刻的ω、κ及对应参数值，然后对上述备选参数进行显著性分析以删选相关参数，最后采用支持向量机等方法根据相关参数进行函数拟合。从上述步骤1和步骤2能够看出，传统策略最大的难题在于对于控制点的需求量特别大，然而目前全球大范围高精度地面控制点获取的难度比较大，难以支持长周期变化模型的估计。
因此，本发明采用具有重叠度的不同轨道的错位关系探测长周期变化模型的方法，该方法不依赖地面控制点，而是对具有重叠度的不同轨道之间的影像进行同名点匹配，根据同名点与长周期变化模型的关系构建函数关系。
a)对具有重叠度的不同轨道之间的影像进行同名点匹配，当前一轨为j，则后一轨为j+1，两轨分别对应的长周期变化模型为：
b)以每轨的影像景为单位，借助于SRTM30米格网数据，将j+1轨上的当前景影像经过模型正算得到地面控制点，然后根据式(17)计算得到两景之间的相对安装值Δω、Δκ，则此值满足如下关系：
c)对具有重叠区域的景分别构建式(21)，将式(21)与相对安装值的方程联立得到方程组，求解方程组得到式(18)的各个绝对值。需要注意的是，由于式(21)为相对模型，故组成的方程组实际为秩亏方程。为了解决这个问题，对当中任意几景影像采用式(17)直接求取外方位元素误差值，误差值代入式(21)中以解决秩亏问题。相对于传统方法，所有轨的数据都需要控制，此方法仅仅需要若干景有控制即可。
d)将上述(c)中的由任意几景影像求取的外方位元素误差值构成备选参数组，对参数组中的参数进行显著性分析，选出精度够高的参数；
采用支持向量机等方法根据精度够高的参数进行函数拟合，得到长周期变化模型的函数表达式；通过长周期变化模型的函数表达式即可很好的评估偏置矩阵的精度，并获得引入偏置角的时间变化率偏置矩阵的精度。
有益效果：
1、本发明公开的一种星载推扫式光学传感器外方位元素定标方法，是通过建立地面控制点库，采用半自动方式辅助点位量测，实现推扫式传感器外方位元素定标。节省时间成本，提升效率；
2、本发明公开的一种星载推扫式光学传感器外方位元素定标方法，通过偏置矩阵补偿卫星轨道和姿态误差，提升外定向精度，实现高精度外方位元素定标，定标精度高，进而提升效率。
3、本发明公开的一种星载推扫式光学传感器外方位元素定标方法相对于传统方法，所有轨的数据都需要控制，但此方法仅仅需要若干景有控制即可。极大的提高了外方位元素定标的效率。
附图说明
图1为本发明步骤1.1控制点数据库建立的系统结构体系图；
图2为本发明步骤1.2几何纠正对控制点数据库调用流程图；
图3为本发明步骤1.4计算旋转角度和比例系数示意图；
图4为本发明步骤1.5.3传统爬山法示意图；
图5为本发明步骤1.5.3改进的爬山法示意图；
图6为本发明步骤2.1轨道误差与姿态误差等效示意图，其中：(a)沿轨向等效性、(b)垂轨向等效性；
图7为本发明实施例1中的河南嵩山和天津检校场的相关影像：其中(a)为河南嵩山检校场1：2000正射影像图，(b)为河南嵩山检校场1：2000数字高程模型图，(c)为天津1：2000正射影像图，(d)为天津1：2000数字高程模型图；
图8为本发明实例中常量偏置矩阵补偿定位残差(FWD(a,b)，BWD(c,d),NAD(e,f)，红色代表垂轨，蓝色代表沿轨)。
具体实施方式
为了更好的说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。
实施例1：
收集了覆盖我国河南嵩山区域、天津区域的1:2000数字正射影像及数字高程模型作为检校控制数据，图7(a)、(b)所示。其中，天津区域覆盖范围约为100km(西东)×50km(南北)，区域内地势平坦，最大高差在30m以内；而河南区域覆盖范围约为50km(西东)×50km(南北)，区域内主要为丘陵地形，最大高差不超过1500m；两个区域正射影像分辨率均优于0.2m，数字高程模型分辨率优于1m；对应的，收集了覆盖河南嵩山的资源三号三线阵影像一景(代表前正后三景，成像于2012年2月3日)和覆盖天津区域的资源三号三线阵影像两景(分别成像于2012年2月28日和2012年5月2日)。
本实施例公开的一种星载推扫式光学传感器外方位元素定标方法具体实现步骤如下：
步骤1，控制点量测。
步骤1.1，控制点影像库的建立。
控制点影像库的建立首先需要在地面应用系统的数据库中创建有关地面控制点列表(TABLE)，然后通过提取地面控制点，将有关控制点的信息存入数据库中。
提取地面控制点流程为：(1)在经扫描纠正的地形图上利用选取控制点的工具选取地面控制点；(2)在原始的卫星遥感影像上找到同名点，并选择其周围的邻域(窗口图像)，该窗口图像的元数据和影像被存入控制点数据库的控制点表中。因此，地面控制点影像库在几何纠正处理过程中逐渐建立和丰富，最终形成完整的、覆盖全球的控制点库。
实现上述地面控制点提取，可采用Client/Server系统结构，如图1所示。
步骤1.2，控制点影像库的调用。
由于控制点影像库的不断更新，可能一个地物在数据库中对应着多个同名控制点的信息。当进行卫星遥感影像几何纠正时，系统为了快速从控制点数据库中检索到最合适的控制点，进行自动匹配，可遵循下述原则:选择与待纠正影像采集时间最接近的控制点。流程图如图2所示。
步骤1.3，相似性测度。
采用相关系数作为匹配测度，其定义为：
式中，
步骤1.4，旋转和比例畸变的消除方法。
首先根据严密成像模型或RPC模型计算控制点影像块北方向与原始影像沿着轨道方向的差异，估计一个旋转角度和比例系数，然后将控制点影像模板绕其中心进行旋转，在该过程中，加入比例系数改正，在旋转后的模板影像当中取一个内接的矩形区域作为新的控制点影像模板(此时，控制点仍然在新的模板中心)。
如图3所示，旋转角为：
比例系数m为：
步骤1.5，控制点匹配效率。
(1)利用严密成像模型或RPC模型确定搜索窗口大小。
减小控制点匹配的搜索范围，不仅可以减小计算量，而且可减少系统在读取影像过程中的开销。一般根据匹配点预测的最大视差来决定搜索窗口的大小。根据严密成像模型或RPC模型可以显著提高控制点的预测精度，减小搜索窗口的大小。
(2)利用严密成像模型或RPC模型确定金字塔匹配的层数。
首先，由于生成金字塔时采用了低通滤波，使得高层上的影像保留了原始影像本身大的整体结构特征，滤掉了局部小的非平稳信息和噪声，从而在高层上减少了匹配的不确定性；其次，高层匹配结果可作为初值向下传递，减少下一层匹配的搜索范围，降低匹配过程中的计算复杂度。一幅影像既含有高频信息又含有低频信息。高频信息在影像中占的比例小，低频信息占的比例大。利用金字塔匹配模式，在较粗的分辨率下获得一个比较接近控制点实际位置的初始值，然后再在原始分辨率下计算影像相关系数。
金字塔匹配中层数确定是一个关键问题，采用根据严密成像模型或RPC模型的预测点位精度来确定金字塔层数。
(3)改进的爬山法搜素策略。
在每一层金字塔的搜索过程中，采用改进的爬山法搜索策略，以减少计算相关系数的次数，并且保证搜索的极值点为整体搜索区域内的极值点。
改进起始点的位置使起始点在主峰的区域内是爬山法能迅速达到最佳极值点的有效保证。如图5所示，改进的搜索策略使在起始点位置沿着X方向就是P0A方向做一次整体搜索，找到X方向的整体峰值点P1，然后沿着P1B方向做整体搜索，搜索到Y方向的整体峰值点P2，该点在主峰区域内的概率比较高。因为该重新搜索的初始位置只要尽可能保证在主峰的区域内即可，因此在X方向搜索和Y方向搜索过程中，可以按照一定的步距进行，没有必要每个点都计算相关系数。然后在重新搜索的初始位置开始，用传统的爬山法做启发式搜索，可以快速准确搜索到主峰S2。
步骤1.6，控制点匹配精度。
按照相关运算的过程，以旋转变换和比例校正后的控制点影像为目标窗口，在原始影像的搜索窗口中取和目标窗口相同大小的一影像进行相关运算，求的代表各窗口中心象素的相关关系的相关系数。根据改进的爬山法搜索策略，最终的搜索中止于其周围邻域点的相关系数小于该点的相关系数的位置。其中，最大相关系数点就是认为是所求的同名点。由于控制点影像模板与原始影像拍摄时间的差异，同名象素的中心点一般并不是真正的同名点。真正的同名点可能在两个方向上偏离象素中心半个象素之内，这就使得相关产生误差。其中，单方向的相关精度为其中
因此，
同理σy＝0.29Δ
因此即二维整象素匹配的理论精度为0.4个象素。
为了把同名点求得更精确一些，可以将最大相关系数位置的邻域点中的相关系数同一个平差函数联系起来，从而将其函数的极值点处作为寻求的同名点，精度将会更好些。
最佳匹配点与周围5×5邻域基本符合图8的椭圆抛物面，因此根据最佳匹配点以及周围5×5的范围内相关系数，用一个椭圆抛物面方程拟合，此时椭圆抛物面的顶点即为亚像元精度的同名点。
取椭圆抛物面方程式的一般形式：
式中，(x,y,z)为5×5邻域内每点相对于中心的位置和与控制点影像块的相关系数；p、q为椭圆抛物面的参数，pq＞0；p、q、x0和y0用间接观测平差求得，当p≠q时，X和Y两个方向的匹配精度不一样，当p＝q时，椭圆抛物面退化为旋转抛物面，两个方向的匹配精度一致，此时椭圆抛物面顶点(x0,y0)作为寻求的亚像元匹配的同名点。由相关系数椭圆抛物面拟合可使相关精度达到0.1个象素的精度(当信噪比较高时)。但是相关精度和信噪比成反比关系，当信噪比较小时，即使采用相关系数椭球抛物面拟合，也不能提高相关精度。
步骤2，最小二乘偏置矩阵计算。
本实例用的是顾及误差时间特性的偏置矩阵模型。
针对姿态系统误差引起的平移、旋转误差，常量偏置矩阵能很好的补偿。但是，由于平台稳定性不够以及陀螺漂移对姿态确定精度的影响，几何定位模型中常常存在着漂移误差。
常量偏置矩阵无法补偿姿态漂移误差，为此，在常量偏置矩阵中引入偏置角的时间变化率：
按照常量偏置矩阵类似方法，利用控制点基于最小二乘原理解求各未知数。
采用河南区域三线阵一景、天津区域三线阵两景构建联合检校模型，利用步骤1高精度匹配算法，从河南区域、天津区域正射影像中获取控制点，分别为前、后、正视影像获取控制点42198、35186、33208个，所有控制点在检校时段影像上均匀、密集分布；检校平差过程中解求三景影像的常量偏置矩阵，而仅解求同一套内方位元素模型参数。
上表中以河南景正视影像为例，对比了偏置矩阵解求前后的几何定位精度。其中，A代表利用卫星发射前实验室测量的相机安装等参数直接定位的精度，而B代表解求常量偏置矩阵后的定位精度。由表中结果可知，由于卫星发射过程中的受力等因素影响，相机安装在轨后发生改变，导致直接定位精度接近1km；而偏置矩阵能够很好地吸收相机安装等外方位元素误差，定位精度提升至2个像素左右(约2.1m×2＝4.2m)。但是，常量偏置矩阵并未能彻底消除该景外方位元素误差；图8中给出了前、后、正视影像的定位残差图，其中(a)、(c)、(e)代表解求常量偏置后的定位残差，而(b)、(d)、(f)代表解求顾及误差时间特性的偏置矩阵后的定位残差。显然，(a)、(c)、(e)上均可看到残差随时间的变化趋势。因此，需采用顾及误差时间特性的偏置矩阵作为资源三号外检校模型。
步骤3，评估偏置矩阵精度。
步骤3.1，偏置补偿矩阵。
目前卫星的定轨精度都是比较高的，所以一般用式(15)中的偏置补偿矩阵RU去表示待消除的外方位元素误差(包括轨道和姿态误差、安装误差)。RU一般定义成式(16)：
步骤3.2，长周期变化模型估计方法。
上述偏置补偿矩阵仅仅能对一景时间范围内的外方位元素误差进行建模。实际上，由于卫星姿态精度受到低频漂移误差和结构变形的而影响，存在一个长周期的变化，其变化模型目前来说仍是未知的。因此，对外方位元素的长周期变化模型进行研究。
设此长周期变化模型表示如下：
其中，θ1,…,θn,t为备选的与长周期变化相关的物理和几何参数，例如太阳角度、地理位置、季节以及时间等已知参数。
本实例采用具有重叠度的不同轨道的错位关系探测长周期变化模型的方法，该方法不依赖地面控制点，而是对具有重叠度的不同轨道之间的影像进行同名点匹配，根据同名点与长周期变化模型的关系构建函数关系。
①对具有重叠度的不同轨道之间的影像进行同名点匹配，假设前一轨为j，则后一轨为j+1，两轨分别对应的长周期变化模型为：
②以每轨的影像景为单位，借助于SRTM30米格网数据，将j+1轨上的当前景影像经过模型正算算到地面上作为控制点，然后根据式(16)计算得到两景之间的相对安装值Δω、Δκ，则此值满足如下关系：
③对具有重叠区域的景分别构建式(20)，联立相对安装值的方程得到方程组，求解方程组得到式(17)的各个绝对值。这里需要注意的是，由于式(20)为相对模型，故组成的方程组实际为秩亏方程。为了解决这个问题，可以对当中某几景影像采用式(16)直接求取外方位元素误差值，代入式(20)中以解决秩亏问题。相对于传统方法需要所有轨都需要控制，此方法仅仅需要若干景有控制即可。
④对上述备选参数进行显著性分析以删选相关参数；
⑤采用支持向量机等方法根据相关参数进行函数拟合，得到模型的函数表达式。
图8中给出了前、后、正视影像的定位残差图，其中(a)、(c)、(e)代表解求常量偏置后的定位残差，而(b)、(d)、(f)代表解求顾及误差时间特性的偏置矩阵后的定位残差。由图可知，采用顾及误差时间特性的偏置矩阵比采用常量偏置矩阵的精度高，能够很好地吸收相机安装等的外方位元素误差。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (5)
1.一种星载推扫式光学传感器外方位元素定标方法，其特征在于：包括如下步骤：
步骤1，控制点量测；
遥感卫星应用需要控制点；采用控制点影像库的卫星遥感影像控制点的半自动选取技术，准确而又快速地量测出控制点，且能够保证控制点影像量测的精度；半自动选取技术主要包括控制点影像库的建立与应用；
步骤1.1，建立控制点数据库；
建立地面控制点影像库首先需要在地面应用系统的数据库中创建有关地面控制点列表(TABLE)，然后通过提取地面控制点，将有关控制点的信息存入数据库中，得到地面控制点影像库，地面控制点影像库表结构如表1；
表1地面控制点影像库表结构
提取地面控制点，建立控制点数据库的流程为：(1)在经扫描纠正的地形图上利用选取控制点的工具选取地面控制点；(2)在原始的卫星遥感影像上找到同名点，并选择同名点周围的邻域，即窗口图像，将所述窗口图像的元数据和影像存入控制点数据库中；重复此操作，控制点数据库在几何纠正处理过程中逐渐建立和丰富，最终形成完整的、覆盖全球的控制点数据库；
若发生变化的地面特征为控制点数据库中控制点所对应的地面区域，则需要及时更新控制点库；
步骤1.2，调用控制点数据库；
由于控制点数据库的不断更新，针对一个地物在数据库中对应着多个同名控制点的信息的情况；当进行卫星遥感影像几何纠正时，系统为了快速从控制点数据库中检索到最合适的控制点，遵循下述原则：选择与待纠正影像采集时间最接近的控制点进行自动匹配；
步骤1.3，相似性测度；
采用相关系数作为步骤1.2影响匹配的匹配测度，相关系数的定义为：
式中，
其中，m、n为模板大小，且N＝m*n，相关系数的取值在[-1，1]，越接近1则影像间相似程度越高；(i，j)为目标区中的像元行列号，(c，r)为搜索区中心的坐标；gi，j为待纠正影响控制点(i，j)的灰度值，g′i，j为影像数据库中与待纠正影像采集时间最接近的控制点的灰度值；gi+r，j+c为待纠正影响控制点(i+r，j+c)的灰度值，g′i+r，j+c为影像数据库中与待纠正影像采集时间最接近的控制点的灰度值；
步骤1.4，消除影像库中影像的旋转和比例畸变；
影像的几何变形主要表现为旋转关系和平移关系的组合，比例关系也涵盖在旋转关系中，而平移关系能够由影像匹配来解决，因此通过旋转、缩放控制点影像块的方式来解决旋转和比例畸变带来的问题；首先根据严密成像模型或RPC模型计算控制点影像的北方向与原始影像沿着轨道方向的差异，估计卫星的旋转角度和比例系数，方法如下：
分别建立原始影像和控制点影像的影像坐标OXY、oxy，在控制点影像的影像坐标上选取两个点P1(x0，y0)、P2(x1，y1)，并使P1和P2的连线与oy轴平行；通过反转换，得出P1点在原始影像的影像坐标上对应的点为P1′(X0，Y0)，P2点在原始影像的影像坐标上对应的点为P2′(X1，Y1′)；过P1′作OY轴的平行线与过P2′作OX轴的平行线交与一点P3′(X2，Y2)；此时线段P1′P2′与线段P2′P3′的夹角称为旋转角，令旋转角为：
比例系数m为：
式中，P′2P′3＝X0-X1；P′2P′3＝Y1-Y0；P1(x0,y0)为控制点影像中心，P2(x1,y1)为控制点影像块下边缘中点；P′1(X0,Y0)和P′2(X1,Y′1)是控制点影像上的对应点利用变换模型投影到原始影像上的投影点，D(P1，P2)和D(P′1，P2′)分别为两点之间的像素距离；
然后将控制点影像模板绕卫星本体坐标的质心进行旋转，在旋转过程中，加入比例系数改正，在旋转后的模板影像当中取一个内接的矩形区域作为新的控制点影像模板，此时，需保证控制点数据库中的控制点仍然在新的模板中心；控制点数据库中的控制点根据新的控制点影像模板生成新的控制点影像；即完成了旋转和比例畸变的消除，生成新的控制点影像；
步骤1.5，采用改进的爬山法搜索策略实现快速控制点匹配，即对步骤1.4生成的新的控制点影像的控制点匹配；
步骤1.5.1：利用严密成像模型或RPC模型，根据匹配点预测的最大视差确定搜索窗口大小；
步骤1.5.2：利用严密成像模型或RPC模型的预测点位精度确定金字塔匹配的层数；
首先，由于生成金字塔时采用低通滤波，使得高层上的影像保留原始影像本身大的整体结构特征，滤掉局部非平稳信息和噪声，从而在高层上减少匹配的不确定性；其次，高层匹配结果能够作为初值向下传递，减少下一层匹配的搜索范围，降低匹配过程中的计算复杂度；一幅影像既含有高频信息又含有低频信息；高频信息在影像中占的比例小，低频信息占的比例大；利用金字塔匹配模式，获得一个接近控制点实际位置的初始值，然后再在原始分辨率下计算影像相关系数，若相关系数达到阈值，确定金字塔匹配的层数；
步骤1.5.3：改进的爬山法搜索策略；
在每一层金字塔的搜索过程中，采用改进的爬山法搜索策略，以减少计算相关系数的次数，并且保证搜索的极值点为整体搜索区域内的极值点；
所述改进的爬山法搜索策略为：在任意一个空间直角坐标系中，预测点P0位于非主峰的区域内，根据爬山法的原理，找到一个峰顶S1，在S1位置沿着X方向就是P0A方向做一次整体搜索，找到X方向的整体峰值点P1，然后沿着X方向就是P1B方向做整体搜索，搜索到Y方向的整体峰值点P2，然后将峰顶S1与整体峰值点P1、整体峰值点P2对比得出最高峰值，最高峰值为整体搜索区域内的极值点；
步骤2，计算最小二乘偏置矩阵；
步骤2.1，建立常量偏置矩阵模型；
设备安装误差与姿轨测量误差等效，因此仅需根据姿轨测量误差特性构建补偿模型；对于视场角较小的高分光学卫星而言，轨道位置误差引起的几何定位误差为平移误差，与俯仰角误差、滚动角误差具有等效性；对沿卫星飞行方向和垂直卫星飞行方向对卫星的轨道位置误差进行分析可知，轨道位置误差与姿态误差的等效性；卫星真实位置与卫星带误差位置之间的距离能够转换成俯仰角误差
将平移误差等效为姿态角误差，采用偏置矩阵Ru补偿姿态误差，修正真实光线指向与带误差光线指向间的偏差，则由式(5)得常量偏置矩阵模型式(6)：
式中，为t时刻GPS相位中心在WGS84坐标系下的位置矢量；为t时刻J2000坐标系相对于WGS84坐标系的转换矩阵；为t时刻body坐标系相对于J2000坐标系的转换矩阵；为camera坐标系相对body坐标系的转换矩阵；(Dx Dy Dz)T为GPS相位中心在本体坐标系的坐标；为相机坐标系原点与本体坐标系原点偏移；代表相机主距；Ψ为CCD阵列沿轨向偏航角；为影像坐标；为主视轴即过主点垂直于CCS线阵的垂点对应的位置；为探元大小；m为比例系数；
Ru为正交旋转矩阵，Ru的作用是使得与之相乘的坐标系绕坐标系本身的y轴、x轴、z轴旋转角度ωu、κu；即：
则式(6)写成：
令：
则：
显然，(Xb Yb Zb)T是由地面点坐标确定的光线本体系下的指向；而(xb xb xb)T是由像方坐标确定的光线本体系下的指向；Ru用于修正两者的偏差从而实现姿态误差补偿；
展开式(9)，变化为：
其中：
对式(11)进行线性化并构建误差方程：
v误差＝A误差x误差-L误差；p误差 (12)
其中，x误差为L误差为根据初值计算的p误差为观测权值，A误差为系数矩阵具体为：
则：
偏置矩阵中待求未知数为三个偏置角，而一个平高控制点只能列两个方程；因此，已知两个控制点即可解求偏置矩阵；
步骤2.2，建立顾及误差时间特性的最小二乘偏置矩阵模型；
针对姿态系统误差引起的平移误差和旋转误差，常量偏置矩阵能很好的补偿；但是，由于平台稳定性不够以及陀螺漂移对姿态确定精度的影响，几何定位模型中常常存在着漂移误差；
常量偏置矩阵无法补偿姿态漂移误差，为此，在常量偏置矩阵中引入偏置角的时间变化率：
将式(14)带入公式(6)以替换Ru，得到最小二乘偏置矩阵模型：
2.如权利要求1所述的一种星载推扫式光学传感器外方位元素定标方法，其特征在于：还可以对步骤1.5求得的控制点进行精度匹配：
采用相邻域点相关系数拟合椭球抛物面来获得该亚像元的二维匹配精度；
二维匹配精度的方法为：以步骤1.4得到的新的控制点影像为目标窗口，在原始影像的搜索窗口中取和目标窗口相同大小的多个影像窗口进行运算，求得代表各个影像窗口中心象素的相关关系的相关系数；根据改进的爬山法搜索策略，对每个影像窗口进行搜索，得到各自的最大相关系数点；最大相关系数点就是所求的待匹配影像控制点的同名点；
当每个影像窗口的同名象素的中心点是真正的同名点时，通过计算同名点即可求得控制点匹配精度；
由于控制点影像模板与原始影像拍摄时间的差异，每个影像窗口的同名象素的中心点不是真正的同名点时；根据偏离情况，按照如下方法计算控制点匹配精度；
情况一、真正的同名点在两个方向上偏离象素中心半个象素之内，则产生了误差，所述误差的大小为象素级；显然，该项误差在两个方向分别均服从内的均匀分布，因此单方向的相关精度为其中
因此，σx＝0.29Δ
同理σy＝0.29Δ
因此即二维整象素匹配的理论精度为0.4个象素；
结合理论精度，将最大相关系数位置的邻域点中的相关系数与平差函数联系起来，从而将平差函数的极值点处作为寻求的同名点；
情况二、其他任何偏离情况均采用下述方法计算；
椭圆抛物面方程拟合式：
式中，(x，y)为5×5邻域内每点相对于中心的位置；z为控制点影像块的相关系数；p、q为椭圆抛物面的参数，pq＞0；p、q、x0和y0用间接观测平差求得，当p≠q时，X和Y两个方向的匹配精度不一样，当p＝q时，椭圆抛物面退化为旋转抛物面，两个方向的匹配精度一致，此时椭圆抛物面顶点(x0,y0)作为寻求的亚像元匹配的同名点；
由式(4)求得的椭圆抛物面的顶点即为亚像元精度的同名点；原因如下：同名点与周围5×5邻域符合椭圆抛物面，因此根据同名点以及周围5×5的范围内相关系数，结合式(4)求得同名点。
3.如权利要求1所述的一种星载推扫式光学传感器外方位元素定标方法，其特征在于：对外方位元素定标结果进行校验，即对步骤2的结果进行校验，方法如下：
评估偏置矩阵的精度，以及引入偏置角的时间变化率偏置矩阵的精度；
步骤3.1，偏置补偿矩阵RU；
考虑到偏置补偿矩阵，式(5)变化为：
RU表示待消除的外方位元素误差；所述外方位元素误差包括轨道误差、姿态误差和安装误差；RU定义成式(17)：
其中，ω、κ分别是相对于X、Y、Z轴旋转的欧拉角；偏置补偿矩阵补偿的是一景短周期范围内的外方位元素误差；
步骤3.2，估计长周期变化模型；
上述偏置补偿矩阵仅仅能对一景时间范围内的外方位元素误差进行建模；实际上，由于卫星姿态精度受到低频漂移误差和结构变形的而影响，存在一个长周期的变化，其变化模型目前来说仍是未知的；目前提高无控定位精度的方法为大区域网平差方法，但由于长周期变化模型的存在，会破坏区域网平差模型中的正态分布，所以不适用于长周期变化；
因此需要对外方位元素的长周期变化模型估计，长周期变化模型表示如下：
其中，θ1,…,θn,t为备选的与长周期变化相关的物理和几何参数，所述物理参数包括：太阳角度、季节以及时间；所述几何参数包括：地理位置；
采用具有重叠度的不同轨道的错位关系探测长周期变化模型的方法，对具有重叠度的不同轨道之间的影像进行同名点匹配，根据同名点与长周期变化模型的关系构建函数关系；
a)对具有重叠度的不同轨道之间的影像进行同名点匹配，当前一轨为j，则后一轨为j+1，两轨分别对应的长周期变化模型为：
b)以每轨的影像景为单位，借助于SRTM30米格网数据，将j+1轨上的当前景影像经过模型正算得到地面控制点，然后根据式(17)计算得到两景之间的相对安装值Δω、Δκ，则此值满足如下关系：
c)对具有重叠区域的景分别构建式(21)，将式(21)与相对安装值的方程联立得到方程组，求解方程组得到式(18)的各个绝对值；需要注意的是，由于式(21)为相对模型，故组成的方程组实际为秩亏方程；为了解决这个问题，对当中任意几景影像采用式(17)直接求取外方位元素误差值，误差值代入式(21)中以解决秩亏问题；
d)将上述(c)中的由任意几景影像求取的外方位元素误差值构成备选参数组，对参数组中的参数进行显著性分析，选出精度高的参数；
采用支持向量机等方法根据精度够高的参数进行函数拟合，得到长周期变化模型的函数表达式；通过长周期变化模型的函数表达式即可评估偏置矩阵的精度，并获得引入偏置角的时间变化率偏置矩阵的精度。
4.如权利要求1所述的一种星载推扫式光学传感器外方位元素定标方法，其特征在于：所述更新控制点库的方法如下；
(1)当地面特征发生明显变化时，及时更新控制点库；
(2)当地面特征发生非明显变化时，通过预先设置阈值，然后将观察日志文件中控制点的匹配相关系数与阈值进行对比，当达到阈值时，及时更新控制点库。
5.如权利要求1所述的一种星载推扫式光学传感器外方位元素定标方法，其特征在于：所述选择与待纠正影像采集时间最接近的控制点进行自动匹配的具体流程如下：从控制点数据库中检索位于待纠正景或启算景影像范围内的控制点经纬度以及控制带影像拍摄条件；若检索的结果中有相同经纬度的控制点，则读取控制点的拍摄时间与待纠正景或启算景影像时间最接近的控制点数据；若检索的结果中无相同经纬度的控制点，则只读取当前控制点数据；根据以上的读取结果判断能否获取读取范围内所有控制点；若可以，则用控制点影像与原始卫星遥感影像匹配；若不可以，则重复本步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910351621.7A CN110111274B (zh) | 2019-04-28 | 2019-04-28 | 一种星载推扫式光学传感器外方位元素定标方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910351621.7A CN110111274B (zh) | 2019-04-28 | 2019-04-28 | 一种星载推扫式光学传感器外方位元素定标方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110111274A true CN110111274A (zh) | 2019-08-09 |
| CN110111274B CN110111274B (zh) | 2020-06-19 |
Family
ID=67487264
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910351621.7A Active CN110111274B (zh) | 2019-04-28 | 2019-04-28 | 一种星载推扫式光学传感器外方位元素定标方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110111274B (zh) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111161186A (zh) * | 2019-12-30 | 2020-05-15 | 国家卫星气象中心(国家空间天气监测预警中心) | 一种推扫式遥感器通道配准方法及装置 |
| CN111272196A (zh) * | 2020-02-29 | 2020-06-12 | 武汉大学 | 特定拍摄条件下的在轨外方位元素自检校方法及系统 |
| CN111637872A (zh) * | 2020-05-14 | 2020-09-08 | 中国人民解放军空军预警学院 | 一种光学微小卫星可见光与红外双波段亚像元成像装置 |
| CN111665218A (zh) * | 2020-05-21 | 2020-09-15 | 武汉大学 | 一种提高二氧化碳差分吸收激光雷达反演精度的方法 |
| CN113096163A (zh) * | 2021-04-16 | 2021-07-09 | 湘潭大学 | 一种无需先验升降轨信息的星载sar影像高精度匹配方法 |
| CN114241064A (zh) * | 2022-02-24 | 2022-03-25 | 中国科学院空天信息创新研究院 | 一种遥感卫星内外方位元素实时几何定标方法 |
| CN117928494A (zh) * | 2024-03-18 | 2024-04-26 | 中国人民解放军战略支援部队航天工程大学 | 一种光学卫星切片影像的几何定位测量方法、系统和设备 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103148870A (zh) * | 2013-03-01 | 2013-06-12 | 国家测绘地理信息局卫星测绘应用中心 | 基于高精度配准控制点的卫星ccd阵列影像几何检校方法 |
| CN103398701A (zh) * | 2013-07-31 | 2013-11-20 | 国家测绘地理信息局卫星测绘应用中心 | 一种基于物方投影面的星载非共线tdi ccd影像拼接方法 |
| CN103927744A (zh) * | 2014-03-27 | 2014-07-16 | 航天东方红卫星有限公司 | 一种基于指向姿态的敏捷卫星观测目标条带分割方法 |
| CN107014399A (zh) * | 2017-04-25 | 2017-08-04 | 国家测绘地理信息局卫星测绘应用中心 | 一种星载光学相机‑激光测距仪组合系统联合检校方法 |
| US10151588B1 (en) * | 2016-09-28 | 2018-12-11 | Near Earth Autonomy, Inc. | Determining position and orientation for aerial vehicle in GNSS-denied situations |
-
2019
- 2019-04-28 CN CN201910351621.7A patent/CN110111274B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103148870A (zh) * | 2013-03-01 | 2013-06-12 | 国家测绘地理信息局卫星测绘应用中心 | 基于高精度配准控制点的卫星ccd阵列影像几何检校方法 |
| CN103398701A (zh) * | 2013-07-31 | 2013-11-20 | 国家测绘地理信息局卫星测绘应用中心 | 一种基于物方投影面的星载非共线tdi ccd影像拼接方法 |
| CN103927744A (zh) * | 2014-03-27 | 2014-07-16 | 航天东方红卫星有限公司 | 一种基于指向姿态的敏捷卫星观测目标条带分割方法 |
| US10151588B1 (en) * | 2016-09-28 | 2018-12-11 | Near Earth Autonomy, Inc. | Determining position and orientation for aerial vehicle in GNSS-denied situations |
| CN107014399A (zh) * | 2017-04-25 | 2017-08-04 | 国家测绘地理信息局卫星测绘应用中心 | 一种星载光学相机‑激光测距仪组合系统联合检校方法 |
Non-Patent Citations (5)
| Title |
|---|
| PETER S 等: "In-flight quality and accuracy of attitude measurements from the CHAMP advanced stellar compass", 《ACTA ASTRONAUTICA》 * |
| XIAOHUATONG 等: "Detection and estimation of ZY-3 three-line array image distortions caused by attitude oscillation", 《ISPRS JOURNAL OF PHOTOGRAMMETRY AND REMOTE SENSING》 * |
| 张过: "缺少控制点的高分辫率卫星遥感影像几何纠正", 《中国优秀博硕士学位论文全文数据库(博士)基础科学辑》 * |
| 李德仁 等: "国产光学卫星影像几何精度研究", 《航天器工程》 * |
| 蒋永华 等: "线阵推扫光学卫星外方位元素自检校方法", 《同济大学学报（自然科学版）》 * |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111161186A (zh) * | 2019-12-30 | 2020-05-15 | 国家卫星气象中心(国家空间天气监测预警中心) | 一种推扫式遥感器通道配准方法及装置 |
| CN111161186B (zh) * | 2019-12-30 | 2023-10-03 | 国家卫星气象中心(国家空间天气监测预警中心) | 一种推扫式遥感器通道配准方法及装置 |
| CN111272196A (zh) * | 2020-02-29 | 2020-06-12 | 武汉大学 | 特定拍摄条件下的在轨外方位元素自检校方法及系统 |
| CN111637872A (zh) * | 2020-05-14 | 2020-09-08 | 中国人民解放军空军预警学院 | 一种光学微小卫星可见光与红外双波段亚像元成像装置 |
| CN111665218A (zh) * | 2020-05-21 | 2020-09-15 | 武汉大学 | 一种提高二氧化碳差分吸收激光雷达反演精度的方法 |
| CN111665218B (zh) * | 2020-05-21 | 2021-07-02 | 武汉大学 | 一种提高二氧化碳差分吸收激光雷达反演精度的方法 |
| CN113096163A (zh) * | 2021-04-16 | 2021-07-09 | 湘潭大学 | 一种无需先验升降轨信息的星载sar影像高精度匹配方法 |
| CN114241064A (zh) * | 2022-02-24 | 2022-03-25 | 中国科学院空天信息创新研究院 | 一种遥感卫星内外方位元素实时几何定标方法 |
| CN117928494A (zh) * | 2024-03-18 | 2024-04-26 | 中国人民解放军战略支援部队航天工程大学 | 一种光学卫星切片影像的几何定位测量方法、系统和设备 |
| CN117928494B (zh) * | 2024-03-18 | 2024-05-24 | 中国人民解放军战略支援部队航天工程大学 | 一种光学卫星切片影像的几何定位测量方法、系统和设备 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110111274B (zh) | 2020-06-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110111274A (zh) | 一种星载推扫式光学传感器外方位元素定标方法 | |
| Li et al. | Rigorous photogrammetric processing of HiRISE stereo imagery for Mars topographic mapping | |
| CN108305237B (zh) | 考虑不同光照成像条件的多立体影像融合制图方法 | |
| CN109903352A (zh) | 一种卫星遥感影像大区域无缝正射影像制作方法 | |
| CN105046251B (zh) | 一种基于环境一号卫星遥感影像的自动正射校正方法 | |
| JP5389964B2 (ja) | 地図情報生成装置 | |
| CN102735216B (zh) | Ccd立体相机三线阵影像数据平差处理方法 | |
| CN107144293A (zh) | 一种视频卫星面阵相机的几何定标方法 | |
| CN102410831B (zh) | 多条带扫描成像模型的设计及定位方法 | |
| CN102968631A (zh) | 山区多光谱遥感卫星影像的自动几何纠正与正射校正方法 | |
| CN106127697A (zh) | 无人机机载成像高光谱几何校正方法 | |
| Wang et al. | Pictometry’s proprietary airborne digital imaging system and its application in 3D city modelling | |
| CN111473802A (zh) | 一种基于线阵推扫的光学传感器内方位元素定标方法 | |
| CN105701830A (zh) | 基于几何模型的lasis波段影像配准方法及系统 | |
| CN110986888A (zh) | 一种航空摄影一体化方法 | |
| Tao et al. | On-orbit geometric calibration of the panchromatic/multispectral camera of the ZY-1 02C satellite based on public geographic data | |
| CN107705272A (zh) | 一种空间影像的高精度几何校正方法 | |
| CN116824079A (zh) | 基于全信息摄影测量的三维实体模型构建方法和装置 | |
| CN105093222A (zh) | 一种sar影像区域网平差连接点自动提取方法 | |
| Whitehead | An integrated approach to determining short-term and long-term patterns of surface change and flow characteristics for a polythermal arctic glacier | |
| CN114777745A (zh) | 一种基于无迹卡尔曼滤波的倾斜取证建模方法 | |
| Grün et al. | Automatic reconstruction and visualization of a complex Buddha Tower of Bayon, Angkor, Cambodia | |
| Reinartz et al. | Accuracy analysis for DEM and orthoimages derived from SPOT HRS stereo data without using GCP | |
| Lee et al. | Automatic building reconstruction with satellite images and digital maps | |
| Cao et al. | Precise sensor orientation of high-resolution satellite imagery with the strip constraint |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |