一种数字化几何定标场和GCP相结合的内方位元素定标方法与流程

本发明涉及一种数字化几何定标场和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＝(a^tpa)^-1(a^tpl)

其中

式中，a为第一中间矩阵，p为第二中间矩阵，l为第三中间矩阵，k为控制点的个数；

步骤(4d)、计算外定标参数xec：

将计算所得的外定标参数xec的值作为外定标参数的当前值

当外定标参数改正数x≤10^-3时，转入步骤(4e)，否则转入步骤(4a)；

步骤(4e)、建立每个内定标参数的误差方程：

viic＝biy-li，权为pi

其中

y＝dxic＝[da0da1da2da3db0db1db2db3]^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＝(b^tpb)^-1(b^tpl)

其中

式中，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δa0dδa1dδa2dδa3dδb0dδb1dδb2dδ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＝(a^tpa)^-1(a^tpl)

其中

式中，a为第一中间矩阵，p为第二中间矩阵，l为第三中间矩阵，k为控制点的个数；

步骤(4d)、计算外定标参数xec：

将计算所得的外定标参数xec的值作为外定标参数的当前值

当外定标参数改正数x≤10^-3时，转入步骤(4e)，否则转入步骤(4a)；

步骤(4e)、建立每个内定标参数的误差方程：

viic＝biy-li，权为pi

其中

y＝dxic＝[da0da1da2da3db0db1db2db3]^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＝(b^tpb)^-1(b^tpl)

其中

式中，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δa0dδa1dδa2dδa3dδb0dδb1dδb2dδ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＝(a^tpa)^-1(a^tpl)(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＝[da0da1da2da3db0db1db2db3]^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＝(b^tpb)^-1(b^tpl)(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δa0dδa1dδa2dδa3dδb0dδb1dδb2dδ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)：

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。