CN104298887B - 一种多片线阵ccd相机的相对辐射定标方法 - Google Patents

Abstract

本发明涉及一种多片线阵CCD相机的相对辐射定标方法，包括以下步骤：建立遥感卫星倾斜布置的多片线阵CCD相机的焦平面坐标系与地理坐标系；根据地球半径以及卫星轨道高度、俯仰角与滚动角，计算成像时刻地理坐标系的坐标原点偏移星下点的地心角大小；建立从地理坐标系到焦平面坐标系的转换关系；在定标模式的成像过程中，依据转换关系计算卫星的偏航角速度，使得多片线阵CCD的所有像元对地面同一区域成像；利用直方图匹配方法建立直方图查找表。本发明建立在严密转换关系的基础之上，将倾斜布置的多片线阵CCD所有像元对同一区域成像，使得所有像元具有相同的输入。

Description

一种多片线阵CCD相机的相对辐射定标方法

技术领域

本发明涉及一种多片线阵CCD相机的相对辐射定标方法，属于遥感卫星定标技术领域。

背景技术

目前，空间遥感CCD相机朝着更高分辨率、大幅宽等方向发展，从而保证更高的空间分辨率和时间分辨率。在轨道高度一定的前提下，空间相机地面像元分辨率的提高导致相机焦距的增长，而焦距的增长和视场角的增大导致需要采用多片线阵CCD拼接来组成较大幅宽空间相机的焦平面。

线阵CCD不同像元的光响应非一致性、读出电路自身及其与CCD耦合的非一致性、不同通道放大电路的非一致性、光学系统的加工精度等原因，导致光学系统和CCD探测器等响应的不一致性，直接表现为同样的辐射亮度输入对各个像元有着不一致性的输出灰度值，从而造成获取的图像出现条带现象。同时，大视场相机焦面由多片线阵CCD拼接而成，不同片CCD由于制造差异，响应度相差较大，导致不同片的图像存在较大的偏差。因此，有必要对多片线阵CCD相机进行相对辐射定标。目前，为减小光学畸变、地球曲率的影响，将多片线阵CCD阵列倾斜布置在焦平面上，用“弧线”形状代替常用的“直线”形状布置方式。例如，法国Pleiades卫星的全色谱段探测器组件由5片线阵CCD组成，地面分辨率0.7m，为减少光学畸变对TDI同步的影响，采取的措施是5片线阵CCD非共线布置，两两相邻两片之间的倾斜角度为0.3°。

目前常用的遥感卫星相对辐射定标方法中，实验室相对辐射定标获取的定标系数，随着卫星发射之后遥感器性能的下降，辐射定标系数会失效；均匀地面场景法需要湖泊、草原、冰川、人工均匀场地等定标场，寻找或建设难度大、成本高；大量卫星图像数据的在轨统计是目前国内采用比较广的一种方法，积累多轨图像数据，但该方法的前提是确保各轨图像数据的稳定性；国外应用广泛的Side-slither定标方法，将卫星偏航角旋转90°，选择雨林、沙漠、冰盖等相对均匀的场景作为定标源，因此，对定标源的要求高，导致对定标区域的选择受限。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供一种多片线阵CCD相机的相对辐射定标方法，建立从地理坐标系到焦平面坐标系的严密转换关系，针对多片线阵CCD在焦平面上倾斜布置的特点、以及光学畸变与地球曲率的影响，建立多片线阵CCD拟合后的成像曲线方程，在此基础之上精确计算偏航角速度，使得倾斜布置的多片线阵CCD所有像元对同一区域成像，该方法无论对星下点成像还是非星下点成像都具备模型严密、求解精度高等优点。

本发明的技术方案是：一种多片线阵CCD相机的相对辐射定标方法，步骤如下：

1)建立遥感卫星多片线阵CCD相机的焦平面坐标系P₀-P₁P₂P₃及其对应的地理坐标系G₀-G₁G₂G₃；其中多片线阵CCD的相邻两片之间以预设的倾斜角度布置在焦平面上；其中，P₀为焦平面坐标系的原点，位于多片线阵CCD的中点；P₁、P₂、P₃三个坐标轴分别平行于卫星本体坐标系S₀-S₁S₂S₃的滚动轴、俯仰轴和偏航轴；G₀为地理坐标系的坐标原点，且G₀是P₀点对应的地面点；G₁轴过G₀且平行于卫星飞行方向；G₃轴指向天顶并过地球中心；G₂轴垂直于G₁G₃构成的平面且符合右手法则；

建立多片线阵CCD拟合后的成像曲线方程p₁＝a_n|p₂|ⁿ+a_n-1|p₂|^n-1+…+a₁|p₂|+a₀，其中，(p₁,p₂)表示多片线阵CCD在焦平面坐标系中的坐标，a_n,a_n-1,…,a₁,a₀为方程系数，n为方程次数，|·|表示取绝对值；

2)建立从地理坐标系G₀-G₁G₂G₃到焦平面坐标系P₀-P₁P₂P₃的转换关系；

21)卫星绕本体坐标系的滚动轴旋转滚动角之后，P₀点对应的地面点从星下点K₀指向点K₁，则从星下点K₀指向点K₁对应的地心角表示为∠K₀OK₁，其中，O为地心；卫星再绕本体坐标系的俯仰轴旋转俯仰角θ之后，P₀点对应的地面点从点K₁指向点G₀，则从点K₁指向点G₀对应的地心角表示为∠K₁OG₀；卫星最后绕本体坐标系的偏航轴旋转偏航角ψ之后，P₀点对应的地面点仍为G₀；根据地球半径R、卫星轨道高度H以及卫星的俯仰角θ、滚动角，计算获得∠K₀OK₁与∠K₁OG₀；

22)从地理坐标系G₀-G₁G₂G₃到焦平面坐标系P₀-P₁P₂P₃的转换关系表示为：

其中，M₁为从地理坐标系G₀-G₁G₂G₃到地固坐标系O-E₁E₂E₃的转换矩阵，具体表示为：

其中h为地面点G₀的地形高度，γ₀为升交点与卫星之间的地心角，i₀为卫星轨道倾角；

M₂为从地固坐标系O-E₁E₂E₃到惯性坐标系O-I₁I₂I₃的转换矩阵，具体表示为：

其中ω为地球自转角速度；

M₃为从惯性坐标系O-I₁I₂I₃到轨道坐标系B₀-B₁B₂B₃的转换矩阵，具体表示为：

其中γ＝γ₀+Ωt，Ω为卫星轨道运行相对地心的角速度；

M₄为从轨道坐标系B₀-B₁B₂B₃到本体坐标系S₀-S₁S₂S₃的转换矩阵，具体表示为：

其中ψ₀、θ₀、分别表示卫星的初始偏航角、俯仰角与滚动角，分别表示卫星的偏航角速度、俯仰角速度、滚动角速度；

M₅为从本体坐标系S₀-S₁S₂S₃到相机坐标系C₀-C₁C₂C₃的转换矩阵，具体表示为：

其中f为相机的焦距；

M₆为从相机坐标系C₀-C₁C₂C₃到焦平面坐标系P₀-P₁P₂P₃的转换矩阵，具体表示为：

3)根据步骤2)中建立的转换关系，计算获得焦平面坐标系P₀-P₁P₂P₃的原点P₀对应的初始偏流角β₀；将卫星偏航角旋转90°-β₀后，进入相对辐射定标模式；

4)根据步骤1)建立的多片线阵CCD拟合后的成像曲线方程，在相对辐射定标模式的成像过程中，根据步骤2)中建立的转换关系计算获得卫星偏航角速度；多片线阵CCD的所有像元在该偏航角速度下对地面同一区域成像，获取图像数据；

5)根据步骤4)获取的图像数据，利用直方图匹配方法获得多片线阵CCD所有像元的直方图查找表。

所述步骤21)中根据地球半径R、卫星轨道高度H以及卫星的俯仰角θ、滚动角计算获得∠K₀OK₁与∠K₁OG₀的具体方法如下：

211)卫星绕本体坐标系的滚动轴旋转滚动角之后，计算获得从星下点K₀指向点K₁对应的地心角

212)卫星绕本体坐标系的俯仰轴旋转俯仰角θ之后，计算从点K₁指向点G₀对应的地心角其中

B₀为卫星轨道坐标系的原点。

所述步骤4)中卫星偏航角速度获得的具体方法为：

41)卫星相对辐射定标模式下，设定偏航角速度的范围为所述偏航角速度不超过卫星姿态机动能力范围；将偏航角速度作为循环变量，以作为偏航角速度的初始值，按照递进步长逐次递加，重复步骤411)至步骤414)，直到偏航角速度到达则循环结束，并转至步骤42)；

411)卫星相对辐射定标成像的起始时刻t₀记为0，t₀时刻对应的焦平面坐标系称为原焦平面坐标系，坐标原点为P₀，t₀时刻对应的地理坐标系称为原地理坐标系，坐标原点为G₀；

412)原焦平面坐标系以线性CCD的像元尺寸d_e作为步长沿着卫星飞行方向进行平移，同时原焦平面坐标系绕着坐标原点P₀以偏航角速度进行旋转，经过时间t₁＝n₁×t_e后，原焦平面坐标系的原点P₀运动到新焦平面坐标系P₀′-P₁′P₂′P₃′的原点P₀′；其中t_e为曝光时间，n₁为原焦平面坐标系移动的像元个数；计算获得P₀′在原焦平面坐标系中的坐标(p₁₀′,p₂₀′,0)，其中

根据公式(1)，计算获得新坐标原点P₀′在原地理坐标系G₀-G₁G₂G₃中的地理坐标，记为G₀′(g₁₀′,g₂₀′,0)；

413)设t₀时刻原焦平面坐标系中的某点B与t₁时刻新焦平面坐标系中的某点A′指向地面同一个区域，根据步骤1)建立的多片线阵CCD拟合后的成像曲线方程，设定B点在原焦平面坐标系中的坐标表示为P_B(p_1B,p_2B,0)，根据公式(1)计算获得B点在原地理坐标系G₀-G₁G₂G₃中对应的地理坐标为G_B(g_1B,g_2B,0)；设定A′点在新焦平面坐标系P₀′-P₁′P₂′P₃′中的坐标表示为P_{A′_new}(p_{1A′_new},p_{2A′_new},0)，经过时间t₁之后，A′点与B点的关系表示为p_{2A′_new}＝(p_2B-p₂₀′/cosβ₀)，其中p_{1A′_new}由步骤1)建立的多片线阵CCD拟合后的成像曲线方程计算得到；根据公式(1)计算A′点在新地理坐标系G₀′-G₁′G₂′G₃′中对应的地理坐标表示为G_{A′_new}(g_{1A′_new},g_{2A′_new},0)；

414)将A′点与B点的地理坐标统一到原地理坐标系中，则A′点与B点在原地理坐标系中的偏移量Δx＝g_{1A′_new}+g₁₀′-g_1B、Δy＝g_{2A′_new}+g₂₀′-g_2B，计算获得平方和Δs＝Δx²+Δy²；

42)遍历步骤41)中获得的所有平方和，并取平方和最小值对应的偏航角速度作为卫星偏航角速度。

所述步骤1)中预设的倾斜角度小于等于3°。

本发明与现有技术相比有益效果为：

(1)本发明在定标成像模式下使得倾斜布置的线阵CCD所有像元对同一个区域进行成像，较目前传统的相对辐射定标方法，建立了相同的输入数据源，即能够在具有同样的辐亮度景物的输入前提下，根据各个像元的输出对其进行相对辐射定标。

(2)本发明提出了基于严密数学模型的倾斜布置多片线阵CCD相机的相对辐射定标方法，建立了从地理坐标系到焦平面坐标系的转换关系，据此计算卫星偏航角速度，同时利用基于直方图匹配的相对辐射定标非线性模型，无论对星下点还是非星下点成像该方法都具备模型严密、精度高等优点。

(3)本发明提出了在焦平面上基于曲线拟合多片线阵CCD的策略，考虑到在轨运行时，由于镜头光学畸变、地球曲率的影响，导致在倾斜布置的多片线阵CCD成像模型一般表现为曲线，该策略符合卫星的实际成像情况，利用严密转换关系以及卫星姿态、轨道等参数将光学系统的变形体现在像方和物方上。

(4)本发明采用试探法计算偏航角速度，避免了解算复杂三角函数方程，在地面获得定标成像过程中的姿态角速度。

附图说明

图1是本发明一种多片线阵CCD相机的相对辐射定标流程图；

图2是本发明在常规成像模式下倾斜布置的5片线阵CCD相机焦平面与卫星飞行方向的关系示意图；

图3是本发明由于卫星滚动角与俯仰角引起的偏离星下点的地心角∠K₀OK₁与∠K₁OG₀示意图；

图4是本发明从地理坐标系到焦平面坐标系的转换示意图；

图5a是本发明从常规成像模式到相对辐射定标模式示意图；

图5b是本发明在相对辐射定标模式下倾斜布置的5片线阵CCD相机焦平面与卫星飞行方向关系的示意图；

图6是本发明对同一区域成像作为条件计算偏航角速度的示意图。

具体实施方式

本发明一种多片线阵CCD相机的相对辐射定标方法，具体步骤如图1所示，该方法由以下步骤实现：

1、建立遥感卫星多片线阵CCD相机的焦平面坐标系P₀-P₁P₂P₃及其对应的地理坐标系G₀-G₁G₂G₃；其中多片线阵CCD的相邻两片之间以预设的倾斜角度布置在焦平面上。

其中，P₀为焦平面坐标系的原点，位于多片线阵CCD的中点；P₁、P₂、P₃三个坐标轴过P₀点，且分别平行于卫星本体坐标系S₀-S₁S₂S₃的滚动轴、俯仰轴和偏航轴；G₀为地理坐标系的坐标原点，且G₀是P₀点对应的地面点；G₁轴过G₀且平行于卫星飞行方向；G₃轴指向天顶并过地球中心；G₂轴垂直于G₁G₃构成的平面且符合右手法则。

由于光学部件设计与加工缺陷的存在，相机镜头会存在一定的光学畸变，为减少光学畸变、地球曲率等对线阵CCD成像质量的影响，遥感卫星的多片线阵CCD一般不是完全共线布置在焦平面上，而是相邻两片之间以一定预设的倾斜偏置角度形成折线布置，一般情况下倾斜偏置角度小于3°。在轨运行时由于光学畸变、地球曲率等因素的影响，导致在倾斜布置的多片线阵CCD成像模型一般表现为曲线，建立多片线阵CCD拟合后的成像曲线方程p₁＝a_n|p₂|ⁿ+a_n-1|p₂|^n-1+…+a₁|p₂|+a₀，其中，(p₁,p₂)表示多片线阵CCD在焦平面坐标系中的坐标，a_n,a_n-1,…,a₁,a₀为方程系数，n为方程次数，一般选择n＝3、4、5，|·|表示取绝对值。

如图2所示为在常规成像模式下倾斜布置的5片线阵CCD相机焦平面与卫星飞行方向关系的示意图，可以看出，在卫星常规成像模式下P₁方向近似沿着卫星飞行方向，其中，相邻两片之间存在0.3°的偏置角度，每片线阵CCD包括6000个像元，本发明的实施例，将利用线阵CCD上的5个点：A(0,1500)，B(7.854,4499.979)，C(23.5618,7499.938)，D(47.123,10499.835)，E(78.539,13499.671)，以多项式方程p₁＝a₃|p₂|³+a₂|p₂|²+a₁|p₂|+a₀来描述线阵CCD的数学模型，利用最小二乘拟合得到方程的系数为：a＝-2.8149×10^-12，b＝4.96×10^-7，c＝-3.192×10^-4，d＝-0.6019。实际卫星在轨定标时，可以通过几何定标得到实际焦平面曲线的数学模型。

2、建立从地理坐标系G₀-G₁G₂G₃到焦平面坐标系P₀-P₁P₂P₃的转换关系，转换过程中卫星姿态转动的次序设定为先滚动、后俯仰、最后偏航。

2.1、卫星绕本体坐标系的滚动轴旋转滚动角之后，如图3所示，P₀点对应的地面点从星下点K₀指向点K₁，则从星下点K₀指向点K₁对应的地心角，即滚动角旋转引起的地心角表示为∠K₀OK₁，其中，O为地心；卫星再绕本体坐标系的俯仰轴旋转俯仰角θ之后，P₀点对应的地面点从点K₁指向点G₀，则从点K₁指向点G₀对应的地心角，即俯仰角旋转引起的地心角表示为∠K₁OG₀；卫星最后绕本体坐标系的偏航轴旋转偏航角ψ之后，P₀点对应的地面点仍为G₀；根据地球半径R、卫星轨道高度H以及卫星的俯仰角θ、滚动角计算获得∠K₀OK₁与∠K₁OG₀。

(1)卫星绕本体坐标系的滚动轴旋转滚动角之后，计算从星下点K₀指向点K₁对应的地心角∠K₀OK₁：

B₀为卫星轨道坐标系的原点，轨道高度B₀K₀＝H，在ΔB₀OK₁中，由正弦定理可得到：

在ΔB₀OK₁中，由余弦定理得到：

(2)卫星绕本体坐标系的俯仰轴旋转俯仰角θ之后，计算从点K₁指向点G₀对应的地心角∠K₁OG₀：

在ΔB₀OG₀中，得到：

由正弦定理、以及三角形三角之和等于180°，得到：

因此，在ΔB₀OG₀中，由余弦定理得到：

在ΔB₀K₁G₀中，由余弦定理得到：

在等腰ΔK₁OG₀中，得到：

2.2、从地理坐标系G₀-G₁G₂G₃到焦平面坐标系P₀-P₁P₂P₃的转换关系表示为：

其中：

M₁为从地理坐标系G₀-G₁G₂G₃到地固坐标系O-E₁E₂E₃的转换矩阵，包括沿G₃轴平移-(R+h)、绕G₂轴旋转∠K₁OG₀、绕G₁轴旋转∠K₀OK₁、绕G₂轴旋转-γ₀、绕G₃轴旋转-i₀，具体表示为：

其中h为地面点G₀的地形高度，γ₀为摄影时刻升交点与卫星之间所对应的地心角，i₀为卫星轨道倾角。

M₂为从地固坐标系O-E₁E₂E₃到惯性坐标系O-I₁I₂I₃的转换矩阵，绕E₂轴旋转-ωt，具体表示为：

其中ω为地球自转角速度。

M₃为从惯性坐标系O-I₁I₂I₃到轨道坐标系B₀-B₁B₂B₃的转换矩阵，绕I₃轴旋转i₀，绕I₂轴旋转γ，沿I₃轴平移(R+H)，具体表示为：

其中γ＝γ₀+Ωt，Ω为卫星轨道运行相对地心的角速度。

M₄为从轨道坐标系B₀-B₁B₂B₃到本体坐标系S₀-S₁S₂S₃的转换矩阵，绕B₁转滚动角，绕B₂转俯仰角θ、绕B₃转偏航角ψ，具体表示为：

其中ψ₀、θ₀、分别表示卫星的初始偏航角、俯仰角、滚动角，分别卫星的偏航角速度、俯仰角速度、滚动角速度。

M₅为从本体坐标系S₀-S₁S₂S₃到相机坐标系C₀-C₁C₂C₃的转换矩阵，缩放比例为f/(H-h)，具体表示为：

其中f为相机的焦距。

M₆为从相机坐标系C₀-C₁C₂C₃到焦平面坐标系P₀-P₁P₂P₃的转换矩阵，沿C₃平移-f，具体表示为：

所述惯性坐标系O-I₁I₂I₃原点O在地心，I₂轴指向北极，I₃轴指向卫星轨道平面与赤道面的交点，I₁轴垂直I₂和I₃形成的平面；地固坐标系O-E₁E₂E₃原点O在地心，E₀E₂轴沿地轴指向北极方向，该坐标系在惯性坐标系内，绕逆时针方向以地球角速度ω自转；卫星轨道坐标系B₀-B₁B₂B₃原点B₀在轨道上，B₁轴指向轨道前进方向，B₃轴指向天顶，B₂轴与轨道面垂直，卫星以角速度Ω作轨道运动；卫星本体坐标系S₀-S₁S₂S₃：坐标原点S₀与B₀重合，卫星无姿态运动时S与B系重合，滚动角，俯仰角θ，偏航角ψ即指S系在B系内的三轴姿态(顺序→θ→ψ)；相机坐标系C₀-C₁C₂C₃：无安装误差时相机坐标系与本体坐标系重合；焦平面坐标系P₀-P₁P₂P₃：C系沿C₃轴平移得到P系。

3、根据步骤2中建立的转换关系，计算获得卫星常规成像模式下焦平面坐标系P₀-P₁P₂P₃原点P₀的初始偏流角β₀；将卫星偏航角旋转90°-β₀后，进入相对辐射定标模式。

如图5a所示，卫星在轨道A位置为常规成像模式，根据公式(1)计算卫星在轨道A位置对应的初始偏流角β₀，然后将卫星偏航角旋转90°-β₀，在轨道位置B卫星进入相对辐射定标成像模式，卫星偏航角旋转90°-β₀而不是90°的目的是消除地球偏转引起的偏流角的影响；图5b是在相对辐射定标模式下倾斜布置的5片线阵CCD相机焦平面与卫星飞行方向关系的示意图,此时的P₂轴与卫星飞行方向的夹角为偏流角β₀。

所述初始偏流角的计算通过公式计算得到，其中，分别为相机焦平面中P₁和P₂方向的像移速度，可通过将公式(1)两边对时间t求微分，即但是，由于微分运算十分复杂，本发明利用差分法(参见：《基于差分法的空间相机像移速度矢量计算》，光学精密工程，2011.5)计算初始偏流角β₀，具体步骤为：首先，在T时刻，根据式(1)计算出焦平面上点P₀(P_{1_P0},P_{2_P0},0)对应的地理坐标(G_{1_P0},G_{2_P0},0)；然后，在T+Δt时刻，利用式(1)计算出同一个地理坐标(G_{1_P0},G_{2_P0},0)在焦平面坐标中的坐标(P_{1_P0}′,P_{2_P0}′,0)，Δt要足够短，例如可取Δt＝0.0002s；最后，利用差分法计算像移速度：

4、根据步骤1建立的多片线阵CCD拟合后的成像曲线方程，在相对辐射定标模式的成像过程中，根据步骤2中建立的转换关系计算获得卫星偏航角速度；多片线阵CCD的所有像元在该偏航角速度下对地面同一区域成像，获取图像数据。

根据所建立的从地理坐标系到焦平面坐标系转换的严密模型，把M₁～M₆矩阵分别代入公式(1)，得到复杂的三角函数方程，将偏航角速度作为循环变量，通过不断循环计算得到满足条件的偏航角速度，具体步骤如下：

4.1、卫星相对辐射定标模式下，设定偏航角速度的可取范围为其中，不应超过卫星姿态机动能力范围；将偏航角速度作为循环变量，以作为偏航角速度的初始值，按照递进步长逐次递加，重复步骤4.1.1至步骤4.1.4，直到偏航角速度到达则循环结束，并转至步骤4.2。

4.1.1、卫星相对辐射定标成像的起始时刻t₀记为0，t₀时刻所对应的焦平面坐标系称为原焦平面坐标系，坐标原点P₀，t₀时刻对应的地理坐标系称为原地理坐标系，坐标原点G₀。

4.1.2、原焦平面坐标系以线性CCD的像元尺寸d_e作为步长沿着卫星飞行方向进行平移，同时原焦平面坐标系绕着坐标原点P₀以偏航角速度进行旋转，经过时间t₁＝n₁×t_e后，原焦平面坐标系的原点P₀运动到新焦平面坐标系P₀′-P₁′P₂′P₃′的原点P₀′；其中t_e为曝光时间，n₁为原焦平面坐标系移动的像元个数；计算获得P₀′在原焦平面坐标系中的坐标(p₁₀′,p₂₀′,0)，其中：

根据公式(1)，计算获得新坐标原点P₀′在原地理坐标系G₀-G₁G₂G₃中的地理坐标，记为G₀′(g₁₀′,g₂₀′,0)。

4.1.3、根据多片线阵CCD的所有像元通过地面同一个区域这一条件，设t₀时刻原焦平面坐标系中的B点与t₁时刻新焦平面坐标系中的A′点指向地面同一个区域，根据步骤1)建立的多片线阵CCD拟合后的成像曲线方程，设定B点在原焦平面坐标系中的坐标表示为P_B(p_1B,p_2B,0)，根据公式(1)计算获得B点在原地理坐标系G₀-G₁G₂G₃中对应的地理坐标为G_B(g_1B,g_2B,0)；设定A′点在新焦平面坐标系P₀′-P₁′P₂′P₃′中的坐标表示为P_{A′_new}(p_{1A′_new},p_{2A′_new},0)，经过时间t₁之后，A′点与B点的关系表示为p_{2A′_new}＝(p_2B-p₂₀′/cosβ₀)，其中p_{1A′_new}由步骤1建立的多片线阵CCD拟合后的成像曲线方程计算得到；根据公式(1)计算A′点在新地理坐标系G₀′-G₁′G₂′G₃′中对应的地理坐标表示为G_{A′_new}(g_{1A′_new},g_{2A′_new},0)。

4.1.4、将A′点与B点的地理坐标统一到原地理坐标系中，则A′点与B点在原地理坐标系中的偏移量Δx＝g_{1A′_new}+g₁₀′-g_1B、Δy＝g_{2A′_new}+g₂₀′-g_2B，计算获得平方和Δs＝Δx²+Δy²。

4.2、遍历步骤4.1中获得的所有平方和，并取平方和最小值对应的偏航角速度作为卫星偏航角速度。

如图6所示，t₀时刻，原焦平面坐标系P₀-P₁P₂P₃、原地理坐标系G₀-G₁G₂G₃，经过时间t₁之后，得到新焦平面坐标系P₀′-P₁′P₂′P₃′，设定t₀时刻原焦平面坐标系中的B点与t₁时刻新焦平面坐标系中的A′点指向地面同一个区域。

在执行定标任务时对计算偏航角速度进行了仿真，设相机焦距f＝2.6m，轨道高度H＝645km，轨道倾角i₀＝97°，t₀时刻卫星在北纬40°，相机曝光时间t_e＝0.3255μs；计算偏航角速度选取的为[-2,2]°/s，递进步长像元尺寸d_e＝10μm，如图6所示，以原焦平面坐标系下的点B为基准，p_2B＝15000，焦平面坐标系的原点沿着卫星飞行方向移动的像元个数n₁＝6000，则根据上述步骤，计算偏航角速度为-0.1064°/s。

5、根据步骤4获取的图像数据，利用直方图匹配方法获得多片线阵CCD所有像元的直方图查找表。

(1)对获得多片线阵CCD的所有像元对同一区域成像的图像数据进行45°旋转，其旋转图像的同一行为所有像元对同一区域的成像。

每个像元的成像数据在输出图像上表现为一列数据，而由于线阵CCD所有像元对同一区域先后依次成像，导致输出图像的45°对角线上的数据为不同像元对同一区域的成像。

(2)计算整个图像的直方图概率密度函数，作为目标期望直方图。

(3)计算每个像元图像的直方图概率密度函数，作为原始直方图。

(4)基于直方图匹配建立原始直方图与目标期望直方图的直方图映射关系，得到该像元的非线性定标映射关系表。

(5)重复步骤(3)和(4)，用同样的方法处理所有的像元，便可以得到所有像元的直方图查找表。

本发明未详细说明部分属本领域技术人员公知常识。

Claims (3)

1.一种多片线阵CCD相机的相对辐射定标方法，其特征在于步骤如下：

1)建立遥感卫星多片线阵CCD相机的焦平面坐标系P₀-P₁P₂P₃及其对应的地理坐标系G₀-G₁G₂G₃；其中多片线阵CCD的相邻两片之间以预设的倾斜角度布置在焦平面上；其中，P₀为焦平面坐标系的原点，位于多片线阵CCD的中点；P₁、P₂、P₃三个坐标轴分别平行于卫星本体坐标系S₀-S₁S₂S₃的滚动轴、俯仰轴和偏航轴；G₀为地理坐标系的坐标原点，且G₀是P₀点对应的地面点；G₁轴过G₀且平行于卫星飞行方向；G₃轴指向天顶并过地球中心；G₂轴垂直于G₁G₃构成的平面且符合右手法则；

建立多片线阵CCD拟合后的成像曲线方程p₁＝a_n|p₂|ⁿ+a_n-1|p₂|^n-1+L+a₁|p₂|+a₀，其中，(p₁,p₂)表示多片线阵CCD在焦平面坐标系中的坐标，a_n,a_n-1,L,a₁,a₀为方程系数，n为方程次数，|·|表示取绝对值；

2)建立从地理坐标系G₀-G₁G₂G₃到焦平面坐标系P₀-P₁P₂P₃的转换关系；

21)卫星绕本体坐标系的滚动轴旋转滚动角之后，P₀点对应的地面点从星下点K₀指向点K₁，则从星下点K₀指向点K₁对应的地心角表示为∠K₀OK₁，其中，O为地心；卫星再绕本体坐标系的俯仰轴旋转俯仰角θ之后，P₀点对应的地面点从点K₁指向点G₀，则从点K₁指向点G₀对应的地心角表示为∠K₁OG₀；卫星最后绕本体坐标系的偏航轴旋转偏航角ψ之后，P₀点对应的地面点仍为G₀；根据地球半径R、卫星轨道高度H以及卫星的俯仰角θ、滚动角计算获得∠K₀OK₁与∠K₁OG₀；

22)从地理坐标系G₀-G₁G₂G₃到焦平面坐标系P₀-P₁P₂P₃的转换关系表示为：

$\left[\begin{array}{c}{P}_1\\ P_2\\ P_3\\ 1\end{array}\right]=M_6{M}_5{M}_4{M}_3{M}_2{M}_1\left[\begin{array}{c}{G}_1\\ G_2\\ 0\\ 1\end{array}\right];---\left(1\right)$

其中，M₁为从地理坐标系G₀-G₁G₂G₃到地固坐标系O-E₁E₂E₃的转换矩阵，具体表示为：

$\begin{array}{c}{M}_1=\left[\begin{array}{cccc}{\mathrm{cosi}}_0& -{\mathrm{sini}}_0& 0& 0\\ {\mathrm{sini}}_0& {\mathrm{cosi}}_0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{\mathrm{cos\γ}}_0& 0& {\mathrm{sin\γ}}_0& 0\\ 0& 1& 0& 0\\ -{\mathrm{sin\γ}}_0& 0& {\mathrm{cos\γ}}_0& 0\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \∠K_0{\mathrm{OK}}_1& \sin \∠K_0{\mathrm{OK}}_1& 0\\ 0& -\sin \∠K_0{\mathrm{OK}}_1& \cos \∠K_0{\mathrm{OK}}_1& 0\\ 0& 0& 0& 1\end{array}\right]\\ \left[\begin{array}{cccc}\cos \∠K_1{\mathrm{OG}}_0& 0& -\sin \∠K_1{\mathrm{OG}}_0& 0\\ 0& 1& 0& 0\\ \sin \∠K_1{\mathrm{OG}}_0& 0& \cos \∠K_1{\mathrm{OG}}_0& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& (R+h)\\ 0& 0& 0& 1\end{array}\right]\end{array};$

其中h为地面点G₀的地形高度，γ₀为升交点与卫星之间的地心角，i₀为卫星轨道倾角；

M₂为从地固坐标系O-E₁E₂E₃到惯性坐标系O-I₁I₂I₃的转换矩阵，具体表示为：

$M_2=\left[\begin{array}{cccc}cos\mathrm{\ω}t& 0& sin\mathrm{\ω}t& 0\\ 0& 1& 0& 0\\ -sin\mathrm{\ω}t& 0& \cos \mathrm{\ω}t& 0\\ 0& 0& 0& 1\end{array}\right];$

其中ω为地球自转角速度；

M₃为从惯性坐标系O-I₁I₂I₃到轨道坐标系B₀-B₁B₂B₃的转换矩阵，具体表示为：

$M_3=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& -(R+H)\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}cos\mathrm{\γ}& 0& -sin\mathrm{\γ}& 0\\ 0& 1& 0& 0\\ sin\mathrm{\γ}& 0& \cos \mathrm{\γ}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{\mathrm{cosi}}_0& {\mathrm{sini}}_0& 0& 0\\ -{\mathrm{sini}}_0& {\mathrm{cosi}}_0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right];$

其中γ＝γ₀+Ωt，Ω为卫星轨道运行相对地心的角速度；

M₄为从轨道坐标系B₀-B₁B₂B₃到本体坐标系S₀-S₁S₂S₃的转换矩阵，具体表示为：

其中ψ₀、θ₀、分别表示卫星的初始偏航角、俯仰角与滚动角，分别表示卫星的偏航角速度、俯仰角速度、滚动角速度；

M₅为从本体坐标系S₀-S₁S₂S₃到相机坐标系C₀-C₁C₂C₃的转换矩阵，具体表示为：

$M_5=\left[\begin{array}{cccc}\frac{f}{H-h}& 0& 0& 0\\ 0& \frac{f}{H-h}& 0& 0\\ 0& 0& \frac{f}{H-h}& 0\\ 0& 0& 0& 1\end{array}\right];$

其中f为相机的焦距；

M₆为从相机坐标系C₀-C₁C₂C₃到焦平面坐标系P₀-P₁P₂P₃的转换矩阵，具体表示为：

$M_6=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& f\\ 0& 0& 0& 1\end{array}\right];$

3)根据步骤2)中建立的转换关系，计算获得焦平面坐标系P₀-P₁P₂P₃的原点P₀对应的初始偏流角β₀；将卫星偏航角旋转90°-β₀后，进入相对辐射定标模式；

4)根据步骤1)建立的多片线阵CCD拟合后的成像曲线方程，在相对辐射定标模式的成像过程中，根据步骤2)中建立的转换关系计算获得卫星偏航角速度；多片线阵CCD的所有像元在该偏航角速度下对地面同一区域成像，获取图像数据；其中根据步骤2)中建立的转换关系计算获得卫星偏航角速度的具体过程如下：

41)卫星相对辐射定标模式下，设定偏航角速度的范围为所述偏航角速度不超过卫星姿态机动能力范围；将偏航角速度作为循环变量，以作为偏航角速度的初始值，按照递进步长逐次递加，重复步骤411)至步骤414)，直到偏航角速度到达则循环结束，并转至步骤42)；

411)卫星相对辐射定标成像的起始时刻t₀记为0，t₀时刻对应的焦平面坐标系称为原焦平面坐标系，坐标原点为P₀，t₀时刻对应的地理坐标系称为原地理坐标系，坐标原点为G₀；

412)原焦平面坐标系以线性CCD的像元尺寸d_e作为步长沿着卫星飞行方向进行平移，同时原焦平面坐标系绕着坐标原点P₀以偏航角速度进行旋转，经过时间t₁＝n₁×t_e后，原焦平面坐标系的原点P₀运动到新焦平面坐标系P₀′-P₁′P₂′P₃′的原点P₀′；其中t_e为曝光时间，n₁为原焦平面坐标系移动的像元个数；计算获得P₀′在原焦平面坐标系中的坐标(p₁₀′,p₂₀′,0)，其中

$\left\{\begin{array}{c}{p_{10}}^{\′}=n_1\×d_e\×{\mathrm{sin\β}}_0\\ {p_{20}}^{\′}=n_1\×d_e\×{\mathrm{cos\β}}_0\end{array}\right.;$

根据公式(1)，计算获得新坐标原点P₀′在原地理坐标系G₀-G₁G₂G₃中的地理坐标，记为G₀′(g₁₀′,g₂₀′,0)；

413)设t₀时刻原焦平面坐标系中的某点B与t₁时刻新焦平面坐标系中的某点A′指向地面同一个区域，根据步骤1)建立的多片线阵CCD拟合后的成像曲线方程，设定B点在原焦平面坐标系中的坐标表示为P_B(p_1B,p_2B,0)，根据公式(1)计算获得B点在原地理坐标系G₀-G₁G₂G₃中对应的地理坐标为G_B(g_1B,g_2B,0)；设定A′点在新焦平面坐标系P₀′-P₁′P₂′P₃′中的坐标表示为P_{A′_new}(p_{1A′_new},p_{2A′_new},0)，经过时间t₁之后，A′点与B点的关系表示为p_{2A′_new}＝(p_2B-p₂₀′/cosβ₀)，其中p_{1A′_new}由步骤1)建立的多片线阵CCD拟合后的成像曲线方程计算得到；根据公式(1)计算A′点在新地理坐标系G₀′-G₁′G₂′G₃′中对应的地理坐标表示为G_{A′_new}(g_{1A′_new},g_{2A′_new},0)；

414)将A′点与B点的地理坐标统一到原地理坐标系中，则A′点与B点在原地理坐标系中的偏移量Δx＝g_{1A′_new}+g₁₀′-g_1B、Δy＝g_{2A′_new}+g₂₀′-g_2B，计算获得平方和Δs＝Δx²+Δy²；

42)遍历步骤41)中获得的所有平方和，并取平方和最小值对应的偏航角速度作为卫星偏航角速度；

5)根据步骤4)获取的图像数据，利用直方图匹配方法获得多片线阵CCD所有像元的直方图查找表。

2.根据权利要求1所述的一种多片线阵CCD相机的相对辐射定标方法，其特征在于：所述步骤21)中根据地球半径R、卫星轨道高度H以及卫星的俯仰角θ、滚动角计算获得∠K₀OK₁与∠K₁OG₀的具体方法如下：

211)卫星绕本体坐标系的滚动轴旋转滚动角之后，计算获得从星下点K₀指向点K₁对应的地心角

212)卫星绕本体坐标系的俯仰轴旋转俯仰角θ之后，计算从点K₁指向点G₀对应的地心角其中 B₀为卫星轨道坐标系的原点。

3.根据权利要求1所述的一种多片线阵CCD相机的相对辐射定标方法，其特征在于：所述步骤1)中预设的倾斜角度小于等于3°。