CN102636159A - 多镜头航天线阵相机***在轨几何自检校方法 - Google Patents

多镜头航天线阵相机***在轨几何自检校方法 Download PDF

Info

Publication number
CN102636159A
CN102636159A CN2012101131294A CN201210113129A CN102636159A CN 102636159 A CN102636159 A CN 102636159A CN 2012101131294 A CN2012101131294 A CN 2012101131294A CN 201210113129 A CN201210113129 A CN 201210113129A CN 102636159 A CN102636159 A CN 102636159A
Authority
CN
China
Prior art keywords
delta
overbar
centerdot
omega
kappa
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
Application number
CN2012101131294A
Other languages
English (en)
Other versions
CN102636159B (zh
Inventor
方勇
江振治
胡海彦
陈虹
王刃
马永社
张丽
杨韫澜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
SURVEYING AND MAPPING INST HEADQUARTERS OF GENERAL STAFF CPLA
Original Assignee
SURVEYING AND MAPPING INST HEADQUARTERS OF GENERAL STAFF CPLA
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by SURVEYING AND MAPPING INST HEADQUARTERS OF GENERAL STAFF CPLA filed Critical SURVEYING AND MAPPING INST HEADQUARTERS OF GENERAL STAFF CPLA
Priority to CN201210113129.4A priority Critical patent/CN102636159B/zh
Publication of CN102636159A publication Critical patent/CN102636159A/zh
Application granted granted Critical
Publication of CN102636159B publication Critical patent/CN102636159B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Processing (AREA)

Abstract

本发明涉及一种多镜头航天线阵相机***在轨几何自检校方法,用于挖掘和提升传感器的量测性能,从而满足航天线阵相机摄影测量定位精度及军事地形图测制的实际需求。其基本思想是首先在检校试验区域布设和测量高精度的地面标志控制点,并获取相机***的地面几何标称参数。然后,获取卫星在轨飞行的轨道数据,并在传感器***所摄对应试验区域影像上量测标志控制点对应的像点坐标;最后基于自建的多镜头线阵传感器成像模型和轨道模型,采用最小二乘光束法自检校平差技术推求影响相机***几何定位精度的参数改正量,从而,提升传感器***的量测性能。实验表明本发明具有较好的可靠性和较高的精度,对提升传感器的对地定位几何精度效果明显。

Description

多镜头航天线阵相机***在轨几何自检校方法
技术领域
本发明涉及一种卫星在轨飞行状态下,其所搭载的多镜头航天线阵传感器几何***性误差标定与改正方法,属于航天摄影测量与遥感测绘技术领域。
背景技术
线阵扫描式传感器按照镜头多少可以分为单镜头和多镜头传感器,为了进行立体交会,通常情况下,单镜头传感器会在焦面上布置多个平行的线阵CCD,而多镜头则会在每个单独的相机的焦面上布置一条(组)线阵CCD。这种结构设计上的差异,决定了传感器的***误差以及在轨检校无论是在几何成像模型上还是在数据处理方法上都有所不同。
近年国产航天摄影测量线阵相机基本都以多镜头设计为主,例如“天绘”及“资源”系列等国产航天摄影测量卫星***均属于多镜头线阵传感器这一类。相机***一旦上天,处于在轨飞行状态下,由于太空中的温湿度、重力、压力等与地面差别极大,这势必引起相机***的地面标称参数发生很大变化,从而直接影响传感器***的对地定位能力。多镜头航天线阵相机***在轨自检校的意义在于其是在传感器实际的航天摄影过程中进行标定,标定的参数更能够满足摄影测量实际作业的需要,利用多余的控制点,可一并对被标定参数的正确性及自检校对传感器摄影定位能力的改善程度做出检验。
在航天线阵传感器的摄影定位研究方面,欧美国家的技术比较领先,无论是传感器硬件制造工艺,还是传感器在轨标定方法研究,都存在明显的技术优势。例如,法国的SPOT系列卫星,美国的IKONOS-2、Quickbird等卫星线阵传感器***都具备高稳定度、高精度对地定位能力。高精度的对地定位能力不仅建立在传感器硬件工艺水平之上、而且卫星在轨飞行运控及传感器在轨检校等技术环节对传感器的高精度定位能力同样重要,在实际应用中,必须将这些技术环节进行有机结合才能达到高精度定位的目的。而传感器在轨检校技术在所有的技术环节中显得尤为重要,因为其直接面向应用,衔接摄影测量作业、并决定测绘产品的质量。出于技术参数保密和军事利益、经济效益的考虑,这些国外的传感器***对商业用户所提供的技术参数均经过加密处理,并且传感器在轨检校技术细节也从未公开。目前国内在这方面的研究几乎为空白,未见相似或同样的方法出现。
发明内容
本发明需要解决技术问题是提供一种能够用于多镜头航天线阵相机***在轨几何自检校方法,挖掘和提升传感器的量测性能,从而满足航天线阵相机摄影测量定位精度及军事地形图测制的实际需求。
为解决上述技术问题,本发明提出一种多镜头航天线阵相机***在轨几何自检校方法,其基本思想是首先在检校试验区域布设和测量高精度的地面标志控制点,并获取相机***的地面几何标称参数。然后,获取卫星在轨飞行的轨道数据(包括GPS位置数据和星敏姿态数据),并在传感器***所摄对应试验区域影像上量测标志控制点对应的像点坐标;最后基于自建的成像模型和轨道模型,采用最小二乘光束法自检校平差(Least Squares Self-calibration Bundle Adiustment)技术推求影响相机***几何定位精度的参数改正量,从而,提升传感器***的量测性能,即提升传感器的对地定位几何精度。
本发明的技术解决方案是:
一种多镜头航天线阵相机***在轨几何自检校方法,其特征在于:该方法包括以下几个步骤:
(1)根据传感器飞行轨迹和影像条带宽度,在检校区内布设高精度标志控制点;
(2)对获取检校区的前、下、后视相机三个影像条带进行控制点对应像点量测;
(3)获取传感器的轨道数据,包括GPS及星敏姿态数据观测量;
(4)根据多镜头航天线阵传感器成像模型和外部定向参数模型,建立观测方程组;
<1>、多镜头航天线阵传感器成像模型
光学***由多个镜头组成的传感器,必须要考虑描述其它镜头相对于下视镜头(平台)的位置与姿态关系,对于镜头j,用dxj,dyj,dzj表示相对于下视镜头的相对位置,用αj,βj,γj表示相对于下视镜头的相对姿态,fj,xpj,ypj为镜头j的焦距和主点坐标,成像模型为:
x - x pj = - f j &CenterDot; s 11 ( X - X C ) + s 21 ( Y - Y C ) + s 31 ( Z - Z C ) - ( m 11 d xj + m 21 d yj + m 31 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj ) + &Delta; xj
y - y pj = - f j &CenterDot; s 12 ( X - X C ) + s 22 ( Y - Y C ) + s 32 ( Z - Z C ) - ( m 12 d xj + m 22 d yj + m 32 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj ) + &Delta; yj - - - ( 1 )
该模型对单镜头与多镜头都是适用的,其中:j=1,2,3,…,n,
[X,Y,Z]:物点在地面坐标***中的坐标;
[XC,YC,ZC]:下视相机投影中心在地面坐标***中的坐标;
[x,y]:像点在相机坐标***中的坐标,理论上x为0或某一常数;
[xp,yp]:像主点在相机坐标***中的坐标;
f:焦距;
从下视相机坐标***到地面坐标***的旋转矩阵,旋角分别为
M(αj,βj,γj)为从下视镜头相机坐标***到镜头j相机坐标***的旋转矩阵;
Figure BDA0000154445880000035
为从镜头j相机坐标***到地面坐标***的旋转矩阵;
并且,Δxj,Δyj是关于镜头畸变、CCD线阵畸变以及内定向参数(改正量)成像***的***误差模型,如下:
&Delta; x j = &Delta; x pj - &Delta;f f x &OverBar; pj + ( k 1 j r j 2 + k 2 j r j 4 ) x &OverBar; pj + p 1 j ( r j 2 + 2 x &OverBar; pj 2 ) + 2 p 2 j x &OverBar; pj y &OverBar; pj + x &OverBar; pj sin &theta;
&Delta; y j = &Delta; y pj - &Delta;f f y &OverBar; pj + ( k 1 j r j 2 + k 2 j r j 4 ) y &OverBar; pj + + 2 p 1 j x &OverBar; pj y &OverBar; pj + p 2 j ( r j 2 + 2 y &OverBar; pj 2 ) + s yj y &OverBar; pj - - - ( 2 )
参数说明:
·每个镜头的主点改正量(位移量)Δxp,Δyp
·k1,k2,p1,p2:每个镜头的对称性径向和切向畸变;
·sy:每个线阵CCD单元y方向上的比例因子;
·每个CCD线阵在焦面上的旋角θ(其产生的畸变主要在x方向,在y方向上的畸变影响可以忽略不计), x &OverBar; pj = x - x pj , y &OverBar; pj = y - y pj , r 2 = x &OverBar; pj 2 + y &OverBar; pj 2 ;
<2>、外部定向参数模型
外部定向参数模型——轨道模型(trajectory model)是用来描述传感器外部定向(EO,external orientation)状态的,其是用关于时间的分段多项式(PPM,PiecewisePolynomial Function)来建模的,可以根据控制点(GCPs)和连接点(TPs)的分布来确定轨道模型的分段数,对于被分割轨道的第i段,该段的首末端点的线阵影像获取时刻记为
Figure BDA0000154445880000046
则定义时间变量
Figure BDA0000154445880000047
为:
t &OverBar; = l - l ini i l fin i - l ini i - - - ( 3 )
其中,
Figure BDA0000154445880000049
为第i段轨道的首末端在对应影像中的线阵行数(首末线阵影像行),这样,对于第i段轨道所对应影像中的第1行线阵影像,其位置与姿态
Figure BDA00001544458800000410
的数学模型可以表示为该行线阵影像的外部定向参数观测量
Figure BDA00001544458800000411
与关于时间的2阶多项式的和,如下式:
X C ( t &OverBar; ) = X instr + X 0 i + X 1 i t &OverBar; + X 2 i t &OverBar; 2
Y C ( t &OverBar; ) = Y instr + Y 0 i + Y 1 i t &OverBar; + Y 2 i t &OverBar; 2
(4)
Z C ( t &OverBar; ) = Z instr + Z 0 i + Z 1 i t &OverBar; + Z 2 i t &OverBar; 2
&omega; C ( t &OverBar; ) = &omega; instr + &omega; 0 i + &omega; 1 i t &OverBar; + &omega; 2 i t &OverBar; 2
Figure BDA0000154445880000055
&kappa; C ( t &OverBar; ) = &kappa; instr + &kappa; 0 i + &kappa; 1 i t &OverBar; + &kappa; 2 i t &OverBar; 2
其中,[X0,X1,X2,…,κ0,κ1,κ2]i是第i段轨道关于GPS偏移量及INS(星敏)轴线偏差、GPS及INS(星敏)***漂移等***误差的18个模型参数;
<3>、观测方程组
下面将详细列出关于上述传感器模型的观测方程组——包括像坐标观测方程、外部定向参数观测方程(连续性条件方程及外部定向参数自身观测方程)、自检校参数观测方程、地面点观测方程,其中方程组中的所有观测量都带有上划线记号,待平差的定向参数无上划线标记:
x &OverBar; + v x = x = x pj - f j &CenterDot; N x D + &Delta; xj - - - ( 5 )
y &OverBar; + v y = y = y pj - f j &CenterDot; N y D + &Delta; yj
——像坐标观测量方程
( X 0 i &OverBar; + X 1 i &OverBar; + X 2 i &OverBar; - X 0 i + 1 &OverBar; ) + v X C 0 = X 0 i + X 1 i + X 2 i - X 0 i + 1
( Y 0 i &OverBar; + Y 1 i &OverBar; + Y 2 i &OverBar; - Y 0 i + 1 &OverBar; ) + v Y C 0 = Y 0 i + Y 1 i + Y 2 i - Y 0 i + 1
( Z 0 i &OverBar; + Z 1 i &OverBar; + Z 2 i &OverBar; - Z 0 i + 1 &OverBar; ) + v Z C 0 = Z 0 i + Z 1 i + Z 2 i - Z 0 i + 1
( &omega; 0 i &OverBar; + &omega; 1 i &OverBar; + &omega; 2 i &OverBar; - &omega; 0 i + 1 &OverBar; ) + v &omega; C 0 = &omega; 0 i + &omega; 1 i + &omega; 2 i - &omega; 0 i + 1
Figure BDA00001544458800000513
( &kappa; 0 i &OverBar; + &kappa; 1 i &OverBar; + &kappa; 2 i &OverBar; - &kappa; 0 i + 1 &OverBar; ) + v &kappa; C 0 = &kappa; 0 i + &kappa; 1 i + &kappa; 2 i - &kappa; 0 i + 1 - - - ( 6 )
——PPM连续性约束条件(0阶)
( X 1 i &OverBar; + 2 X 2 i &OverBar; - X 1 i + 1 &OverBar; ) + v X C 1 = X 1 i + 2 X 2 i - X 1 i + 1
( Y 1 i &OverBar; + 2 Y 2 i &OverBar; - Y 1 i + 1 &OverBar; ) + v Y C 1 = Y 1 i + 2 Y 2 i - Y 1 i + 1
( Z 1 i &OverBar; + 2 Z 2 i &OverBar; - Z 1 i + 1 &OverBar; ) + v Z C 1 = Z 1 i + 2 Z 2 i - Z 1 i + 1
( &omega; 1 i &OverBar; + 2 &omega; 2 i &OverBar; - &omega; 1 i + 1 &OverBar; ) + v &omega; C 1 = &omega; 1 i + 2 &omega; 2 i - &omega; 1 i + 1
Figure BDA0000154445880000065
( &kappa; 1 i &OverBar; + 2 &kappa; 2 i &OverBar; - &kappa; 1 i + 1 &OverBar; ) + v &kappa; C 1 = &kappa; 1 i + 2 &kappa; 2 i - &kappa; 1 i + 1 - - - ( 7 )
——PPM连续性约束条件(1阶)
( X 2 i &OverBar; - X 2 i + 1 &OverBar; ) + v X C 2 = X 2 i - X 2 i + 1
( Y 2 i &OverBar; - Y 2 i + 1 &OverBar; ) + v Y C 2 = Y 2 i - Y 2 i + 1
( Z 2 i &OverBar; - Z 2 i + 1 &OverBar; ) + v Z C 2 = Z 2 i - Z 2 i + 1
( &omega; 2 i &OverBar; - &omega; 2 i + 1 &OverBar; ) + v &omega; C 2 = &omega; 2 i - &omega; 2 i + 1
Figure BDA00001544458800000611
( &kappa; 2 i &OverBar; - &kappa; 2 i + 1 &OverBar; ) + v &kappa; C 2 = &kappa; 2 i - &kappa; 2 i + 1 - - - ( 8 )
——PPM连续性约束条件(2阶)
X 0 &OverBar; + v X 0 = X 0 Y 0 &OverBar; + v Y 0 = Y 0 Z 0 &OverBar; + v Z 0 = Z 0
X 1 &OverBar; + v X 1 = X 1 Y 1 &OverBar; + v Y 1 = Y 1 Z 1 &OverBar; + v Z 1 = Z 1 - - - ( 9 )
X 2 &OverBar; + v X 2 = X 2 Y 2 &OverBar; + v Y 2 = Y 2 Z 2 &OverBar; + v Z 2 = Z 2
&omega; 0 &OverBar; + v &omega; 0 = &omega; 0
Figure BDA00001544458800000623
&kappa; 0 &OverBar; + v &kappa; 0 = &kappa; 0
&omega; 1 &OverBar; + v &omega; 1 = &omega; 1
Figure BDA00001544458800000626
&kappa; 1 &OverBar; + v &kappa; 1 = &kappa; 1
&omega; 2 &OverBar; + v &omega; 2 = &omega; 2
Figure BDA00001544458800000629
&kappa; 2 &OverBar; + v &kappa; 2 = &kappa; 2
——外部定向参数自身观测方程
&Delta; x P &OverBar; + v &Delta; x P = &Delta; x P
&Delta; y P &OverBar; + v &Delta; y P = &Delta; y P
&Delta;f &OverBar; + v &Delta;f = &Delta;f
k 1 &OverBar; + v k 1 = k 1
k 2 &OverBar; + v k 2 = k 2
p 1 &OverBar; + v p 1 = p 1
p 2 &OverBar; + v p 2 = p 2
(10)
s y &OverBar; + v s y = s y
&theta; &OverBar; + v &theta; = &theta;
——自检校参数观测方程
X &OverBar; + v X = X
Y &OverBar; + v Y = Y - - - ( 11 )
Z &OverBar; + v Z = Z
——地面点观测方程
(5)利用最小二乘法进行光束法自检校平差数据处理,求解几何***误差检校参数估值——首先进行线性化、系数阵构造、法化后线性***求解;
(6)对所解得参数进行显著性检验,并进行数据分析,包括内部精度、RMS、相关系数和粗差探测计算。
所述的利用最小二乘法进行光束法自检校平差数据处理,求解几何***误差检校参数估值——首先进行线性化、系数阵构造、法化后线性***求解如下:
<1>、线性化方程组(偏导数)
对成像模型与轨道模型关于未知数求偏导,使之线性化,以便求解:
&delta;x &delta; X 0 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 0 &delta;x &delta; X 1 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 1 &delta;x &delta; X 2 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 2
&delta;x &delta; Y 0 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 0 &delta;x &delta; Y 1 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 1 &delta;x &delta; Y 2 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 2
&delta;x &delta; Z 0 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Z 0 &delta;x &delta; Z 1 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Y 1 &delta;x &delta; Z 2 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Y 2
&delta;x &delta; &omega; 0 = &delta;x &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 0 &delta;x &delta; &omega; 1 = &delta;x &delta;&omega; C &CenterDot; &delta; &omega; C &delta; &omega; 1 &delta;x &delta; &omega; 2 = &delta;x &delta;&omega; C &CenterDot; &delta;&omega; C &delta;&omega; 2
Figure BDA00001544458800000813
Figure BDA00001544458800000814
&delta;x &delta; &kappa; 0 = &delta;x &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 0 &delta;x &delta; &kappa; 1 = &delta;x &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 1 &delta;x &delta; &kappa; 2 = &delta;x &delta; &kappa; C &CenterDot; &delta;&kappa; C &delta; &kappa; 2 - - - ( 12 )
——成像模型中x关于外部定向参数的偏导数
&delta;y &delta; X 0 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 0 &delta;y &delta; X 1 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 1 &delta;y &delta; X 2 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 2
&delta;y &delta; Y 0 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 0 &delta;y &delta; Y 1 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 1 &delta;y &delta; Y 2 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 2
&delta;y &delta; Z 0 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Z 0 &delta;y &delta; Z 1 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Y 1 &delta;y &delta; Z 2 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Y 2
&delta;y &delta; &omega; 0 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 0 &delta;y &delta; &omega; 1 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 1 &delta;y &delta; &omega; 2 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 2
Figure BDA00001544458800000831
Figure BDA00001544458800000832
&delta;y &delta; &kappa; 0 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 0 &delta;y &delta; &kappa; 1 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 1 &delta;y &delta; &kappa; 2 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 2 - - - ( 13 )
——成像模型中y关于外部定向参数的偏导数
其中,
Figure BDA0000154445880000091
Figure BDA0000154445880000092
Figure BDA0000154445880000093
并且,
&PartialD; x &PartialD; X C = - f D 2 ( s 13 N x - s 11 D ) = a 1,1
&PartialD; x &PartialD; Y C = - f D 2 ( s 23 N x - s 21 D ) = a 1,2
&PartialD; x &PartialD; Z C = - f D 2 ( s 33 N x - s 31 D ) = a 1,3
&PartialD; x &PartialD; &omega; C = f D { [ ( X - X C ) &delta; s 13 &PartialD; &omega; C + ( Y - Y C ) &delta; s 23 &PartialD; &omega; C - ( Z - Z C ) &delta; s 33 &PartialD; &omega; C ] N x D - ( X - X 0 ) &delta; s 11 &PartialD; &omega; C - ( Y - Y 0 ) &delta; s 21 &PartialD; &omega; C - ( Z - Z 0 ) &delta; s 31 &PartialD; &omega; C } = a 1,4
Figure BDA0000154445880000098
&PartialD; x &PartialD; &kappa; C = f D { [ ( X - X C ) &delta;s 13 &delta;&kappa; C + ( Y - Y C ) &delta;s 33 &PartialD; &kappa; C ] N x D - ( X - X 0 ) &delta;s 11 &delta; &kappa; C - ( Y - Y 0 ) &delta;s 21 &PartialD; &kappa; C - ( Z - Z 0 ) &delta;s 31 &PartialD; &kappa; C } = a 1,6
&PartialD; y &PartialD; X C = - f N 2 ( r 13 Z y - r 13 N ) = a 2,1
&PartialD; y &PartialD; Y C = - f N 2 ( r 23 Z y - r 22 N ) = a 2 , 2
&PartialD; y &PartialD; Z C = - f N 2 ( r 33 Z y - r 32 N ) = a 2,3
&PartialD; y &PartialD; &omega; C = f N { [ ( Y - Y 0 ) r 33 - ( Z - Z 0 ) r 23 ] Z y N - ( Y - Y 0 ) r 32 - ( Z - Z 0 ) r 22 } = a 2,7
Figure BDA00001544458800000914
&PartialD; y &PartialD; &kappa; C = - f N Z x = a 2,9
系数ai,j指的是在矩阵AGCP和ATP对应的元素(AGCP、ATP定义见下文);
&PartialD; x &PartialD; X = - f D 2 ( Ds 11 - N x s 13 ) = b 1,1 &PartialD; y &PartialD; X = - f D 2 ( Ds 12 - N y s 13 ) = b 2,1 - - - ( 14 )
&PartialD; x &PartialD; Y = - f D 2 ( Ds 31 - N x s 23 ) = b 1,2 &PartialD; y &PartialD; Y = - f D 2 ( Ds 22 - N y s 23 ) = b 2 , 2
&PartialD; x &PartialD; Z = - f D 2 ( Ds 31 - N x s 33 ) = b 1,3 &PartialD; y &PartialD; Z = - f D 2 ( Ds 32 - N y s 33 ) = b 2,3
——成像模型关于GCPs及TPs的偏导数
系数bi,j指的是在矩阵BGCP和BTP对应的元素(BGCP、BTP定义见下文);
&PartialD; x &PartialD; &theta; = y &OverBar; p cos &theta; = s 1,9 &PartialD; y &PartialD; &theta; = 0 = s 2,9 - - - ( 15 )
——成像模型关于自检校未知数的偏导数
除了θ外,其它参数均为线性函数,系数si,j指的是在矩阵S中对应的元素(S定义见下文),参数Δxp,Δyp,ΔxC,ΔyC,c,k1,k2,p1,p2和sy的si,j值见下表:
Figure BDA0000154445880000109
<2>、系数阵构造
对于单轨道单条带区域影像,系数矩阵记为A,平差模型中的解向量x包含:
◆xEO:PPM参数(18×N_S)
◆xTP:TPs地面点坐标(3×N_TP)
◆xGCPP:GCP地面点坐标(3×N_GCP)
◆xSC:自检校参数(单镜头光学***:5+3×N_LINES;多镜头光学***:8×N_LINES)
未知数向量x及系数矩阵形式为:
x = x EO x TP x GCP x SC ; A = A GCP S GCP A TP B TP S TP C 0 C 1 C 2 F E S
观测方程组即为:
其中,矩阵Pxx是一类观测量的权重矩阵,其可进一步构造出所有观测量的整体权重对角矩阵,如下:
P ll = P GCP P TP P C 0 P C 1 P C 2 P EO P E P S
记n为观测量个数,u为未知数个数,则未知数向量x,系数阵A及权重矩阵Pll的维数分别为[u,1],[n,u]以及[n,n],法方程矩阵N为一个维数是[u×u]的对称系数阵;
<3>、线性***求解
至此,得线性***
Nx=z                  (17)
需求解未知数向量x,法方程矩阵N是一个块状矩阵,每一个单元块都有着不同的结构,如:对角阵、超对角阵、零矩阵以及稀疏阵,该矩阵较为复杂,它根据建模中的参数个数不同而有着不同的子矩阵数量,采用LU因式分解方法进行解算,即法方程矩阵被分解为:
N=LU                  (18)
其中,L为下三角阵(元素αij仅位于矩阵对角线及以下位置),U为上三角阵(元素βij仅位于矩阵对角线及以上位置),对于一个[n,n]大小的矩阵,(18)式变为:
n 11 n 12 &CenterDot; &CenterDot; &CenterDot; n 1 n n 21 n 22 &CenterDot; &CenterDot; &CenterDot; n 2 n &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; n n 1 n n 2 &CenterDot; &CenterDot; &CenterDot; n nn = &alpha; 11 0 &CenterDot; &CenterDot; &CenterDot; 0 &alpha; 21 &alpha; 22 &CenterDot; &CenterDot; &CenterDot; 0 &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &alpha; n 1 &alpha; n 2 &CenterDot; &CenterDot; &CenterDot; &alpha; nn &beta; 11 &beta; 12 &CenterDot; &CenterDot; &CenterDot; &beta; 1 n 0 &beta; 22 &CenterDot; &CenterDot; &CenterDot; &beta; 2 n &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; 0 0 &CenterDot; &CenterDot; &CenterDot; &beta; nn - - - ( 19 )
使用这样的分解式,线性集变换为:
Nx=(LU)x=L(Ux)=z    (20)
这样,可以分两步解算,首先,解算
Ly=z                  (21)
然后,解算
Ux=y                  (22)
之所以这样进行解算,是因为解算三角系数阵线性***是比较简单的,方程(21)可以进行向前替代解算,如下:
y 1 = z 1 &alpha; 11 - - - ( 23 )
y i = 1 &alpha; ii + [ z i - &Sigma; j = 1 i - 1 &alpha; ij y j , i = 2,3 , &CenterDot; &CenterDot; &CenterDot; n
同样,方程(28)可以进行向后替代解算,如下:
x n = y n &beta; nn - - - ( 24 ) .
x i = 1 &beta; ii + [ y i - &Sigma; j = i + 1 n &beta; ij x j , i = n - 1 , n - 2 , &CenterDot; &CenterDot; &CenterDot; 1
对所解得参数进行显著性检验,并进行数据分析,包括内部精度、RMS、相关系数和粗差探测计算如下:
<1>、显著性检验
下面给出基本的显著性检验方法,
当附加参数接近正交或正交时,可使用数理统计中的t分布,对所求得的参数,逐个进行显著性检验,t分布是按下式定义的变量:
t = &zeta; &eta;
其中ζ为标准正态分布N(0,1);η定义为
Figure BDA0000154445880000136
变量,υ是χ2的自由度,此时统计假设为:
Figure BDA0000154445880000137
Figure BDA0000154445880000138
为第i个附加参数的估值,取
&zeta; = a ^ i - E ( a ^ i ) &sigma; 0 q ii ~ N ( 0,1 )
&eta; = ( n - &mu; ) s 0 2 &sigma; 0 2 / ( n - u ) ~ &chi; 2 v , v = n - u
其中为单位权方差,
Figure BDA00001544458800001312
为由平差运算中得出的单位权中误差的平方,其期望值为
Figure BDA0000154445880000141
qii取自平差中未知数协因数阵Q的对角元素,由于假设
Figure BDA0000154445880000142
因此得出分布 t = &zeta; &eta; = a ^ i &sigma; 0 q ii = &sigma; 0 2 s 0 2 = a ^ i s 0 q ii , 当给定显著水平α,可由分布表查出临界值tα,若t<tα,则原假设成立,即该参数不显著,可在下次迭代平差中剔除,当附加参数间的相关较大时,一维的t-检验会导致错误的结论,由于t的相关往往仅出现在某一组的附加参数中,所以应该把这一组参数放在一起,使用多维检验的方法,此时,原假设为:
Figure BDA0000154445880000145
为在一起检验的k个附加参数,取 &xi; &prime; = [ a ^ - E ( a ^ ) ] T Q a ^ a ^ - 1 [ a ^ - E ( a ^ ) ] &sigma; 0 2 / k ~ &chi; 2 / v , v = k &eta; &prime; = ( n - u ) s 0 2 &sigma; 0 2 ( n - u ) ~ &chi; 2 / v , v = n - u , 因此得出F的统计量为:
Figure BDA0000154445880000147
根据两个自由度v(即k和(n-u))和假定显著水平α,可查F分布表,假如H0为真,则去除每个参数
Figure BDA0000154445880000148
当平差模型加入附加参数后,可能会由于附加参数之间或附加参数与其它未知参数之间的强相关而使得解的精度和可靠度恶化,为此应求出整个未知数的协因数阵,进而求出相关系数矩阵,进行逐一检查,以确定附加参数的可测定性检验、附加参数的外部可靠性检验及附加参数的相关性检验;
<2>、内部精度
Sigma naught(记
Figure BDA0000154445880000149
)为平差后单位权标准差,由***残差向量v计算而得,
&sigma; ^ 0 = v T pv n - u
定义Qxx
Qxx=N-1    (25)
则协方差矩阵Kxx计算公式为
K xx = &sigma; ^ 0 Q xx
Kxx为对称矩阵,对角元素为
&sigma; ^ kk = &sigma; ^ 0 q kk
其描述的是某一未知数xk的平均方差,在(i,j)位置上的非对角线元素描述的是未知数xi,xj之间的协方差,GCPs及TPs的平面及高程精度可分开按如下公式计算而得:
&sigma; x 2 = tr ( K xx x ) g &sigma; y 2 = tr ( K xx y ) g &sigma; z 2 = tr ( K xx z ) g
其中
σx,σy和σz:X、Y及Z方向上的标准差;
Figure BDA0000154445880000156
GCPs(或TPs)的X,Y,Z坐标的Kxx
g:点数;
<3>、RMS计算
已知地面点坐标的TPs可以作为检查点使用(CPs),方法是将地面坐标平差所得估值
Figure BDA0000154445880000157
与相应的“正确”坐标值(准真值)[Xcorr,Ycorr,Zcorr]相比较,RMSE按下式计算
RMSE x 2 = &Sigma; i = 1 N _ CP ( X ^ - X corr ) 2 N _ CP RMSE y 2 = &Sigma; i = 1 N _ CP ( Y ^ - Y corr ) 2 N _ CP RMSE z 2 = &Sigma; i = 1 N _ CP ( Z ^ - Z corr ) 2 N _ CP
其中N_CP为CPs点个数;
<4>、相关系数计算
需要对附加参数的确定性进行研究,以避免过度参数化问题的产生,基于此,需计算未知数之间的相关系数ρij
&rho; ij = q ij q ii q jj
其中,qij为N-1中所定义矩阵(i,j)位置上的元素,ρij的绝对值越接近1即说明xi,xj之间的存在强相关性,如果存在强相关参数,则应使得这两者之一以伪观测量的形式作为常值出现;
<5>、粗差探测
由于公式(1)中的“真值残差”不可得,所以使用Grün方法进行摄影测量区域平差像点坐标观测量粗差探测,该方法基于Baarda数据挖掘技术,即对于每个观测量i,值wi按以下公式计算:
w i = - v i &sigma; ^ v i
其中
Figure BDA0000154445880000162
Figure BDA0000154445880000163
为矩阵Qvv中的第i个元素,Qvv定义如下:
Q vv = P ll - 1 - A ( A T P ll A ) - 1 A T
如果零假设
H 0 i v : E ( v i ) = 0
为真,则wi服从τ分布,如果***残差足够大,则τ分布可用更易于处理的学生分布进行替代,粗差探测一般在平差结束后进行,人工剔除被检测出的观测量粗差后方可进行新一轮的平差处理。
本发明可以完成相机的在轨几何检校,有望对影响传感器定位精度的不可用、不可知或已变化的相关标称参数等进行推求,近而改正影响相机定位能力的几何***误差参数,从而挖掘和提升传感器的量测性能。本发明有望增强国内多镜头航天线阵相机***在轨检校这一技术短板,为提升国产航天线阵传感器的对地定位能力提供技术支撑。
附图说明
图1是实施本发明的多镜头航天三线阵相机立体影像获取示意图,1,2,3分别指代前视、下视和后视相机。
图2是实施本发明的相机j与下视相机的几何位置相对关系图。
图3是实施本发明的分段多项式轨道模型示例图,s1,si,sn分别代表第1段,第i段及第n段卫星平台飞行轨道。
图4是关于GCP和TP的系数矩阵A结构,S为对应轨道分段,O为对应观测方程。
图5是系数矩阵AGCP和ATP中的单元块,大小为2行×18列。
图6是关于GCP和TP的系数矩阵B结构。
图7是系数矩阵BGCP和BTP中的单元块,大小为2行×3列。
图8是单镜头、多镜头推扫式传感器光学***系数矩阵SGCP结构(line代表CCD线阵,lens代表镜头)。
图9是SGCP和STP中的矩阵块,大小分别为2行×3列,2行×6列。
图10是矩阵C0结构。
图11是矩阵C0中的单元块。
图12是矩阵C1结构。
图13是矩阵C1中的单元块。
图14是矩阵C2中的单元块。
图15是矩阵C2中的单元块。
图16是最终的系数阵(右)及未知数向量(左)。
具体实施方式
下面以利用多镜头航天三线阵相机***获取的同一条带上的三个影像带为例,如图1所示,说明本发明提出的多镜头航天线阵相机***在轨几何自检校方法的实施方式。
图1为多镜头航天三线阵相机***在飞行中获取下视、前视和后视影像的成像过程。一种多镜头航天线阵相机***在轨几何自检校方法的具体实施步骤如下:
(1)根据传感器飞行轨迹和影像条带宽度,在所获取的影像对应的地面试验区域内布设高精度控制点。
控制点的布设可以在卫星过顶之前,根据卫星的星下点轨迹,在地面试验区内进行高精度标志控制点布设;也可以在卫星过顶之后,根据影像特征判点后进行相应的控制测量(全野外量测)。前者的精度高,但成本也高。后者的精度略低,但布控灵活,外业成本低。无论是何种布控形式,为了检校参数的全面性和可靠性,最终的控制点对应像点量测精度应保证在1个像元左右——标志控制点的像点量测精度应在一个像元内,全野外量测的控制点对应像点量测精度不应超过2个像元。另外,控制点的密度和数量也应予以考虑。
(2)对获取检校区的前、下、后视相机三个影像条带进行控制点同名像点量测。
关于像点坐标的量测,可进行像点的高精度影像匹配量测,也可以人工立体观测。目前可用的商用全数字摄影测量像点量测工具很多,比较通用的有Eardas、PCI、ENVI及Ermaper等。
(3)获取传感器的轨道数据,包括GPS位置观测量及星敏的姿态观测量。
试验区域对应的传感器轨道GPS位置及星敏的姿态观测数据是连同对应区域的影像数据由传感器拥有者一同提供的。
(4)根据多镜头航天线阵传感器成像模型和外部定向参数模型,建立观测方程组。
(5)利用最小二乘法进行光束法自检校平差数据处理,求解几何***误差检校参数估值——首先进行线性化、系数阵构造、法化后线性***求解。
(6)对所解得参数进行显著性检验,并进行数据分析,包括内部精度、RMS、相关系数和粗差探测计算。
由于(4)~(6)三个步骤为本方法的核心步骤,下面就其所涉及的理论和数据处理方法详述如下:
第(4)步中涉及三个数学模型,并以所建立的三个模型为基础,建立观测方程组如下:
①多镜头航天线阵传感器成像模型
对于光学***由多个镜头组成的传感器,需要考虑描述其它镜头相对于下视镜头(平台)的位置与姿态关系。例如,对于镜头j,用dxj,dyj,dzj表示相对于下视镜头的相对位置,用αj,βj,γj表示相对于下视镜头的相对姿态(图2)。fj,xpj,ypj为镜头j的焦距和主点坐标。成像模型为:
x - x pj = - f j &CenterDot; s 11 ( X - X C ) + s 21 ( Y - Y C ) + s 31 ( Z - Z C ) - ( m 11 d xj + m 21 d yj + m 31 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj )
y - y pj = - f j &CenterDot; s 12 ( X - X C ) + s 22 ( Y - Y C ) + s 32 ( Z - Z C ) - ( m 12 d xj + m 22 d yj + m 32 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj ) - - - ( 2 )
其中,
[X,Y,Z]:物点在地面坐标***中的坐标;
[XC,YC,ZC]:下视相机投影中心在地面坐标***中的坐标;
[x,y]:像点在相机坐标***中的坐标,理论上x为0或某一常数;
[xp,yp]:像主点在相机坐标***中的坐标;
f:焦距;
从下视相机坐标***到地面坐标***的旋转矩阵,旋角分别为
Figure BDA0000154445880000195
M(αj,βj,γj)为从下视镜头相机坐标***到镜头j相机坐标***的旋转矩阵;
Figure BDA0000154445880000196
为从镜头j相机坐标***到地面坐标***的旋转矩阵。
可以看出,该模型对单镜头与多镜头都是适用的(j=1,2,3,…,n)。
考虑到相机自检校,上述模型可进一步扩展为:
x - x pj = - f j &CenterDot; s 11 ( X - X C ) + s 21 ( Y - Y C ) + s 31 ( Z - Z C ) - ( m 11 d xj + m 21 d yj + m 31 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj ) + &Delta; xj
y - y pj = - f j &CenterDot; s 12 ( X - X C ) + s 22 ( Y - Y C ) + s 32 ( Z - Z C ) - ( m 12 d xj + m 22 d yj + m 32 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj ) + &Delta; yj
(3)
其中,Δxj,Δyj是关于镜头畸变、CCD线阵畸变以及内定向参数(改正量)等成像***的***误差模型。
②外部定向参数模型及轨道观测量的引入
外部定向参数模型——轨道模型(trajectory model)是用来描述传感器外部定向(EO,external orientation)状态的,其是用关于时间的分段多项式(PPM,PiecewisePolynomial Function)来建模的。如果多项式分段数为1,则模型即为一个抛物线函数,图3。
可以根据控制点(GCPs)和连接点(TPs)的分布来确定轨道模型的分段数。对于被分割轨道的第i段(该段的首末端点的线阵影像获取时刻记为
Figure BDA0000154445880000201
),则定义时间变量
Figure BDA0000154445880000202
为:
t &OverBar; = t - t ini i t fin i - t ini i - - - ( 4 )
其中t为该段轨道上某线阵影像的获取时刻。以2阶多项式为例(时间
Figure BDA0000154445880000204
为因变量),则在分段多项式的每一段,传感器外部定向参数
Figure BDA0000154445880000205
模型为:
X C ( t &OverBar; ) = X 0 i + X 1 i t &OverBar; + X 2 i t &OverBar; 2
Y C ( t &OverBar; ) = Y 0 i + Y 1 i t &OverBar; + Y 2 i t &OverBar; 2 - - - ( 5 )
Z C ( t &OverBar; ) = Z 0 i + Z 1 i t &OverBar; + Z 2 i t &OverBar; 2
&omega; C ( t &OverBar; ) = &omega; 0 i + &omega; 1 i t &OverBar; + &omega; 2 i t &OverBar; 2
Figure BDA00001544458800002010
&kappa; C ( t &OverBar; ) = &kappa; 0 i + &kappa; 1 i t &OverBar; + &kappa; 2 i t &OverBar; 2
[X0,X1,X2,…,κ0,κ1,κ2]i为第i段外部定向模型系数。在多项式的相邻两段连接点处,需要加入多项式的连续性及光滑性条件,即多项式的0阶、1阶和2阶导数在各相邻分段连接点处相等。
得到每个线阵成像时刻的位置和姿态观测量后,可以用线阵扫描行数代替时间,公式(4)可以被改写为:
t &OverBar; = l - l ini i l fin i - l ini i - - - ( 6 )
其中,
Figure BDA0000154445880000211
为第i段轨道的首末端在对应影像中的线阵行数(首末线阵影像行)。这样,对于第i段轨道所对应影像中的第1行线阵影像,其位置与姿态
Figure BDA0000154445880000212
的数学模型可以表示为该行线阵影像的外部定向参数观测量与关于时间的2阶多项式的和,如下式:
X C ( t &OverBar; ) = X instr + X 0 i + X 1 i t &OverBar; + X 2 i t &OverBar; 2
Y C ( t &OverBar; ) = Y instr + Y 0 i + Y 1 i t &OverBar; + Y 2 i t &OverBar; 2 - - - ( 7 )
Z C ( t &OverBar; ) = Z instr + Z 0 i + Z 1 i t &OverBar; + Z 2 i t &OverBar; 2
&omega; C ( t &OverBar; ) = &omega; instr + &omega; 0 i + &omega; 1 i t &OverBar; + &omega; 2 i t &OverBar; 2
Figure BDA0000154445880000219
&kappa; C ( t &OverBar; ) = &kappa; instr + &kappa; 0 i + &kappa; 1 i t &OverBar; + &kappa; 2 i t &OverBar; 2
其中,[X0,X1,X2,…,κ0,κ1,κ2]i是第i段轨道关于GPS偏移量及INS(星敏)轴线偏差、GPS及INS(星敏)***漂移等***误差的18个模型参数。具体讲,
Figure BDA00001544458800002111
常数项用来补偿成像***与GPS之间的偏移量(shifts)及与INS(星敏)的轴线偏差(drifts),一次项和二次项
Figure BDA00001544458800002113
是除此之外的GPS及INS(星敏)观测量中的***漂移等误差项参数。还有,为了保证满足多项式的连续性及光滑性约束条件,在轨道各段之间的连接点处,要用模型的0阶、1阶及2阶连续性及光滑性条件进行约束,即连接点处,相邻轨道段的多项式函数值、1阶及2阶导数值要求相等。
③相机成像***的像差模型
CCD及物镜畸变主要包含在公式(3)中的Δx,Δy像差模型中,对于多镜头线阵CCD传感器,现假设传感器仅有N_LENS个镜头,每个镜头的焦面上安置一组CCD线阵,总共N_LINES个CCD线阵。需要说明的是,在这种多镜头传感器结构中,考虑到数据量和传输率的问题,往往每个镜头的焦面上只安置一个CCD线阵(即N_LINES=N_LENS),这样,则需要对N_LENS个镜头畸变及对应的N_LINES(或N_LENS)个CCD线阵焦面上的位移进行建模。附加参数说明如下:
·每个镜头的主点改正量(位移量)Δxp,Δyp
·k1,k2,p1,p2:每个镜头的对称性径向和切向畸变;
·sy:每个线阵CCD单元y方向上的比例因子;
·每个CCD线阵在焦面上的旋角θ(其产生的畸变主要在x方向,在y方向上的畸变影响可以忽略不计)。
这样,总共的附加参数个数为N_LENS×8,对各项误差进行累加即为单个镜头j的内定向参数误差,公式为:
&Delta; x j = &Delta; x pj - &Delta;f f x &OverBar; pj + ( k 1 j r j 2 + k 2 j r j 4 ) x &OverBar; pj + p 1 j ( r j 2 + 2 x &OverBar; pj 2 ) + 2 p 2 j x &OverBar; pj y &OverBar; pj + x &OverBar; pj sin &theta;
&Delta; y j = &Delta; y pj - &Delta;f f y &OverBar; pj + ( k 1 j r j 2 + k 2 j r j 4 ) y &OverBar; pj + + 2 p 1 j x &OverBar; pj y &OverBar; pj + p 2 j ( r j 2 + 2 y &OverBar; pj 2 ) + s yj y &OverBar; pj - - - ( 8 )
其中, x &OverBar; pj = x - x pj , y &OverBar; pj = y - y pj , r 2 = x &OverBar; pj 2 + y &OverBar; pj 2 .
①②③中的数学模型以共线方程为基础,对其进行扩展以适应于传感器的内部定向建模及外部定向建模。内部定向模型全面考虑了镜头畸变、CCD线阵在焦面内的移位与旋转。外部定向模型为关于时间或影像扫描行数的2阶分段多项式,据此,可对轨道进行适当的分段,每一段轨道的轨迹为一抛物线函数,而且,在航空飞机或直升机上的推扫式传感器必须使用GPS、INS观测设备,这些GPS、INS观测量可被整合到平差过程中去,并用***性误差加以改正。轨道分段数要视轨道长度以及光滑性要求等具体情况而定。通常,卫星轨道分为2段即可;但对于航空飞机或直升机平台而言,由于轨道平滑性差,所以轨道分段数要多得多。在提供适量的地面控制点前提下,采用最小二乘技术进行数据平差处理以得到***可能的最佳估值解。观测量(诸如像坐标、地面控制点,未知数的伪观测量)的权重根据相应观测量的精度而定。
④建立观测方程组
下面将详细列出关于上述线阵传感器模型的观测方程组。此处,方程组中的所有观测量都带有上划线记号,待平差的定向参数无上划线标记。
·像坐标观测量
此处的像坐标观测量是指控制点及连接点的像坐标观测值。像坐标观测方程用共线方程(3)来描述像坐标观测量
Figure BDA0000154445880000226
及其对应的改正量vx,vy,传感器内部、外部定向参数以及对应物点地面坐标之间的关系:
x &OverBar; + v x = x = x pj - f j &CenterDot; N x D + &Delta; xj - - - ( 9 )
y &OverBar; + v y = y = y pj - f j &CenterDot; N y D + &Delta; yj
由于该观测方程组是非线性方程组,所以在平差处理中需要进行线性化,观测量的权重与其观测精度相关,后面详述。
·外部定向参数观测量
外部定向参数观测量分为两类,第一类是关于PPM函数的连续性约束条件的观测量。第二类观测量是带一定权重的PPM参数初始值观测量。
★连续性约束条件观测量
轨道模型的多项式0阶、1阶及2阶导数值连续性约束条件观测量为权值很高的伪观测值。连续性约束条件观测方程组为:
( X 0 i &OverBar; + X 1 i &OverBar; + X 2 i &OverBar; - X 0 i + 1 &OverBar; ) + v X C 0 = X 0 i + X 1 i + X 2 i - X 0 i + 1
( Y 0 i &OverBar; + Y 1 i &OverBar; + Y 2 i &OverBar; - Y 0 i + 1 &OverBar; ) + v Y C 0 = Y 0 i + Y 1 i + Y 2 i - Y 0 i + 1
( Z 0 i &OverBar; + Z 1 i &OverBar; + Z 2 i &OverBar; - Z 0 i + 1 &OverBar; ) + v Z C 0 = Z 0 i + Z 1 i + Z 2 i - Z 0 i + 1
( &omega; 0 i &OverBar; + &omega; 1 i &OverBar; + &omega; 2 i &OverBar; - &omega; 0 i + 1 &OverBar; ) + v &omega; C 0 = &omega; 0 i + &omega; 1 i + &omega; 2 i - &omega; 0 i + 1
Figure BDA0000154445880000237
( &kappa; 0 i &OverBar; + &kappa; 1 i &OverBar; + &kappa; 2 i &OverBar; - &kappa; 0 i + 1 &OverBar; ) + v &kappa; C 0 = &kappa; 0 i + &kappa; 1 i + &kappa; 2 i - &kappa; 0 i + 1 - - - ( 10 )
1阶连续性约束条件观测方程组
( X 1 i &OverBar; + 2 X 2 i &OverBar; - X 1 i + 1 &OverBar; ) + v X C 1 = X 1 i + 2 X 2 i - X 1 i + 1
( Y 1 i &OverBar; + 2 Y 2 i &OverBar; - Y 1 i + 1 &OverBar; ) + v Y C 1 = Y 1 i + 2 Y 2 i - Y 1 i + 1
( Z 1 i &OverBar; + 2 Z 2 i &OverBar; - Z 1 i + 1 &OverBar; ) + v Z C 1 = Z 1 i + 2 Z 2 i - Z 1 i + 1
( &omega; 1 i &OverBar; + 2 &omega; 2 i &OverBar; - &omega; 1 i + 1 &OverBar; ) + v &omega; C 1 = &omega; 1 i + 2 &omega; 2 i - &omega; 1 i + 1
Figure BDA00001544458800002313
( &kappa; 1 i &OverBar; + 2 &kappa; 2 i &OverBar; - &kappa; 1 i + 1 &OverBar; ) + v &kappa; C 1 = &kappa; 1 i + 2 &kappa; 2 i - &kappa; 1 i + 1 - - - ( 11 )
2阶连续性约束条件观测方程组
( X 2 i &OverBar; - X 2 i + 1 &OverBar; ) + v X C 2 = X 2 i - X 2 i + 1
( Y 2 i &OverBar; - Y 2 i + 1 &OverBar; ) + v Y C 2 = Y 2 i - Y 2 i + 1
( Z 2 i &OverBar; - Z 2 i + 1 &OverBar; ) + v Z C 2 = Z 2 i - Z 2 i + 1
( &omega; 2 i &OverBar; - &omega; 2 i + 1 &OverBar; ) + v &omega; C 2 = &omega; 2 i - &omega; 2 i + 1
Figure BDA0000154445880000245
( &kappa; 2 i &OverBar; - &kappa; 2 i + 1 &OverBar; ) + v &kappa; C 2 = &kappa; 2 i - &kappa; 2 i + 1 - - - ( 12 )
相应的伪观测量的权值应设置为很高(~1020)。
★外部定向参数观测量
PPM参数的初始值需看作伪观测值来处理。对于每一段多项式,可以写出如下外部定向参数方程式:
X 0 &OverBar; + v X 0 = X 0 Y 0 &OverBar; + v Y 0 = Y 0 Z 0 &OverBar; + v Z 0 = Z 0
X 1 &OverBar; + v X 1 = X 1 Y 1 &OverBar; + v Y 1 = Y 1 Z 1 &OverBar; + v Z 1 = Z 1 - - - ( 13 )
X 2 &OverBar; + v X 2 = X 2 Y 2 &OverBar; + v Y 2 = Y 2 Z 2 &OverBar; + v Z 2 = Z 2
&omega; 0 &OverBar; + v &omega; 0 = &omega; 0
Figure BDA00001544458800002417
&kappa; 0 &OverBar; + v &kappa; 0 = &kappa; 0
&omega; 1 &OverBar; + v &omega; 1 = &omega; 1
Figure BDA00001544458800002420
&kappa; 1 &OverBar; + v &kappa; 1 = &kappa; 1
&omega; 2 &OverBar; + v &omega; 2 = &omega; 2
Figure BDA00001544458800002423
&kappa; 2 &OverBar; + v &kappa; 2 = &kappa; 2
如果外部定向参数函数模型的最高阶为1阶或0阶,那么此时2阶或1阶参数初始值设置为0,对应的伪观测量的权重值可设置为一个非常大的数(~1020)。
·自检校参数
如果实验室已经给出内部定向参数的检测值,那么我们可以把它们当作带权观测量来处理。例如,对于镜头j的内部定向参数观测方程为:
&Delta; x P &OverBar; + v &Delta; x P = &Delta; x P
&Delta; y P &OverBar; + v &Delta; y P = &Delta; y P
&Delta;f &OverBar; + v &Delta;f = &Delta;f
k 1 &OverBar; + v k 1 = k 1
k 2 &OverBar; + v k 2 = k 2
p 1 &OverBar; + v p 1 = p 1
p 2 &OverBar; + v p 2 = p 2
s y &OverBar; + v s y = s y - - - ( 14 )
&theta; &OverBar; + v &theta; = &theta;
权重可设置为实验室检测结果的协方差的倒数。
·地面点
在更高级或要求更高的平差处理过程中,地面点坐标并不作为固定值(真值,无误差值)来处理。实际上,其也是经过测量而得,所以也应当作为观测量来对待(精度或权重由其标准差而定)。所以地面点坐标观测量观测方程可写为:
X &OverBar; + v X = X
(15)
Y &OverBar; + v Y = Y
Z &OverBar; + v Z = Z
第(5)步中包括线性化、给定未知数初值、系数阵构造、线性***求解等环节,如下:
①线性化(对第4步中的三个模型——成像模型和轨道模型及像差模型求偏导)
参考外部定向模型(5)和(7),成像模型(3)中x关于外部定向模型参数的偏导数如下:
&delta;x &delta; X 0 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 0 &delta;x &delta; X 1 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 1 &delta;x &delta; X 2 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 2
&delta;x &delta; Y 0 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 0 &delta;x &delta; Y 1 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 1 &delta;x &delta; Y 2 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 2
&delta;x &delta; Z 0 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Z 0 &delta;x &delta; Z 1 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Y 1 &delta;x &delta; Z 2 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Y 2
&delta;x &delta; &omega; 0 = &delta;x &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 0 &delta;x &delta; &omega; 1 = &delta;x &delta;&omega; C &CenterDot; &delta; &omega; C &delta; &omega; 1 &delta;x &delta; &omega; 2 = &delta;x &delta;&omega; C &CenterDot; &delta;&omega; C &delta;&omega; 2
Figure BDA00001544458800002613
Figure BDA00001544458800002614
Figure BDA00001544458800002615
&delta;x &delta; &kappa; 0 = &delta;x &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 0 &delta;x &delta; &kappa; 1 = &delta;x &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 1 &delta;x &delta; &kappa; 2 = &delta;x &delta; &kappa; C &CenterDot; &delta;&kappa; C &delta; &kappa; 2
y关于外部定向模型参数的偏导数同理可得:
&delta;y &delta; X 0 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 0 &delta;y &delta; X 1 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 1 &delta;y &delta; X 2 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 2
&delta;y &delta; Y 0 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 0 &delta;y &delta; Y 1 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 1 &delta;y &delta; Y 2 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 2
&delta;y &delta; Z 0 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Z 0 &delta;y &delta; Z 1 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Y 1 &delta;y &delta; Z 2 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Y 2
&delta;y &delta; &omega; 0 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 0 &delta;y &delta; &omega; 1 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 1 &delta;y &delta; &omega; 2 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 2
Figure BDA00001544458800002631
Figure BDA00001544458800002632
Figure BDA00001544458800002633
&delta;y &delta; &kappa; 0 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 0 &delta;y &delta; &kappa; 1 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 1 &delta;y &delta; &kappa; 2 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 2
其中,
Figure BDA00001544458800002637
Figure BDA00001544458800002638
关于传感器位置和姿态的偏导数为:
&PartialD; x &PartialD; X C = - f D 2 ( s 13 N x - s 11 D ) = a 1,1
&PartialD; x &PartialD; Y C = - f D 2 ( s 23 N x - s 21 D ) = a 1,2
&PartialD; x &PartialD; Z C = - f D 2 ( s 33 N x - s 31 D ) = a 1,3
&PartialD; x &PartialD; &omega; C = f D { [ ( X - X C ) &delta; s 13 &PartialD; &omega; C + ( Y - Y C ) &delta; s 23 &PartialD; &omega; C - ( Z - Z C ) &delta; s 33 &PartialD; &omega; C ] N x D - ( X - X 0 ) &delta; s 11 &PartialD; &omega; C - ( Y - Y 0 ) &delta; s 21 &PartialD; &omega; C - ( Z - Z 0 ) &delta; s 31 &PartialD; &omega; C } = a 1,4
Figure BDA0000154445880000275
&PartialD; x &PartialD; &kappa; C = f D { [ ( X - X C ) &delta;s 13 &delta;&kappa; C + ( Y - Y C ) &delta;s 33 &PartialD; &kappa; C ] N x D - ( X - X 0 ) &delta;s 11 &delta; &kappa; C - ( Y - Y 0 ) &delta;s 21 &PartialD; &kappa; C - ( Z - Z 0 ) &delta;s 31 &PartialD; &kappa; C } = a 1,6
&PartialD; y &PartialD; X C = - f N 2 ( r 13 Z y - r 13 N ) = a 2,1
&PartialD; y &PartialD; Y C = - f N 2 ( r 23 Z y - r 22 N ) = a 2 , 2
&PartialD; y &PartialD; Z C = - f N 2 ( r 33 Z y - r 32 N ) = a 2,3
&PartialD; y &PartialD; &omega; C = f N { [ ( Y - Y 0 ) r 33 - ( Z - Z 0 ) r 23 ] Z y N - ( Y - Y 0 ) r 32 - ( Z - Z 0 ) r 22 } = a 2,7
Figure BDA00001544458800002711
&PartialD; y &PartialD; &kappa; C = - f N Z x = a 2,9 - - - ( 16 )
系数ai,j指的是在矩阵AGCP和ATP对应的元素(AGCP、ATP定义见下文)。在公式(3)中对控制点(GCPs)及连接点(TPs)求偏导,有,
&PartialD; x &PartialD; X = - f D 2 ( Ds 11 - N x s 13 ) = b 1,1 &PartialD; y &PartialD; X = - f D 2 ( Ds 12 - N y s 13 ) = b 2,1
&PartialD; x &PartialD; Y = - f D 2 ( Ds 31 - N x s 23 ) = b 1,2 &PartialD; y &PartialD; Y = - f D 2 ( Ds 22 - N y s 23 ) = b 2 , 2 - - - ( 17 )
&PartialD; x &PartialD; Z = - f D 2 ( Ds 31 - N x s 33 ) = b 1,3 &PartialD; y &PartialD; Z = - f D 2 ( Ds 32 - N y s 33 ) = b 2,3
系数bi,j指的是在矩阵BGCP和BTP对应的元素(BGCP、BTP定义见下文)。在自检校未知数中,除了θ外,其它参数均为线性函数,根据公式(8)可求得关于参数θ的偏导数为:
&PartialD; x &PartialD; &theta; = y &OverBar; p cos &theta; = s 1,9 &PartialD; y &PartialD; &theta; = 0 = s 2,9 - - - ( 18 )
系数si,j指的是在矩阵S中对应的元素(S定义见下文)。参数Δxp,Δyp,ΔxC,ΔyC,c,k1,k2,p1,p2和sy的si,j值在表1中列出。
表1x,y关于Δxp,Δyp,ΔxC,ΔyC,c,k1,k2,p1,p2和sy参数在矩阵S中相应位置的系数
Figure BDA0000154445880000281
②给定未知数初值
为了计算偏导数系数以进行平差迭代,必须给定未知数初值。
传感器外都定向模型参数初值由平台轨道观测量而得。
可以采用离散数学的方法计算每段轨道的PPM参数[X0,X1,X2,…,κ0,κ1,κ2]初值。以第i段轨道的(该段的首末端点的线阵影像获取时刻记为)函数XC(t)为例,根据起点时刻
Figure BDA0000154445880000283
中点时刻
Figure BDA0000154445880000284
和末点时刻
Figure BDA0000154445880000285
函数值XC可由三次样条内插而得,即一阶导数及二阶导数计算公式如下:
d 1 = X ( t fin i ) - X ( t ini i ) t fin i - t ini i
d 2 = 2 &CenterDot; X ( t mid i ) - X ( t ini i ) - X ( t fin i ) ( t fin i - t ini i ) 2
由于
X C ( t &OverBar; ) = X 0 + X 1 t &OverBar; + X 2 t &OverBar; 2
&delta; X C &delta; t &OverBar; = X 1 + 2 X 2 t &OverBar;
&delta; 2 X C &delta; t &OverBar; 2 = 2 X 2
Figure BDA0000154445880000294
X0=Xini,则有
Xfin=Xini+X1+X2
&delta; X C &delta; t &OverBar; = X 1 + 2 X 2
&delta; 2 X C &delta; t &OverBar; 2 = 2 X 2
可设
d 1 = &delta; X C &delta; t &OverBar;
d 2 = &delta; 2 X C &delta; t &OverBar; 2
对于每段轨道,可计算X2
X 2 = d 2 X &prime; &prime; ,
X1=d1-2X2
同样的方法可以用来估算分段多项式模型
Figure BDA00001544458800002910
的系数初值。
地面控制点、连接点的地面坐标的近似值可以由外方位元素的初值进行前方交会计算而得。
还有,自检校参数初值可以设置为实验室检校值,若未进行实验室检校,自检校参数初值可设置为0。
③系数阵构造
这里主要说明系数阵的分块构造,并进行简单的图示。矩阵分块所涉及的维数用以下符号记:N_GCP表示GCPs个数,N_TP表示TPs个数,N_S表示参与平差的分段多项式的分段数,N_LINES表示CCD线阵数,N_LENS表示镜头数。
·关于GCPs、TPs(像点)的共线方程
线性化后,关于GCPs及TPs的系数阵是由共线方程对未知数求偏导后,进而代入未知数初值计算而得。
对于控制点而言,所生成的矩阵有:
★AGCP,包含由公式(16)求得的关于外方位未知数的偏导数;
★BGCP,包含由公式(17)求得的关于地面点坐标的偏导数;
★SGCP,包含由公式(18)求得的自检校参数未知数的偏导数。
★对于连接点而言,所生成的矩阵有:
★ATP,包含由公式(16)求得的关于外方位未知数的偏导数;
★BTP,包含由公式(17)求得的关于地面点坐标的偏导数;
★STP,包含由公式(18)求得的自检校参数未知数的偏导数。
AGCP,ATP为块状矩阵(图4),每块的维数为(2,18)(图5)。AGCP的维数为[2×N_LINES×N_GCP,18×N_S],ATP的维数为[2×N_LINES×N_TP,18×N_S]。
BGCP,BTP为块状对角矩阵(图6、图7),每块的维数为(2,3)。BGCP的维数为[3×N_GCP,3×N_GCP],BTP的维数为[3×N_TP,3×N_TP]。
矩阵SGCP和STP有着相同的结构,只是针对单镜头和多镜头由于参数个数不同,矩阵结构存在细微的差别。对于单镜头光学***SGCP和STP的维数为[2×N_LINES×N_GCP,5+3×N_LINES]和[2×N_LINES×N_TP,5+3×N_LINES];对于多镜头光学***SGCP和STP的维数为[2×N_LINES×N_GCP,8×N_LINES]和[2×N_LINES×N_TP,8×N_LINES]。图8为单镜头、多镜头推扫式传感器光学***系数矩阵SGCP结构,STP结构与之类似。图9为相应的矩阵块。
·轨道模型约束条件
描述多项式连续性约束条件的系数矩阵记为Cn(n=0,1,2,根据偏导数阶数而定)。可根据公式(10)、(11)和(12)进行构造,矩阵维数为[6×(N_S-1),18×N_S]。
矩阵C0所描述的函数连续性约束条件见图10、图11,一阶约束条件系数阵记为C1(图12、图13),二阶约束条件系数阵记为C2(图14、图15)。
·PPM参数初值
可以伪观测量的形式对PPM参数未知数初值进行赋值,此处,用以描述PPM参数伪观测量的系数矩阵记为F。矩阵F为单位对角阵,维数与PPM参数个数相等,即[18×N_S,18×N_S]。
·APs初值
描述自检校参数(Δxp,Δyp,k1,k2,p1,p2,sy和θ)伪观测量矩阵为一单位对角阵,记S,维数[dim_S,dim_S]。dim_S等于N_LINES×8(单镜头);dim_S等于5+3×N_LINES(多镜头)。
·GCP坐标初值
如同上述情况一样,控制点地面坐标观测方程组的系数阵也是以伪观测量的形式被引入平差模型中来的,记E,为一[3×N_GCP,3×N_GCP]对角矩阵。
·最终的系数阵
对于单轨道单条带区域影像,最终系数矩阵记为A,平差模型中的解向量x包含:
★xEO:PPM参数(18×N_S)
★xTP:TPs地面点坐标(3×N_TP)
★xGCP:GCP地面点坐标(3×N_GCP)
★xSC:自检校参数(单镜头光学***:5+3×N_LINES;多镜头光学***:8×N_LINES)
公式(19)、(20)及图16示意了未知数向量x及系数矩阵A。
x = x EO x TP x GCP x SC - - - ( 19 )
A = A GCP S GCP A TP B TP S TP C 0 C 1 C 2 F E S - - - ( 20 )
最终的观测方程组为
Figure BDA0000154445880000321
矩阵Pxx是一类观测量的权重矩阵,其可进一步构造出所有观测量的整体权重对角矩阵。
P ll = P GCP P TP P C 0 P C 1 P C 2 P EO P E P S - - - ( 22 )
记n为观测量个数,u为未知数个数,则未知数向量x,系数阵A及权重矩阵Pll的维数分别为[u,1],[n,u]以及[n,n]。
最终的法方程矩阵N为一个维数是[u×u]的对称系数阵。
④线性***求解
至此,得线性***
Nx=z    (23)
需解得未知数向量x。
法方程矩阵N是一个块状矩阵,每一个单元块都有着不同的结构,如:对角阵、超对角阵、零矩阵以及稀疏阵。该矩阵较为复杂,它根据建模中的参数个数不同而有着不同的子矩阵数量。为了提高存储及计算效率,需要针对这样的矩阵找到合适的算法。问题的关键是如何最大程度的限制非零元素的规模。
最小二乘解算一般都是采用LU因式分解,Cholesky因式分解以及QR因式分解等方法。在这三种方法中,矩阵会被分解为两个具有使得解算简化的一定特性的两个矩阵。
此处,采用LU因式分解方法进行解算。即法方程矩阵被分解为:
N=LU                  (24)
其中,L为下三角阵(元素αij仅位于矩阵对角线及以下位置),U为上三角阵(元素βij仅位于矩阵对角线及以上位置)。对于一个[n,n]大小的矩阵,(24)式变为:
n 11 n 12 &CenterDot; &CenterDot; &CenterDot; n 1 n n 21 n 22 &CenterDot; &CenterDot; &CenterDot; n 2 n &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; n n 1 n n 2 &CenterDot; &CenterDot; &CenterDot; n nn = &alpha; 11 0 &CenterDot; &CenterDot; &CenterDot; 0 &alpha; 21 &alpha; 22 &CenterDot; &CenterDot; &CenterDot; 0 &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &alpha; n 1 &alpha; n 2 &CenterDot; &CenterDot; &CenterDot; &alpha; nn &beta; 11 &beta; 12 &CenterDot; &CenterDot; &CenterDot; &beta; 1 n 0 &beta; 22 &CenterDot; &CenterDot; &CenterDot; &beta; 2 n &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; 0 0 &CenterDot; &CenterDot; &CenterDot; &beta; nn - - - ( 25 )
使用这样的分解式,线性集变换为:
Nx=(LU)x=L(Ux)=z    (26)
这样,可以分两步解算,首先,解算
Ly=z                  (27)
然后,解算
Ux=y                  (28)
之所以这样进行解算,是因为解算三角系数阵线性***是比较简单的。方程(27)可以进行向前替代解算,如下:
y 1 = z 1 &alpha; 11 - - - ( 29 )
y i = 1 &alpha; ii + [ z i - &Sigma; j = 1 i - 1 &alpha; ij y j , i = 2,3 , &CenterDot; &CenterDot; &CenterDot; n
同样,方程(28)可以进行向后替代解算,如下:
x n = y n &beta; nn - - - ( 30 )
x i = 1 &beta; ii + [ y i - &Sigma; j = i + 1 n &beta; ij x j , i = n - 1 , n - 2 , &CenterDot; &CenterDot; &CenterDot; 1
第(6)步为整个方法的最后一步,即对所解得参数进行显著性检验,进行数据分析,包括内部精度、RMS、相关系数、粗差探测等计算,最后给出标定参数最终结果。
当平差完成后,地面点坐标的内部精度、外部精度应予以分析。同时也应该研究未知数之间的相关性,以便移除过度参数化的参数。此外,还应执行粗差探测,等工作,详述如下:
·附加参数的显著性检验
数据处理的一个重要任务就是要判别所建立的模型是否有效,并从影响随机变量的诸变量中判别哪些变量的影响是显著的,哪些是不显著的。最后,利用最终所得到的经验公式进行预测和控制。
在进行数据平差处理时,由于附加参数所形成的扩展了的法方程式,将不一定仍能保持其解的稳定性,这是与附加参数的形式、个数及摄影条件有关的。在引进附加参数时,可以首先由一组比较复杂的参数开始。然后使用数理统计的方法,消除那些不能以足够精度加以确定的和不能足够保证可以区分的参数。或者相反,也可以根据经验,首先取用少量的认为最基本的一组参数。有了这些可靠的基础之后,再适当增补一些参数,并检查其显著性。下面给出基本的显著性检验方法。
当附加参数接近正交或正交时,可使用数理统计中的t分布,对所求得的参数,逐个进行显著性检验。
t分布是按下式定义的变量:
t = &zeta; &eta;
其中ζ为标准正态分布N(0,1);η定义为
Figure BDA0000154445880000342
变量,υ是χ2的自由度。此时统计假设为:
Figure BDA0000154445880000343
Figure BDA0000154445880000344
为第i个附加参数的估值。取
&zeta; = a ^ i - E ( a ^ i ) &sigma; 0 q ii ~ N ( 0,1 )
&eta; = ( n - &mu; ) s 0 2 &sigma; 0 2 / ( n - u ) ~ &chi; 2 v , v = n - u
其中
Figure BDA0000154445880000347
为单位权方差。
Figure BDA0000154445880000348
为由平差运算中得出的单位权中误差的平方,其期望值为
Figure BDA0000154445880000351
qii取自平差中未知数协因数阵Q的对角元素。由于假设
Figure BDA0000154445880000352
因此得出分布 t = &zeta; &eta; = a ^ i &sigma; 0 q ii = &sigma; 0 2 s 0 2 = a ^ i s 0 q ii , 当给定显著水平α,可由分布表查出临界值tα。若t<tα,则原假设成立,即该参数不显著,可在下次迭代平差中剔除。
当附加参数间的相关较大时,一维的t-检验会导致错误的结论。由于t的相关往往仅出现在某一组的附加参数中,所以应该把这一组参数放在一起,使用多维检验的方法。此时,原假设为:
为在一起检验的k个附加参数。
&xi; &prime; = [ a ^ - E ( a ^ ) ] T Q a ^ a ^ - 1 [ a ^ - E ( a ^ ) ] &sigma; 0 2 / k ~ &chi; 2 / v , v = k &eta; &prime; = ( n - u ) s 0 2 &sigma; 0 2 ( n - u ) ~ &chi; 2 / v , v = n - u , 因此得出F的统计量为:
根据两个自由度v(即k和(n-u))和假定显著水平α,可查F分布表。假如H0为真,则去除每个参数
Figure BDA0000154445880000358
当平差模型加入附加参数后,可能会由于附加参数之间或附加参数与其它未知参数之间的强相关而使得解的精度和可靠度恶化。为此应求出整个未知数的协因数阵,进而求出相关系数矩阵,进行逐一检查,以确定附加参数的可测定性检验、附加参数的外部可靠性检验及附加参数的相关性检验。
·内部精度
Sigma naught(记
Figure BDA0000154445880000359
)为平差后单位权标准差,由***残差向量v计算而得,
&sigma; ^ 0 = v T pv n - u
定义Qxx
Qxx=N-1    (31)
则协方差矩阵Kxx计算公式为
K xx = &sigma; ^ 0 Q xx
Kxx为对称矩阵,对角元素为
&sigma; ^ kk = &sigma; ^ 0 q kk
其描述的是某一未知数xk的平均方差。在(i,j)位置上的非对角线元素描述的是未知数xi,xj之间的协方差。
GCPs及TPs的平面及高程精度可分开按如下公式计算而得:
&sigma; x 2 = tr ( K xx x ) g &sigma; y 2 = tr ( K xx y ) g &sigma; z 2 = tr ( K xx z ) g
其中
σx,σy和σz:X、Y及Z方向上的标准差;
Figure BDA0000154445880000366
GCPs(或TPs)的X,Y,Z坐标的Kxx
g:点数。
·RMS计算
已知地面点坐标的TPs可以作为检查点使用(CPs)。方法是将地面坐标平差所得估值
Figure BDA0000154445880000367
与相应的“正确”坐标值(准真值)[Xcorr,Ycorr,Zcorr]相比较,RMSE按下式计算
RMSE x 2 = &Sigma; i = 1 N _ CP ( X ^ - X corr ) 2 N _ CP RMSE y 2 = &Sigma; i = 1 N _ CP ( Y ^ - Y corr ) 2 N _ CP RMSE z 2 = &Sigma; i = 1 N _ CP ( Z ^ - Z corr ) 2 N _ CP
其中N_CP为CPs点个数。
·相关系数计算
需要对附加参数的确定性进行研究,以避免过度参数化问题的产生。基于此,需计算未知数之间的相关系数ρij
&rho; ij = q ij q ii q jj
其中,qij为N-1中所定义矩阵(i,j)位置上的元素。ρij的绝对值越接近1即说明xi,xj之间的存在强相关性。如果存在强相关参数,则应使得这两者之一以伪观测量的形式作为常值出现。
·粗差探测
由于公式(3)中的“真值残差”不可得,所以使用Grün方法进行摄影测量区域平差像点坐标观测量粗差探测(Grün,1982)。该方法基于Baarda数据挖掘技术,即对于每个观测量i,值wi按以下公式计算:
w i = - v i &sigma; ^ v i
其中
Figure BDA0000154445880000372
Figure BDA0000154445880000373
为矩阵Qvv中的第i个元素,Qvv定义如下:
Q vv = P ll - 1 - A ( A T P ll A ) - 1 A T
如果零假设
H 0 i v : E ( v i ) = 0
为真,则wi服从τ分布。如果***残差足够大,则τ分布可用更易于处理的学生分布进行替代。粗差探测一般在平差结束后进行,人工剔除被检测出的观测量粗差后方可进行新一轮的平差处理。
本发明中所设计的传感器模型,对航空、航天线阵传感器均有很好的适用性及通用性,根本的原因在于外部定向所采用的函数有很好的适应性,其同时考虑了传感器自检校成像***的***误差建模。在建模过程中,主导思想是建立一个可适用于不同特性推扫式传感器的成像数学模型,并且根据不同的检校内容,仅估计感兴趣的参数。例如,可以根据平台的类型——卫星、飞机或直升机,以及轨道长度,还有地面控制点及连接点(每段轨道内至少有一个连接点)的分布,而确定轨道分段数。由于使用了伪观测量,从而使得对仅感兴趣的参数进行估计成为可能。实际上,任何参数都可以常数形式对待,只要赋予其很高的权重,当作伪观测量看待即可。这种灵活性可用来克服过度参数化对平差过程的影响。为了选择适当的兴趣参数,对参数间的相关性分析及显著性检验时必不可少的。
由于不同的研究内容待求解的参数个数很可能不同,所以平差所要求的最少地面控制点个数也不同。通常,待求参数个数主要取决于轨道PPM的分段数、PPM的阶数以及自检校参数个数等。考虑较为简单的情况:卫星平台上搭载双线阵传感器,轨道模型使用2阶多项式,并且无连接点,则未知数个数为36。由于关于PPM的连续性约束条件有18个方程,并且每个GCP给出4个观测方程,所以最少应有5个控制点参与平差。为了提高多余观测,通常建议给出6个GCP,当然,GCPs分布会影响到所要求的控制点数量。
由于所建立的模型有很好的适应性,所以,新相机很容易被整合到该模型中。实做中,可以从网上或文献中可以得到所需的传感器信息(线阵个数,镜头个数,焦距,每线阵的探测器件个数,视场角等)。如果传感器的内部参数不可知,则模型所假设的前提是CCD线阵在焦面上相互平行、并垂直于飞行方向,用自检校技术估计***误差。关于外部定向,在航天情况下,用于计算PPM参数初值的卫星位置及速度向量均包含在星厉数据中,如果卫星星厉数据不可用,则用开普勒参数估计名义上的轨道。按照本方法所做的大量试验表明,所发明的方法对多镜头航天线阵相机***在轨几何自检校适用性很好,对传感器的量测性能的挖掘和提升效果明显。

Claims (3)

1.一种多镜头航天线阵相机***在轨几何自检校方法,其特征在于:该方法包括以下几个步骤:
(1)根据传感器飞行轨迹和影像条带宽度,在检校区内布设高精度标志控制点;
(2)对获取检校区的前、下、后视相机三个影像条带进行控制点对应像点量测;
(3)获取传感器的轨道数据,包括GPS及星敏姿态数据观测量;
(4)根据多镜头航天线阵传感器成像模型和外部定向参数模型,建立观测方程组:
<1>、多镜头航天线阵传感器成像模型
光学***由多个镜头组成的传感器,必须要考虑描述其它镜头相对于下视镜头(平台)的位置与姿态关系,对于镜头j,用dxj,dyj,dzj表示相对于下视镜头的相对位置,用αj,βj,γj表示相对于下视镜头的相对姿态,fj,xpj,ypj为镜头j的焦距和主点坐标,成像模型为:
x - x pj = - f j &CenterDot; s 11 ( X - X C ) + s 21 ( Y - Y C ) + s 31 ( Z - Z C ) - ( m 11 d xj + m 21 d yj + m 31 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj ) + &Delta; xj
y - y pj = - f j &CenterDot; s 12 ( X - X C ) + s 22 ( Y - Y C ) + s 32 ( Z - Z C ) - ( m 12 d xj + m 22 d yj + m 32 d zj ) s 13 ( X - X C ) + s 23 ( Y - Y C ) + s 33 ( Z - Z C ) - ( m 13 d xj + m 23 d yj + m 33 d zj ) + &Delta; yj - - - ( 1 )
该模型对单镜头与多镜头都是适用的,其中:j=1,2,3,…,n,
[X,Y,Z]:物点在地面坐标***中的坐标;
[XC,YC,ZC]:下视相机投影中心在地面坐标***中的坐标;
[x,y]:像点在相机坐标***中的坐标,理论上x为0或某一常数;
[xp,yp]:像主点在相机坐标***中的坐标;
f:焦距;
Figure FDA0000154445870000021
从下视相机坐标***到地面坐标***的旋转矩阵,旋角分别为
Figure FDA0000154445870000022
M(αj,βj,γj)为从下视镜头相机坐标***到镜头j相机坐标***的旋转矩阵;
Figure FDA0000154445870000023
为从镜头j相机坐标***到地面坐标***的旋转矩阵;
并且,Δxj,Δyj是关于镜头畸变、CCD线阵畸变以及内定向参数(改正量)成像***的***误差模型,如下:
&Delta; x j = &Delta; x pj - &Delta;f f x &OverBar; pj + ( k 1 j r j 2 + k 2 j r j 4 ) x &OverBar; pj + p 1 j ( r j 2 + 2 x &OverBar; pj 2 ) + 2 p 2 j x &OverBar; pj y &OverBar; pj + x &OverBar; pj sin &theta;
&Delta; y j = &Delta; y pj - &Delta;f f y &OverBar; pj + ( k 1 j r j 2 + k 2 j r j 4 ) y &OverBar; pj + + 2 p 1 j x &OverBar; pj y &OverBar; pj + p 2 j ( r j 2 + 2 y &OverBar; pj 2 ) + s yj y &OverBar; pj - - - ( 2 )
参数说明:
·每个镜头的主点改正量(位移量)Δxpj,Δypj
·k1j,k2j,p1j,p2j:每个镜头的对称性径向和切向畸变;
·syj:每个线阵CCD单元y方向上的比例因子;
·每个CCD线阵在焦面上的旋角θ(其产生的畸变主要在x方向,在y方向上的畸变影响可以忽略不计), x &OverBar; pj = x - x pj , y &OverBar; pj = y - y pj , r 2 = x &OverBar; pj 2 + y &OverBar; pj 2 ;
<2>、外部定向参数模型
外部定向参数模型——轨道模型(trajectory model)是用来描述传感器外部定向(EO,external orientation)状态的,其是用关于时间的分段多项式(PPM,PiecewisePolynomial Function)来建模的,可以根据控制点(GCPs)和连接点(TPs)的分布来确定轨道模型的分段数,对于被分割轨道的第i段,该段的首末端点的线阵影像获取时刻记为
Figure FDA0000154445870000029
则定义时间变量
Figure FDA00001544458700000210
为:
t &OverBar; = l - l ini i l fin i - l ini i - - - ( 3 )
其中,
Figure FDA0000154445870000031
为第i段轨道的首末端在对应影像中的线阵行数(首末线阵影像行),这样,对于第i段轨道所对应影像中的第1行线阵影像,其位置与姿态
Figure FDA0000154445870000032
的数学模型可以表示为该行线阵影像的外部定向参数观测量
Figure FDA0000154445870000033
与关于时间
Figure FDA0000154445870000034
的2阶多项式的和,如下式:
X C ( t &OverBar; ) = X instr + X 0 i + X 1 i t &OverBar; + X 2 i t &OverBar; 2
Y C ( t &OverBar; ) = Y instr + Y 0 i + Y 1 i t &OverBar; + Y 2 i t &OverBar; 2 - - - ( 4 )
Z C ( t &OverBar; ) = Z instr + Z 0 i + Z 1 i t &OverBar; + Z 2 i t &OverBar; 2
&omega; C ( t &OverBar; ) = &omega; instr + &omega; 0 i + &omega; 1 i t &OverBar; + &omega; 2 i t &OverBar; 2
Figure FDA0000154445870000039
&kappa; C ( t &OverBar; ) = &kappa; instr + &kappa; 0 i + &kappa; 1 i t &OverBar; + &kappa; 2 i t &OverBar; 2
其中,[X0,X1,X2,…,κ0,κ1,κ2]i是第i段轨道关于GPS偏移量及INS(星敏)轴线偏差、GPS及INS(星敏)***漂移等***误差的18个模型参数;
<3>、观测方程组
下面将详细列出关于上述传感器模型的观测方程组——包括像坐标观测方程、外部定向参数观测方程(连续性条件方程及外部定向参数自身观测方程)、自检校参数观测方程、地面点观测方程,其中方程组中的所有观测量都带有上划线记号,待平差的定向参数无上划线标记:
x &OverBar; + v x = x = x pj - f j &CenterDot; N x D + &Delta; xj
y &OverBar; + v y = y = y pj - f j &CenterDot; N y D + &Delta; yj - - - ( 5 )
——像坐标观测量方程
( X 0 i &OverBar; + X 1 i &OverBar; + X 2 i &OverBar; - X 0 i + 1 &OverBar; ) + v X C 0 = X 0 i + X 1 i + X 2 i - X 0 i + 1
( Y 0 i &OverBar; + Y 1 i &OverBar; + Y 2 i &OverBar; - Y 0 i + 1 &OverBar; ) + v Y C 0 = Y 0 i + Y 1 i + Y 2 i - Y 0 i + 1
( Z 0 i &OverBar; + Z 1 i &OverBar; + Z 2 i &OverBar; - Z 0 i + 1 &OverBar; ) + v Z C 0 = Z 0 i + Z 1 i + Z 2 i - Z 0 i + 1
( &omega; 0 i &OverBar; + &omega; 1 i &OverBar; + &omega; 2 i &OverBar; - &omega; 0 i + 1 &OverBar; ) + v &omega; C 0 = &omega; 0 i + &omega; 1 i + &omega; 2 i - &omega; 0 i + 1
( &kappa; 0 i &OverBar; + &kappa; 1 i &OverBar; + &kappa; 2 i &OverBar; - &kappa; 0 i + 1 &OverBar; ) + v &kappa; C 0 = &kappa; 0 i + &kappa; 1 i + &kappa; 2 i - &kappa; 0 i + 1 - - - ( 6 )
——PPM连续性约束条件(0阶)
( X 1 i &OverBar; + 2 X 2 i &OverBar; - X 1 i + 1 &OverBar; ) + v X C 1 = X 1 i + 2 X 2 i - X 1 i + 1
( Y 1 i &OverBar; + 2 Y 2 i &OverBar; - Y 1 i + 1 &OverBar; ) + v Y C 1 = Y 1 i + 2 Y 2 i - Y 1 i + 1
( Z 1 i &OverBar; + 2 Z 2 i &OverBar; - Z 1 i + 1 &OverBar; ) + v Z C 1 = Z 1 i + 2 Z 2 i - Z 1 i + 1
( &omega; 1 i &OverBar; + 2 &omega; 2 i &OverBar; - &omega; 1 i + 1 &OverBar; ) + v &omega; C 1 = &omega; 1 i + 2 &omega; 2 i - &omega; 1 i + 1
Figure FDA00001544458700000411
( &kappa; 1 i &OverBar; + 2 &kappa; 2 i &OverBar; - &kappa; 1 i + 1 &OverBar; ) + v &kappa; C 1 = &kappa; 1 i + 2 &kappa; 2 i - &kappa; 1 i + 1 - - - ( 7 )
——PPM连续性约束条件(1阶)
( X 2 i &OverBar; - X 2 i + 1 &OverBar; ) + v X C 2 = X 2 i - X 2 i + 1
( Y 2 i &OverBar; - Y 2 i + 1 &OverBar; ) + v Y C 2 = Y 2 i - Y 2 i + 1
( Z 2 i &OverBar; - Z 2 i + 1 &OverBar; ) + v Z C 2 = Z 2 i - Z 2 i + 1
( &omega; 2 i &OverBar; - &omega; 2 i + 1 &OverBar; ) + v &omega; C 2 = &omega; 2 i - &omega; 2 i + 1
Figure FDA00001544458700000417
( &kappa; 2 i &OverBar; - &kappa; 2 i + 1 &OverBar; ) + v &kappa; C 2 = &kappa; 2 i - &kappa; 2 i + 1 - - - ( 8 )
——PPM连续性约束条件(2阶)
X 0 &OverBar; + v X 0 = X 0 Y 0 &OverBar; + v Y 0 = Y 0 Z 0 &OverBar; + v Z 0 = Z 0
X 1 &OverBar; + v X 1 = X 1 Y 1 &OverBar; + v Y 1 = Y 1 Z 1 &OverBar; + v Z 1 = Z 1 - - - ( 9 )
X 2 &OverBar; + v X 2 = X 2 Y 2 &OverBar; + v Y 2 = Y 2 Z 2 &OverBar; + v Z 2 = Z 2
&omega; 0 &OverBar; + v &omega; 0 = &omega; 0
Figure FDA00001544458700000511
&kappa; 0 &OverBar; + v &kappa; 0 = &kappa; 0
&omega; 1 &OverBar; + v &omega; 1 = &omega; 1
Figure FDA00001544458700000514
&kappa; 1 &OverBar; + v &kappa; 1 = &kappa; 1
&omega; 2 &OverBar; + v &omega; 2 = &omega; 2
Figure FDA00001544458700000517
&kappa; 2 &OverBar; + v &kappa; 2 = &kappa; 2
——外部定向参数自身观测方程
&Delta; x P &OverBar; + v &Delta; x P = &Delta; x P
&Delta; y P &OverBar; + v &Delta; y P = &Delta; y P
&Delta;f &OverBar; + v &Delta;f = &Delta;f
k 1 &OverBar; + v k 1 = k 1
k 2 &OverBar; + v k 2 = k 2
p 1 &OverBar; + v p 1 = p 1
p 2 &OverBar; + v p 2 = p 2
s y &OverBar; + v s y = s y - - - ( 10 )
&theta; &OverBar; + v &theta; = &theta;
——自检校参数观测方程
X &OverBar; + v X = X
Y &OverBar; + v Y = Y - - - ( 11 )
Z &OverBar; + v Z = Z
——地面点观测方程
(5)利用最小二乘法进行光束法自检校平差数据处理,求解几何***误差检校参数估值——首先进行线性化、系数阵构造、法化后线性***求解;
(6)对所解得参数进行显著性检验,并进行数据分析,包括内部精度、RMS、相关系数和粗差探测计算。
2.如权利要求1所述的多镜头航天线阵相机***在轨几何自检校方法,其特征在于:所述的利用最小二乘法进行光束法自检校平差数据处理,求解几何***误差检校参数估值——首先进行线性化、系数阵构造、法化后线性***求解如下:
<1>、线性化方程组(偏导数)
对成像模型与轨道模型关于未知数求偏导,使之线性化,以便求解:
&delta;x &delta; X 0 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 0 &delta;x &delta; X 1 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 1 &delta;x &delta; X 2 = &delta;x &delta; X C &CenterDot; &delta; X C &delta; X 2
&delta;x &delta; Y 0 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 0 &delta;x &delta; Y 1 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 1 &delta;x &delta; Y 2 = &delta;x &delta; Y C &CenterDot; &delta; Y C &delta; Y 2
&delta;x &delta; Z 0 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Z 0 &delta;x &delta; Z 1 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Y 1 &delta;x &delta; Z 2 = &delta;x &delta; Z C &CenterDot; &delta; Z C &delta; Y 2
(12)
&delta;x &delta; &omega; 0 = &delta;x &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 0 &delta;x &delta; &omega; 1 = &delta;x &delta;&omega; C &CenterDot; &delta; &omega; C &delta; &omega; 1 &delta;x &delta; &omega; 2 = &delta;x &delta;&omega; C &CenterDot; &delta;&omega; C &delta;&omega; 2
Figure FDA00001544458700000614
Figure FDA00001544458700000615
&delta;x &delta; &kappa; 0 = &delta;x &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 0 &delta;x &delta; &kappa; 1 = &delta;x &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 1 &delta;x &delta; &kappa; 2 = &delta;x &delta; &kappa; C &CenterDot; &delta;&kappa; C &delta; &kappa; 2
——成像模型中x关于外部定向参数的偏导数
&delta;y &delta; X 0 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 0 &delta;y &delta; X 1 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 1 &delta;y &delta; X 2 = &delta;y &delta; X C &CenterDot; &delta; X C &delta; X 2
&delta;y &delta; Y 0 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 0 &delta;y &delta; Y 1 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 1 &delta;y &delta; Y 2 = &delta;y &delta; Y C &CenterDot; &delta; Y C &delta; Y 2
&delta;y &delta; Z 0 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Z 0 &delta;y &delta; Z 1 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Y 1 &delta;y &delta; Z 2 = &delta;y &delta; Z C &CenterDot; &delta; Z C &delta; Y 2
(13)
&delta;y &delta; &omega; 0 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 0 &delta;y &delta; &omega; 1 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 1 &delta;y &delta; &omega; 2 = &delta;y &delta; &omega; C &CenterDot; &delta; &omega; C &delta; &omega; 2
Figure FDA00001544458700000632
Figure FDA00001544458700000633
&delta;y &delta; &kappa; 0 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 0 &delta;y &delta; &kappa; 1 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 1 &delta;y &delta; &kappa; 2 = &delta;y &delta; &kappa; C &CenterDot; &delta; &kappa; C &delta; &kappa; 2
——成像模型中y关于外部定向参数的偏导数
其中,
Figure FDA0000154445870000071
Figure FDA0000154445870000072
并且,
&PartialD; x &PartialD; X C = - f D 2 ( s 13 N x - s 11 D ) = a 1,1
&PartialD; x &PartialD; Y C = - f D 2 ( s 23 N x - s 21 D ) = a 1,2
&PartialD; x &PartialD; Z C = - f D 2 ( s 33 N x - s 31 D ) = a 1,3
&PartialD; x &PartialD; &omega; C = f D { [ ( X - X C ) &delta;s 13 &PartialD; &omega; C + ( Y - Y C ) &delta;s 23 &PartialD; &omega; C - ( Z - Z C ) &delta;s 33 &PartialD; &omega; C ] N x D - ( X - X 0 ) &delta;s 11 &PartialD; &omega; C - ( Y - Y 0 ) &delta;s 21 &PartialD; &omega; C - ( Z - Z 0 ) &delta;s 31 &PartialD; &omega; C } = a 1,4
Figure FDA0000154445870000078
&PartialD; x &PartialD; &kappa; C = f D { [ ( X - X C ) &delta;s 13 &PartialD; &kappa; C + ( Y - Y C ) &delta; s 33 &PartialD; &kappa; C - ( Z - Z C ) &delta;s 33 &PartialD; &kappa; C ] N x D - ( X - X 0 ) &delta;s 11 &PartialD; &kappa; C - ( Y - Y 0 ) &delta;s 21 &PartialD; &kappa; C - ( Z - Z 0 ) &delta;s 31 &delta;&kappa; C } = a 1,6
&PartialD; y &PartialD; X C = - f N 2 ( r 13 Z y - r 13 N ) = a 2,1
&PartialD; y &PartialD; Y C = - f N 2 ( r 23 Z y - r 22 N ) = a 2,2
&PartialD; y &PartialD; Z C = - f N 2 ( r 33 Z y - r 32 N ) = a 2,3
&PartialD; y &PartialD; &omega; C = f N { [ ( Y - Y 0 ) r 33 - ( Z - Z 0 ) r 23 ] Z y N - ( Y - Y 0 ) r 32 - ( Z - Z 0 ) r 22 } = a 2,7
Figure FDA00001544458700000714
&PartialD; y &PartialD; &kappa; C = - f N Z x = a 2,9
系数ai,j指的是在矩阵AGCP和ATP对应的元素(AGCP、ATP定义见下文);
&PartialD; x &PartialD; X = - f D 2 ( Ds 11 - N x s 13 ) = b 1,1 &PartialD; y &PartialD; X = - f D 2 ( Ds 12 - N y s 13 ) = b 2,1 - - - ( 14 )
&PartialD; x &PartialD; Y = - f D 2 ( Ds 21 - N x s 23 ) = b 1,2 &PartialD; y &PartialD; Y = - f D 2 ( Ds 22 - N y s 23 ) = b 2,2
&PartialD; x &PartialD; Z = - f D 2 ( Ds 31 - N x s 33 ) = b 1,3 &PartialD; y &PartialD; Z = - f D 2 ( Ds 32 - N y s 33 ) = b 2,3
——成像模型关于GCPs及TPs的偏导数
系数bi,j指的是在矩阵BGCP和BTP对应的元素(BGCP、BTP定义见下文);
&PartialD; x &PartialD; &theta; = y &OverBar; p cos &theta; = s 1,9 &PartialD; y &PartialD; &theta; = 0 = s 2,9 - - - ( 15 )
——成像模型关于自检校未知数的偏导数
除了θ外,其它参数均为线性函数,系数si,j指的是在矩阵S中对应的元素(S定义见下文),参数Δxp,Δyp,ΔxC,ΔyC,c,k1,k2,p1,p2和sy的si,j值见下表:
Figure FDA0000154445870000083
<2>、系数阵构造
对于单轨道单条带区域影像,系数矩阵记为A,平差模型中的解向量x包含:
◆xEO:PPM参数(18×N_S)
◆xTP:TPs地面点坐标(3×N_TP)
◆xGCP:GCP地面点坐标(3×N_GCP)
◆xSC:自检校参数(单镜头光学***:5+3×N_LINES;多镜头光学***:8×N_LINES)
未知数向量x及系数矩阵形式为:
x = x EO x TP x GCP x SC ; A = A GCP S GCP A TP B TP S TP C 0 C 1 C 2 F E S
观测方程组即为:
Figure FDA0000154445870000093
其中,矩阵Pxx是一类观测量的权重矩阵,其可进一步构造出所有观测量的整体权重对角矩阵,如下:
P ll = P GCP P TP P C 0 P C 1 P C 2 P EO P E P S
记n为观测量个数,u为未知数个数,则未知数向量x,系数阵A及权重矩阵Pll的维数分别为[u,1],[n,u]以及[n,n],法方程矩阵N为一个维数是[u×u]的对称系数阵;
<3>、线性***求解
至此,得线性***
Nx=z                  (17)
需求解未知数向量x,法方程矩阵N是一个块状矩阵,每一个单元块都有着不同的结构,如:对角阵、超对角阵、零矩阵以及稀疏阵,该矩阵较为复杂,它根据建模中的参数个数不同而有着不同的子矩阵数量,采用LU因式分解方法进行解算,即法方程矩阵被分解为:
N=LU                  (18)
其中,L为下三角阵(元素αij仅位于矩阵对角线及以下位置),U为上三角阵(元素βij仅位于矩阵对角线及以上位置),对于一个[n,n]大小的矩阵,(18)式变为:
n 11 n 12 &CenterDot; &CenterDot; &CenterDot; n 1 n n 21 n 22 &CenterDot; &CenterDot; &CenterDot; n 2 n &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; n n 1 n n 2 &CenterDot; &CenterDot; &CenterDot; n nn = &alpha; 11 0 &CenterDot; &CenterDot; &CenterDot; 0 &alpha; 21 &alpha; 22 &CenterDot; &CenterDot; &CenterDot; 0 &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &alpha; n 1 &alpha; n 2 &CenterDot; &CenterDot; &CenterDot; &alpha; nn &beta; 11 &beta; 12 &CenterDot; &CenterDot; &CenterDot; &beta; 1 n 0 &beta; 22 &CenterDot; &CenterDot; &CenterDot; &beta; 2 n &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; &CenterDot; 0 0 &CenterDot; &CenterDot; &CenterDot; &beta; nn - - - ( 19 )
使用这样的分解式,线性集变换为:
Nx=(LU)x=L(Ux)=z    (20)
这样,可以分两步解算,首先,解算
Ly=z                  (21)
然后,解算
Ux=y                  (22)
之所以这样进行解算,是因为解算三角系数阵线性***是比较简单的,方程(21)可以进行向前替代解算,如下:
y 1 = z 1 &alpha; 11 - - - ( 23 )
y i = 1 &alpha; ii + [ z i - &Sigma; j = 1 i - 1 &alpha; ij y j , i = 2,3 , &CenterDot; &CenterDot; &CenterDot; n
同样,方程(28)可以进行向后替代解算,如下:
x n = y n &beta; nn - - - ( 24 ) .
x i = 1 &beta; ii + [ y i - &Sigma; j = i + 1 n &beta; ij x j , i = n - 1 , n - 2 , &CenterDot; &CenterDot; &CenterDot; 1
3.如权利要求1所述的多镜头航天线阵相机***在轨几何自检校方法,其特征在于:对所解得参数进行显著性检验,并进行数据分析,包括内部精度、RMS、相关系数和粗差探测计算如下:
<1>、显著性检验
下面给出基本的显著性检验方法,
当附加参数接近正交或正交时,可使用数理统计中的t分布,对所求得的参数,逐个进行显著性检验,t分布是按下式定义的变量:
t = &zeta; &eta;
其中ζ为标准正态分布N(0,1);η定义为变量,υ是χ2的自由度,此时统计假设为:
Figure FDA0000154445870000115
Figure FDA0000154445870000116
为第i个附加参数的估值,取
&zeta; = a ^ i - E ( a ^ i ) &sigma; 0 q ii ~ N ( 0,1 )
&eta; = ( n - &mu; ) s 0 2 &sigma; 0 2 / ( n - u ) ~ &chi; 2 v , v = n - u
其中
Figure FDA0000154445870000119
为单位权方差,为由平差运算中得出的单位权中误差的平方,其期望值为
Figure FDA00001544458700001111
qii取自平差中未知数协因数阵Q的对角元素,由于假设
Figure FDA00001544458700001112
因此得出分布 t = &zeta; &eta; = a ^ i &sigma; 0 q ii = &sigma; 0 2 s 0 2 = a ^ i s 0 q ii , 当给定显著水平α,可由分布表查出临界值tα,若t<tα,则原假设成立,即该参数不显著,可在下次迭代平差中剔除,当附加参数间的相关较大时,一维的t-检验会导致错误的结论,由于t的相关往往仅出现在某一组的附加参数中,所以应该把这一组参数放在一起,使用多维检验的方法,此时,原假设为:
Figure FDA0000154445870000121
Figure FDA0000154445870000122
为在一起检验的k个附加参数,取 &xi; &prime; = [ a ^ - E ( a ^ ) ] T Q a ^ a ^ - 1 [ a ^ - E ( a ^ ) ] &sigma; 0 2 / k ~ &chi; 2 / v , v = k &eta; &prime; = ( n - u ) s 0 2 &sigma; 0 2 ( n - u ) ~ &chi; 2 / v , v = n - u , 因此得出F的统计量为:
根据两个自由度v(即k和(n-u))和假定显著水平α,可查F分布表,假如H0为真,则去除每个参数
Figure FDA0000154445870000125
当平差模型加入附加参数后,可能会由于附加参数之间或附加参数与其它未知参数之间的强相关而使得解的精度和可靠度恶化,为此应求出整个未知数的协因数阵,进而求出相关系数矩阵,进行逐一检查,以确定附加参数的可测定性检验、附加参数的外部可靠性检验及附加参数的相关性检验;
<2>、内部精度
Sigma naught(记)为平差后单位权标准差,由***残差向量v计算而得,
&sigma; ^ 0 = v T pv n - u
定义Qxx
Qxx=N-1    (25)
则协方差矩阵Kxx计算公式为
K xx = &sigma; ^ 0 Q xx
Kxx为对称矩阵,对角元素为
&sigma; ^ kk = &sigma; ^ 0 q kk
其描述的是某一未知数xk的平均方差,在(i,j)位置上的非对角线元素描述的是未知数xi,xj之间的协方差,GCPs及TPs的平面及高程精度可分开按如下公式计算而得:
&sigma; x 2 = tr ( K xx x ) g &sigma; y 2 = tr ( K xx y ) g &sigma; z 2 = tr ( K xx z ) g
其中
σx,σy和σz:X、Y及Z方向上的标准差;
Figure FDA0000154445870000134
GCPs(或TPs)的X,Y,Z坐标的Kxx
g:点数;
<3>、RMS计算
已知地面点坐标的TPs可以作为检查点使用(CPs),方法是将地面坐标平差所得估值与相应的“正确”坐标值(准真值)[Xcorr,Ycorr,Zcorr]相比较,RMSE按下式计算
RMSE x 2 = &Sigma; i = 1 N _ CP ( X ^ - X corr ) 2 N _ CP RMSE y 2 = &Sigma; i = 1 N _ CP ( Y ^ - Y corr ) 2 N _ CP RMSE z 2 = &Sigma; i = 1 N _ CP ( Z ^ - Z corr ) 2 N _ CP
其中N_CP为CPs点个数;
<4>、相关系数计算
需要对附加参数的确定性进行研究,以避免过度参数化问题的产生,基于此,需计算未知数之间的相关系数ρij
&rho; ij = q ij q ii q jj
其中,qij为N-1中所定义矩阵(i,j)位置上的元素,ρij的绝对值越接近1即说明xi,xj之间的存在强相关性,如果存在强相关参数,则应使得这两者之一以伪观测量的形式作为常值出现;
<5>、粗差探测
由于公式(1)中的“真值残差”不可得,所以使用Grün方法进行摄影测量区域平差像点坐标观测量粗差探测,该方法基于Baarda数据挖掘技术,即对于每个观测量i,值wi按以下公式计算:
w i = - v i &sigma; ^ v i
其中
Figure FDA0000154445870000142
Figure FDA0000154445870000143
为矩阵Qvv中的第i个元素,Qvv定义如下:
Q vv = P ll - 1 - A ( A T P ll A ) - 1 A T
如果零假设
H 0 i v : E ( v i ) = 0
为真,则wi服从τ分布,如果***残差足够大,则τ分布可用更易于处理的学生分布进行替代,粗差探测一般在平差结束后进行,人工剔除被检测出的观测量粗差后方可进行新一轮的平差处理。
CN201210113129.4A 2012-04-18 2012-04-18 多镜头航天线阵相机***在轨几何自检校方法 Expired - Fee Related CN102636159B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210113129.4A CN102636159B (zh) 2012-04-18 2012-04-18 多镜头航天线阵相机***在轨几何自检校方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210113129.4A CN102636159B (zh) 2012-04-18 2012-04-18 多镜头航天线阵相机***在轨几何自检校方法

Publications (2)

Publication Number Publication Date
CN102636159A true CN102636159A (zh) 2012-08-15
CN102636159B CN102636159B (zh) 2014-05-07

Family

ID=46620649

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210113129.4A Expired - Fee Related CN102636159B (zh) 2012-04-18 2012-04-18 多镜头航天线阵相机***在轨几何自检校方法

Country Status (1)

Country Link
CN (1) CN102636159B (zh)

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103115627A (zh) * 2013-01-21 2013-05-22 武汉大学 遥感卫星线阵传感器多轨联合在轨几何检校方法
CN103148870A (zh) * 2013-03-01 2013-06-12 国家测绘地理信息局卫星测绘应用中心 基于高精度配准控制点的卫星ccd阵列影像几何检校方法
CN103679673A (zh) * 2013-11-22 2014-03-26 中国资源卫星应用中心 一种宽视场线阵ccd影像几何畸变模拟方法
CN103822644A (zh) * 2014-02-14 2014-05-28 北京航天控制仪器研究所 一种三维激光成像***的相机标定方法
CN104729529A (zh) * 2013-12-24 2015-06-24 北京市测绘设计研究院 地形图测量***误差判断的方法和***
CN104807477A (zh) * 2015-04-24 2015-07-29 国家测绘地理信息局卫星测绘应用中心 一种基于靶标控制点的卫星ccd阵列影像几何检校方法
CN105931248A (zh) * 2016-05-06 2016-09-07 西安航天天绘数据技术有限公司 卫星影像定位的方法和装置
CN106500729A (zh) * 2016-11-29 2017-03-15 武汉大学 一种无需控制信息的智能手机自检校方法
CN106643669A (zh) * 2016-11-22 2017-05-10 北京空间机电研究所 一种多镜头多探测器航空相机单中心投影转换方法
CN106885585A (zh) * 2016-12-30 2017-06-23 国家测绘地理信息局卫星测绘应用中心 一种基于光束法平差的星载摄影测量***一体化检校方法
CN107991558A (zh) * 2017-11-23 2018-05-04 国网福建省电力有限公司 基于t分布检验法的数字校准方法
CN108303117A (zh) * 2017-01-12 2018-07-20 中国农业大学 一种基于后方交会测量的云镜摄***参数测量方法及***
CN109724625A (zh) * 2019-01-22 2019-05-07 中国人民解放军61540部队 一种光学复合大面阵测绘相机的像差改正方法
CN111272196A (zh) * 2020-02-29 2020-06-12 武汉大学 特定拍摄条件下的在轨外方位元素自检校方法及***
CN112070843A (zh) * 2020-08-04 2020-12-11 北京空间机电研究所 一种空间相机几何参数在轨标定方法
CN112634371A (zh) * 2019-09-24 2021-04-09 北京百度网讯科技有限公司 用于输出信息、标定相机的方法和装置
WO2021185219A1 (zh) * 2020-03-16 2021-09-23 左忠斌 一种用于太空领域的3d采集和尺寸测量方法
CN117606447A (zh) * 2023-09-25 2024-02-27 中国人民解放军61540部队 视场分割型航空面阵相机的试验场检校方法和装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101644570A (zh) * 2009-09-17 2010-02-10 北京空间机电研究所 航天三线阵ccd相机视主点在轨监测方法
CN101858755A (zh) * 2010-06-01 2010-10-13 北京控制工程研究所 一种星敏感器的标定方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101644570A (zh) * 2009-09-17 2010-02-10 北京空间机电研究所 航天三线阵ccd相机视主点在轨监测方法
CN101858755A (zh) * 2010-06-01 2010-10-13 北京控制工程研究所 一种星敏感器的标定方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
涂辛茹等: "机载三线阵传感器ADS40的几何验校", 《测绘学报》 *

Cited By (29)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103115627A (zh) * 2013-01-21 2013-05-22 武汉大学 遥感卫星线阵传感器多轨联合在轨几何检校方法
CN103148870A (zh) * 2013-03-01 2013-06-12 国家测绘地理信息局卫星测绘应用中心 基于高精度配准控制点的卫星ccd阵列影像几何检校方法
CN103148870B (zh) * 2013-03-01 2015-07-08 国家测绘地理信息局卫星测绘应用中心 基于高精度配准控制点的卫星ccd阵列影像几何检校方法
CN103679673B (zh) * 2013-11-22 2016-06-08 中国资源卫星应用中心 一种宽视场线阵ccd影像几何畸变模拟方法
CN103679673A (zh) * 2013-11-22 2014-03-26 中国资源卫星应用中心 一种宽视场线阵ccd影像几何畸变模拟方法
CN104729529A (zh) * 2013-12-24 2015-06-24 北京市测绘设计研究院 地形图测量***误差判断的方法和***
CN104729529B (zh) * 2013-12-24 2017-11-10 北京市测绘设计研究院 地形图测量***误差判断的方法和***
CN103822644B (zh) * 2014-02-14 2016-08-17 北京航天控制仪器研究所 一种三维激光成像***的相机标定方法
CN103822644A (zh) * 2014-02-14 2014-05-28 北京航天控制仪器研究所 一种三维激光成像***的相机标定方法
CN104807477A (zh) * 2015-04-24 2015-07-29 国家测绘地理信息局卫星测绘应用中心 一种基于靶标控制点的卫星ccd阵列影像几何检校方法
CN105931248A (zh) * 2016-05-06 2016-09-07 西安航天天绘数据技术有限公司 卫星影像定位的方法和装置
CN106643669A (zh) * 2016-11-22 2017-05-10 北京空间机电研究所 一种多镜头多探测器航空相机单中心投影转换方法
CN106643669B (zh) * 2016-11-22 2018-10-19 北京空间机电研究所 一种多镜头多探测器航空相机单中心投影转换方法
CN106500729A (zh) * 2016-11-29 2017-03-15 武汉大学 一种无需控制信息的智能手机自检校方法
CN106500729B (zh) * 2016-11-29 2019-05-17 武汉大学 一种无需控制信息的智能手机自检校方法
CN106885585B (zh) * 2016-12-30 2020-01-21 自然资源部国土卫星遥感应用中心 一种基于光束法平差的星载摄影测量***一体化检校方法
CN106885585A (zh) * 2016-12-30 2017-06-23 国家测绘地理信息局卫星测绘应用中心 一种基于光束法平差的星载摄影测量***一体化检校方法
CN108303117A (zh) * 2017-01-12 2018-07-20 中国农业大学 一种基于后方交会测量的云镜摄***参数测量方法及***
CN108303117B (zh) * 2017-01-12 2020-06-02 中国农业大学 一种基于后方交会测量的云镜摄***参数测量方法及***
CN107991558A (zh) * 2017-11-23 2018-05-04 国网福建省电力有限公司 基于t分布检验法的数字校准方法
CN109724625A (zh) * 2019-01-22 2019-05-07 中国人民解放军61540部队 一种光学复合大面阵测绘相机的像差改正方法
CN109724625B (zh) * 2019-01-22 2021-05-04 中国人民解放军61540部队 一种光学复合大面阵测绘相机的像差改正方法
CN112634371A (zh) * 2019-09-24 2021-04-09 北京百度网讯科技有限公司 用于输出信息、标定相机的方法和装置
CN112634371B (zh) * 2019-09-24 2023-12-15 阿波罗智联(北京)科技有限公司 用于输出信息、标定相机的方法和装置
CN111272196A (zh) * 2020-02-29 2020-06-12 武汉大学 特定拍摄条件下的在轨外方位元素自检校方法及***
WO2021185219A1 (zh) * 2020-03-16 2021-09-23 左忠斌 一种用于太空领域的3d采集和尺寸测量方法
CN112070843A (zh) * 2020-08-04 2020-12-11 北京空间机电研究所 一种空间相机几何参数在轨标定方法
CN112070843B (zh) * 2020-08-04 2024-03-15 北京空间机电研究所 一种空间相机几何参数在轨标定方法
CN117606447A (zh) * 2023-09-25 2024-02-27 中国人民解放军61540部队 视场分割型航空面阵相机的试验场检校方法和装置

Also Published As

Publication number Publication date
CN102636159B (zh) 2014-05-07

Similar Documents

Publication Publication Date Title
CN102636159B (zh) 多镜头航天线阵相机***在轨几何自检校方法
CN104931022A (zh) 基于星载激光测高数据的卫星影像立体区域网平差方法
CN103557841B (zh) 一种提高多相机合成影像摄影测量精度的方法
Leprince et al. In-flight CCD distortion calibration for pushbroom satellites based on subpixel correlation
CN101750619B (zh) 自检校pos直接对地目标定位方法
CN107527328B (zh) 一种兼顾精度与速度的无人机影像几何处理方法
CN102693542A (zh) 一种影像特征匹配方法
KR20070096370A (ko) Los벡터 조정모델을 이용한 영상의 기하보정 방법 및 그장치
CN102778224A (zh) 一种基于极坐标参数化的航空摄影测量光束法平差的方法
Wang et al. Planar block adjustment and orthorectification of ZY-3 satellite images
CN106526593A (zh) 基于sar严密成像模型的子像素级角反射器自动定位方法
CN113514829A (zh) 面向InSAR的初始DSM的区域网平差方法
Mostafa et al. Airborne direct georeferencing of frame imagery: An error budget
CN102721410A (zh) 一种基于gps/imu定位定向技术的海岛空中三角测量方法
CN102538764A (zh) 一种复合式像对立体定位方法
Kedzierski et al. Detection of gross errors in the elements of exterior orientation of low-cost UAV images
Mostafa Boresight calibration of integrated inertial/camera systems
Rokhmana et al. Cadastral surveys with non-metric camera using UAV: a feasibility study
Zhang et al. Fusion of ascending and descending polarimetric SAR data for color orthophoto generation
Takaku et al. An Overview of Geometric Calibration and DSM Generation for ALOS-3 Optical Imageries
Mitishita et al. Approach for improving the integrated sensor orientation
Seo et al. KOMPSAT-3A direct georeferencing mode and geometric Calibration/Validation
Haase et al. Bundle adjustment of spaceborne double-camera push-broom imagers and its application to LROC NAC imagery
Deltsidis et al. Orthorectification of World View 2 stereo pair using a new rigorous orientation model
CN117788281B (zh) 低pos精度的机载线阵高光谱影像区域影像拼接方法及***

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20140507

Termination date: 20170418