2
$\begingroup$

I am a beginner in the field of aerospace. Recently, I wanted to use the extended Kalman filter to improve the orbit.
My idea is to use TLE data as observations values, the CowellPropagator in python's third library polistro.twobody.propagation module as the physical model to output the estimated value. Then I will make a data fusion between the two values.
However, there seems to be something wrong with the propagator I designed.
For the following two TLE data:

0 STARLETTE
1  7646U 75010A   22193.48472273 -.00000146  00000-0 -15590-5 0  9998
2  7646  49.8266  11.4131 0205800   0.2699 359.8273 13.82312865395926

0 STARLETTE
1 07646U 75010A   22193.84620530 -.00000147  00000-0 -21717-5 0  9998
2 07646  49.8266   9.9864 0205799   1.4641 358.6814 13.82312857395994

I use the following python code to convert and get their observation time, position and velocity under GCRS

2022-07-12 11:38:00.043871  7048467.162860288   1385737.9412575748  -15526.126134106882 -930.2068177030378  4762.392389570778   5755.6164340768755
2022-07-12 20:18:32.137920  7080903.488568736   1209362.4334778101  -16177.630768791158 -813.6988128053417  4783.958444397058   5755.325126774668

Then I want to use the propagator designed by myself to propagate from the first observation time to the second

A = 0.723456
m = 47
A_over_m = A / m
C_D = 2.0
kernel = SPK.open('..\sample_data\de421.bsp')
gravity_model = GGM03S('..\sample_data\GGM03S.txt')
df_d = sw.sw_daily(update=True)
f107a = df_d.loc["2022-07-12"].f107_81ctr_obs
ap = df_d.loc["2022-07-12"].Apavg
f107 = df_d.loc["2022-07-11"].f107_81ctr_obs
def f_accel(t0, u_, k):
    du_kep = func_twobody(t0, u_, k)
    # J2_perturbation
    ax, ay, az = J2_perturbation(
        t0, u_, k, J2=Earth.J2.value, R=Earth.R.to(u.km).value
    )
    ax_J3, ay_J3, az_J3 = J3_perturbation(
        t0, u_, k, J3=Earth.J3.value, R=Earth.R.to(u.km).value
    )
  
    # atmospheric_drag
    time_gap = datetime.timedelta(seconds=t0)
    now_time = row_epoch + time_gap
    # a = gravity_model.AccelHarmonic(u_[:3]*1000, now_time, 10, 10)
    lat, lon, alt = pymap3d.eci2geodetic(u_[0] * 1000, u_[1] * 1000, u_[2] * 1000, now_time)
    densities, temperatures = msise_model(now_time, alt, lat, lon, f107a, f107, ap)
    a_x, a_y, a_z = atmospheric_drag(
        t0, u_, k, C_D, A_over_m, densities[5] * 1e12
    )

    # Sun perturbation
    fd, jr = jday(now_time.year, now_time.month, now_time.day, now_time.hour, now_time.minute, now_time.second)
    position_Sun = kernel[0, 10].compute(fd+jr)
    position_Sun -= kernel[0, 3].compute(fd+jr)
    position_Sun -= kernel[3, 399].compute(fd + jr)
    delta_r = position_Sun - u_[:3]
    Sun_k = Sun.k.value / 1000000000
    a_x1, a_y1, a_z1 = (
            Sun_k * delta_r / norm(delta_r) ** 3
            - Sun_k * position_Sun / norm(position_Sun) ** 3
    )

    # Moon perturbation
    position_Moon = kernel[3, 301].compute(fd+jr)
    delta_r = position_Moon - u_[:3]
    Moon_k = Moon.k.value / 1000000000
    # print('norm(delta_r) ** 3:', norm(delta_r) ** 3)
    a_x2, a_y2, a_z2 = (
            Moon_k * delta_r / norm(delta_r) ** 3
            - Moon_k * position_Moon / norm(position_Moon) ** 3
    )

    du_ad = np.array([0, 0, 0, ax + a_x + a_x2 + a_x1 + ax_J3, ay + a_y + a_y2 + a_y1 + ay_J3, az + a_z + a_z2 + a_z1 + az_J3])
    return du_kep + du_ad
t2 = datetime.datetime.strptime("2022-07-12 11:38:00", "%Y-%m-%d %H:%M:%S")
t3 = datetime.datetime.strptime("2022-07-12 20:18:32", "%Y-%m-%d %H:%M:%S")
elapsed_seconds = (t3 - t2).total_seconds()
Y = [7048467.162860288, 1385737.9412575748, -15526.126134106882, -930.2068177030378, 4762.392389570778, 5755.6164340768755]
orbit = Orbit.from_vectors(Earth, Y[:3] * u.m, Y[3:6] * u.m/u.s)
orbit = orbit.propagate(elapsed_seconds * u.s, method=CowellPropagator(f=f_accel))
r, v = orbit.rv()
r_m = r.to(u.m).to_value()
v_ms = v.to(u.m / u.s).to_value()

The results obtained are quite different from the observed values:

r_x  7079656.958443244 
r_y  1216598.9037586488 
r_z  -7390.772128518734 
v_x  -825.2388349140476 
v_y  4781.665127771911 
v_z  5755.610689348027

So my question is:

  • If I must use numerical integration instead of SGP4, how can I improve?
  • I have also used this propagator on starlink data, but the results are accurate. Is the difference between the two caused by orbital height?
$\endgroup$
1
  • 3
    $\begingroup$ NORAD has said that their propagation model is the only thing that should be used with TLE's. This is because they fit the TLE data to that model, and even more accurate propagation methods will produce worse results. $\endgroup$ Oct 5, 2022 at 13:40

1 Answer 1

8
$\begingroup$

This is a bad idea. Two line elements must be propagated with the SGP4 algorithm. While some of the two line element names are the same as classical Keplerian element names, they are not the same. The intent of TLEs and SGP4 is to provide an approximate Cartesian state (position and velocity) that accounts for some perturbations.

Do not treat those Cartesian states as observations. The error can be well into tens of kilometers, sometimes more. The intent is to enable

  • Quick calculation of satellite passes, where an error in the hundreds of kilometers is fine,
  • Quick calculation of acquisition time and pointing for an omnidirectional antenna, where an error of about a hundred kilometers is fine,
  • Quick calculation of acquisition time and pointing for a tracking antenna, where an error of a few tens of kilometers is fine,
  • Quick calculation of possible collisions, where a several kilometer error is fine.

The point is "quick calculation". The SGP4 algorithm inherently values speed of calculation over accuracy. Do not use TLEs and the SGP4 algorithm as "observations".

$\endgroup$
1
  • $\begingroup$ Thank you for your reply! So I want to know, where can I get the observations of resident space object (such as the satellite starlette I use)? $\endgroup$ Oct 6, 2022 at 3:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.