# Automatic Task Planning and Its On-Orbit Verification of Agile Remote Sensing Satellite

Authors: Zhang Yahang; Yang Haiyue; Yang Mengfei; Yang Ruochu; Leng Shuhang
Journal: International Journal of Aerospace Engineering
 (2023)

---

## Abstract

An on-board autonomous task planning system is designed and implemented in this article, aiming at the problem that the current remote sensing satellite needs complex instruction support to perform tasks and depends on the ground system too much. The complex earth observation task description and injection decomposition modules are designed in the system. No more than 128 ordered points describe the curve area task or irregular polygon area task, and the complex task is decomposed into several strips according to the satellite imaging width. Then, the task’s orbit, attitude, mobility, energy, and time constraints are calculated through the modules of on-board orbit prediction, agile attitude calculation, track planning, energy prediction, and visible arc calculation. Finally, the on-board autonomous task planning and execution are completed through task solution space search and metatask command generation modules. The on-orbit flight verification is carried out on the high-resolution multimode (GFDM) satellite. The results show that this paper’s on-board autonomous task planning system can complete complex task injection and autonomous planning and finally execute.


---

## Body

## 1. Introduction

Space-based remote sensing systems play a critical role in scenarios of flood and forest disaster prevention, land survey, and other emerging fields. Ground operation and management are becoming more complex due to the increasing trends in the complexity of observation requirements and the number of satellites in constellation or swarm [1–4]. Hence, the traditional task planning mode is hard to apply [5, 6]. In addition, current satellite management and control procedure relies on the ground telemetry and telecontrol (TM & TC) station, which leads to a longer time duration from the input of task requirement to remote sensing product generation. In such a centralized mode, the task planning system has flaws in image information timeliness, emergency task response, multisatellite collaborative scheduling, etc. Therefore, it is urgent to develop on-board intelligence and autonomy, including on-board task planning, real-time evaluation, decision-making, multisatellite, multipayload cooperation, and other technical research.Most current task planning systems are applied in the ground computing environment; they do not consider the computing and storage resources of on-board computers [7–14]. Correspondingly, most currently developed satellites have minimal computing ability and memory space, which is impossible to plan complex tasks on orbit. So, the algorithms have to be redesigned to fit the ability of on-board computers. To overcome shortage of computational resource, Zhang et al. designed a satellite and ground station hybrid method to plan remote sensing tasks [15]. An autonomous task plan framework that could run on-board was proposed to solve agile satellite imaging task planning problems [16, 17]. On JWST satellite, the event-driven system was established to deal on-board tasks [18]. Here, high-resolution multimode satellite (GFDM) is taken as an attempt to experiment with on-orbit task and attitude planning, its processing capacity is no more than 50 MHz, and SRAM is no more than 6 MByte. To suit the GFDM satellite, the on-board task planning system was modified to ensure computational efficiency in limited on-board resources and complete planning in small storage space.This paper designed an on-board automatic task planning system and conducted a set of on-orbit flight verification, aiming at the new requirements of task planning for agile remote sensing satellites. The second section will give the overall design, planning algorithm design, and processing flow design of the on-board task planning system. The third section will detail the on-orbit flight verification based on the GFDM satellite. Finally, the conclusion and prospect of this paper are given in the Section4.


## 2. Design of On-Board Task Planning System

The existing agile satellite takes a single strip or a spot area as ametatask to complete remote sensing [19]; here, the metatask refers to the task that can be imaged or executed by switching on and off once by the payload [20]. The on-board task planning system designed in this paper is compatible with these metatasks; the architecture is shown in Figure 1. Except for metatasks, for a better description, the complex task is defined in this paper as the task that can only be imaged entirely requiring several contiguous strips, such as observation tasks of curve strip or polygonal area. The on-board task planning system of agile remote sensing satellites considers the satellite orbit, imaging width, and agile ability. The metatask is the basic observation task executed by the on-board task planning system, and the complex task needs to be decomposed into multiple metatasks. The task planning system will decompose the curve strips or irregular polygons into a series of strips or spots according to the comprehensive constraints of satellite imaging width, subsatellite tack, memory size, power, and acceleration of angle velocity.Figure 1
On-board task planning architecture.


### 2.1. Overall Design

The ground operation and control center receives users’ observation requirements and injects multiple complex tasks or metatasks into satellite. The on-board task planning system autonomously completes task planning, instruction generation, and other operations. The system architecture is shown in Figure1. After receiving and decomposing the complex tasks, the satellite implements task planning according to the real-time on-orbit state telemetry. A task planning process is divided into 8 subprocedures, that is, from requirement gathering to instruction generating, including the following:Step 1.
Gather requirements. The users send tasks to the ground station by means of polygons, curve strips, gaze, and tracking tasks and then forward to satellites through the satellite-ground link.Step 2.
Examine the validity of the task. Check task time, attributes, and types to determine whether the task is valid or not.Step 3.
Decompose complex tasks. The on-board task planning system decomposes and processes complex tasks according to the subsatellite track vector, generating the metatask.Step 4.
Planning execution window. Select imaging orbital arc.Step 5.
Calculate external constraints. Calculate the constraints of power [21], storage, attitude, and so on.Step 6.
Plan the attitude and attitude maneuver trajectory.Step 7.
Generate metatasks.Step 8.
Generate instructions.


### 2.2. Complex Task Description and Decomposition

#### 2.2.1. Task Description

Complex tasks include irregular curve strips, irregular areas, point group target, etc. Tasks can be described by the following essential attributes:(1)
Task ID: the unique ID of the task, which is used to retrieve, delete, or modify the task(2)
Task type: it includes point target task, nontrack strip task, curve strip task, and region task(3)
Task priority: it includes high, medium, and low priorities. Task priority is mainly used to ensure the execution of tasks with high priority after multiple tasks are mutually exclusive(4)
Payload type: it can be divided into visible light, multispectral, infrared, and other types according to the type of remote sensor(5)
Time range of observation task: the observation and data return of the task must be completed within the time range of the observation task. Some emergency observation tasks, such as flood and other disaster monitoring, will lead to the loss of due value of observation data beyond a specific time range(6)
Solar elevation angle: solar elevation angles for visible light imaging tasks(7)
Cloud cover: the prediction of cloudiness over imaging target and the probability of covering to target


#### 2.2.2. Decompose Tasks

Wu et al. presented a novel splitting algorithm to region target for satellite swarm [8], which is the first comprehensive solution for satellite task planning. Taking the method in the paper as a reference, a region target is handled by strips paralleled to subsatellite track.(1) Decomposition of Polygonal Area Task. It is necessary to decompose the polygonal areas according to the trend of the subsatellite track since a single strip cannot cover a complex task area, as shown in Figure 2.Figure 2
A polygonal area complex task decomposition.The first thing is to examine whether all area points are within the visibility corridor boundaries related to the subsatellite track and satellite maneuver before decomposing. The area would be decomposed according to the boundaries’ width if it passed the examination. Otherwise, this area becomes unobservable this time, so there is no need to decompose.The visibility corridor range can be calculated by(1)βmin=αmin−θ2,βmax=αmax+θ2,where βmin and βmax are the satellite’s minimum and maximum observation angle, αmin,αmax is the satellite maneuver range, and θ is the FOV (field of view).The polygon area is decomposed based on the widthk after obtaining the maneuver range, and the dynamic decomposition process is as follows:Step 1.
The irregular polygon is initialized into a convex polygonP by convex hull algorithm, and the subtask (metatask) set is initialized, M=Θ.Step 2.
Find the left most point positionx1 and the right-most point position x2 in the irregular area, set the left direction to negative and the right direction as positive, and let x=x1.Step 3.
Obtain the stripl, and let M=M∪l. Strip l is obtained by finding the intersected figure of polygon P and straight lines μ1 and μ2, where μ1 is obtained by drawing a line through x1 paralleled to the subsatellite track, and μ2 is obtained by translating μ1 by width k.Step 4.
Update the current location,x=x+k; the updated result is shown in Figure 3.Step 5.
Check whetherx is less than x2; if true, turn to Step 2.Step 6.
OutputM, end decomposition.Figure 3
Obtaining the strip.(2) Decomposition of Curve Strip Task. The curve strip task is described by point set P in which Pi1<i≤n is the characteristic point of a curve and is checked according to agile ability of satellite and imaging width k. Only the curve whose characteristic points set P are all in observable corridor can be decomposed. The procedure is as follows:Step 1.
SelectingPi as a start point of vector Vij and selecting point Pji<j≤n, one by one, as end point of vector Vij, then, the vector Vij=<Pi,Pj> is obtained.Step 2.
According to imaging widthk and vector Vij, the rectangle B whose length L is variable can be obtained.Step 3.
SelectingPj in i<j≤n, to get m_max end point, then to form vector Vijm_max=<Pi,Pjm_max>. Further, strip Bi that can cover all characteristic points between Pi and Pj can be obtained.Step 4.
Usingi=i+m_max to update i, then repeating Step 1 to Step 3, until i+m_max>n.The decomposition procedure of the curve strip task is shown in Figure4.Figure 4
Decomposition of the curve strip task.


### 2.3. Task Execution Time Window Calculation Algorithm

Visible window can be calculated based on the satellite’s maximum maneuver angleαmax and FOV θ, as is shown in Figure 5.Figure 5
Scheme of visible window calculation.A brief process of the visible window algorithm is as follows:Step 1.
Initialize the satellite observation rangeδ=α+θ and the execution window set T.Step 2.
The starting execution windowT1s=t1s,t′1s and the ending execution window T1e=t1e,t′1e for metatask m1s1,e1∈M are calculated according to the satellite orbit prediction and δ, where s1 is the starting point and e1 is the ending point.Step 3.
LetT=T∪T1.Step 4.
For∀mi∈M, repeat Step 2 and Step 3.Step 5.
Final setT will be the prediction of execution time window for every metatask decomposed from this complex task.


### 2.4. Task Planning

Mm1,m2,⋯,mn is recorded as the complex task M after decomposing (briefly introduced in Section 2.2.2), where n is the number of the meta tasks.For meta-taskmi∈M,
(i)
sloni,slati is longitude and altitude at starting point si(ii)
eloni,elati is longitude and altitude at ending point ei(iii)
li is strip lengthThe execution time is finally determined by the attitude and track planning, satisfying multiple constraints including attitude angle, angle velocity, and execution time window.


#### 2.4.1. Attitude Planning

The attitude representation method is based on 1-2-3 Euler angular rotation sequence. The roll angle isφ, pitch angle is θ, and yaw angle is ψ. The attitude planning algorithm is as follows:Step 1.
Calculate the observation point position at the execution time in J2000 coordinate system.r⟶p,ECI=x1p,y1p,z1p,vx1p,vy1p,vz1p; the starting point s1 position at t1s in the J2000 coordinate system can be calculated by its longitude and latitude, where the execution time window T1=t1s,t′1e for metatask m1s1,e1 is obtained in Section 2.3.Step 2.
Calculate the position and velocity of the satellite at the execution time in the J2000 coordinate system.Orbital data att1s can be obtained from the satellite GPS receiver, including a: semimajor axis; e: eccentricity; i: orbit inclination angle; Ω: right ascension of ascending node; w: argument of perigee; and M: mean anomaly.According to spacecraft dynamics equation [22, 23], r⟶SAT,ECI=x1s,y1s,z1s,vx1s,vy1s,vz1s, the position and velocity in the J2000 coordinate system can be obtained by the following formulas, after calculating the true anomaly f:
(2)x=rcosucosΩ−sinucosisinΩ,y=rcosusinΩ+sinucosicosΩ,z=rsinusini,vx=vsinγcosucosΩ−sinucosisinΩ+cosγ−sinucosΩ−cosucosisinΩ,vy=vsinγcosusinΩ+sinucosicosΩ+cosγ−sinusinΩ+cosucosicosΩ,vz=vsinγsinusini+cosγcosusini,r=a1−e21+ecosf,v=μ2r−1a,u=ω+f,sinγ=esinf1+e2+2ecosf,where r is the geocentric distance of the satellite; v is the satellite velocity; u is the argument of latitude of the satellite; γ is the track angle; μ is the geocentric gravitational constant, and μ=398600.44km3/s2.Step 3.
Calculate the satellite attitude parameters.Roll angleφ and pitch angle θ are obtained by the following formulas:
(3)θ=arcsinrx,orbitrorbit,rx,orbit>0,ry,orbit>0,θ∈0,π2,−arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit<0,θ∈−π2,0,π2−arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit>0,θ∈π2,π,−π2+arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit>0,θ∈−π,−π2,θ=arcsinrx,orbitrorbit,where Δr⟶ECI=r⟶SAT,ECI−r⟶p,ECI is the vector difference between the imaging point and the satellite in the J2000 coordinate system and Δr⟶orbit=rx,orbit,ry,orbit,rz,orbit is the orbit coordinate system converted from Δr⟶ECI.And roll angular velocityφ̇ and the pitch angular velocity θ̇ can be calculated using quadratic interpolation method, to improve the calculating speed and suit the limited on-board computational resources.Then, obtain the yaw angleψ according to φ, θ, φ̇, and θ̇. The yaw angular velocity ψ̇ is calculated in the same way as obtaining φ̇ and θ̇.Step 4.
Calculate the attitude parameters at the ending point.t1e is obtained by the formula
(4)t1e=t1s+l1vobs,where vobs is the push-broom speed of the satellite and li is the strip length of a metatask.Ift1e∉t1s,t′1e, which means the execution window T1 cannot cover t1e; this metatask has no solution. Otherwise, φ1e,θ1e,ψ1e,φ̇1e,θ̇1e,ψ̇1e can be obtained by repeating Step 1–Step 3 according to the orbit data at t1e.


#### 2.4.2. Track Planning

The satellite needs to switch attitude through track planning, fromφi−1e,θi−1e,ψi−1e,φ̇i−1e,θ̇i−1e,ψ̇i−1e at ti−1e to φis,θis,ψis,φ̇is,θ̇is,ψ̇is at tis between executing metatasks mi−1 and mi, with constraints of
(i)
maximum switching timeΔti−1,i=tis−ti−1e(ii)
maximum angular velocityvlim(iii)
maximum angular accelerationalimOther than that, the variation curves of angleθt, angular velocity θ̇t, and angular acceleration θ¨t in every axis (yaw, pitch, and roll) of the satellite need to be smooth. This is to ensure that the remote sensing satellite can perform tasks without jitter.The specific algorithm is as follows:(i)
Δti−1,i=τ, where τ is the minimum track planning time (we set τ=5 in this paper)(ii)
Calculate the satellite attitude parametersφi,tmps,θi,tmps,ψi,tmps,φ̇i,tmps,θ̇i,tmps,ψ̇i,tmps at ti,tmps=ti−1e+Δti−1,i(iii)
Obtain the accelerationφ¨i,tmps,θ¨i,tmps,ψ¨i,tmps using the quadratic interpolation methodStep 1.
InitializeΔti−1,i, and calculate the attitude parameters.Step 2.
Obtain the attitude parameters at every control period time (0.1 second in this paper) inΔti−1,i, based on the 5th-order track planning method.The formulas of this method are shown below:(5)θt=a0+a1t+a2t2+a3t3+a4t4+a5t5,θ̇t=a1+2a2t+3a3t2+4a4t3+5a5t4,θ¨t=2a2+6a3t+12a4t2+20a5t3,a0=θ0,a1=θ̇0,a2=θ¨02,a3=20θf−20θ0−8θ̇f+12θ̇0tf−3θ¨0‐θ¨ftf22tf2,a4=30θ0−30θf−14θ̇f+16θ̇0tf−3θ¨0‐2θ¨ftf22tf2,a5=12θf−12θ0−6θ̇f+6θ̇0tf−θ¨0‐θ¨ftf22tf2,where θ0,θ̇0,θ¨0 is the track planning starting angle, angular velocity, and angular acceleration at ti−1e and θf,θ̇f,θ¨f is the track planning ending angle, angular velocity, and angular acceleration at ti,tmpe.Step 3.
Calculate the compound angular velocityv and compound angular acceleration a using 1-2-3 Euler angle rotation kinematical equation, which can obtain maximum compound angular velocity vmax and maximum compound angular acceleration.Step 4.
Letℓv=vmax/vlim and ℓa=amax/alim. If ℓv≤1 and ℓa≤1, let tis=ti,tmps, track planning end. Otherwise, let ℓmax=maxℓv,ℓa, update Δti−1,i=Δti−1,i×ℓmax, and repeat Step 1–Step 3.


#### 2.4.3. Analyzing Time Complexity

The proposed method calculates angle relations between each orbital points and target vertexes that compose strip, region, and curve. Orbital points are orbital position in ECI (Earth Central of Inertial) coordinate system, and the points are sorted by time. The count of orbital points isnt second; for example, there are nt=24h×60min×60sec=86400 seconds in a day. If the count of vertexes is mv, the time complexity of the proposed algorithm is Omvnt.


## 2.1. Overall Design

The ground operation and control center receives users’ observation requirements and injects multiple complex tasks or metatasks into satellite. The on-board task planning system autonomously completes task planning, instruction generation, and other operations. The system architecture is shown in Figure1. After receiving and decomposing the complex tasks, the satellite implements task planning according to the real-time on-orbit state telemetry. A task planning process is divided into 8 subprocedures, that is, from requirement gathering to instruction generating, including the following:Step 1.
Gather requirements. The users send tasks to the ground station by means of polygons, curve strips, gaze, and tracking tasks and then forward to satellites through the satellite-ground link.Step 2.
Examine the validity of the task. Check task time, attributes, and types to determine whether the task is valid or not.Step 3.
Decompose complex tasks. The on-board task planning system decomposes and processes complex tasks according to the subsatellite track vector, generating the metatask.Step 4.
Planning execution window. Select imaging orbital arc.Step 5.
Calculate external constraints. Calculate the constraints of power [21], storage, attitude, and so on.Step 6.
Plan the attitude and attitude maneuver trajectory.Step 7.
Generate metatasks.Step 8.
Generate instructions.


## 2.2. Complex Task Description and Decomposition

### 2.2.1. Task Description

Complex tasks include irregular curve strips, irregular areas, point group target, etc. Tasks can be described by the following essential attributes:(1)
Task ID: the unique ID of the task, which is used to retrieve, delete, or modify the task(2)
Task type: it includes point target task, nontrack strip task, curve strip task, and region task(3)
Task priority: it includes high, medium, and low priorities. Task priority is mainly used to ensure the execution of tasks with high priority after multiple tasks are mutually exclusive(4)
Payload type: it can be divided into visible light, multispectral, infrared, and other types according to the type of remote sensor(5)
Time range of observation task: the observation and data return of the task must be completed within the time range of the observation task. Some emergency observation tasks, such as flood and other disaster monitoring, will lead to the loss of due value of observation data beyond a specific time range(6)
Solar elevation angle: solar elevation angles for visible light imaging tasks(7)
Cloud cover: the prediction of cloudiness over imaging target and the probability of covering to target


### 2.2.2. Decompose Tasks

Wu et al. presented a novel splitting algorithm to region target for satellite swarm [8], which is the first comprehensive solution for satellite task planning. Taking the method in the paper as a reference, a region target is handled by strips paralleled to subsatellite track.(1) Decomposition of Polygonal Area Task. It is necessary to decompose the polygonal areas according to the trend of the subsatellite track since a single strip cannot cover a complex task area, as shown in Figure 2.Figure 2
A polygonal area complex task decomposition.The first thing is to examine whether all area points are within the visibility corridor boundaries related to the subsatellite track and satellite maneuver before decomposing. The area would be decomposed according to the boundaries’ width if it passed the examination. Otherwise, this area becomes unobservable this time, so there is no need to decompose.The visibility corridor range can be calculated by(1)βmin=αmin−θ2,βmax=αmax+θ2,where βmin and βmax are the satellite’s minimum and maximum observation angle, αmin,αmax is the satellite maneuver range, and θ is the FOV (field of view).The polygon area is decomposed based on the widthk after obtaining the maneuver range, and the dynamic decomposition process is as follows:Step 1.
The irregular polygon is initialized into a convex polygonP by convex hull algorithm, and the subtask (metatask) set is initialized, M=Θ.Step 2.
Find the left most point positionx1 and the right-most point position x2 in the irregular area, set the left direction to negative and the right direction as positive, and let x=x1.Step 3.
Obtain the stripl, and let M=M∪l. Strip l is obtained by finding the intersected figure of polygon P and straight lines μ1 and μ2, where μ1 is obtained by drawing a line through x1 paralleled to the subsatellite track, and μ2 is obtained by translating μ1 by width k.Step 4.
Update the current location,x=x+k; the updated result is shown in Figure 3.Step 5.
Check whetherx is less than x2; if true, turn to Step 2.Step 6.
OutputM, end decomposition.Figure 3
Obtaining the strip.(2) Decomposition of Curve Strip Task. The curve strip task is described by point set P in which Pi1<i≤n is the characteristic point of a curve and is checked according to agile ability of satellite and imaging width k. Only the curve whose characteristic points set P are all in observable corridor can be decomposed. The procedure is as follows:Step 1.
SelectingPi as a start point of vector Vij and selecting point Pji<j≤n, one by one, as end point of vector Vij, then, the vector Vij=<Pi,Pj> is obtained.Step 2.
According to imaging widthk and vector Vij, the rectangle B whose length L is variable can be obtained.Step 3.
SelectingPj in i<j≤n, to get m_max end point, then to form vector Vijm_max=<Pi,Pjm_max>. Further, strip Bi that can cover all characteristic points between Pi and Pj can be obtained.Step 4.
Usingi=i+m_max to update i, then repeating Step 1 to Step 3, until i+m_max>n.The decomposition procedure of the curve strip task is shown in Figure4.Figure 4
Decomposition of the curve strip task.


## 2.2.1. Task Description

Complex tasks include irregular curve strips, irregular areas, point group target, etc. Tasks can be described by the following essential attributes:(1)
Task ID: the unique ID of the task, which is used to retrieve, delete, or modify the task(2)
Task type: it includes point target task, nontrack strip task, curve strip task, and region task(3)
Task priority: it includes high, medium, and low priorities. Task priority is mainly used to ensure the execution of tasks with high priority after multiple tasks are mutually exclusive(4)
Payload type: it can be divided into visible light, multispectral, infrared, and other types according to the type of remote sensor(5)
Time range of observation task: the observation and data return of the task must be completed within the time range of the observation task. Some emergency observation tasks, such as flood and other disaster monitoring, will lead to the loss of due value of observation data beyond a specific time range(6)
Solar elevation angle: solar elevation angles for visible light imaging tasks(7)
Cloud cover: the prediction of cloudiness over imaging target and the probability of covering to target


## 2.2.2. Decompose Tasks

Wu et al. presented a novel splitting algorithm to region target for satellite swarm [8], which is the first comprehensive solution for satellite task planning. Taking the method in the paper as a reference, a region target is handled by strips paralleled to subsatellite track.(1) Decomposition of Polygonal Area Task. It is necessary to decompose the polygonal areas according to the trend of the subsatellite track since a single strip cannot cover a complex task area, as shown in Figure 2.Figure 2
A polygonal area complex task decomposition.The first thing is to examine whether all area points are within the visibility corridor boundaries related to the subsatellite track and satellite maneuver before decomposing. The area would be decomposed according to the boundaries’ width if it passed the examination. Otherwise, this area becomes unobservable this time, so there is no need to decompose.The visibility corridor range can be calculated by(1)βmin=αmin−θ2,βmax=αmax+θ2,where βmin and βmax are the satellite’s minimum and maximum observation angle, αmin,αmax is the satellite maneuver range, and θ is the FOV (field of view).The polygon area is decomposed based on the widthk after obtaining the maneuver range, and the dynamic decomposition process is as follows:Step 1.
The irregular polygon is initialized into a convex polygonP by convex hull algorithm, and the subtask (metatask) set is initialized, M=Θ.Step 2.
Find the left most point positionx1 and the right-most point position x2 in the irregular area, set the left direction to negative and the right direction as positive, and let x=x1.Step 3.
Obtain the stripl, and let M=M∪l. Strip l is obtained by finding the intersected figure of polygon P and straight lines μ1 and μ2, where μ1 is obtained by drawing a line through x1 paralleled to the subsatellite track, and μ2 is obtained by translating μ1 by width k.Step 4.
Update the current location,x=x+k; the updated result is shown in Figure 3.Step 5.
Check whetherx is less than x2; if true, turn to Step 2.Step 6.
OutputM, end decomposition.Figure 3
Obtaining the strip.(2) Decomposition of Curve Strip Task. The curve strip task is described by point set P in which Pi1<i≤n is the characteristic point of a curve and is checked according to agile ability of satellite and imaging width k. Only the curve whose characteristic points set P are all in observable corridor can be decomposed. The procedure is as follows:Step 1.
SelectingPi as a start point of vector Vij and selecting point Pji<j≤n, one by one, as end point of vector Vij, then, the vector Vij=<Pi,Pj> is obtained.Step 2.
According to imaging widthk and vector Vij, the rectangle B whose length L is variable can be obtained.Step 3.
SelectingPj in i<j≤n, to get m_max end point, then to form vector Vijm_max=<Pi,Pjm_max>. Further, strip Bi that can cover all characteristic points between Pi and Pj can be obtained.Step 4.
Usingi=i+m_max to update i, then repeating Step 1 to Step 3, until i+m_max>n.The decomposition procedure of the curve strip task is shown in Figure4.Figure 4
Decomposition of the curve strip task.


## 2.3. Task Execution Time Window Calculation Algorithm

Visible window can be calculated based on the satellite’s maximum maneuver angleαmax and FOV θ, as is shown in Figure 5.Figure 5
Scheme of visible window calculation.A brief process of the visible window algorithm is as follows:Step 1.
Initialize the satellite observation rangeδ=α+θ and the execution window set T.Step 2.
The starting execution windowT1s=t1s,t′1s and the ending execution window T1e=t1e,t′1e for metatask m1s1,e1∈M are calculated according to the satellite orbit prediction and δ, where s1 is the starting point and e1 is the ending point.Step 3.
LetT=T∪T1.Step 4.
For∀mi∈M, repeat Step 2 and Step 3.Step 5.
Final setT will be the prediction of execution time window for every metatask decomposed from this complex task.


## 2.4. Task Planning

Mm1,m2,⋯,mn is recorded as the complex task M after decomposing (briefly introduced in Section 2.2.2), where n is the number of the meta tasks.For meta-taskmi∈M,
(i)
sloni,slati is longitude and altitude at starting point si(ii)
eloni,elati is longitude and altitude at ending point ei(iii)
li is strip lengthThe execution time is finally determined by the attitude and track planning, satisfying multiple constraints including attitude angle, angle velocity, and execution time window.


### 2.4.1. Attitude Planning

The attitude representation method is based on 1-2-3 Euler angular rotation sequence. The roll angle isφ, pitch angle is θ, and yaw angle is ψ. The attitude planning algorithm is as follows:Step 1.
Calculate the observation point position at the execution time in J2000 coordinate system.r⟶p,ECI=x1p,y1p,z1p,vx1p,vy1p,vz1p; the starting point s1 position at t1s in the J2000 coordinate system can be calculated by its longitude and latitude, where the execution time window T1=t1s,t′1e for metatask m1s1,e1 is obtained in Section 2.3.Step 2.
Calculate the position and velocity of the satellite at the execution time in the J2000 coordinate system.Orbital data att1s can be obtained from the satellite GPS receiver, including a: semimajor axis; e: eccentricity; i: orbit inclination angle; Ω: right ascension of ascending node; w: argument of perigee; and M: mean anomaly.According to spacecraft dynamics equation [22, 23], r⟶SAT,ECI=x1s,y1s,z1s,vx1s,vy1s,vz1s, the position and velocity in the J2000 coordinate system can be obtained by the following formulas, after calculating the true anomaly f:
(2)x=rcosucosΩ−sinucosisinΩ,y=rcosusinΩ+sinucosicosΩ,z=rsinusini,vx=vsinγcosucosΩ−sinucosisinΩ+cosγ−sinucosΩ−cosucosisinΩ,vy=vsinγcosusinΩ+sinucosicosΩ+cosγ−sinusinΩ+cosucosicosΩ,vz=vsinγsinusini+cosγcosusini,r=a1−e21+ecosf,v=μ2r−1a,u=ω+f,sinγ=esinf1+e2+2ecosf,where r is the geocentric distance of the satellite; v is the satellite velocity; u is the argument of latitude of the satellite; γ is the track angle; μ is the geocentric gravitational constant, and μ=398600.44km3/s2.Step 3.
Calculate the satellite attitude parameters.Roll angleφ and pitch angle θ are obtained by the following formulas:
(3)θ=arcsinrx,orbitrorbit,rx,orbit>0,ry,orbit>0,θ∈0,π2,−arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit<0,θ∈−π2,0,π2−arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit>0,θ∈π2,π,−π2+arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit>0,θ∈−π,−π2,θ=arcsinrx,orbitrorbit,where Δr⟶ECI=r⟶SAT,ECI−r⟶p,ECI is the vector difference between the imaging point and the satellite in the J2000 coordinate system and Δr⟶orbit=rx,orbit,ry,orbit,rz,orbit is the orbit coordinate system converted from Δr⟶ECI.And roll angular velocityφ̇ and the pitch angular velocity θ̇ can be calculated using quadratic interpolation method, to improve the calculating speed and suit the limited on-board computational resources.Then, obtain the yaw angleψ according to φ, θ, φ̇, and θ̇. The yaw angular velocity ψ̇ is calculated in the same way as obtaining φ̇ and θ̇.Step 4.
Calculate the attitude parameters at the ending point.t1e is obtained by the formula
(4)t1e=t1s+l1vobs,where vobs is the push-broom speed of the satellite and li is the strip length of a metatask.Ift1e∉t1s,t′1e, which means the execution window T1 cannot cover t1e; this metatask has no solution. Otherwise, φ1e,θ1e,ψ1e,φ̇1e,θ̇1e,ψ̇1e can be obtained by repeating Step 1–Step 3 according to the orbit data at t1e.


### 2.4.2. Track Planning

The satellite needs to switch attitude through track planning, fromφi−1e,θi−1e,ψi−1e,φ̇i−1e,θ̇i−1e,ψ̇i−1e at ti−1e to φis,θis,ψis,φ̇is,θ̇is,ψ̇is at tis between executing metatasks mi−1 and mi, with constraints of
(i)
maximum switching timeΔti−1,i=tis−ti−1e(ii)
maximum angular velocityvlim(iii)
maximum angular accelerationalimOther than that, the variation curves of angleθt, angular velocity θ̇t, and angular acceleration θ¨t in every axis (yaw, pitch, and roll) of the satellite need to be smooth. This is to ensure that the remote sensing satellite can perform tasks without jitter.The specific algorithm is as follows:(i)
Δti−1,i=τ, where τ is the minimum track planning time (we set τ=5 in this paper)(ii)
Calculate the satellite attitude parametersφi,tmps,θi,tmps,ψi,tmps,φ̇i,tmps,θ̇i,tmps,ψ̇i,tmps at ti,tmps=ti−1e+Δti−1,i(iii)
Obtain the accelerationφ¨i,tmps,θ¨i,tmps,ψ¨i,tmps using the quadratic interpolation methodStep 1.
InitializeΔti−1,i, and calculate the attitude parameters.Step 2.
Obtain the attitude parameters at every control period time (0.1 second in this paper) inΔti−1,i, based on the 5th-order track planning method.The formulas of this method are shown below:(5)θt=a0+a1t+a2t2+a3t3+a4t4+a5t5,θ̇t=a1+2a2t+3a3t2+4a4t3+5a5t4,θ¨t=2a2+6a3t+12a4t2+20a5t3,a0=θ0,a1=θ̇0,a2=θ¨02,a3=20θf−20θ0−8θ̇f+12θ̇0tf−3θ¨0‐θ¨ftf22tf2,a4=30θ0−30θf−14θ̇f+16θ̇0tf−3θ¨0‐2θ¨ftf22tf2,a5=12θf−12θ0−6θ̇f+6θ̇0tf−θ¨0‐θ¨ftf22tf2,where θ0,θ̇0,θ¨0 is the track planning starting angle, angular velocity, and angular acceleration at ti−1e and θf,θ̇f,θ¨f is the track planning ending angle, angular velocity, and angular acceleration at ti,tmpe.Step 3.
Calculate the compound angular velocityv and compound angular acceleration a using 1-2-3 Euler angle rotation kinematical equation, which can obtain maximum compound angular velocity vmax and maximum compound angular acceleration.Step 4.
Letℓv=vmax/vlim and ℓa=amax/alim. If ℓv≤1 and ℓa≤1, let tis=ti,tmps, track planning end. Otherwise, let ℓmax=maxℓv,ℓa, update Δti−1,i=Δti−1,i×ℓmax, and repeat Step 1–Step 3.


### 2.4.3. Analyzing Time Complexity

The proposed method calculates angle relations between each orbital points and target vertexes that compose strip, region, and curve. Orbital points are orbital position in ECI (Earth Central of Inertial) coordinate system, and the points are sorted by time. The count of orbital points isnt second; for example, there are nt=24h×60min×60sec=86400 seconds in a day. If the count of vertexes is mv, the time complexity of the proposed algorithm is Omvnt.


## 2.4.1. Attitude Planning

The attitude representation method is based on 1-2-3 Euler angular rotation sequence. The roll angle isφ, pitch angle is θ, and yaw angle is ψ. The attitude planning algorithm is as follows:Step 1.
Calculate the observation point position at the execution time in J2000 coordinate system.r⟶p,ECI=x1p,y1p,z1p,vx1p,vy1p,vz1p; the starting point s1 position at t1s in the J2000 coordinate system can be calculated by its longitude and latitude, where the execution time window T1=t1s,t′1e for metatask m1s1,e1 is obtained in Section 2.3.Step 2.
Calculate the position and velocity of the satellite at the execution time in the J2000 coordinate system.Orbital data att1s can be obtained from the satellite GPS receiver, including a: semimajor axis; e: eccentricity; i: orbit inclination angle; Ω: right ascension of ascending node; w: argument of perigee; and M: mean anomaly.According to spacecraft dynamics equation [22, 23], r⟶SAT,ECI=x1s,y1s,z1s,vx1s,vy1s,vz1s, the position and velocity in the J2000 coordinate system can be obtained by the following formulas, after calculating the true anomaly f:
(2)x=rcosucosΩ−sinucosisinΩ,y=rcosusinΩ+sinucosicosΩ,z=rsinusini,vx=vsinγcosucosΩ−sinucosisinΩ+cosγ−sinucosΩ−cosucosisinΩ,vy=vsinγcosusinΩ+sinucosicosΩ+cosγ−sinusinΩ+cosucosicosΩ,vz=vsinγsinusini+cosγcosusini,r=a1−e21+ecosf,v=μ2r−1a,u=ω+f,sinγ=esinf1+e2+2ecosf,where r is the geocentric distance of the satellite; v is the satellite velocity; u is the argument of latitude of the satellite; γ is the track angle; μ is the geocentric gravitational constant, and μ=398600.44km3/s2.Step 3.
Calculate the satellite attitude parameters.Roll angleφ and pitch angle θ are obtained by the following formulas:
(3)θ=arcsinrx,orbitrorbit,rx,orbit>0,ry,orbit>0,θ∈0,π2,−arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit<0,θ∈−π2,0,π2−arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit>0,θ∈π2,π,−π2+arcsinrx,orbitrorbit,rx,orbit<0,ry,orbit>0,θ∈−π,−π2,θ=arcsinrx,orbitrorbit,where Δr⟶ECI=r⟶SAT,ECI−r⟶p,ECI is the vector difference between the imaging point and the satellite in the J2000 coordinate system and Δr⟶orbit=rx,orbit,ry,orbit,rz,orbit is the orbit coordinate system converted from Δr⟶ECI.And roll angular velocityφ̇ and the pitch angular velocity θ̇ can be calculated using quadratic interpolation method, to improve the calculating speed and suit the limited on-board computational resources.Then, obtain the yaw angleψ according to φ, θ, φ̇, and θ̇. The yaw angular velocity ψ̇ is calculated in the same way as obtaining φ̇ and θ̇.Step 4.
Calculate the attitude parameters at the ending point.t1e is obtained by the formula
(4)t1e=t1s+l1vobs,where vobs is the push-broom speed of the satellite and li is the strip length of a metatask.Ift1e∉t1s,t′1e, which means the execution window T1 cannot cover t1e; this metatask has no solution. Otherwise, φ1e,θ1e,ψ1e,φ̇1e,θ̇1e,ψ̇1e can be obtained by repeating Step 1–Step 3 according to the orbit data at t1e.


## 2.4.2. Track Planning

The satellite needs to switch attitude through track planning, fromφi−1e,θi−1e,ψi−1e,φ̇i−1e,θ̇i−1e,ψ̇i−1e at ti−1e to φis,θis,ψis,φ̇is,θ̇is,ψ̇is at tis between executing metatasks mi−1 and mi, with constraints of
(i)
maximum switching timeΔti−1,i=tis−ti−1e(ii)
maximum angular velocityvlim(iii)
maximum angular accelerationalimOther than that, the variation curves of angleθt, angular velocity θ̇t, and angular acceleration θ¨t in every axis (yaw, pitch, and roll) of the satellite need to be smooth. This is to ensure that the remote sensing satellite can perform tasks without jitter.The specific algorithm is as follows:(i)
Δti−1,i=τ, where τ is the minimum track planning time (we set τ=5 in this paper)(ii)
Calculate the satellite attitude parametersφi,tmps,θi,tmps,ψi,tmps,φ̇i,tmps,θ̇i,tmps,ψ̇i,tmps at ti,tmps=ti−1e+Δti−1,i(iii)
Obtain the accelerationφ¨i,tmps,θ¨i,tmps,ψ¨i,tmps using the quadratic interpolation methodStep 1.
InitializeΔti−1,i, and calculate the attitude parameters.Step 2.
Obtain the attitude parameters at every control period time (0.1 second in this paper) inΔti−1,i, based on the 5th-order track planning method.The formulas of this method are shown below:(5)θt=a0+a1t+a2t2+a3t3+a4t4+a5t5,θ̇t=a1+2a2t+3a3t2+4a4t3+5a5t4,θ¨t=2a2+6a3t+12a4t2+20a5t3,a0=θ0,a1=θ̇0,a2=θ¨02,a3=20θf−20θ0−8θ̇f+12θ̇0tf−3θ¨0‐θ¨ftf22tf2,a4=30θ0−30θf−14θ̇f+16θ̇0tf−3θ¨0‐2θ¨ftf22tf2,a5=12θf−12θ0−6θ̇f+6θ̇0tf−θ¨0‐θ¨ftf22tf2,where θ0,θ̇0,θ¨0 is the track planning starting angle, angular velocity, and angular acceleration at ti−1e and θf,θ̇f,θ¨f is the track planning ending angle, angular velocity, and angular acceleration at ti,tmpe.Step 3.
Calculate the compound angular velocityv and compound angular acceleration a using 1-2-3 Euler angle rotation kinematical equation, which can obtain maximum compound angular velocity vmax and maximum compound angular acceleration.Step 4.
Letℓv=vmax/vlim and ℓa=amax/alim. If ℓv≤1 and ℓa≤1, let tis=ti,tmps, track planning end. Otherwise, let ℓmax=maxℓv,ℓa, update Δti−1,i=Δti−1,i×ℓmax, and repeat Step 1–Step 3.


## 2.4.3. Analyzing Time Complexity

The proposed method calculates angle relations between each orbital points and target vertexes that compose strip, region, and curve. Orbital points are orbital position in ECI (Earth Central of Inertial) coordinate system, and the points are sorted by time. The count of orbital points isnt second; for example, there are nt=24h×60min×60sec=86400 seconds in a day. If the count of vertexes is mv, the time complexity of the proposed algorithm is Omvnt.


## 3. On-Orbit Verification Results and Analysis

An on-orbit flight test on satellite (GFDM) with a single-core BM3803 processor (main frequency: 50 MHz) has been taken to verify the system and method proposed in this paper. Each satellite subsystem executes the generated instruction by receiving it from the on-board routing system using the 1553B avionics bus.This verification was carried out on July 5th, 2020. This verification was carried out using 1/6 of the computing resources for automatic task planning. The result turns out that the proposed method can accomplish complex task, including decomposition and automatic planning.A complex task consisting of a non-along-track strip task and a curve area task is taken as an example in this paper to introduce the on-orbit verification briefly. The specific process and results of this verification are as follows:(1)
Orbit data in verificationGFDM satellite takes the trajectory between two descending nodes as one orbit. Data output interval is set as 1 second. In on-orbit verification, the orbit data is between 2020-7-5 11:40:00 (UTCG) and 2020-7-5 19:40:00 (BJT), semimajor axis 7020.634788 km, eccentricity 0.000381, inclination angle 97.905112°, right ascension of ascending node 260.231848°, argument of perigee 155.560876°, and mean anomaly 229.827296°.(2)
Task injection in verificationData of the two subtasks (a curve area task and a non-along-track strip task) are as follows:(1)
Geographic position of curve area task: <-51.7°, 4.3°>; <-51.2°, 3.9°>; <-50.7°, 2.2°>; <-50.5°, 1.9°>; and <-50.0°, 1.7°>(2)
Geographic position of non-along-track strip task: <-58.4°, -33°>; <-58.4°, -35°>The generated instruction package taking the curve area task as an example is shown in Figure6.
(3)
On-board task decomposition in verificationFigure 6
Instruction package on curve complex task.The command to start the on-board autonomous mission planning is sent at 19:50:11, July 5th, 2020, Beijing time. After about 117 seconds, the on-board mission decomposition and planning results are received at 19:52:08, July 5th, 2020, Beijing time, and transmitted to the ground center in the form of telemetry package. As shown in Figure6, the curve task is decomposed into three imaging metatasks. Taking metatask 1 as an example, its planned start time is 173887001 (0x0a5d4e19) seconds of on-board time. The corresponding Beijing time is 2020-7-5 21:55:48. The time coding method is similar to most typical embedded system, like IoT (Internet of Things) [24, 25].The command to boot the on-board autonomous task planning is at BJT 2020-7-5 19:50:11. And the results of on-board task decomposition and planning are at BJT 2020-7-5 19:52:08. The instruction package is shown in Figure7.
(4)
Attitude data after planningFigure 7
Results of metatasks planned on board for the curve area task.The results of attitude planning and track planning are transmitted down to the ground station in the form of instruction package, as shown in Figure8.
(5)
Curve strip task verificationFigure 8
On-board attitude planning and trajectory planning of curve complex tasks.To verify the curve strip task, a curve task with 5 characteristic points is planned on board. Five points by longitude and latitude with degree are <-51.7, 4.3>; <-51.2, 3.9>; <-50.7, 2.2>; <-50.5, 1.9>; and <-50.0, 1.7>. The curve and its partial zoomed position are shown in Figure9 in which Figure 9(a) is the global view. The curve strip task is planned and executed by GFDM satellite, and its partial image (corresponding to Figure 9(b)) is shown in Figure 10.
(6)
Non-along-track strip task verificationFigure 9
The curve strip task shown in map software.Figure 10
Imaging result of Figure9(b) by on-board mission planning.A non-along-track strip with start point at <-58.4, -33> and end point at <-58.4, -35> is injected to on-orbit satellite. The non-along-track strip is lined in Figure11(a), and its partial is zoomed into Figure 11(b).Figure 11
Non-along-track strip task shown in map software.Non-along-track strip task was successfully planned and executed on-board on July 5th, 2020, Beijing time, and its result of image at Buenos Aires is shown in Figure12.Figure 12
Imaging result of Figure11(b) by on-board mission planning.


## 4. Conclusion

The results from the experiment of the GFDM show that the automatic on-board task planning system can bring the following capabilities’ improvement:(1)
The on-orbit task planning is feasible. The automatic level of satellite is improved; meanwhile, the dependence of satellites on the ground has been decreased(2)
Facilitated customer manipulation. Customers can focus on task-related matters and do not need to consider the satellite’s complex parameters and knowledge(3)
Fully utilized satellite capability. The task execution is optimized based on the real-time constraints on the satellite, including energy, storage, and orbit



---

*Source: 8923088-2023-02-28.xml*