CN102322850B - 一种基于成像质量预估的任务预处理方法 - Google Patents

Abstract

一种基于成像质量预估的任务预处理方法，首先根据卫星轨道数据、姿态机动能力和地面目标位置信息，计算卫星对地面目标的可见时间窗口，从任务列表中剔除无可见时间窗口及时间窗口在地影区的任务。然后分析确定影响图像质量的影响因素，构造成像质量确定函数，继续从任务列表中剔除不满足成像质量要求的任务。对于任务列表中剩余的每一个任务，对与各任务对应的可见时间窗口进行裁剪，得到任务的执行时间窗口，同时根据任务执行过程中各时间点卫星指向目标顶点的姿态角计算任务的观测持续时间。当任务的观测持续时间在任务执行时间窗口之内时，该任务为可执行任务，所有可执行的任务将会提交给卫星任务规划系统进行任务规划与调度。

Description

一种基于成像质量预估的任务预处理方法

技术领域

本发明属于卫星任务规划与调度领域，涉及一种卫星任务规划过程中的任务预处理方法。

背景技术

快速姿态机动成像卫星借助快速姿态机动能力，可实现多种复杂成像模式。与采用星下点成像的传统对地观测卫星相比，快速姿态机动能力大大增加了卫星对目标的观测机会，因而具有更强的观测能力。快速姿态机动成像卫星的每一种成像模式都伴随着多个姿态机动、相机开关机等操作，这些操作形成一个前后连贯的控制指令序列。由于指令繁多，无法保证指令编排和上注的可靠性和指令执行的实时性问题，因此必须建立一套任务规划与调度系统，完成大批量观测任务的自动化分析与处理。

对任务的预处理是快速姿态机动成像卫星任务规划的关键环节。任务预处理方法在多篇公开文献中已有描述，如刘雄在2006年第23卷7期《计算机仿真》上发表的《卫星全球普查任务规划系统预处理模块的开发》，李曦在其硕士毕业论文《多星区域观测任务的效率优化方法研究》，以及李菊芳在其博士毕业论文《航天侦察多星多地面站任务规划问题研究》中，都提到了对观测任务进行预处理。这些文献都是采用STK软件进行任务预处理的，其基本步骤为在STK软件中构造仿真场景，建立卫星、地面目标和传感器等仿真实体，设置卫星轨道、传感器视场等参数，即可通过传感器与目标的可见时间得到时间窗口数据。这种方法实现简单，在卫星任务规划方法研究和软件产品中得到了广泛应用。但是，基于STK软件进行任务预处理，会使得最终的优化结果受到STK软件的影响，STK软件的误差也不可避免的会带入到最终的优化结果中，并且这种方法没有考虑不同观测条件下图像质量的变化，如卫星侧摆角或俯仰角过大时图像的分辨率下降等问题。

针对上述问题，也有部分文献在时间窗口计算过程中引入了对图像质量因素的考虑，如宗建建在2005年第26卷6期《计算机工程与设计》上发表的“卫星任务规划系统时间窗口模块的设计与实现”一文，将时间窗口长度、分辨率和目标太阳高度角作为时间窗口计算的三项主要约束来确定时间窗口，当这三项约束比预设临界值差时剔除时间窗口。实际上，卫星图像质量的影响因素众多，根据卫星设计和任务需求的不同，评定方法也具有多样化，仅仅挑取其中有限的几项参数作为图像质量估计的指标难以满足不同的需求，如有的观测目标图像质量侧重于几何分辨率，而有的侧重于信噪比，也有可能引入其它指标。此外，快速姿态机动成像卫星由于具备大角度姿态机动能力，在同一段可见时间区间采用不同姿态角对同一目标进行观测时成像质量存在很大差异，因此该方法对快速姿态机动成像卫星不适用。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种通用性强的基于成像质量预估的卫星任务预处理方法。

本发明的技术解决方案是：一种基于成像质量预估的任务预处理方法，步骤如下：

(1)根据相机幅宽将任务区域条带划分；

(2)根据卫星的最大俯仰角和最大滚动角以及任务条带信息确定卫星对各条带的可见时间区间[t₀，t_n]，剔除无可见时间区间及时间区间在地影区的任务；每个条带的可见时间区间[t₀，t_n]由卫星对该条带四个顶点的可见时间区间[t₁，t₂]_k，k＝1，2，3，4求交集并剔除在地影区的时间区间后得到；卫星对条带单个顶点的可见时间区间由卫星的最大俯仰角和最大滚动角确定；

(3)将调制传递函数MTF、推扫方向的几何分辨率、线阵方向的几何分辨率、信噪比和幅宽作为影响图像质量的主要因素，在步骤(2)确定的可见时间区间[t₀，t_n]内将图像质量表示为时间的一元函数R(t)，

R(t)＝η_MTFR_MTF+η_GSDxR_GSDx(t)+η_GSDyR_GSDy(t)+η_snrR_snr(t)+η_breathR_breath(t)，

式中η_MTF+η_GSDx+η_GSDy+η_snr+η_breath＝1，R_MTF，R_GSDx(t)，R_GSDy(t)，R_snr(t)，R_breath(t)分别为调制传递函数MTF、推扫方向的几何分辨率、线阵方向的几何分辨率、信噪比和幅宽与时间的一元函数关系；

(4)采用步骤(3)得到的一元函数表达式分别计算在步骤(2)中确定的可见时间区间[t₀，t_n]的两个端点t₀和t_n以及区间中间点

的图像质量R₀＝R(t₀)，R_n＝R(t_n)，

(5)记最低许可成像质量为R_p，将R_p与步骤(4)得到的R₀，R_n，R_max进行比较，确定满足成像质量要求的可见时间区间，剔除不满足成像质量要求的任务；

(6)采用弦截法计算每个任务条带的观测持续时间T_last，T_last为从该条带最先观测的顶点开始到该条带最后观测的顶点结束所持续的时间；当任务条带的观测持续时间T_last小于步骤(5)的任务观测窗口时间区间长度时，该任务为可执行任务，将所有可执行任务提交给任务规划系统进行任务规划与调度。

本发明与现有技术相比的优点在于：

(1)本发明方法在根据用户的成像质量要求确定的时间窗口基础上，对任务进行了初步筛选，剔除了“无可见时间窗口”、“不满足成像质量要求”、“时间窗口在地影区”、“不满足观测持续时间”的任务，对立体成像任务，根据用户要求的立体成像角度对其可执行性进行了判断，达到了削减任务规划问题的决策数量及可行域的目的，降低了问题的复杂度；

(2)传统任务规划问题中，都将目标简化为点目标，设定统一的较小的观测持续时间以保证能获取指定的地面目标图像；快速姿态机动成像卫星由于具备多种复杂成像模式，面临的目标大小存在很大差异，并且任务之间姿态机动频繁，耗费的时间各不相同，设定统一的观测持续时间，将影响快速姿态机动成像卫星发挥其高效能的优势。本发明方法根据观测目标的具体位置，通过确定任务对应星下点的时间，得到任务的观测开始时间和结束时间，计算出较精确的任务观测持续时间，以满足快速姿态机动成像卫星任务规划的需要；

(3)本发明方法将成像质量作为时间窗口计算的依据，对影响图像质量的主要动态指标进行了综合考虑，对于不同地表属性和观测需求，能够通过调整权重、地表反射率等参数，将不同观测目标区别对待，适应不同观测目标时间窗口裁剪要求多样化的需求，能够满足快速姿态机动成像卫星任务规划的需要，并且成像质量计算函数与时间窗口裁剪方法解耦，具有良好的通用性和灵活性；

(4)本发明方法采用弦截法作为求解任务的观测开始时间和结束时间，以及最低许可成像质量方程的核心算法，不仅能够得到任意精度的数值解，而且迭代次数较少，计算速度快，能够满足大规模优化问题对于算法时间复杂度的约束。

附图说明

图1为本发明方法的流程框图；

图2为本发明可见时间窗口计算流程图；

图3为本发明幅宽计算公式推导示意图；

图4为本发明实施例中任务3在可见时间窗口[t₀，t_n]内观测俯仰角的变化情况示意图；

图5为本发明实施例中任务3在可见时间窗口[t₀，t_n]内观测滚动角的变化情况示意图；

图6为本发明实施例中任务3在可见时间窗口[t₀，t_n]内太阳高度角的变化情况示意图；

图7为本发明实施例中任务3在可见时间窗口[t₀，t_n]内成像质量的变化情况示意图。

具体实施方式

如图1所示，本发明快速姿态机动成像卫星批量观测任务预处理方法的流程包括：1.任务区域条带划分；2.计算卫星的位置及速度；3.计算卫星对每个任务条带的可见时间区间；4.构造成像质量确定函数；5.判断时间窗口的存在性；6.采用弦截法裁剪时间窗口；7.确定任务的观测持续时间；8.输出各观测目标的信息。下面分别进行详细介绍：

一、任务区域条带划分

通常，观测任务区域通过多个顶点的经纬度描述，将这些顶点依次连接即得到任务区域。本发明中，采用如下的方法将任务区域划分为平行于卫星轨道的条带，便于卫星实施观测：

1.从卫星星下线出发，以星下点相机幅宽(卫星侧摆角为0时的相机幅宽)为宽度，做星下线的平行线，直至覆盖任务区域；星下线即星下点(卫星位置点在地球表面上的垂直投影点)的集合；

2.从任务区域的每个顶点做卫星星下线的垂线，并计算垂足之间的距离，记垂足之间距离最长的两个垂足分别为L1和L2，与L1对应的任务区域顶点为D1，与L2对应的任务区域顶点为D2；

3.分别连接L1和D1，L2和D2，得到与步骤1中星下线平行线的交点，构成四边形的4个交点即组成一个任务条带，由此对任务区域进行划分；

4.选择能够完全覆盖任务区域的一个或多个任务条带，作为任务规划与调度的基本单元。

二、计算卫星的位置及速度

采用对轨道动力学方程数值积分求解的方法预报限定时间段内卫星在J2000惯性坐标系下的轨道位置和速度。

根据卫星的轨道根数，能够推算出任务规划初始时刻J2000惯性坐标系下的轨道位置R_sat和速度V_sat，再采用Cowell方法求解轨道动力学方程(选用高斯型摄动运动方程)，得到限定时间段内卫星在J2000惯性坐标系下的轨道位置R_sat和速度V_sat。高斯型摄动运动方程及Cowell方法在国防工业出版社出版的《航天器轨道理论》(刘林著，2000年)一书中有详细的说明。J2000惯性坐标系定义见参考文献“地球卫星运动中坐标系附加摄动与参考系选择问题”(《空间科学学报》2008年第28卷第2期，作者刘林、汤靖师)。

三、计算卫星对每个任务条带的可见时间区间

对每个条带，根据上述第二部分得到的轨道位置R_sat和速度V_sat，计算各离散时刻点卫星指向各条带顶点的姿态角，再根据卫星姿态机动范围对各离散时刻点进行遍历，得到各条带的可见时间窗口，最后根据目标太阳高度角的计算结果，剔除在地影区的观测窗口。下面以一个条带的计算为例进行说明。

1.根据J2000坐标系下的轨道位置R_sat和速度V_sat，计算限定时间段内各时刻点卫星指向条带各顶点的姿态角。下面仅以一个点的计算为例进行说明。

已知卫星的轨道位置R_sat、速度V_sat，地面目标点的大地经纬度

及协调世界时UTC时间t。首先根据目标点的大地经纬度，计算出t时刻目标点在J2000惯性坐标系下的位置矢量R_T，f(t)，然后根据R_T，f(t)与卫星t时刻的位置矢量R_sat，得到t时刻卫星指向该目标点的姿态角。具体步骤如下：

将地面目标点大地经纬度转化为地心经纬度

计算公式为：

λ_c＝λ_d，

其中

表示地球扁率，然后计算目标点地心距：

R_e＝6378.140km，为地球赤道半径。

根据UTC时间计算地固坐标系到J2000惯性坐标系的转换矩阵R_if(t)，计算方法在国防工业出版社出版的《航天器轨道理论》(刘林著，2000年)中有详细描述。通过坐标变换，得到目标点在J2000惯性坐标系下的位置矢量：

R_x(α)、R_y(α)、R_z(α)分别表示绕x、y、z轴旋转的基元变换矩阵：

$R_x\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& \sin \mathrm{\α}\\ 0& -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right]$

$R_y\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& 0& -\sin \mathrm{\α}\\ 0& 1& 0\\ \sin \mathrm{\α}& 0& \cos \mathrm{\α}\end{array}\right]$

$R_z\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 1\end{array}\right]$

然后计算J2000惯性坐标系下卫星指向地面目标点的矢量：

R_f(t)＝R_T，f(t)-R_sat

将矢量R_f(t)由J2000惯性坐标系变换到卫星轨道坐标系：

$R_o\left(t\right)=R_{\mathrm{oi}}{R}_f\left(t\right)=\left(\begin{array}{c}{v}_x\left(t\right)\\ v_y\left(t\right)\\ v_z\left(t\right)\end{array}\right)$

其中，R_oi表示J2000惯性坐标系到卫星轨道坐标系的转换矩阵。上述地固坐标系定义、卫星轨道坐标系定义及J2000惯性坐标系到卫星轨道坐标系的变换推导见北京航空航天大学出版社出版的《卫星轨道姿态动力学与控制》(章仁为编著，1998年)。

偏航角为0时，根据姿态欧拉角之间几何关系，得到卫星对目标的观测姿态角

(转序为312)：

$\left[\begin{array}{c}\mathrm{yaw}\left(t\right)\\ \mathrm{roll}\left(t\right)\\ \mathrm{pitch}\left(t\right)\end{array}\right]=\left[\begin{array}{c}0\\ \arcsin \left(\frac{v_x\left(t\right)}{\left|R_o\left(t\right)\right|}\right)\\ -\arctan \left(\frac{v_y\left(t\right)}{v_z\left(t\right)}\right)\end{array}\right]$

上式中yaw(t)、roll(t)和pitch(t)分别表示卫星指向目标点的姿态角中时间t与偏航角、滚动角和俯仰角的对应关系。

2.根据卫星姿态机动范围和步骤1得到的各时刻点卫星指向条带各顶点的姿态角，计算卫星对条带各顶点的可见时间区间[t₁，t₂]_k，k＝1，2，3，4。

卫星对目标的观测受限于卫星的姿态机动能力，因此只有当卫星指向目标的姿态在卫星的姿态机动范围内，才能执行观测任务。对于第k个顶点，

T表示限定的时间范围，如果对应姿态角[yaw roll pitch]满足|roll|≤roll_max，|pitch|≤pitch_max，其中roll_max、pitch_max表示卫星最大滚动角和最大俯仰角，则t∈[t₁，t₂]_k，即[t₁，t₂]_k为满足上述条件的t的集合。计算流程如图2所示。

3.根据卫星对条带各顶点的可见时间区间[t₁，t₂]_k，k＝1，2，3，4，计算卫星对条带的可见时间区间[t₀，t_n]。

对卫星对各顶点的可见时间区间[t₁，t₂]_k，k＝1，2，3，4求交集，即得到卫星对条带的可见时间区间[t₀，t_n]。

4.根据地面太阳高度角计算结果，剔除在地影区的时间区间。

通常情况下步骤3计算得到的可见区间[t₀，t_n]存在多个解，其中部分可见区间在地影区内，无法满足光学相机的成像条件，应予以剔除。

对步骤3得到的可见区间[t₀，t_n]，从其中任选一点，计算该时刻点观测目标的太阳高度角ε，如果ε＜0，则该点在地影区内，予以剔除，否则保留。

太阳高度角与时间有关，计算过程如下：

(1)根据公历时间(格里高利历)计算儒略日和儒略世纪数

设给出公历时间的年、月、日(含日的小数部分)分别为Y、M、D，则对应的儒略日为：

$\mathrm{JD}=D-32075+[1461\×(Y+4800+[\frac{M-14}{12}\left]\right)\÷4]$

$+[367\×(M-2-[\frac{M-14}{12}]\×12)\÷12]$

$-[3\×[(Y+4900+[\frac{M-14}{12}\left]\right)\÷100]\÷4]-0.5$

式中[X]表示取X的整数部分，小数部分省略。

从1950年1月1日12时起算的约简儒略日为MJD＝JD-2433283修正儒略世纪数为MJC＝MJD/36525

(2)计算太阳相对于地球的轨道参数的变化

长半轴a_h＝1.49597927×10⁸；

偏心率e_h＝1.67301085×10^-2-4.1926×10^-5MJC-1.26×10^-7MJC²；

轨道倾角(即黄赤交角)

i_h＝23.4457888616-1.30141669×10^-2MJC-9.445×10^-7MJC²+5.000×10^-7MJC³；近地点幅角ω_h＝1.67301085×10^-2-4.1926×10^-5MJC-1.26×10^-7MJC²；

平近点角

M_h＝358.000682-0.9856002623MJD-1.550000×10^-4MJC²-3.3333×10^-6MJC³；

升交点赤经Ω_h＝0；

(3)采用迭代法求解方程E_h-e_hsinE_h＝M_h，得到偏近点角E_h，迭代公式为E_h(k+1)＝e_hsinE_h(k)+M_h或

初值E_h(0)＝M_h；

(4)由E_h和e_h求真近点角θ_h

$\mathrm{tg}\frac{{\mathrm{\θ}}_h}{2}=\sqrt{\frac{1-e_h}{1+e_h}}\mathrm{tg}\frac{E_h}{2};$

(5)根据公式u_h＝ω_h+θ_h计算纬度幅角(即黄经)；

(6)根据公式tgα_h＝cosi_htgu_h计算太阳赤经α_h；

(7)根据公式sinδ_h＝sini_hsinu_h计算太阳赤纬；

(8)根据公式α_G＝α_G0+6.3003881t计算格林成治赤经；

(9)α_G0和α_G以弧度(rad)为单位，t为从该年初始时刻(按世界时)起算的时间，以平太阳日为单位。α_G0值从当年天文历书中查到；

(10)计算地面目标太阳高度角

λ、

为观测目标经纬度。

上述方法的推导详见宇航出版社出版的《航天器飞行动力学原理》(肖业伦编著，1995年)。

四、构造成像质量确定函数

成像质量主要由推扫和线阵两个方向的几何分辨率、调制传递函数MTF、信噪比、幅宽等参数决定。其中MTF主要和相机设计制造等静态因素及卫星姿态稳定度、颤振等有关，因此对于观测时间点的选择没有实际的影响，在相机成像参数设置保持不变时为常值。信噪比主要由地面目标反射率和太阳高度角决定，对于指定区域，反射率是一定的，太阳高度角越大，信噪比越高；在相机设计参数与卫星轨道高度确定的前提下，幅宽与观测滚动角和俯仰角有关，滚动角和俯仰角越大，幅宽越大，由于大幅宽与高分辨率存在矛盾，因此幅宽不作为成像质量计算的主要指标。推扫分辨率主要受俯仰角影响，俯仰角越大，推扫分辨率越差。线阵分辨率主要受滚动角影响，滚动角越大，线阵分辨率越差。因此影响成像质量的动态因素主要包括滚动角、俯仰角和太阳高度角。根据第三部分的内容，在卫星、目标确定的前提下，观测时间是滚动角、俯仰角和太阳高度角的唯一决定因素，因此成像质量可表示为时间的一元函数。

本发明中，成像质量综合评定由MTF、推扫方向(GSDx)与线阵方向(GSDy)的几何分辨率、信噪比(snr)和幅宽(breath)加权求和得到。首先根据成像条件分别计算各项指标，然后将各项指标的计算结果转化为单项指标的成像质量，再将各项指标。加权，得到综合评定结果。计算过程如下：

1.计算几何分辨率

已知卫星轨道高度H、姿态角[yaw roll pitch]、相机像元尺寸d、焦距f、视场角FOV、地球半径Re，卫星推扫方向几何分辨率GSD_x和线阵方向几何分辨率GSD_y为

${\mathrm{GSD}}_x=H\·\sec \mathrm{\β}\·\frac{d}{f}{\sec }^2\mathrm{\γ}$

${\mathrm{GSD}}_y=H\·\mathrm{se}{c}^2\mathrm{\β}\·\frac{d}{f}\sec \mathrm{\γ}$

其中β＝roll，

表示将相机摄影点投影到推扫方向和线阵方向后的等效侧摆角和等效俯仰角。

2.计算成像幅宽

幅宽计算方法如图3所示，图中O表示地球球心，S表示卫星，S′表示星下点，+X表示卫星的前进方向，即相机推扫方向，+Y表示相机线阵方向，假设地球为理想球体，XY平面过星下点与地球球面相切。A、B两点表示相机视场角边缘与XY平面的交点，A′、B′表示相机视场角边缘与地球表面交点，则A′、B′两点间的球面距离为卫星姿态机动后的幅宽。F点为相机光轴与XY平面的交点，F′为相机光轴与地球表面的交点。

记卫星光轴指向角为α，A、B点的姿态机动指向角分别为α₁、α₂，根据直角三角形的性质，有

$\mathrm{\α}=\arctan \left(\frac{\sqrt{{\sin }^2\mathrm{roll}+{\tan }^2\mathrm{pitch}}}{\cos \mathrm{roll}}\right)$

${\mathrm{\α}}_1=\arctan \left(\frac{\sqrt{{\sin }^2(\mathrm{roll}-\frac{\mathrm{FOV}}{2})+{\tan }^2\mathrm{pitch}}}{\cos (\mathrm{rpll}-\frac{\mathrm{FOV}}{2})}\right)$

${\mathrm{\α}}_2=\arctan \left(\frac{\sqrt{{\sin }^2(\mathrm{roll}+\frac{\mathrm{FOV}}{2})+{\tan }^2\mathrm{pitch}}}{\cos (\mathrm{rpll}+\frac{\mathrm{FOV}}{2})}\right)$

根据三角运算关系，在三角形SOF′中可以得到：

于是有

${\mathrm{SF}}^{\′}=(R_e+H)\cos \mathrm{\α}-\sqrt{{R_e}^2-{(R_e+H)}^2{\sin }^2\mathrm{\α}}$

同理可得

${\mathrm{SA}}^{\′}=(R_e+H)\cos {\mathrm{\α}}_1-\sqrt{{R_e}^2-{(R_e+H)}^2{\sin }^2{\mathrm{\α}}_1}$

${\mathrm{SB}}^{\′}=(R_e+H)\cos {\mathrm{\α}}_2-\sqrt{{R_e}^2-{(R_e+H)}^2{\sin }^2{\mathrm{\α}}_2}$

由图3可知∠A′SB′＝∠ASB＝FOV，在三角形A′SB′中采用余弦定理计算线段A′B′的长度：

$A^{\′}{B}^{\′}=\sqrt{{\mathrm{SA}}^{\′2}+{\mathrm{SB}}^{\′2}-{2\mathrm{SA}}^{\′}\·{\mathrm{SB}}^{\′}\·\cos \mathrm{FOV}}$

由线段A′B′的长度计算对应的地心角

$\mathrm{\φ}=2\arcsin \left(\frac{A^{\′}{B}^{\′}}{{2R}_e}\right)$

幅宽为A′、B′两点间的球面距离L_A′B′＝φ·R_e。

3.计算信噪比

当CCD相机TDI积分级数选定后，信噪比由太阳高度角和地表反射率决定。信噪比计算采用构造不同条件下信噪比数表，再对数表进行二元插值的方法。信噪比数表的构造采用6S(Second Simulation of the Satellite Signal in theSolar Spectrum)软件计算得到，孙吉福、吴雁林在2003年第24卷第2期的《航天返回与遥感》上发表的“HJ-1A卫星CCD相机入射光谱辐射亮度的分析”对该方法有比较具体的介绍。通过6S软件，得到太阳高度角为10°、20°、30°、40°、50°、60°、70°、80°、90°，地表反射率为7％、10％、15％、20％、26％、40％、50％、60％、70％、80％、90％、100％时分别对应的信噪比，在此基础上采用二元函数插值计算得到任意太阳高度角和任意地表反射率对应的信噪比，二元函数插值在北京航空航天大学出版社出版的《数值分析》(颜庆津编著，2000年)一书中有详细的说明，方法如下：

已知信噪比数表snr(r_i，ε_j)，i＝0，1，…，11，r_i＝0.07，0.15，…，1.00，表示地表反射率，j＝0，1，…，8，ε_j＝10，20，…，90，表示太阳高度角。若任意地表反射率和太阳高度角(r，ε)满足：

r_i-1≤r＜r_i，ε_j-1≤ε＜ε_j

则对应的信噪比计算公式为：

$\mathrm{snr}(r,\mathrm{\ϵ})=\frac{(r-r_i)(\mathrm{\ϵ}-{\mathrm{\ϵ}}_j)}{(r_{i-1}-r_i)({\mathrm{\ϵ}}_{j-1}-{\mathrm{\ϵ}}_j)}\mathrm{snr}(r_{i-1},{\mathrm{\ϵ}}_{j-1})+\frac{(r-r_i)(\mathrm{\ϵ}-{\mathrm{\ϵ}}_{j-1})}{(r_{i-1}-r_i)({\mathrm{\ϵ}}_j-{\mathrm{\ϵ}}_{j-1})}\mathrm{snr}(r_{i-1},{\mathrm{\ϵ}}_j)$

$+\frac{(r-r_{i-1})(\mathrm{\ϵ}-{\mathrm{\ϵ}}_j)}{(r_i-r_{i-1})({\mathrm{\ϵ}}_{j-1}-{\mathrm{\ϵ}}_j)}\mathrm{snr}(r_i,{\mathrm{\ϵ}}_{j-1})+\frac{(r-r_{i-1})(\mathrm{\ϵ}-{\mathrm{\ϵ}}_{j-1})}{(r_i-r_{i-1})({\mathrm{\ϵ}}_j-{\mathrm{\ϵ}}_{j-1})}\mathrm{snr}(r_i,{\mathrm{\ϵ}}_j)$

上述算法中地表反射率由地貌特征决定，对于特定目标，地表反射率为常值。

4.计算各项指标单项质量等级

根据各项指标计算单项质量等级，其中MTF、几何分辨率和幅宽采用线性内插的方法计算，信噪比以标准成像条件(太阳高度角70°，地表反射率26％，85分)为参考计算。

$R_{\mathrm{MTF}}=\frac{\mathrm{MTF}-{\mathrm{MTF}}_{\mathrm{worst}}}{{\mathrm{MTF}}_{\mathrm{best}}-{\mathrm{MTF}}_{\mathrm{worst}}}\·100$

其中，MTF_worst、MTF_best分别表示最小、最大MTF。

$R_{\mathrm{GSD}}=\frac{\mathrm{GSD}-{\mathrm{GSD}}_{\mathrm{worst}}}{{\mathrm{GSD}}_{\mathrm{best}}-{\mathrm{GSD}}_{\mathrm{worst}}}\·100$

GSD_worst、GSD_best分别表示最差、最好分辨率。

$R_{\mathrm{breath}}=\frac{L-L_{\mathrm{worst}}}{L_{\mathrm{best}}-L_{\mathrm{worst}}}\·100$

L_worst、L_best分别表示最小、最大幅宽。

$R_{\mathrm{snr}}=\left\{\begin{array}{cc}85+\frac{\mathrm{snr}(r,\mathrm{\ϵ})-\mathrm{snr}\left(\mathrm{0.26,70}\right)}{\mathrm{snr}\left(\mathrm{1.00,90}\right)-\mathrm{snr}\left(\mathrm{0.07,10}\right)}\·15,& \mathrm{snr}(r,\mathrm{\ϵ})>85\\ \frac{\mathrm{snr}(r,\mathrm{\ϵ})-\mathrm{snr}\left(\mathrm{0.07,10}\right)}{\mathrm{snr}\left(\mathrm{1.00,90}\right)-\mathrm{snr}\left(\mathrm{0.07,10}\right)}\·15,& \mathrm{snr}(r,\mathrm{\ϵ})\≤85\end{array}\right.$

5.计算成像质量综合等级

成像质量等级R＝R(t)表达式如下：

R(t)＝∑η_iR_i(t)＝η_MTFR_MTF+η_GSDxR_GSDx(t)+η_GSDyR_GSDy(t)+η_snrR_snr(t)+η_breathR_breath(t)η_MTF、η_GSDx、η_GSDy、η_snr和η_breath分别表示单项指标的权重，单项指标的权重，满足∑η＝1。根据用户需求不同，各项指标的权重可以调整。综合评分为各项指标分数的加权和。除MTF为常值外，其它指标评分均为时间的函数。

五、判断时间窗口的存在性

采用第四部分得到的一元函数表达式R＝R(t)分别计算在第三部分中确定的可见时间区间[t₀，t_n]的两个端点t₀和t_n以及区间中间点的图像质量R₀＝R(t₀)，R_n＝R(t_n)，

记最低许可成像质量为R_p，将R_p与R₀，R_n，R_max进行比较，确定任务的满足成像质量要求的可见时间区间[t₁，t₂]是否存在。

光学卫星成像时，通常采用星下点成像质量最好，因此R_mid约等于该可见时间区间的最高成像质量等级，在

区间，R＝R(t)可近似为单调递增函数，在

区间，R＝R(t)可近似为单调递减函数。

(1)如果R_mid＜R_p，由函数的单调性可知，则时间区间[t₁，t₂]不存在，将本目标标记为“不满足成像质量要求”；

(2)如果R_mid＞R_p，R₀≥R_p，R_n≥R_p，由函数的单调性可知，可见时间区间完全满足成像质量要求，t₁＝t₀，t₂＝t_n；

(3)如果R_mid＞R_p，R₀＜R_p，R_n≥R_p，由连续函数介值定理和函数的单调性可知，

t₂＝t_n，采用弦截法求t₁；

(4)如果R_mid＞R_p，R₀≥R_p，R_n＜R_p，由连续函数介值定理和函数的单调性可知，t₁＝t₀，

采用弦截法求t₂；

(5)如果R_mid＞R_p，R₀＜R_p，R_n＜R_p，由连续函数介值定理，

采用弦截法求t₁，t₂。

六、采用弦截法裁剪时间窗口

根据最低许可成像质量R_p，在满足成像质量要求的可实传区间[t₁，t₂]内采用弦截法求解方程R_p＝R(t_p)，得到第五部分步骤(3)、(5)中的时间窗口起始点t₁，步骤(4)、(5)中的时间窗口结束点t₂。

弦截法(也称割线法)是求解复杂非线性方程常用的一种数值解法，其优点在于收敛速度较快，收敛速度的阶至少为1.618，且避免了牛顿法需计算函数导数的不足。弦截法在北京航空航天大学出版社出版的《数值分析》(颜庆津编著，2000年)一书中有详细的说明。

弦截法求解方程R_p＝R(t_p)的步骤如下：

当时，且满足R_p＝R(t_p)。令x_-1＝t₀，对于k＝0，1，…，M执行

(1)计算 $x_{k+1}=x_k-\frac{R\left(x_k\right)\·(x_k-x_{k-1})}{R\left(x_k\right)-R\left(x_{k-1}\right)};$

(2)若|x_k+1-x_k|＜σ，取t_p≈x_k+1，得到时间窗口起点，否则转(1)。其中M表示最大迭代次数，σ为许可的计算误差。

令

x₀＝t_n，同法可得到R_p＝R(t_p)在区间上的另一个根，即时间窗口结束点。

对于立体成像任务，根据任务的观测窗口时间区间[t₁，t₂]、用户要求的成像俯仰角roll_i(i表示要求的立体成像次数)，以及第三部分步骤2计算得到的各时间点指向目标的姿态角，判断是否能实现立体成像。根据卫星向目标运动时俯仰角的单调性，如果roll₁≤roll_i≤roll₂，其中roll₁、roll₂分别是t₁、t₂时刻的俯仰角，则能按要求的角度完成立体成像任务，否则，将本目标标记为“不满足立体成像要求”。

七、计算每个任务条带的观测持续时间

条带的观测持续时间，即从开始观测到结束观测的持续时间，取决于条带的长短。对一个多条带任务，由于所有条带的长度相等，因此只需任选一个条带计算，即能得到所有条带的观测持续时间。对每个条带，首先根据第三部分的步骤1得到的卫星指向各顶点的姿态角，再采用弦截法计算条带的开始观测时间和结束观测时间，得到条带的观测持续时间，最后判断条带的可见时间窗口是否能包含条带的观测持续时间，如果不能，则该任务不能执行。下面以一个条带为例进行说明：

1.根据第三部分的步骤1得到的各时刻点卫星指向条带各顶点的姿态角，在第三部分的步骤2得到的[t₁，t₂]_k，k＝1，2，3，4内，采用弦截法计算俯仰角为0时，卫星指向条带各顶点的时间。

下面以顶点k为例，说明采用弦截法求解方程Pitch(t_pk)＝0的步骤：

卫星对目标的可见时间由卫星的姿态机动能力决定，卫星在轨运动过程中，对固定目标的观测俯仰角由正的最大值变化为0，再由0变为负的最大值，也即在顶点k的可见时间区间[t₁，t₂]_k上，卫星在第一点t₁的俯仰角为roll_max，在最后一点t₂的俯仰角为-roll_max，因此区间[t₁，t₂]_k内必有卫星俯仰角为0的时刻点，即当Pitch(t₁)·Pitch(t₂)＜0时，

且满足Pitch(t_pk)＝0。

令x_-1＝t₁，x₀＝t_pk，对于i＝0，1，…，M执行

(1)计算 $x_{i+1}=x_i-\frac{\mathrm{Pitch}\left(x_i\right)\·(x_i-x_{i-1})}{\mathrm{Pitch}\left(x_i\right)-\mathrm{Pitch}\left(x_{i-1}\right)};$

(2)若|x_i+1-x_i|＜σ，取t_pk≈x_i+1，得到俯仰角为0时，卫星指向目标顶点的时间，否则转(1)。

其中M表示最大迭代次数，σ为许可的计算误差。

2.由于卫星对条带的观测必定开始于某个顶点，也结束于某个顶点，因此根据卫星对各顶点的可见时间，能够确定卫星对条带的开始观测时间和结束观测时间。

将t_pk按照先后顺序排序，记最早的时间为T_start，最晚的时间为T_end，分别表示条带的开始观测时间以及结束观测时间。

3.根据下式计算条带的持续观测时间

T_last＝T_end-T_start

T_last是根据条带各顶点的可见时间区间计算得到的，没有考虑地影的影响，因此还需要判断剔除了在地影区的时间区间后，第三部分得到的条带的可见区间[t₀，t_n]是否能够完成对条带的观测：

如果T_last＞t_n-t₀，本条带的可见时间窗口比观测持续时间短，整个观测任务无法完成。

八、输出所有目标的时间窗口信息，对于被标记的任务，输出不能完成的原因，将未被标记的任务提交给任务规划系统，进行任务规划与调度。

实施例

考虑一颗运行于太阳同步圆轨道的快速姿态机动成像卫星，历元时刻2009年7月26日00:00:00.000UTC瞬根数为半长轴7051.2km，轨道倾角97.3087°，升交点赤经249.784°，纬度辐角0°。卫星姿态机动范围为俯仰方向和滚动方向±45°。目标的地理经纬度信息如表1所示，每个目标的四个顶点的经度和纬度按次序一一对应。

表1目标地理经纬度

本例共包含6个观测任务，任务1、2、5、6是常规任务，最低许可成像质量分别为55、55、65和55。任务3和任务4是立体成像任务，任务3要求2次立体成像，俯仰角分别为-30°和30°，最低许可成像质量55；任务4要求3次立体成像，俯仰角分别为-30°，0°和30°，最低许可成像质量70。

卫星光学相机采用全色和多光谱两种CCD，全色像元尺寸10μm，多光谱像元尺寸40μm，相机焦距10m，视场角1.07°，对应的星下点幅宽为13km，MTF最大为0.0902，最小0.0818。全色CCD TDI积分级数设置为24时，MTF为0.0895。地表反射率50％。信噪比数表如下：

表2相机信噪比数表

对这批任务的预处理步骤如下：

(1)任务区域条带划分

根据相机幅宽、任务区域的经纬度及卫星星下点，对任务区域进行条带划分，划分结果如表3所示。

表3任务区域条带划分结果

(2)计算可见时间窗口

根据卫星轨道根数，推算2009年7月26日00:00:00.000UTC至2009年7月27日00:00:00.000UTC时间段内卫星在J2000坐标系下的轨道位置、速度，然后计算各整秒时刻卫星对各任务的观测角度。根据卫星姿态机动范围，得到卫星与任务的可见窗口。

任务1和任务2的观测角度超出了卫星的最大姿态机动范围，因此任务1和任务2被标记为“无可见时间窗口”，其他任务的计算结果如表4所示。

表4任务的可见时间窗口

任务序号	可见开始时刻t0	可见结束时刻tn
			3	02:15:33	02:19:24
4	00:47:34	00:50:56
			5	12:27:25	12:30:53
6	00:37:16	00:39:45

除了任务1和任务2“无可见时间窗口”之外，其余任务均处于阳照区。

(3)确定成像质量的主要影响因素

由于MTF与观测时间点的选择无关，因此不是影响成像质量的主要因素。

在任务3的可见时间窗口02:15:33至02:19:24之间，俯仰角和滚动角的计算结果如图4、图5所示，俯仰角变化范围为±45°，滚动角变化范围为20.32°～25.73°，因此俯仰角变化是影响成像质量的主要因素。信噪比由太阳高度角决定，太阳高度角计算结果如图6所示，在整个可见窗口内，太阳高度角的变化范围不超过1°，因此不是影响成像质量的主要因素。

任务4～任务6的分析结果与任务3相同。

(4)计算成像质量

由于观测姿态角和太阳高度角均可表示为时间的一元函数，因此成像质量也可表示为时间的一元函数R(t)。

由于大幅宽和高分辨率存在矛盾，通过限制姿态角保证高分辨率时，必然使幅宽减小。通常成像质量确定方法将高分辨率作为一项更重要的指标。本实施例中将幅宽在成像质量中的权重设为0，即不予考虑。

各项指标的权重分配如下：

指标	权重
		MTF	0.125
信噪比	0.375
		全色推扫方向分辨率	0.125
全色线阵方向分辨率	0.125
		多光谱推扫方向分辨率	0.125
多光谱线阵方向分辨率	0.125
		成像幅宽	0

即

R(t)＝0.125R_MTF+0.125R_全色GSDx(t)+0.125R_全色GSDy(t)

+0.125R_{多光谱GSDx}(t)+0.125R_{多光谱GSDy}(t)+0.375R_snr(t)

由于任务数量众多，仅以任务3的综合计算结果为代表说明。

根据各可见时刻点各项指标的单项质量等级，得到任务3的成像质量综合等级，图7为任务3在可见时间窗口内成像质量随观测时间的变化曲线。

(5)判断时间窗口[t₁，t₂]的存在性

由于任务数量众多，仅以任务3的计算过程为代表说明。

对任务3，以t₀＝2009-7-2600:00:00为参考时间，由t₀＝8133秒(即t₀＝2009-07-2602:15:33)，t_n＝8364秒(即t_n＝2009-07-2602:19:24)，

得到R(t₀)＝52.85，

R(t_n)＝50.08。

由于R_mid＞R_p，R₀＜R_p，R_n＜R_p，由连续函数介值定理，t₁∈(t₀，t_m)，t₂∈(t_m，t_n)，需要采用弦截法求t₁，t₂。

其他任务的计算结果如下：

任务4：

R(t₀)＝56.32，

R(t_n)＝55.16，由于R_mid＞R_p，R₀＜R_p，R_n＜R_p，由连续函数介值定理，t₁∈(t₀，t_m)，t₂∈(t_m，t_n)，需要采用弦截法求t₁，t₂。

任务5：

R(t₀)＝27.93，

R(t_n)＝28.38，由于R_mid＜R_p，由函数的单调性可知，满足成像质量要求的时间区间[t₁，t₂]不存在，本任务被标记为“不满足成像质量要求”。

任务6：

R(t₀)＝52.49，

R(t_n)＝52.35，由于R_mid＞R_p，R₀＜R_p，R_n＜R_p，由连续函数介值定理，t₁∈(t₀，t_m)，t₂∈(t_m，t_n)，需要采用弦截法求t₁，t₂。

(6)采用弦截法裁减时间窗口[t₁，t₂]

对于可见时间窗口内存在满足成像质量要求时段的任务3、4、6，采用弦截法分别计算t₁，t₂。计算结果如表5所示。

表5满足成像质量的可见区间

任务编号	可见开始时刻t1	可见结束时刻t2
			3	02:15:48.44	02:18:52.95
4	00:48:21.69	00:50:05.77
			6	00:37:31.07	00:39:28.15

对于立体成像根据要求的成像角度和在满足成像质量的可见时间窗口内指向目标的姿态角，判断是否能实现立体成像。对于任务3，卫星在[t₁，t₂]上的俯仰角变化范围为40.88°～-37.36°，能满足要求的30°和-30°的立体成像角度；对于任务4，卫星在[t₁，t₂]上的俯仰角变化范围为27.98°～-26.48°，不能满足30°、0°和-30°的立体成像角度，将任务4标记为“不满足立体成像要求”。

(7)确定观测持续时间

根据计算，任务3的观测持续时间为T_last＝1.93秒，任务6的观测持续时间为T_last＝2.31秒，均远远小于各自满足成像质量要求的可见窗口长度，因此不存在“不满足观测持续时间”的任务。

(8)将预处理结果(如表6所示)输出给任务规划系统，进行任务规划与调度。

表6任务的预处理结果

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

Claims (1)

1.一种基于成像质量预估的任务预处理方法，其特征在于步骤如下：

(1)根据相机幅宽将任务区域条带划分；

(2)根据卫星的最大俯仰角和最大滚动角以及任务条带信息确定卫星对各条带的可见时间区间[t₀，t_n]，剔除无可见时间区间及时间区间在地影区的任务；每个条带的可见时间区间[t₀,t_n]由卫星对该条带四个顶点的可见时间区间[t₁，t₂]_k，k＝1，2，3，4求交集并剔除在地影区的时间区间后得到；卫星对条带单个顶点的可见时间区间由卫星的最大俯仰角和最大滚动角确定；

(3)将调制传递函数MTF、推扫方向的几何分辨率、线阵方向的几何分辨率、信噪比和幅宽作为影响图像质量的主要因素，在步骤(2)确定的可见时间区间[t₀，t_n]内将图像质量表示为时间的一元函数R(t)，

R(t)＝η_MTFR_MTF+η_GSDxR_GSDx(t)+η_GSDyR_GSDy(t)+η_snrR_snr(t)+η_breathR_breath(t)，

式中η_MTF+η_GSDx+η_GSDy+η_snr+η_breath＝1，R_MTF，R_GSDx(t)，R_GSDy(t)，R_snr(t)，R_breath(t)分别为调制传递函数MTF、推扫方向的几何分辨率、线阵方向的几何分辨率、信噪比和幅宽与时间的一元函数关系；

(4)采用步骤(3)得到的一元函数表达式分别计算在步骤(2)中确定的可见时间区间[t₀，t_n]的两个端点t₀和t_n以及区间中间点

的图像质量R₀＝R(t₀)，R_n＝R(t_n)， $R_{\mathrm{mid}}=R\left(\frac{t_0+t_n}{2}\right);$

(5)记最低许可成像质量为R_p，将R_p与步骤(4)得到的R₀，R_n，R_mid进行比较，确定满足成像质量要求的可见时间区间，剔除不满足成像质量要求的任务；

(6)采用弦截法计算每个任务条带的观测持续时间T_last，T_last为从该条带最先观测的顶点开始到该条带最后观测的顶点结束所持续的时间；当任务条带的观测持续时间T_last小于步骤(5)的可见时间区间长度时，该任务为可执行任务，将所有可执行任务提交给任务规划系统进行任务规划与调度。

Images