CN109143295A - 一种数字化几何定标场和gcp相结合的内方位元素定标方法 - Google Patents
一种数字化几何定标场和gcp相结合的内方位元素定标方法 Download PDFInfo
- Publication number
- CN109143295A CN109143295A CN201811270813.7A CN201811270813A CN109143295A CN 109143295 A CN109143295 A CN 109143295A CN 201811270813 A CN201811270813 A CN 201811270813A CN 109143295 A CN109143295 A CN 109143295A
- Authority
- CN
- China
- Prior art keywords
- calibration
- gcp
- parameter
- internal
- coordinate system
- 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 51
- 238000005457 optimization Methods 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 64
- 238000012937 correction Methods 0.000 claims description 46
- 239000000523 sample Substances 0.000 claims description 44
- 238000005259 measurement Methods 0.000 claims description 25
- 238000005096 rolling process Methods 0.000 claims description 9
- 238000003860 storage Methods 0.000 claims description 5
- 230000003044 adaptive effect Effects 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 3
- 238000012880 independent component analysis Methods 0.000 claims description 3
- 238000004590 computer program Methods 0.000 claims description 2
- 238000010276 construction Methods 0.000 abstract description 6
- 238000012545 processing Methods 0.000 abstract description 6
- 230000000694 effects Effects 0.000 abstract description 3
- 230000008685 targeting Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 230000032683 aging Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000011524 similarity measure Methods 0.000 description 1
- 238000002187 spin decoupling employing ultra-broadband-inversion sequences generated via simulated annealing Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/42—Determining position
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Length Measuring Devices By Optical Means (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
一种数字化几何定标场和GCP相结合的内方位元素定标方法，包括高精度外业地面控制点布设、采集及人工内业立体环境刺点方案、利用数字化几何定标场数据采用密集匹配的方法进行控制点的自动量测、探元指向角相机内定校模型构建及模型参数求解策略和方法、相机内定标参数优化模型构建及模型参数求解。经实践验证表明该内定标方法可行、有效，结果稳定、可靠、内定标精度高，目前已工程化应用到商遥卫星地面数据处理系统中，运行稳定，取得了良好的工程化应用效果。
Description
技术领域
本发明涉及一种数字化几何定标场和GCP相结合的内方位元素定标方法，属于遥感影像在轨几何检校与处理技术领域。
背景技术
虽然卫星上的相机在卫星发射前会进行严格的实验室内标定，然而由于发射过程中的震动、材料放气、在轨运行时成像条件的改变以及器件的老化等因素的影响，使得相机内方位元素定标参数发生改变，因此需要对卫星进行在轨内方位元素定标，卫星在轨内方位元素定标是光学遥感卫星实现高精度几何定位的关键环节，直接影响卫星影像的内部几何精度和波段配准精度。
由于高分辨卫星地面像元分辨率(GSD)越来越高(注：高景卫星的GSD为0.5m、WorldViewⅢ更是达到了0.3m)、以及内方位元素定标精度要求越来越高(注：通常要求CCD探元的指向角标定精度优于0.3像元)，从而对数字化几何定标场数据(包括数字正射影像数据(DOM数据)和数字高程影像数据(DEM数据))的内部几何精度要求极高(注：DOM数据内部几何精度优于0.1m、DEM数据高程精度优于0.5m)，然而，目前已有的数字化几何检校场数据的内部几何精度为0.3m左右，采用传统的内方位元素定标方法根本无法满足要求，如果采用更好精度的数字化几何定标场数据将耗费巨大的经费，因此，必须针对GSD优于0.5米的高分辨率卫星的高精度内方位元素定标的特殊要求，采用适当的方法进行高精度内方位元素定标。
发明内容
本发明要解决的技术问题是：克服现有技术的不足，提供了一种数字化几何定标场和GCP相结合的内方位元素定标方法，包括高精度外业地面控制点布设、采集及人工内业立体环境刺点方案、利用数字化几何定标场数据采用密集匹配的方法进行控制点的自动量测、探元指向角相机内定校模型构建及模型参数求解策略和方法、相机内定标参数优化模型构建及模型参数求解。经实践验证表明该内定标方法可行、有效，结果稳定、可靠、内定标精度高，目前已工程化应用到商遥卫星地面数据处理系统中，运行稳定，取得了良好的工程化应用效果。
本发明目的通过以下技术方案予以实现：
一种数字化几何定标场和GCP相结合的内方位元素定标方法，包括如下步骤：
步骤一、在待定标影像的某一范围内进行外业地面控制点(GCP)采集并制成点之计，根据所述点之计在立体环境下生成像控点文件；
步骤二、将数字化几何定标场数据与步骤一中的待定标影像进行控制点自动量测，获取控制点坐标；
步骤三，利用卫星的姿态、轨道、行时和卫星的实验室定标参数，建立基于一元三次曲线的探元指向角在轨几何检校模型；
步骤四、利用所述步骤二中控制点坐标和所述步骤三中的探元指向角在轨几何检校模型，迭代求解相机内标定参数；
步骤五，利用所述步骤一中生成的像控点文件，对所述步骤四中相机内定标参数进行优化。
上述数字化几何定标场和GCP相结合的内方位元素定标方法，所述步骤二中控制点的数量不少于所述步骤一中待定标影像的某一范围的像素点的万分之五。
上述数字化几何定标场和GCP相结合的内方位元素定标方法，所述步骤二中，采用灰度匹配和线特征匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点初步自动量测，然后采用最小二乘匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点精确自动量测。
上述数字化几何定标场和GCP相结合的内方位元素定标方法，所述步骤三中基于一元三次曲线的探元指向角在轨几何检校模型为：
其中
式中，(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标；λ为比例系数，分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；(t)表示当前参数是一个随时间变化的量；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表探元列号；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数。
上述数字化几何定标场和GCP相结合的内方位元素定标方法，所述步骤四中利用步骤二中控制点坐标和步骤三中的探元指向角在轨几何检校模型，迭代求解相机内标定参数的具体方法包括如下步骤：
步骤(4a)、计算沿轨指向角残差Fx(XEC,XIC)与垂轨指向角残差Fy(XEC,XIC)：
式中，XEC为外定标参数；XIC为内定标参数；a1、b1、c1、a2、b2、c2、a3、b3、c3分别代表相机安装矩阵的9个元素；矢量为物方矢量U，代表从相机投影中心到物方点的矢量在本体坐标系下的坐标；
步骤(4b)、建立每个外定标参数的误差方程：
ViEC＝AiX-Li，权为Pi，
其中
式中，ViEC为外定标参数的误差；Li是基于内外定标参数当前值的误差向量，为外定标参数的当前值；为内定标参数的当前值；i为序号；Ai是外定标参数误差方程式的系数矩阵；X代表外定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；λ为比例系数，Δpitch(t)、Δroll(t)、Δyaw(t)分别为随时间自适应的俯仰、翻滚以及偏航方向夹角变化量；
步骤(4c)、计算外定标参数改正数X：
X＝(ATPA)-1(ATPL)
其中
式中，A为第一中间矩阵，P为第二中间矩阵，L为第三中间矩阵，K为控制点的个数；
步骤(4d)、计算外定标参数XEC：
将计算所得的外定标参数XEC的值作为外定标参数的当前值
当外定标参数改正数X≤10-3时，转入步骤(4e)，否则转入步骤(4a)；
步骤(4e)、建立每个内定标参数的误差方程：
ViIC＝BiY-Li，权为Pi
其中
Y＝dXIC＝[dA0 dA1 dA2 dA3 dB0 dB1 dB2 dB3]T
式中，ViIC为内定标参数的误差；Li是基于内外定标参数当前值的误差向量；Bi是内定标参数误差方程式的系数矩阵；Y代表内定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
步骤(4f)、计算外定标参数改正数Y：
Y＝(BTPB)-1(BTPL)
其中
式中，B为第四中间矩阵，P为第二中间矩阵，L为第三中间矩阵，K为控制点的个数；
步骤(4g)、计算内定标参数XIC：
将计算所得的内定标参数XIC的值作为内定标参数的当前值
当内定标参数改正数Y≤10-3时，转入步骤(4h)，否则转入步骤(4a)；
步骤(4h)、利用步骤(4d)中计算的外定标参数XEC和步骤(4g)中计算内定标参数XIC，更新相机在轨几何检校参数文件。
上述数字化几何定标场和GCP相结合的内方位元素定标方法，所述利用步骤一中生成的像控点文件，对步骤四中相机内定标参数进行优化的方法包括如下步骤：
步骤(5a)、建立每个内定标参数优化的误差方程：
ViGCP＝BiGCPYGCP-LiGCP权为WiGCP
其中
YGCP＝dΔXIC＝[dΔA0 dΔA1 dΔA2 dΔA3 dΔB0 dΔB1 dΔB2 dΔB3]T
式中，ViGCP为内定标参数优化的误差，LiGCP是基于步骤四中内定标参数XIC和外定标参数XEC当前值的误差向量；BiGCP是误差方程式的系数矩阵；YGCP代表内定标参数优化改正数dΔXIC，d代表改正数符号；WiGCP是当前GCP的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
步骤(5b)、计算内定标参数优化改正数YGCP：
YGCP＝(BGCP TWGCPBGCP)-1(BGCP TWGCPLGCP)
其中
式中，BGCP为第五中间矩阵，WGCP为第六中间矩阵，LGCP为第七中间矩阵；
步骤(5c)、根据步骤四中计算的相机内定标参数XIC的当前值和步骤(5b)计算的内定标参数优化改正数YGCP，确定优化后的相机内定标参数XIC：
一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述方法的步骤。
本发明相比于现有技术具有如下有益效果：
一种数字化几何定标场和GCP相结合的内方位元素定标方法，包括高精度外业地面控制点布设、采集及人工内业立体环境刺点方案、利用数字化几何定标场数据采用密集匹配的方法进行控制点的自动量测、探元指向角相机内定校模型构建及模型参数求解策略和方法、相机内定标参数优化模型构建及模型参数求解；该方法首先，采用一元三次曲线探元指向角模型作为相机内定校模型，以及利用数字化几何定标场数据采用密集匹配的方法进行控制点的自动量测成果，然后，采用分步、分级迭代解求策略,解算相机内定标参数，最后，利用在待定标影像上某1000～2000行范围内相对均匀分布的外业地面控制点(GCP)采集和人工内业立体环境刺点成果，采用最小二乘平差策略对相机内定标参数进行优化，从而进一步提高相机内定标精度。经中国首颗0.5米分辨率商业遥感卫星及星座——高景01/02/03/04星在轨内定标的实践表明该内定标方法可行、有效，结果稳定、可靠、内定标精度高，目前已工程化应用到商遥卫星地面数据处理系统中，运行稳定，取得了良好的工程化应用效果。
附图说明
图1为本发明的步骤流程图；
图2为本发明实施例的影响控制点的示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明的实施方式作进一步详细描述。
一种数字化几何定标场和GCP相结合的内方位元素定标方法，如图1所示，包括如下步骤：
步骤101、在待定标影像的某一范围内进行外业地面控制点(GCP)采集并制成点之计，根据所述点之计在立体环境下生成像控点文件。
步骤102、将数字化几何定标场数据与步骤101中的待定标影像进行控制点自动量测，获取控制点坐标；控制点的数量不少于所述步骤一中待定标影像的某一范围的像素点的万分之五。采用灰度匹配和线特征匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点初步自动量测，然后采用最小二乘匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点精确自动量测。
步骤103，利用卫星的姿态、轨道、行时和卫星的实验室定标参数，建立基于一元三次曲线的探元指向角在轨几何检校模型。
基于一元三次曲线的探元指向角在轨几何检校模型为：
其中
式中，(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标；λ为比例系数，分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；(t)表示当前参数是一个随时间变化的量；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表探元列号；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数。
步骤104、利用所述步骤102中控制点坐标和所述步骤103中的探元指向角在轨几何检校模型，迭代求解相机内标定参数。其具体步骤如下：
步骤(4a)、计算沿轨指向角残差Fx(XEC,XIC)与垂轨指向角残差Fy(XEC,XIC)：
式中，XEC为外定标参数；XIC为内定标参数；a1、b1、c1、a2、b2、c2、a3、b3、c3分别代表相机安装矩阵的9个元素；矢量为物方矢量U，代表从相机投影中心到物方点的矢量在本体坐标系下的坐标；
步骤(4b)、建立每个外定标参数的误差方程：
ViEC＝AiX-Li，权为Pi，
其中
式中，ViEC为外定标参数的误差；Li是基于内外定标参数当前值的误差向量，为外定标参数的当前值；为内定标参数的当前值；i为序号；Ai是外定标参数误差方程式的系数矩阵；X代表外定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；λ为比例系数，Δpitch(t)、Δroll(t)、Δyaw(t)分别为随时间自适应的俯仰、翻滚以及偏航方向夹角变化量；
步骤(4c)、计算外定标参数改正数X：
X＝(ATPA)-1(ATPL)
其中
式中，A为第一中间矩阵，P为第二中间矩阵，L为第三中间矩阵，K为控制点的个数；
步骤(4d)、计算外定标参数XEC：
将计算所得的外定标参数XEC的值作为外定标参数的当前值
当外定标参数改正数X≤10-3时，转入步骤(4e)，否则转入步骤(4a)；
步骤(4e)、建立每个内定标参数的误差方程：
ViIC＝BiY-Li，权为Pi
其中
Y＝dXIC＝[dA0 dA1 dA2 dA3 dB0 dB1 dB2 dB3]T
式中，ViIC为内定标参数的误差；Li是基于内外定标参数当前值的误差向量；Bi是内定标参数误差方程式的系数矩阵；Y代表内定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
步骤(4f)、计算外定标参数改正数Y：
Y＝(BTPB)-1(BTPL)
其中
式中，B为第四中间矩阵，P为第二中间矩阵，L为第三中间矩阵，K为控制点的个数；
步骤(4g)、计算内定标参数XIC：
将计算所得的内定标参数XIC的值作为内定标参数的当前值
当内定标参数改正数Y≤10-3时，转入步骤(4h)，否则转入步骤(4a)；
步骤(4h)、利用步骤(4d)中计算的外定标参数XEC和步骤(4g)中计算内定标参数XIC，更新相机在轨几何检校参数文件。
步骤105，利用所述步骤101中生成的像控点文件，对所述步骤104中相机内定标参数进行优化。其具体步骤如下：
步骤(5a)、建立每个内定标参数优化的误差方程：
ViGCP＝BiGCPYGCP-LiGCP权为WiGCP
其中
YGCP＝dΔXIC＝[dΔA0 dΔA1 dΔA2 dΔA3 dΔB0 dΔB1 dΔB2 dΔB3]T
式中，ViGCP为内定标参数优化的误差，LiGCP是基于步骤四中内定标参数XIC和外定标参数XEC当前值的误差向量；BiGCP是误差方程式的系数矩阵；YGCP代表内定标参数优化改正数dΔXIC，d代表改正数符号；WiGCP是当前GCP的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
步骤(5b)、计算内定标参数优化改正数YGCP：
YGCP＝(BGCP TWGCPBGCP)-1(BGCP TWGCPLGCP)
其中
式中，BGCP为第五中间矩阵，WGCP为第六中间矩阵，LGCP为第七中间矩阵；
步骤(5c)、根据步骤四中计算的相机内定标参数XIC的当前值和步骤(5b)计算的内定标参数优化改正数YGCP，确定优化后的相机内定标参数XIC：
实施例：
一种数字化几何定标场和GCP相结合的内方位元素定标方法，包括如下步骤：
步骤1、在待定标影像上某1000～2000行范围内，进行外业控制点(GCP)采集并制成点之计，人工内业立体环境刺点，从而获取一定数量、相对均匀分布的像控点文件，通常位于影像中部左右、控制点数量约20个左右、影像两端需要有控制点，其分布如图2所示。
步骤2、利用数字化几何定标场数据(包括DOM数据和DSM数据)，采用密集匹配的方法进行控制点的自动量测，从而获取大量、密集的控制点坐标：首先，在最上层金字塔进行特征点提取(如：SUSAN算子、Harris算子、Forstner算子等)，其次，采用灰度相似性测度使匹配精度达到整像素，然后，利用匹配结果进行几何模型构建方法，消除视差，直到所有金字塔均匹配完成，最后，采用最小二乘匹配，使匹配精度达到亚像素级。
步骤3、利用卫星下传的姿态、轨道、行时等辅助数据以及实验室定标参数，构建基于一元三次曲线的探元指向角在轨几何检校模型：
式中，(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标；λ为比例系数，分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；(t)表示当前参数是一个随时间变化的量；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表探元列号；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数。
在以上几何定标模型中，待定标参数分为外定标参数XEC和内定标参数XIC，外定标参数
步骤4、利用步骤2自动匹配的密集控制点坐标和步骤3的探元指向角在轨几何检校模型，采用分步、分级迭代解求策略，求解内定标参数；实现方式为，首先基于最小二乘平差解算外定标参数，恢复相机坐标系在空间中的姿态；然后在此基础上，基于最小二乘平差解算相机内定标参数，确定相机CCD各探元在相机坐标系下的指向角。
步骤4.1、利用步骤2自动匹配的密集控制点信息，设在待定标影像上自动量测了N个均匀分布的地面控制点，各地面控制点对应的物方点和像方点分别记为GCPGi和GCPMi，物方点GCPGi的WGS84地心直角坐标为(Xi,Yi,Zi)，像方点GCPMi的影像坐标为(si,li)；
步骤4.2、令式(1)中：
式(1)转化为式(2)：
上式中，矢量为物方矢量U，代表从相机投影中心到物方点的矢量在本体坐标系下的坐标；a1,b1,c1a2,b2,c2a3,b3,c3分别代表相机安装矩阵的9个元素；Fx(XEC,XIC)、Fy(XEC,XIC)分别为沿轨指向角残差与垂轨指向角残差；
步骤4.3、对外定标参数XEC和内定标参数XIC赋初值
步骤4.4、将内定标参数XIC的当前值视为真值，将外定标参数XEC视为待求的未知参数，将内定标参数XIC和外定标参数XEC的当前值代入式(2)，对每个定向点，对式(2)进行线性化处理，建立误差方程(3)：
ViEC＝AiX-Li，权为Pi (3)
其中
ViEC为外定标参数的误差；Li是基于内外定标参数当前值的误差向量，为外定标参数的当前值；为内定标参数的当前值；i为序号；Ai是外定标参数误差方程式的系数矩阵；X代表外定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；λ为比例系数，Δpitch(t)、Δroll(t)、Δyaw(t)分别为随时间自适应的俯仰、翻滚以及偏航方向夹角变化量；
按式(4)计算法方程系数矩阵，
上式中，矩阵矩阵矩阵
利用最小二乘平差计算X，如式(5)，
X＝(ATPA)-1(ATPL) (5)
利用式(6)更新外定标参数XE的当前值，当外定标参数改正数X≤10-3时，转入步骤4.5，否则转入步骤4.1；
需要说明的是，除了采用外定标参数改正数X≤10-3作为判定依据，也可以采用外定标参数的误差ViEC作为判定依据。即当外定标参数的误差ViEC不大于预设值时，转入步骤4.5，否则转入步骤4.1。
步骤4.5、解算内定标参数，将步骤4.4所得外定标参数XEC的当前值视为真值，而内定标参数XI则视为待求的未知参数，将内定标参数XIC和外定标参数XEC的当前值代入公式(2)，对每个定向点，对式(2)进行线性化处理，建立误差方程(7)，
ViIC＝BiY-Li权为Pi (7)
其中，
Y＝dXIC＝[dA0 dA1 dA2 dA3 dB0 dB1 dB2 dB3]T
式中，ViIC为内定标参数的误差；Li是基于内外定标参数当前值的误差向量；Bi是内定标参数误差方程式的系数矩阵；Y代表内定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
按式(8)计算法方程系数矩阵；
上式中
利用最小二乘平差计算Y，如式(9)；
Y＝(BTPB)-1(BTPL) (9)
利用式(10)更新内定标参数XI的当前值，当内定标参数改正数Y≤10-3时，转入步骤4.6，否则转入步骤4.1；
需要说明的是，除了采用内定标参数改正数Y≤10-3作为判定依据，也可以采用内定标参数的误差ViIC作为判定依据。即当外定标参数的误差ViEC不大于预设值时，转入步骤4.6，否则转入步骤4.1。
步骤4.6、根据步骤4.4和步骤4.5所得外定标参数的当前值XEC和内定标参数XIC的当前值，更新相机在轨几何检校参数文件。
步骤5、利用步骤1获取的高精度GCP数据，对步骤4所得的相机内定标参数进行优化。
将步骤4所得的相机内定标参数XIC和外定标参数XEC的当前值(XEC,XIC)为真值，对于步骤1获取的每个高精度GCP，将式(2)进行线性化处理，建立误差方程(11)，
ViGCP＝BiGCPYGCP-LiGCP权为WiGCP (11)
其中
YGCP＝dΔXIC＝[dΔA0 dΔA1 dΔA2 dΔA3 dΔB0 dΔB1 dΔB2 dΔB3]T
式中，ViGCP为内定标参数优化的误差，LiGCP是基于步骤4中内定标参数XIC和外定标参数XEC当前值的误差向量；BiGCP是误差方程式的系数矩阵；YGCP代表内定标参数优化改正数dΔXIC，d代表改正数符号；WiGCP是当前GCP的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
按式(12)计算法方程系数矩阵；
上式中
利用最小二乘平差计算YGCP，如式(13)；
YGCP＝(BGCP TWGCPBGCP)-1(BGCP TWGCPLGCP) (13)
根据步骤4计算的相机内定标参数XIC和步骤5所得的相机内定标参数优化参数，更新相机在轨几何检校参数文件，如式(14)：
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。
Claims (10)
1.一种数字化几何定标场和GCP相结合的内方位元素定标方法，其特征在于：包括如下步骤：
步骤一、在待定标影像的某一范围内进行外业地面控制点(GCP)采集并制成点之计，根据所述点之计在立体环境下生成像控点文件；
步骤二、将数字化几何定标场数据与步骤一中的待定标影像进行控制点自动量测，获取控制点坐标；
步骤三，利用卫星的姿态、轨道、行时和卫星的实验室定标参数，建立基于一元三次曲线的探元指向角在轨几何检校模型；
步骤四、利用所述步骤二中控制点坐标和所述步骤三中的探元指向角在轨几何检校模型，迭代求解相机内标定参数；
步骤五，利用所述步骤一中生成的像控点文件，对所述步骤四中相机内定标参数进行优化。
2.根据权利要求1所述的一种数字化几何定标场和GCP相结合的内方位元素定标方法，其特征在于：所述步骤二中控制点的数量不少于所述步骤一中待定标影像的某一范围的像素点的万分之五。
3.根据权利要求1所述的一种数字化几何定标场和GCP相结合的内方位元素定标方法，其特征在于：所述步骤二中，采用灰度匹配和线特征匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点初步自动量测，然后采用最小二乘匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点精确自动量测。
4.根据权利要求1所述的一种数字化几何定标场和GCP相结合的内方位元素定标方法，其特征在于：所述步骤三中基于一元三次曲线的探元指向角在轨几何检校模型为：
其中
式中，(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标；λ为比例系数，分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；(t)表示当前参数是一个随时间变化的量；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表探元列号；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数。
5.根据权利要求4所述的一种数字化几何定标场和GCP相结合的内方位元素定标方法，其特征在于：所述步骤四中利用步骤二中控制点坐标和步骤三中的探元指向角在轨几何检校模型，迭代求解相机内标定参数的具体方法包括如下步骤：
步骤(4a)、计算沿轨指向角残差Fx(XEC,XIC)与垂轨指向角残差Fy(XEC,XIC)：
式中，XEC为外定标参数；XIC为内定标参数；a1、b1、c1、a2、b2、c2、a3、b3、c3分别代表相机安装矩阵的9个元素；矢量为物方矢量U，代表从相机投影中心到物方点的矢量在本体坐标系下的坐标；
步骤(4b)、建立每个外定标参数的误差方程：
ViEC＝AiX-Li，权为Pi，
其中
式中，ViEC为外定标参数的误差；Li是基于内外定标参数当前值的误差向量，为外定标参数的当前值；为内定标参数的当前值；i为序号；Ai是外定标参数误差方程式的系数矩阵；X代表外定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；λ为比例系数，Δpitch(t)、Δroll(t)、Δyaw(t)分别为随时间自适应的俯仰、翻滚以及偏航方向夹角变化量；
步骤(4c)、计算外定标参数改正数X：
X＝(ATPA)-1(ATPL)
其中
式中，A为第一中间矩阵，P为第二中间矩阵，L为第三中间矩阵，K为控制点的个数；
步骤(4d)、计算外定标参数XEC：
将计算所得的外定标参数XEC的值作为外定标参数的当前值
当外定标参数改正数X≤10-3时，转入步骤(4e)，否则转入步骤(4a)；
步骤(4e)、建立每个内定标参数的误差方程：
ViIC＝BiY-Li，权为Pi
其中
Y＝dXIC＝[dA0 dA1 dA2 dA3 dB0 dB1 dB2 dB3]T
式中，ViIC为内定标参数的误差；Li是基于内外定标参数当前值的误差向量；Bi是内定标参数误差方程式的系数矩阵；Y代表内定标参数改正数，d代表改正数符号；Pi是当前定向点RPi的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
步骤(4f)、计算外定标参数改正数Y：
Y＝(BTPB)-1(BTPL)
其中
式中，B为第四中间矩阵，P为第二中间矩阵，L为第三中间矩阵，K为控制点的个数；
步骤(4g)、计算内定标参数XIC：
将计算所得的内定标参数XIC的值作为内定标参数的当前值
当内定标参数改正数Y≤10-3时，转入步骤(4h)，否则转入步骤(4a)；
步骤(4h)、利用步骤(4d)中计算的外定标参数XEC和步骤(4g)中计算内定标参数XIC，更新相机在轨几何检校参数文件。
6.根据权利要求4所述的一种数字化几何定标场和GCP相结合的内方位元素定标方法，其特征在于：所述利用步骤一中生成的像控点文件，对步骤四中相机内定标参数进行优化的方法包括如下步骤：
步骤(5a)、建立每个内定标参数优化的误差方程：
ViGCP＝BiGCPYGCP-LiGCP 权为WiGCP
其中
YGCP＝dΔXIC＝[dΔA0 dΔA1 dΔA2 dΔA3 dΔB0 dΔB1 dΔB2 dΔB3]T
式中，ViGCP为内定标参数优化的误差，LiGCP是基于步骤四中内定标参数XIC和外定标参数XEC当前值的误差向量；BiGCP是误差方程式的系数矩阵；YGCP代表内定标参数优化改正数dΔXIC，d代表改正数符号；WiGCP是当前GCP的像点量测精度对应的权；Fxi和Fyi分别为沿轨指向角残差Fx(XEC,XIC)、垂轨指向角残差Fy(XEC,XIC)的函数模型；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数；
步骤(5b)、计算内定标参数优化改正数YGCP：
YGCP＝(BGCP TWGCPBGCP)-1(BGCP TWGCPLGCP)
其中
式中，BGCP为第五中间矩阵，WGCP为第六中间矩阵，LGCP为第七中间矩阵；
步骤(5c)、根据步骤四中计算的相机内定标参数XIC的当前值和步骤(5b)计算的内定标参数优化改正数YGCP，确定优化后的相机内定标参数XIC：
7.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如下步骤：
步骤一、在待定标影像的某一范围内进行外业地面控制点(GCP)采集并制成点之计，根据所述点之计在立体环境下生成像控点文件；
步骤二、将数字化几何定标场数据与步骤一中的待定标影像进行控制点自动量测，获取控制点坐标；
步骤三，利用卫星的姿态、轨道、行时和卫星的实验室定标参数，建立基于一元三次曲线的探元指向角在轨几何检校模型；
步骤四、利用所述步骤二中控制点坐标和所述步骤三中的探元指向角在轨几何检校模型，迭代求解相机内标定参数；
步骤五，利用所述步骤一中生成的像控点文件，对所述步骤四中相机内定标参数进行优化。
8.根据权利要求7所述的一种计算机可读存储介质，其特征在于：所述步骤二中控制点的数量不少于所述步骤一中待定标影像的某一范围的像素点的万分之五。
9.根据权利要求7所述的一种计算机可读存储介质，其特征在于：所述步骤二中，采用灰度匹配和线特征匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点初步自动量测，然后采用最小二乘匹配方法对数字化几何定标场数据与步骤一中的影像进行控制点精确自动量测。
10.根据权利要求7所述的一种计算机可读存储介质，其特征在于：所述步骤三中基于一元三次曲线的探元指向角在轨几何检校模型为：
其中
式中，(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标；λ为比例系数，分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；(t)表示当前参数是一个随时间变化的量；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表探元列号；pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；A0、A1、A2、A3、B0、B1、B2、B3均一元三次曲线探元指向角内定标模型的系数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811270813.7A CN109143295B (zh) | 2018-10-29 | 2018-10-29 | 一种数字化几何定标场和gcp相结合的内方位元素定标方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811270813.7A CN109143295B (zh) | 2018-10-29 | 2018-10-29 | 一种数字化几何定标场和gcp相结合的内方位元素定标方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109143295A true CN109143295A (zh) | 2019-01-04 |
| CN109143295B CN109143295B (zh) | 2021-04-02 |
Family
ID=64806383
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811270813.7A Active CN109143295B (zh) | 2018-10-29 | 2018-10-29 | 一种数字化几何定标场和gcp相结合的内方位元素定标方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109143295B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114241064A (zh) * | 2022-02-24 | 2022-03-25 | 中国科学院空天信息创新研究院 | 一种遥感卫星内外方位元素实时几何定标方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102901519A (zh) * | 2012-11-02 | 2013-01-30 | 武汉大学 | 一种基于探元指向角光学推扫卫星在轨分步几何定标方法 |
| CN103673995A (zh) * | 2013-11-29 | 2014-03-26 | 航天恒星科技有限公司 | 一种线阵推扫式相机在轨光学畸变参数标定方法 |
| CN106873004A (zh) * | 2016-12-21 | 2017-06-20 | 中国资源卫星应用中心 | 基于太阳高度角自适应的高轨面阵相机在轨几何检校方法 |
| US20170307763A1 (en) * | 2016-04-26 | 2017-10-26 | Uber Technologies, Inc. | Road registration differential gps |
| US20180024562A1 (en) * | 2016-07-21 | 2018-01-25 | Mobileye Vision Technologies Ltd. | Localizing vehicle navigation using lane measurements |
-
2018
- 2018-10-29 CN CN201811270813.7A patent/CN109143295B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102901519A (zh) * | 2012-11-02 | 2013-01-30 | 武汉大学 | 一种基于探元指向角光学推扫卫星在轨分步几何定标方法 |
| CN103673995A (zh) * | 2013-11-29 | 2014-03-26 | 航天恒星科技有限公司 | 一种线阵推扫式相机在轨光学畸变参数标定方法 |
| US20170307763A1 (en) * | 2016-04-26 | 2017-10-26 | Uber Technologies, Inc. | Road registration differential gps |
| US20180024562A1 (en) * | 2016-07-21 | 2018-01-25 | Mobileye Vision Technologies Ltd. | Localizing vehicle navigation using lane measurements |
| CN106873004A (zh) * | 2016-12-21 | 2017-06-20 | 中国资源卫星应用中心 | 基于太阳高度角自适应的高轨面阵相机在轨几何检校方法 |
Non-Patent Citations (2)
| Title |
|---|
| YONGJUN ZHANG ET AL.: "On-Orbit Geometric Calibration of ZY-3 Three-Line Array Imagery With Multistrip Data Sets", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 * |
| 孟伟灿等: "线阵推扫式相机高精度在轨几何标定", 《武汉大学学报·信息科学版》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114241064A (zh) * | 2022-02-24 | 2022-03-25 | 中国科学院空天信息创新研究院 | 一种遥感卫星内外方位元素实时几何定标方法 |
| CN114241064B (zh) * | 2022-02-24 | 2022-05-17 | 中国科学院空天信息创新研究院 | 一种遥感卫星内外方位元素实时几何定标方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109143295B (zh) | 2021-04-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109903352B (zh) | 一种卫星遥感影像大区域无缝正射影像制作方法 | |
| KR101965965B1 (ko) | 위성영상과 제공 rpc로부터 제작된 수치표고모델의 자동 오차보정 방법 | |
| CN107644435B (zh) | 顾及姿态校正的敏捷光学卫星无场地几何标定方法及系统 | |
| CN104897175B (zh) | 多相机光学推扫卫星在轨几何定标方法及系统 | |
| CN112017224B (zh) | Sar数据区域网平差处理方法和系统 | |
| CN107314763B (zh) | 一种基于制约函数非线性估计的卫星影像区域网平差方法 | |
| CN108830889B (zh) | 基于全局几何约束的遥感影像与基准影像的匹配方法 | |
| CN102901519B (zh) | 一种基于探元指向角光学推扫卫星在轨分步几何定标方法 | |
| CN108226982B (zh) | 单线阵卫星激光联合高精度定位处理方法 | |
| CN104807477B (zh) | 一种基于靶标控制点的卫星ccd阵列影像几何检校方法 | |
| CN111044037B (zh) | 一种光学卫星影像的几何定位方法及装置 | |
| CN107564057B (zh) | 顾及大气折光校正的高轨面阵光学卫星在轨几何标定方法 | |
| CN106873004B (zh) | 基于太阳高度角自适应的高轨面阵相机在轨几何检校方法 | |
| CN113538595B (zh) | 利用激光测高数据辅助提升遥感立体影像几何精度的方法 | |
| CN109506627B (zh) | 一种基于rpc模型的卫星影像系统误差修正方法 | |
| CN111486864A (zh) | 基于立体正八边结构的多源传感器联合标定方法 | |
| CN111508028A (zh) | 光学立体测绘卫星相机的自主在轨几何定标方法及系统 | |
| CN109188483B (zh) | 一种时序化高精度外方位元素自动定标方法 | |
| CN105571598B (zh) | 一种卫星激光高度计足印相机姿态的测定方法 | |
| CN114838740A (zh) | 一种考虑不同经纬度区域的卫星图像几何定标方法 | |
| CN109143295B (zh) | 一种数字化几何定标场和gcp相结合的内方位元素定标方法 | |
| CN110992429A (zh) | 一种单一场景大视场相机外参数标定方法 | |
| CN111611525B (zh) | 基于物方匹配高程偏差迭代修正的遥感数据高程解算方法 | |
| CN105510901A (zh) | 基于多定标场的光学卫星影像时变误差定标方法及系统 | |
| CN111275773A (zh) | 无场几何定标方法及系统 |
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 |