CN117830175B

CN117830175B - 一种任意取向条件下图像几何扭曲自定标方法

Info

Publication number: CN117830175B
Application number: CN202410244558.8A
Authority: CN
Inventors: 彭青玉; 陆星; 方雪清; 王仕涛
Original assignee: Jinan University
Current assignee: Jinan University
Priority date: 2024-03-05
Filing date: 2024-03-05
Publication date: 2024-06-07
Anticipated expiration: 2044-03-05
Also published as: CN117830175A

Abstract

本发明公开了一种任意取向条件下图像几何扭曲自定标方法，涉及天体测量技术领域，包括：基于抖动观测获取的图像得到星象像素坐标(x ₀ ,y ₀)；获取参考星表中的天球坐标()，并得到()和(x ₀ ,y ₀)的匹配列表；通过心射投影，得到天球坐标对应的标准坐标理论值；将星象的像素坐标减去当前扭曲，得到观测位置像素坐标，结合标准坐标理论值，拟合出四常数模型的参数，计算出每个星象对应的标准坐标观测值，进一步得出每个星象在标准坐标下的残差；多图像联合，求解当前的扭曲增量；若扭曲增量大于阈值则根据扭曲增量更新当前扭曲，进行迭代，直至小于阈值输出扭曲模型。本发明可求解任意取向条件下的高精度几何扭曲，提高天体测量精度。

Description

一种任意取向条件下图像几何扭曲自定标方法

技术领域

本发明涉及天体测量技术领域，特别是涉及一种任意取向条件下图像几何扭曲自定标方法。

背景技术

光学系统不可避免地存在几何扭曲（GD），对于天体测量观测获得的图像，几何扭曲主要取决于光学系统的结构和成像设备，如哈勃望远镜，在视场的角落几何扭曲可达约5个像素。对于高精度天体测量，求解扭曲模型十分必要。尽管Gaia卫星的极限星等达到了21Gmag，随着深空探测的不断推进，依然有很多地面大口径望远镜以及空间望远镜可以观测到比Gaia星表更暗的目标。如Pan-STARRS的观测图像上存在大量暗弱恒星无法从Gaia星表中的到匹配。可以预料我国空间站望远镜项目将探测到大量暗弱恒星，建立一个不依赖星表精度的扭曲求解方法，有助于提高天体测量精度。

同时，由于不同的科学任务，望远镜的取向（即图像x轴与天球赤经的夹角）并非总是保持与天球一致，开发一种对于任意取向均可适用的高精度的求解扭曲的方法十分重要。

发明内容

本发明的目的是提供一种任意取向条件下图像几何扭曲自定标方法，可以求解任意取向条件下的高精度几何扭曲，从而提高天体测量精度。

为实现上述目的，本发明提供了如下方案：

本发明提供了一种任意取向条件下图像几何扭曲自定标方法，该方法包括以下步骤：

S1，对抖动观测获取的图像进行预处理，并获取图像中的星象的原始像素位置(x ₀ , y ₀)；

S2，下载星象对应天区的星表，并将星表中参考星的天球坐标()与图像中的星象的原始像素位置(x ₀ , y ₀)匹配，获得(/>)和(x ₀ , y ₀)一一对应的匹配列表；

S3，基于匹配列表，搜索图像中最靠近中心位置的星象将其对应的天球坐标作为切点坐标（A, D），并通过切点坐标（A, D）进行心射投影换算出所有星象对应的标准坐标理论值（），初始化当前扭曲GD(x ₀ , y ₀)=(0,0)；

S4，将星象的原始像素位置减去当前扭曲GD(x ₀ , y ₀)，得到观测位置坐标(x,y) ，即(x,y) = (x ₀ , y ₀) - GD(x ₀ , y ₀)，结合标准坐标理论值（），通过公式（2）拟合，得出a、b、c、d四个常数；

式中，，/>是像素坐标与标准坐标的比例系数，/>是像素坐标与标准坐标的取向角，c和d分别是像素坐标与标准坐标的平移因子；通常情况下，地基望远镜可以通过调整CCD，使得取向角/>，而空间望远镜则可以在空间自由翻滚，取向角不固定，而本发明则对任意取向角均适用；

S5，通过a、b、c、d四常数，计算出每个星象对应的标准坐标观测值()，并与投影得到的标准坐标理论值（/>）相减，得出每个星象的残差(/>)；

S6，通过步骤S4获取的a、b、c、d和步骤S5获得的()以及切点坐标中的D值，进行多图像联合求解得到扭曲增量（/>）；

S7，图像划分，将所有图像中观测到的星象置于同一像素坐标下，将图像的视场划分成若干个等面积的子区域，根据星象的原始像素位置(x ₀ , y ₀)，计算落在每个子区域的星象对应的扭曲增量的平均值作为该子区域中心位置的几何扭曲增量；

S8，通过每个星象的原始像素位置(x ₀ , y ₀)对S7的几何扭曲增量进行双线性插值，并更新星象对应的扭曲，判断步骤S7得到的每个子区域中心位置的几何扭曲增量是否小于等于阈值，若小于等于阈值则执行步骤S9，若大于阈值，则将更新得到的星象对应的扭曲作为当前扭曲，执行步骤S4；

S9，输出扭曲模型。将图像进行S7相同的划分，将落在对应区域的星象对应的几何扭曲均值作为该区域中心的扭曲进行输出，即为扭曲模型。

进一步地，所述步骤S1中，预处理，具体包括：边缘裁剪、平场以及背景扣除，其中，背景扣除包括一维中值扣除和二维中值滤波扣除。

进一步地，所述步骤S1中，获取图像中的星象的原始像素位置(x ₀ , y ₀)，具体包括：

对抖动观测获取的图像进行点源检测，通过定心算法获取图像中星象的原始像素位置(x ₀ , y ₀)。

进一步地，所述步骤S3中，通过切点坐标（A, D）进行心射投影换算出所有星象对应的标准坐标理论值（），具体包括：

通过公式（1）将参考星的天球坐标()转换得到标准坐标理论值（/>）：

进一步地，所述步骤S5中，通过a、b、c、d四常数，计算出对于每个星象对应的标准坐标观测值()，并与投影得到的标准坐标理论值（/>）相减，得出每个星象的残差(/>)，具体包括：

；

进一步地，所述步骤S6中，通过步骤S4获取的a、b、c、d和步骤S5获得的()以及切点坐标中的D值，进行多图像联合求解得到扭曲增量（/>），具体包括：

通过公式（4）求解，

式中，i，j表示图像序号，N为同一参考星出现的图像帧数，、，/>，/>分别为第i和j帧图像的切点赤纬坐标，/>，/>分别为第i和j帧图像的常数，/>分别表示同一参考星在第i和j帧中换算到标准坐标下的残差。

进一步地，所述步骤S7中，将视场划分成若干个等面积的子区域，并将落在每个子区域内的所有星象对应的扭曲取平均作为该子区域中心的几何扭曲增量。

进一步地，所述步骤S8中，通过每个星象的原始像素位置(x ₀ , y ₀)对S7的几何扭曲增量进行双线性插值，并更新星象对应的扭曲，具体包括：

通过每个星象的原始像素位置(x ₀ , y ₀)对S7的几何扭曲增量进行双线性插值，得到每个星象的对应扭曲增量；

更新星象对应的扭曲：

GD’(x ₀ , y ₀)=GD(x ₀ , y ₀)+对应扭曲增量（5）

其中，GD’(x ₀ , y ₀)表示更新得到的星象对应的扭曲，GD(x ₀ , y ₀)表示更新前的扭曲。

进一步地，所述步骤S9，输出扭曲模型，具体包括：

S9输出的扭曲模型为视场划分后每个子区域中心坐标对应的几何扭曲，扭曲值为落在对应区域的所有星象对应的扭曲的平均值。具体地，对图像按需分割成若干个子区域，并对子区域中对应的所有星象扭曲取均值作为该区域中心的扭曲值，输出多个子区域中心的扭曲，每个区域的扭曲分为x方向和y方向，因此共计2*区域数个数值，即为扭曲模型。对该模型进行插值方式即可得到图像任意位置的扭曲，为了更方便使用扭曲模型。

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供的任意取向条件下图像几何扭曲自定标方法，可以求解任意取向条件下的高精度几何扭曲，从而提高天体测量精度。该方法不依赖星表精度，不限制取向角，不仅适用于地面望远镜，对空间望远镜系统同样适用，同时，使用本方法，也可以精确测量多芯片拼接CCD的芯片间距，这对于我国即将实施的中国空间站望远镜项目十分重要。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的任意取向条件下图像几何扭曲自定标方法的流程图；

图2为本发明实施例十字抖动观测模式；

图3为本发明实施例井字抖动观测模式；

图4为本发明实施例典型的CCD扭曲示意图;

图5为本发明实施例心射投影的示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提供的任意取向条件下图像几何扭曲自定标方法，是一种在天体测量领域中对于任意取向情况下均可求解几何扭曲的图像处理方法。

如图1所示，本发明提供的任意取向条件下图像几何扭曲自定标方法，可以由任意取向条件下图像几何扭曲自定标系统执行，该系统集成了各个功能模块，用于执行各个步骤的算法，该方法包括以下步骤：

S1，对抖动观测获取的图像进行预处理，并获取图像中的星象的原始像素位置(x ₀ , y ₀)；示例地，预处理，具体包括：边缘裁剪、平场以及背景扣除，其中，背景扣除包括一维中值扣除和二维中值滤波扣除；本领域技术人员可根据实际需要采用相关预处理操作；

获取图像中的星象的原始像素位置(x ₀ , y ₀)，具体包括：

其中，通过公式（1）将天球坐标()转换得到标准坐标理论值（/>）：

式中，A和D分别是图像中心对应的天球切点坐标的赤经和赤纬，和/>分别是参考星的赤经和赤纬，/>和/>分别为参考星通过心射投影转换到切平面标准坐标下的坐标值。其中，标准坐标就是将天球投影到图像中心对应的天球坐标为切点的切平面上的直角坐标，标准坐标也称为理想坐标，对应的平面也叫做理想平面。心射投影如图5所示，理想坐标在上面，图像坐标在下面。

式中，，/>是像素坐标与标准坐标的比例系数，/>是像素坐标与标准坐标的取向角，c和d分别是像素坐标与标准坐标的平移因子。通常情况下，通过调整CCD，使得取向角/>，而本发明则对任意取向角均适用。

S5，通通过a、b、c、d四常数，计算出每个星象对应的标准坐标观测值()，即将(x, y)带入公式（2）进行计算，并与投影得到的标准坐标理论值（/>）相减，得出每个星象在理想坐标下的残差(/>)，具体包括：

S6，通过步骤S4获取的a、b、c、d和步骤S5获得的()以及切点坐标中的D值（也就是切点的赤纬值），进行多图像联合求解得到扭曲增量（/>）；具体地，

通过公式（4）求解，

式中，i，j表示图像序号，N为同一参考星出现的图像帧数，、，/>，/>分别为第i和j帧图像的切点赤纬坐标，/>、/>分别为第i和j帧图像的常数，/>分别表示同一参考星在第i和j帧中换算到标准坐标下的残差。

S7，图像划分，将所有图像中观测到的星象置于同一像素坐标下，将视场划分成若干个等面积的子区域，根据星象原始像素位置(x ₀ , y ₀)，计算落在每个子区域的星象的几何扭曲增量的平均值作为该子区域中心位置的几何扭曲增量；该步骤的目的是：每次求解GD前，需要对整个图像进行格点划分，在每个格点内进行多颗星象扭曲的叠加以便压缩测量误差。在每颗星象应用GD时，需要从相邻格点处的扭曲进一步通过插值法获得星象位置处的具体扭曲值，因此进行步骤S8。

其中，置于同一像素坐标下，指的是将图像中的星象按照原始坐标放置到一个图像框架内。

S8，通过每个星象的原始坐标(x ₀ , y ₀)对S7的扭曲增量进行双线性插值，并更新星象对应的扭曲，判断步骤S7得到的扭曲增量是否小于等于阈值，若小于等于阈值则执行步骤S9输出扭曲模型，若大于阈值，则将更新得到的星象对应的扭曲作为当前扭曲，执行步骤S4；具体包括：

更新星象对应的扭曲：

GD’(x ₀ , y ₀)=GD(x ₀ , y ₀)+对应扭曲增量（5）

每次求解GD’前，需要对整个图像进行格点划分，在每个格点内进行多颗星象扭曲的叠加以便压缩测量误差。在每颗星象应用GD’时，需要从相邻格点处的扭曲通过插值法得出星象位置处的具体扭曲值。

其中，所述步骤S1中，取向和典型的抖动观测模式示意图（方框表示CCD图像范围，黑色十字为望远镜指向，虚线为赤经方向），除了如图2和图3所示十字和井字两种抖动观测模式，也可以是类似的抖动方式，如米字抖动。

对于抖动观测获取的图像进行预处理，包括边缘裁剪（对于部分CCD图像边缘不响应的情况下），平场校正（有平场数据时），背景扣除（背景扣除包括一维中值扣除或者二维中值滤波扣除）。预处理结果对扭曲精度有细微的影响。

获取图像中星象的坐标列表(x ₀ , y ₀)，具体包括：

对抖动观测获取的图像进行点源检测，通过定心算法获取图像中星象的坐标列表(x ₀ , y ₀)。示例地，定心算法可以采用二维高斯定心，可选的算法有修正矩、一维高斯等。

本发明中，图像扭曲定义为GD = GD(x,y)，每次计算GD都是以星象原始像素位置为准。典型的CCD扭曲如图4所示，仅与图像的x，y坐标有关，扭曲增量表示（），单位为像素。

其中，所述步骤S2中，下载星象对应天区的星表，并获取星表中天球坐标()，星象和星表中的参考星匹配，获得(/>)和(x ₀ , y ₀)一一对应的列表。

假设星象误差为()，则有

其中，(x ₀ , y ₀)为通过定心算法得到的星象中心的像素坐标。和/>是由星表通过切点坐标（A,D）进行心射投影换算而来的理论值。一般选取图像最靠近中心的星象对应的恒星的赤经和赤纬作为切点坐标即可。(x ₀ , y ₀) 、（/>）、（A,D）均为已知量，扭曲(dx, dy)初始化为0，参数a~d通过拟合获得。

式（6）通过换算可得：

简化可得

其中；

令

对于两帧有重叠星象的图像i、j，重叠星象的星表误差，/>是一样的，因此有

移项得

计算提取得

由于同一个参考星出现在多帧图像中，如N帧，对于第i帧图像的该星象中心位置，联合其他出现该参考星的图像，则有

由于参考星在N帧图像中均出现且分布每个图像的不同区域（基于抖动观测），因此、/>四项为0；从而得出前边所述公式（4）：

通过步骤S7，将所有图像中观测到的星象置于同一像素坐标下，将图像划分成若干个等面积的区域，根据星象原始像素位置(x ₀ , y ₀)，计算落在每个区域的星象的几何扭曲增量的平均值作为该区域中心位置的几何扭曲的增量。

通过步骤S8，更新扭曲并进入迭代过程，其中，如果是第一次迭代，GD从(0,0)变化为（），如果不是第一次迭代的话，GD’=GD +（/>），注意，此处的GD实际上对应的是多个星象的位置的扭曲，如果一共出现了M个星象，则M个星象所在的位置的扭曲增量（/>）的数值将会通过对S7中获得的划分后区域中心扭曲增量进行双线性插值的方式被重新计算出来，从而获得一个误差更小的值，而不是直接采用原始值。

步骤S8中，还需判断S7获取的每个划分后区域中心的扭曲是否小于阈值，小于则输出扭曲模型，大于则进行S4进行迭代。由于扭曲的变化，a、b、c、d在步骤4重新拟合，会有微小变化。

示例地，一种实施方式，所述步骤S9，输出扭曲模型，具体包括：

对图像按需分割成多个等面积的子区域，并对落在每个子区域的星象扭曲取均值作为该区域中心的扭曲值，每个子区域均有x和y两个方向扭曲值，即为扭曲模型，通过插值方式即可得到图像任意位置的扭曲。

以1900*1900像素的图像为例，一般是分成19*19个100*100像素的子区域，将所有抖动观测获取的图像中星象原始位置（x₀,y₀）落在对应子区域的星象的扭曲进行平滑（取平均），作为对应子区域中心的扭曲值，实际输出的内容为19*19个子区域中心点的扭曲，共计19*19*2个数值即为扭曲模型，在应用时，通过双线性插值方式即可得到图像任意位置的扭曲。更进一步，为了减少参数数量和运算量，可以通过对当前的扭曲模型（即19*19*2个扭曲值）进行多项式拟合的方式，将拟合后的多项式作为扭曲函数，

其中，扭曲函数定义如下：

式中，、/>分别表示x和y方向的扭曲，i、j分别表示阶数，/>、/>分别表示x和y方向扭曲拟合的参数，各15个，共计30个参数。

本发明可以求解任意取向条件下的高精度几何扭曲，且不依赖星表精度，从而提高天体测量精度。该方法不仅适用于地面望远镜，对空间望远镜系统也适用，这对于我国即将实施的中国空间站望远镜项目（CSST）十分重要。

本发明还提供一种电子装置，包括一个或多个处理器、存储器；

一个或多个应用程序，其中所述一个或多个应用程序被存储在所述存储器中并被配置为由所述一个或多个处理器执行，所述一个或多个程序配置用于执行如上述的任意取向条件下图像几何扭曲自定标方法。

本发明还提供一种非暂态计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如上述的任意取向条件下图像几何扭曲自定标方法。

当然，本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关硬件(如处理器，控制器等)来完成，所述的程序可存储于一计算机可读取的存储介质中，所述程序在执行时可包括如上述各方法实施例的流程。其中所述的存储介质可为存储器、磁碟、光盘等。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims

1.一种任意取向条件下图像几何扭曲自定标方法，其特征在于，包括以下步骤：

S1，对抖动观测获取的图像进行预处理，并获取图像中的星象的原始像素位置(x₀，y₀)；

S2，下载星象对应天区的星表，并将星表中参考星的天球坐标(α，δ)与图像中的星象的原始像素位置(x₀，y₀)匹配，获得(α，δ)和(x₀，y₀)一一对应的匹配列裹；

S3，基于匹配列表，搜索图像中最靠近中心位置的星象，并将其对应的天球坐标作为切点坐标(A，D)，并通过切点坐标(A，D)进行心射投影换算出所有星象对应的标准坐标理论值(ξ_c，η_c)，初始化当前扭曲GD(x₀，y₀)＝(0，0)；

S4，将星象的原始像素位置(x₀，y₀)减去当前扭曲GD(x₀，y₀)，得到观测位置坐标(x，y)，即(x，y)＝(x₀，y₀)-GD(x₀，y₀)，结合标准坐标理论值(ξ_c，η_c)，通过公式(2)拟合，得出a、b、c、d四个常数；

式中，ρ是像素坐标与标准坐标的比例系数，/>是像素坐标与标准坐标的取向角，c和d分别是像素坐标与标准坐标的平移因子；

S5，通过a、b、c、d四个常数，计算出每个星象对应的标准坐标观测值(ξ，η)，并与投影得到的标准坐标理论值(ξ_c，η_c)相减，得出每个星象的残差(Δξ，Δη)；

S6，通过步骤S4获取的a、b、c、d四个常数和步骤S5获得的残差(Δξ，Δη)以及切点坐标中的D值，进行多图像联合求解得到扭曲增量(dx_i，dy_i)，具体包括：

通过公式(4)求解，

式中，i，j表示图像序号，N为同一参考星出现的图像帧数， D_i、D_j分别为第i帧和第j帧图像的切点赤纬坐标，a_i、b_i为第i帧图像的常数，(Δξ_i，Δη_i)、(Δξ_j，Δη_j)分别表示同一参考星在第i帧和第j帧图像中换算到标准坐标下的残差；

S7，图像划分，将所有图像中观测到的星象置于同一像素坐标下，将图像的视场划分成若干个等面积的子区域，根据星象的原始像素位置(x₀，y₀)，计算落在每个子区域的星象对应的扭曲增量(dx_i，dy_i)的平均值，并将其作为该子区域中心位置的几何扭曲增量；

S8，通过每个星象的原始像素位置(x₀，y₀)对S7的几何扭曲增量进行双线性插值，并更新星象对应的扭曲，判断步骤S7得到的每个子区域中心位置的几何扭曲增量是否小于等于阈值，若小于等于阈值则执行步骤S9，若大于阈值，则将更新得到的星象对应的扭曲作为当前扭曲，执行步骤S4；

S9，输出扭曲模型。

2.根据权利要求1所述的任意取向条件下图像几何扭曲自定标方法，其特征在于，所述步骤S1中，所述预处理，具体包括：边缘裁剪、平场扣除以及背景扣除，其中，背景扣除包括一维中值滤波扣除和二维中值滤波扣除。

3.根据权利要求1所述的任意取向条件下图像几何扭曲自定标方法，其特征在于，所述步骤S1中，获取图像中的星象的原始像素位置(x₀，y₀)，具体包括：

对抖动观测获取的图像进行点源目标检测，通过定心算法获取图像中星象的原始像素位置(x₀，y₀)。

4.根据权利要求1所述的任意取向条件下图像几何扭曲自定标方法，其特征在于，所述步骤S3中，通过切点坐标(A，D)进行心射投影换算出所有星象对应的标准坐标理论值(ξ_c，η_c)，具体包括：

通过公式(1)将参考星的天球坐标(α，δ)转换得到标准坐标理论值(ξ_c，η_c)：

5.根据权利要求1所述的任意取向条件下图像几何扭曲自定标方法，其特征在于，所述步骤S5中，通过a、b、c、d四个常数，计算出对于每个星象对应的标准坐标观测值(ξ，η)，并与投影得到的标准坐标理论值(ξ_c，η_c)相减，得出每个星象的残差(Δξ，Δη)，具体包括：

6.根据权利要求1所述的任意取向条件下图像几何扭曲自定标方法，其特征在于，所述步骤S8中，通过每个星象的原始像素位置(x₀，y₀)对S7的几何扭曲增量进行双线性插值，并更新星象对应的扭曲，具体包括：

通过每个星象的原始像素位置(x₀，y₀)对S7的几何扭曲增量进行双线性插值，得到每个星象的对应扭曲增量；

更新星象对应的扭曲：

GD′(x₀，y₀)＝GD(x₀，y₀)+对应扭曲增量 (5)

其中，GD′(x₀，y₀)表示更新得到的星象对应的扭曲，GD(x₀，y₀)表示更新前的扭曲。

7.根据权利要求1所述的任意取向条件下图像几何扭曲自定标方法，其特征在于，所述步骤S9中，输出扭曲模型，具体包括：

输出的扭曲模型为图像的视场划分后每个子区域中心位置对应的几何扭曲，扭曲值为落在对应子区域的所有星象对应的扭曲的平均值。