针对国产面阵摆扫成像系统获取的影像生产正射影像的方法
阅读说明：本技术 针对国产面阵摆扫成像系统获取的影像生产正射影像的方法 (Method for producing orthoimage aiming at image acquired by domestic area array swinging imaging system ) 是由 段延松 赵新博 张祖勋 柯涛 于 2021-08-24 设计创作，主要内容包括：本发明涉及一种针对国产面阵摆扫成像系统获取的影像生产正射影像的方法。首先选择合理的局部坐标系作为基准坐标系,将导航数据信息转换到基准坐标系中,然后在DEM辅助下对所有影像进行同名点自动提取,形成影像连接点网,之后使用连接点开展基于DEM的区域网平差,解算每张影像的外方位参数,最后基于最优化理论生产目标区域的整体最优正射影像。本发明引入GPS/IMU参数作为初值,使用DEM数据辅助平差,解决了传统摄影测量无法直接进行平差解算的问题。本发明可以更快速、更准确地完成影像拼接,生产合格的正射影像产品,填补国产面阵摆扫成像系统数据处理的空白。(The invention relates to a method for producing an orthoimage aiming at an image acquired by a domestic area array sweep imaging system. Firstly, a reasonable local coordinate system is selected as a reference coordinate system, navigation data information is converted into the reference coordinate system, then homonymous points are automatically extracted from all images under the assistance of a DEM (digital elevation model) to form an image connection point network, then, the connection points are used for carrying out adjustment of the regional network based on the DEM, the external orientation parameter of each image is solved, and finally, the overall optimal orthoscopic image of a target region is produced based on an optimization theory. The invention introduces GPS/IMU parameters as initial values, uses DEM data to assist adjustment, and solves the problem that the adjustment can not be directly calculated in the traditional photogrammetry. The method can complete image splicing more quickly and accurately, produce qualified orthoscopic image products, and fill the blank of data processing of a domestic area array scanning imaging system.)
技术领域
本发明属于遥感测绘技术领域，特别是涉及一种针对国产面阵摆扫成像系统获取的影像生产正射影像的方法。
背景技术
面阵摆扫成像系统是一种快速获取地面影像信息的设备，与传统固定框幅成像系统和多视倾斜成像系统相比，数据获取效率高，整体成本低，具有很强的实用性。早在2008年，国际上就已经推出商用系统，典型代表是以色列的A3面阵摆扫成像系统和日本的SamplX面阵摆扫成像系统。这些系统已经形成商业化产品，对外提供了完整解决方案，方案包含了从系统航拍安装到最终测绘产品(包括空中三角测量、数字高程模型、正射影像、数字线划图)的全流程生产过程操作说明，但是在具体技术细节上，如空中三角测量处理过程和技术细节等方面没有提供相关信息。厂商在设备说明中也仅提供了如相机分辨率、相机焦距、影像格式等设备输出成果参数，对于更核心的参数如成像系统POS系统工作原理、参数等都没有提供，只提供了对系统获取的影像信息开展生产的软件。经过国内相关
技术领域
的科研院所10年的努力，成功研制出国产面阵摆扫成像系统。国产面阵摆扫成像系统由成像硬件和处理软件组成，成像硬件主要由镜头、成像CCD面板、GPS测位设备、IMU测姿设备、摆扫驱动设备、存储设备等核心部件组成，如图1所示。处理软件主要包含影像信号解码、GPS数据处理、IMU数据处理以及测绘产品生产模块。在成像原理上，面阵摆扫成像系统与传统成像系统存在较大差异，其最大特点是连续摆扫成像，摆扫周期内，任意两张影像可近似理解为同心成像，可无缝拼接，重叠区域为同光线成像，不能进行立体定位。此外，国产面阵摆扫成像系统所使用的IMU受自主技术限制，精度只能做到5～10角秒，无法实现影像精准定位定姿态。因此，面阵摆扫成像系统的处理软件需要充分考虑以上特点。此外，面阵摆扫成像系统的每张影像都具有GPS信息和IMU信息，可为处理软件提供重要的信息。
目前现有处理技术没有充分考虑摆扫影像的特点，特别是同心成像无法进行立体定位的特性，因此现有系统无法生产合格产品。本发明提出的方法针对国产面阵摆扫成像系统获取影像连续摆扫成像、影像重叠小、IMU精度不高等特点，引入DEM数据辅助处理，在处理效率和成果精度之间寻找平衡，可以更快速、更准确地完成影像拼接，生产合格的正射影像产品，有效填补国产面阵摆扫成像系统数据处理的空白。
发明内容
本发明针对现有技术的不足，提供一种针对国产面阵摆扫成像系统获取的影像生产正射影像的方法。首先选择合理的局部坐标系作为基准坐标系，将导航数据信息转换到基准坐标系中，然后在DEM辅助下对所有影像进行同名点自动提取，形成影像连接点网，之后使用连接点开展基于DEM的区域网平差，解算每张影像的外方位参数，最后基于最优化理论生产目标区域的整体最优正射影像。
为了达到上述目的，本发明提供的技术方案是一种针对国产面阵摆扫成像系统获取的影像生产正射影像的方法，包括以下步骤：
步骤1，建立局部坐标系作为统一的基准坐标系，将导航数据信息转换到基准坐标系中；
步骤2，在概略DEM辅助下进行影像同名点自动提取；
步骤2.1，在概略DEM辅助下，根据摆扫成像周期性和先后顺序获取影像间的相邻关系；
步骤2.2，根据步骤2.1获取的影像相邻关系，采用基于两类膨胀进行立体影像密集匹配的方法开展同名点匹配；
步骤3，基于DEM辅助进行区域网平差，将步骤1转换处理后的导航数据信息作为平差初始值，解算每张影像的外方位元素参数；
步骤3.1，利用所有影像的POS信息计算区域网范围，根据区域网范围从全球SRTM数据中获取对应的DEM数据；
步骤3.2，利用外方位角元素计算每张影像的旋转矩阵，使用步骤1得到的角元素作为计算的初始值；
步骤3.3，以单束同名光线为观测单元，分别计算各连接像点光线与DEM的交点，得到对应的物方坐标；
步骤3.4，计算同名像点对应的物方点间的差异大小，如果大于角元素精度推导的理论物方定位精度α，则将同名像点标记为粗差点，不参与平差；否则，取同名像点对应的物方坐标的平均值作为真值，参与平差；
步骤3.5，逐点构建误差方程，解算外方位元素；
步骤3.6，重复步骤3.2至步骤3.5直到影像的每个角元素的变化不超过10-6弧度且所有像点改正数的平方和小于设定的阈值λ；
步骤4，整体最优正射影像生产；
步骤4.1，计算每张原始影像四个角点对应的地面坐标；
步骤4.2，选择正射影像对应的最佳原始影像进行反投影；
步骤4.3，反投影法生产正射影像。
而且，所述步骤1中局部坐标系采用定点切平面坐标系统，即TPC坐标系(TangentPlane Coordinate System)，坐标定义如下：在WGS84坐标框架下，选择飞行区域中几何中心点经度L0，纬度B0作为平面原点，以地球椭球高为0作椭球切平面，定义过原点纬度线的东方向为X坐标轴，过原点经度线北方向为Y坐标轴，过原点背向地心的垂线为Z坐标轴，如图4所示。
TPC坐标系与像空间坐标系之间存在转换关系如下：
即从像空间坐标系→载荷本体坐标系→IMU本体坐标系→导航坐标系→地心地固坐标系→TPC坐标系，各坐标系定义及转换如图5所示。
各坐标系间旋转矩阵计算公式如下：
式中，为TPC坐标系下的影像姿态旋转矩阵，是影像的三个外方位角元素，为待求解参数；是地心地固坐标系与TPC坐标系之间的转换矩阵，由TPC坐标原点的经度L0，纬度B0定义；是导航坐标系与地心地固坐标系之间的转换矩阵，由载荷所在位置的经度Lz、纬度Bz定义；是导航坐标系和IMU本体坐标系之间的旋转矩阵，(Φ,Θ,Ψ)是IMU记录的俯仰、横滚、偏航数值；是载荷本体坐标系与IMU本体坐标系之间的转换矩阵，又称为设备安置矩阵，由IMU安置在载荷上的ex、ey、ez三个角定义；是像空间坐标系与载荷本体坐标系之间的转换矩阵。
求得旋转矩阵后，反解可得到影像在TPC坐标系中的姿态角，即：
ω＝arcsin(-b3) (9)
GPS点位置信息可以使用地球坐标定义模型无损地转换到任意需要的坐标系中，从原始的GPS坐标(L,B,H)转换到局部切平面坐标系中的TPC坐标经过如下变换：
首先将原始的GPS数据转换到地心坐标系中：
Xc＝(N+H)cos B cos L (11)
Yc＝(N+H)cos B sin L (12)
Zc＝[N(1-e2)+H]sin B (13)
式中，N为卯酉圈的半径，a为地球椭球的长半轴，b为地球椭球的短半轴，(Xc,Yc,Zc)为地心坐标系的坐标。
然后将局部切平面坐标系原点坐标(L0,B0,0)带入式(11)～(13)，求得局部切平面坐标系原点在地心坐标系中的坐标，记为从地心系坐标到TPC坐标系的计算公式为：
式中，是地心地固坐标系与TPC坐标系之间的转换矩阵，由式(2)计算可得；(XTPC,YTPC,ZTPC)为TPC坐标系的坐标。
即为转换到基准坐标系下的导航信息。
而且，所述步骤2中根据摆扫成像周期性和先后顺序获取影像间的相邻关系是先将摆扫一个周期得到的N张影像分为一组，每张影像根据其在摆扫周期中的次序将它们编号为1至N，顺序连续的影像必为邻接影像。各组之间的相邻关系可根据GPS和IMU信息，将获取的影像投影到全球DEM表示的地面上，然后根据每张影像的地面范围求影像重叠情况，也可根据重叠区域反算出原始影像的重叠区域，以此确认各组间影像的相邻关系。
而且，所述步骤3.3中以单束光线为观测单元，分别计算各连接点与DEM的交点，得到对应的物方坐标，计算公式如下：
式中，(Xg,Yg)为交点的物方坐标(即TCP坐标)；zdem为交点的DEM高程；Xs、Ys、Zs是影像的外方位线元素，步骤1得到的(XTPC,YTPC,ZTPC)可作为外方位线元素的初始值；(x,y)是像点坐标；f是相机焦距。
而且，所述步骤3.4中理论物方定位精度α的计算方式如下：
式中，H为摄影位置到地面的高度，ε为成像系统的IMU角元素精度，θ为主光轴与重力方向的夹角。
而且，所述步骤3.5中逐点构建误差方程为：
式中，(x,y)是像点坐标，(x)、(y)是用未知数初值获取的像点坐标，vx、vy是误差变量，表示求偏导，d代表误差改正，求解出外方位元素的改正数之后即可更新每张影像的外方位元素。
而且，所述步骤4.3中反投影法生产正射影像的方法如下：设正射影像上任一点P的坐标为(X′,Y′)，由正射影像左下角图廓点地面坐标(X0,Y0)与正射影像比例尺分母M，计算P点所对应的地面坐标(X,Y)如下式所示：
X＝X0+M×X′ (22)
Y＝Y0+M×Y′ (23)
计算正射影像上点P对应原始影像上相应像点坐标p(x,y)为：
式中，X0、Y0为所有原始影像的角点对应的地面坐标中的最小值，Z为P点的高程，由DEM内插求得；ai,bi,ci(i＝1,2,3)是影像的旋转矩阵元素，Xs、Ys、Zs是影像的线元素，f是原始影像的焦距。
由于所得的像点坐标不一定正好落在像元素中心，为此需进行灰度内插，本发明采用双线性内插的方法求得像点p的灰度值g(x,y)，将像点的灰度值赋给正射影像上对应的像素点P(X′,Y′)，即：
G(X′,Y′)＝g(x,y) (26)
依次对每个像元进行上述运算，即能获得正射影像。
与现有技术相比，本发明具有如下优点：
(1)国产面阵摆扫成像系统连续摆扫成像，利用其成像周期内的影像连续即相邻的特点，使用GPS/IMU与DEM数据求交获取影像邻接矩阵和重叠区域，针对性开展影像匹配，快速提取影像连接点。
(2)国产面阵摆扫成像系统的POS精度欠佳，用于直接地理定位精度有限，不能直接生产正射影像，仍需要平差处理；但该系统的影像之间重叠度小、交会条件弱，本发明引入GPS/IMU参数作为初值，使用DEM数据辅助平差解决了传统摄影测量无法直接进行平差解算的问题。
(3)本发明可以更快速、更准确地完成影像拼接，生产合格的正射影像产品，填补国产面阵摆扫成像系统数据处理的空白。
附图说明
图1为本发明实施例国产面阵摆扫成像系统核心部件示意图。
图2为本发明实施例技术流程图。
图3为本发明实施例IMU旋转角的国际标准定义示意图。
图4为本发明实施例中局部切平面坐标系(TPC坐标系)定义示意图。
图5为本发明实施例中IMU坐标系转换为TPC坐标系的示意图。
图6为本发明实施例国产面阵摆扫成像系统摆扫成像的GPS/IMU参数转换结果示意图。
图7为本发明实施例影像邻接矩阵示意图。
图8是本发明实施例的正射影像数据成果。
图9是本发明实施例的正射影像数据拼接部位的局部放大图。
具体实施方式
本发明提供涉及一种针对国产面阵摆扫成像系统获取的影像生产正射影像的方法，首先选择合理的局部坐标系作为基准坐标系，将导航数据信息转换到基准坐标系中，然后在DEM辅助下对所有影像进行同名点自动提取，形成影像连接点网，之后使用连接点开展基于DEM的区域网平差，解算每张影像的外方位参数，最后基于最优化理论生产目标区域的整体最优正射影像。
下面结合附图和实施例对本发明的技术方案作进一步说明。
如图2所示，本发明实施例的流程包括以下步骤：
步骤1，建立局部坐标系作为统一的基准坐标系，将导航数据信息转换到基准坐标系中。
国产面阵摆扫成像系统安装有GPS设备和IMU设备，GPS设备记录飞行器位置信息，包含经度、纬度和高程，IMU记录飞行器姿态信息，包含俯仰、横滚、偏航3个角度值，IMU坐标系按国际标准进行定义，如图3所示。
局部坐标系采用定点切平面坐标系统，即TPC坐标系(Tangent Plane CoordinateSystem)，坐标定义如下：在WGS84坐标框架下，选择飞行区域中几何中心点经度L0，纬度B0作为平面原点，以地球椭球高为0作椭球切平面，定义过原点纬度线的东方向为X坐标轴，过原点经度线北方向为Y坐标轴，过原点背向地心的垂线为Z坐标轴，如图4所示。
TPC坐标系与像空间坐标系之间存在转换关系如下：
即从像空间坐标系→载荷本体坐标系→IMU本体坐标系→导航坐标系→地心地固坐标系→TPC坐标系，各坐标系定义及转换如图5所示。
各坐标系间旋转矩阵计算公式如下：
式中，为TPC坐标系下的影像姿态旋转矩阵，是影像的三个外方位角元素，为待求解参数；是地心地固坐标系与TPC坐标系之间的转换矩阵，由TPC坐标原点的经度L0，纬度B0定义；是导航坐标系与地心地固坐标系之间的转换矩阵，由载荷所在位置的经度Lz、纬度Bz定义；是导航坐标系和IMU本体坐标系之间的旋转矩阵，(Φ,Θ,Ψ)是IMU记录的俯仰、横滚、偏航数值；是载荷本体坐标系与IMU本体坐标系之间的转换矩阵，又称为设备安置矩阵，由IMU安置在载荷上的ex、ey、ez三个角定义；是像空间坐标系与载荷本体坐标系之间的转换矩阵。
求得旋转矩阵后，反解可得到影像在TPC坐标系中的姿态角，即：
ω＝arcsin(-b3) (9)
GPS点位置信息可以使用地球坐标定义模型无损地转换到任意需要的坐标系中，从原始的GPS坐标(L,B,H)转换到局部切平面坐标系中的TPC坐标经过如下变换：
首先将原始的GPS数据转换到地心坐标系中：
Xc＝(N+H)cos B cos L (11)
Yc＝(N+H)cos B sin L (12)
Zc＝[N(1-e2)+H]sin B (13)
式中，N为卯酉圈的半径，a为地球椭球的长半轴，b为地球椭球的短半轴，(Xc,Yc,Zc)为地心坐标系的坐标。
然后将局部切平面坐标系原点坐标(L0,B0,0)带入式(11)～(13)，求得局部切平面坐标系原点在地心坐标系中的坐标，记为从地心系坐标到TPC坐标系的计算公式为：
式中，是地心地固坐标系与TPC坐标系之间的转换矩阵，由式(2)计算可得；(XTPC,YTPC,ZTPC)为TPC坐标系的坐标。
即为转换到基准坐标系下的导航信息。
步骤2，在概略DEM辅助下进行影像同名点自动提取。
针对国产面阵摆扫成像系统获取的影像，可以充分利用其周期摆扫规律和GPS/IMU参数，配合全球DEM数据(主要是SRTM数据)，实现快速有效的影像匹配。面阵摆扫成像系统的影像和其GPS和IMU参数投影到全球DEM数据处理原理和效果示意图如图6所示。
步骤2.1，在概略DEM辅助下，根据摆扫成像周期性和先后顺序获取影像间的相邻关系。
首先，将摆扫一个周期得到的N(本实施例中，N取36)张影像分为一组，每张影像根据其在摆扫周期中的次序将它们编号为1至N，顺序连续的影像必为邻接影像。各组之间的相邻关系可根据GPS和IMU信息，将获取的影像投影到全球DEM表示的地面上，然后根据每张影像的地面范围求影像重叠情况，也可根据重叠区域反算出原始影像的重叠区域，以此确认各组间影像的相邻关系，最终确认的邻接矩阵如图7所示。
步骤2.2，根据步骤2.1获取的影像相邻关系，采用基于两类膨胀进行立体影像密集匹配的方法开展同名点匹配。
步骤3，基于DEM辅助进行区域网平差，将步骤1转换处理后的导航数据信息作为平差初始值，解算每张影像的外方位元素参数。
步骤3.1，利用所有影像的POS信息计算区域网范围，根据区域网范围从全球SRTM数据中获取对应的DEM数据。
步骤3.2，利用外方位角元素计算每张影像的旋转矩阵，使用步骤1得到的角元素作为计算的初始值。
外方位角元素计算影像的旋转矩阵如下式所示：
步骤3.3，以单束同名光线为观测单元，分别计算各连接像点光线与DEM的交点，得到对应的物方坐标，计算公式如下：
式中，(Xg,Yg)为交点的物方坐标(即TCP坐标)；zdem为交点的DEM高程；Xs、Ys、Zs是影像的外方位线元素，步骤1得到的(XTPC,YTPC,ZTPC)可作为外方位线元素的初始值；(x,y)是像点坐标；f是相机焦距。
步骤3.4，计算同名像点对应的物方点间的差异大小，如果大于角元素精度推导的理论物方定位精度α，则将同名像点标记为粗差点，不参与平差；否则，取同名像点对应的物方坐标的平均值作为真值，参与平差。
设成像系统的IMU角元素精度为ε，主光轴与重力方向的夹角为θ，则理论物方定位精度α为：
式中，H为摄影位置到地面的高度。
步骤3.5，逐点构建误差方程，解算未知数改正数，误差方程为：
式中，(x,y)是像点坐标，(x)、(y)是用未知数初值获取的像点坐标，vx、vy是误差变量，表示求偏导，d代表误差改正，求解出外方位元素的改正数之后即可更新每张影像的外方位元素。
步骤3.6，重复步骤3.2至步骤3.5直到影像的每个角元素的变化不超过10-6弧度且所有像点改正数的平方和小于设定的阈值λ。
步骤4，整体最优正射影像生产。
步骤4.1，计算每张原始影像四个角点对应的地面坐标。
利用公式(18)和(19)计算每张原始影像四个角点对应的地面坐标构成四边形，并求得每个四边形的重心。
步骤4.2，选择正射影像对应的最佳原始影像进行反投影。
设正射影像上任一点P对应的地面坐标为(X,Y)，选择与点P(X,Y)距离最近的四边形重心对应的原始影像作为最佳原始影像进行反投影。
步骤4.3，反投影法生产正射影像。
设正射影像上任一点P的坐标为(X′,Y′)，由正射影像左下角图廓点地面坐标(X0,Y0)与正射影像比例尺分母M，计算P点所对应的地面坐标(X,Y)如下式所示：
X＝X0+M×X′ (22)
Y＝Y0+M×Y′ (23)
计算正射影像上点P对应原始影像上相应像点坐标p(x,y)为：
式中，X0、Y0为所有原始影像的角点对应的地面坐标中的最小值，Z为P点的高程，由DEM内插求得；ai,bi,ci(i＝1,2,3)是影像的旋转矩阵元素，Xs、Ys、Zs是影像的线元素，f是原始影像的焦距。
由于所得的像点坐标不一定正好落在像元素中心，为此需进行灰度内插，本发明采用双线性内插的方法求得像点p的灰度值g(x,y)，将像点的灰度值赋给正射影像上对应的像素点P(X′,Y′)，即：
G(X′,Y′)＝g(x,y) (26)
依次对每个像元进行上述运算，即能获得正射影像。
采用步骤4.1-步骤4.3生产整体正射影像，每个像素仅计算一次，避免了冗余计算，大幅提高了处理效率，由于每个像素都选择了最优原始影像，最终得到的整体正射影像也是最优化的。
具体实施时，以上流程可采用计算机软件技术实现自动运行流程。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。