¿Es porque la Tierra no es exactamente una esfera y también debido a la distribución desigual de la masa?

¡Sí!

¿Es un efecto atmosférico?

No.

El término grande J2 o achatado de la Tierra siempre está tratando de perturbar las órbitas de los satélites, y el efecto es más fuerte en LEO.

En esta respuesta , por ejemplo, muestro cómo J2 avanza constantemente el argumento del perigeo de una órbita elíptica a pesar de que tiene una inclinación cero. Esto se debe a que esta fuerza varía en 1/r^4 en comparación con el término del monopolo terrestre que varía en 1/r^2.

Primero busqué una expresión simple para el período de las oscilaciones de excentricidad similar a la ecuación para el cambio de período debido a J2 en esta respuesta , pero no pude encontrar algo simple. Parece ser una dependencia más complicada de la inclinación, el semieje mayor y la excentricidad.

Entonces, en cambio, pensé en hacer una integración numérica directa como lo hice aquí y aquí, pero sentí que sería difícil probar en la respuesta que las oscilaciones no eran solo un efecto numérico.

Así que decidí "engañar" al propagador SPG4 que se ejecuta en Skyfield para propagar las órbitas de los satélites, ya que es bastante confiable y estable para las órbitas "normales", aunque no es extremadamente preciso.

Generé "TLE falsos" ( ¡no intente esto en casa! ) Basados ​​en el TLE del 1 de abril de 2017 de TianGong-1, poniendo a cero los términos de arrastre y varios otros parámetros, dejando solo el tamaño de la órbita (expresado en movimiento medio en revoluciones por día), inclinación y excentricidad. Propagué la órbita durante 120 minutos a partir de la medianoche UTC durante 100 días. Para cada día, registro las altitudes máximas y mínimas. Desde que puse a cero la resistencia, la altitud permanece más o menos constante.

Los gráficos son para inclinaciones de 0 a 90 grados en pasos de 5 grados, y se ven totalmente diferentes entre sí y, aunque son periódicas, las formas no son realmente sinusoidales .

nota: En la pregunta la gráfica muestra una periodicidad de unos 50 días. Tiangong-1 tenía una inclinación de unos 42,8 grados. En el gráfico siguiente, el período a 40 grados es inferior a 50 días y a 45 es superior a 50. Así que al menos podemos confirmar que el cálculo reproduce el comportamiento.

Dado que el período de este LEO es inferior a 225 minutos, SGP4 no utiliza las correcciones de "espacio profundo" de SDP4 para la gravedad del Sol y la Luna. Entonces, el único efecto que queda es la desviación del campo gravitatorio de la Tierra de la forma esférica.

nota: esta ilustración usa SGP4 para propagar un satélite en LEO durante +/- 250 días sobre la época de TLE y esa no es la forma correcta de usar TLE y SGP4 con precisión. Lo he hecho aquí solo para mostrar grandes cambios en la órbita, pero para cálculos cuidadosos, apéguese a los tiempos cercanos a la época del TLE . Puede leer esta excelente respuesta para obtener más información al respecto.

def make_fake_TLE(incdegs, ecc, altkm):
    "don't try this at home!!!"
    a          = (6378+altkm) * 1000.
    revsperday = 24.*3600/(twopi * np.sqrt(a**3 / GMe))

    eccstr = '{:07d}'.format(int(1E+07 * ecc))
    mmstr  = '{:011.8f}'.format(revsperday)
    incstr = '{:08.4f}'.format(incdegs)

    L1 = "1 00000x xxxxxx   17091.00000000 +.00000000 +00000-0 +00000-3 0  0000"
    L2 = "2 00000 " + incstr + " 000.0000 " + eccstr + " 000.0000 000.0000 " + mmstr + "000000"

    return (L1, L2)

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos, EarthSatellite

# Real TLE for TianGong-1 for April 1, 2017
# "1 37820U 11053A   17091.81919743 +.00033376 +00000-0 +21646-3 0  9994"
# "2 37820 042.7591 306.4244 0019428 259.0955 165.0115 15.75087641315914"

twopi = 2.*np.pi
GMe   = 3.986E+14  # earth standard gravitational parameter
Re_km = 6378.

days    = np.arange(500) - 250
minutes = np.arange(120)

delta_inc = 5. # degrees
n_inc     = 19 
incdegs = delta_inc * np.arange(n_inc) 

### NOTE!  The first time you run, set calculate to True! It takes several minutes to
# propagate all of the orbits, then saves the result as a  .npy file
# Then, set calculate to False. Now you can play with plotting without
# doing any more calculations, so it's much much faster.
# https://matplotlib.org/gallery/lines_bars_and_markers/fill_between_demo.html

calculate = True
fname     = 'incgroups19x-250to+250'

if calculate:
    load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
    ts   = load.timescale()

    data    = load('de421.bsp')
    earth   = data['earth']
    ts      = load.timescale()
    incgroups = []
    for incdeg in incdegs: 
        rmins, rmaxs = [], []
        for day in days:
            if not day%100:
                print day,
            time    = ts.utc(2017, 4, day, 0, minutes)  # starting April 1, 2017
            L1, L2  = make_fake_TLE(incdeg, 0.0019, 350.)
            TG1     = EarthSatellite(L1, L2)
            TG1pos  = TG1.at(time).position.km
            r       = np.sqrt((TG1pos**2).sum(axis=0))
            rmins.append(r.min())
            rmaxs.append(r.max())
        rmins = np.array(rmins)
        rmaxs = np.array(rmaxs)
        incgroups.append([rmins, rmaxs])

    incgroups = np.array(incgroups)

    if fname:
        q = np.stack(incgroups)
        np.save(fname, q)
else:
    if fname[-4].lower() != '.npy':
        fname += '.npy'
    incgroups = np.load(fname)

groups = incgroups.copy()

# post-process it
dh = 40.
for i, thing in enumerate(groups):
    thing += dh*i
baseline = groups[0].mean()
groups   = delta_inc * (groups-baseline)/dh


if True:
    N         = len(groups)
    colors    = plt.rcParams['axes.prop_cycle'].by_key()['color']
    for i, (rmin, rmax) in enumerate(groups):
        c  = colors[i%len(colors)]
        dh = 40. * i
        plt.fill_between(days, rmin, rmax,
                         color=None, facecolor=c, linewidth=0)
    plt.xlim(-250, 250)
    plt.ylim(-3, 94)
    plt.xlabel('days', fontsize=16)
    plt.title('relative apo-peri envelopes', fontsize=16)
    plt.ylabel('inclination (degs)', fontsize=16)
    plt.show()

if True:
    colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
    fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
    for i, (rmin, rmax) in enumerate(incgroups):
        c = colors[i%len(colors)]
        ax1.plot(days, rmin-Re_km, c)
        ax1.plot(days, rmax-Re_km, c)
        ax2.fill_between(days, rmin-Re_km, rmax-Re_km, color='b', facecolor=c)
    ax1.set_xlim(-250, 250)
    ax2.set_xlim(-250, 250)
    plt.show()

Después de la excelente respuesta del uhoh, no hay necesidad de otra respuesta, pero con una simulación podemos ver el efecto del achatamiento de la Tierra y de la atmósfera.

La simulación incluye las aceleraciones newtoniana y relativista de todos los planetas, Sol y Luna.
El campo de gravedad de la Tierra se modela con el modelo de gravedad SGG-UGM-1 (calculado utilizando la anomalía de gravedad derivada de EGM2008 y los datos de observación de GOCE) truncado al grado y orden 15 (para ahorrar el tiempo de ejecución, manteniendo una buena precisión en comparación con el total). modelo).
Para el cálculo de la densidad del aire, utilizo el modelo NRLMSISE-00 junto con un archivo de datos actualizado para los índices solar y geomagnético. Los índices reales se pueden encontrar aquí: www.celestrak.com/spacedata/SW-All.txt.

El primer paso consiste en determinar el mejor coeficiente balístico para minimizar un parámetro de simulación en particular. Después de 36 minutos, el programa encuentra un coeficiente balístico de aproximadamente$111\,kg/m^2$(no es fijo, porque el coeficiente de arrastre varía con la composición del aire).

Ahora se puede iniciar la simulación:
1) con un TLE y el propagador SGP4 del CSpOC calcular el estado inicial (posición y velocidad) del satélite para la época TLE;
2) propagar ese estado inicial con un propagador especialmente diseñado (mi propagador se basa en el integrador 8(5,3) Dormand-Prince);
3) cuando la altitud del satélite cae por debajo de 70 km, detener la simulación.

Aquí está el resultado obtenido con el TLE 17080.91114878:

El gráfico muestra que mi modelo reproduce con precisión el gráfico de Wikipedia.

¿Por qué los gráficos son inestables? ¿Es porque la Tierra no es exactamente una esfera y también debido a la distribución desigual de la masa?

Puedes responder a esa pregunta simulando la Tierra como si fuera una esfera; aquí está el resultado:

Ya no se muestra el bamboleo (vemos una tasa de descomposición irregular porque la densidad del aire varía según la actividad solar).

¿Por qué los gráficos son inestables? ¿Es un efecto atmosférico?

Para responder a esa pregunta solo necesitamos desactivar la atmósfera; aquí está el resultado:

El bamboleo es exactamente el mismo que el modelo completo.

En comparación con otros gráficos de la órbita, la oscilación parece ser bastante exagerada.

teorías

No soy un experto en esto, pero mirando un poco, tengo algunas teorías y encontré algunas oscilaciones de altitudes y excentricidades para las estaciones espaciales.

Precesión nodal

Creo que se debe a la precesión nodal de la estación, el período es de unos dos meses para la ISS, creo que sería similar para Tiangong. Encontré esto mientras miraba esta respuesta . La precesión nodal es la inclinación de la estación que oscila debido a que la tierra no es una esfera.

Entonces esperaría que la altitud promedio de la estación descienda con inclinaciones bajas y aumente con inclinaciones altas. Sin embargo, eso no explica completamente por qué Apogee y Perigee se oponen entre sí en este caso.

Eso está mal. Gracias a @MAH por atrapar ese. Estúpido de mí, lo busqué / lo recordé mal.

resonancia lunar

También podría deberse a algún tipo de influencia lunar con esos dos meses, pero solo una suposición.

Diferentes oscilaciones

resonancia de excentricidad de 14 días

En los foros de nasaspaceflight encontraron una oscilación más lenta, que podría estar relacionada.

Media órbita

Otra oscilación, mucho más pequeña, que se puede ver en la imagen de abajo en la parte superior derecha, se debe en realidad a los oblatos de la tierra. La altitud disminuye cada vez más sobre el ecuador que es más grueso que el resto de la tierra. fuente