Ballisztikus pálya előrejelzése Kálmán-szűrővel


A történet

Néhány hete (2018. novemberében) újabb összecsapások voltak az Izraeli és a Palesztin erők között.  A Palesztin erők több száz rakétát lőttek ki Izraeli területre, ezek nagy részét a Vaskupola nevű védelmi rendszer hatástalanította. Általában ez nem követte a politikai eseményeket. Ez most kivétel lett. Az érdeklődésemet az keltette fel, hogy megtudtam azt, hogy az Izraeli védelmi rendszer nem próbál minden egyes támadó rakétát megsemmisíteni.  Hogy mért nem? Ennek nagyon prózai oka van: pénz.  A palesztin Kasszám-rakéták kb. 500-800 dollárba kerülnek darabonként, ezzel szemben az elfogásukra használt izraeli Tamir-ok kb. 45 000 dollárba.  Bár az Izraeli Védelmi erőknek biztosan nagyobb a költségvetése mint a Hamaszé, de azért ez elég nagy különbség. A költséghatékonyság jegyében a Vaskupola először megbecsüli, hogy a támadó rakéta lakott, vagy egyéb kiemelten védett területen fog-e becsapódni és csak akkor indítja az ellencsapást, ha erre a kérdésre igen válasz születik.

Ez a tény felébresztette a kíváncsiságomat. Először is, hogy lehet megbecsülni a várható pályát?  Másodszor: mennyire bizonytalanok azok a radar adatok, amik ehhez a feladathoz rendelkezésre állnak? Ezért utána néztem a ballisztikus pályán mozgó testek pályamodelljének.  Erről készítettem is egy bejegyzést: Ballisztikus pálya state-space modellje.

Megpróbáltam szimulálni egy ilyen feladatot, és megoldást találni rá.  Így érkeztünk el a mai bejegyzéshez és a Kálmán szűrőhöz. Bár az eredeti szűrő lineáris jelekre lett kidolgozva, van nem-lineáris megvalósítása is, és tudom, hogy ezeket használják például levegő-levegő rakéták irányítására. Szóval kézenfekvőnek tűnt, hogy megnézzem alkalmas-e a becsapódási pont előrejelzésére.

A szimuláció során feltettem, hogy az Izraeli erők:

  • ismerik az alapvető meteorológiai körülményeket, mint például légsűrűség
  • beazonosítják az ellenség által használt rakétát a kiégés előtt1, és ismerik ennek a rakétatípusnak az alapadatait
  • a rendszer (EL/M-2084) zaja normál eloszlást követ2, ennek a zajnak a standard eloszlása 125m
  • a támadó rakéta üzemanyagának kiégése után 10 másodperc áll rendelkezésre annak eldöntésere, hogy elindítsák-e a lelövést.

A szimulációnál csak a már említett bejegyzés Ferde hajítás légellenállással stat-space modelljét használtam.  Bár a modell korántsem teljes, de a megoldás során használt Kálmán szűrőben lényegében bármilyen modellt be lehet helyettesíteni egyetlen függvény átírásával. Ennek megfelelően jobb modell esetén könnyen továbbfejleszthető a kód. A Kasszám-rakétákról szinte semmi se elérhető nyilvánosan, szóval a modellben egy BM-21 Grad által kilőtt rakéta hozzávetőleges adatait használtam.  Ezzel az eszközzel szintén rendelkezik a Hamasz, és a Vaskupola első éles tesztjét is ilyen típussal végezték.  A feltételezett rakéta adatai:

  • Kaliber: 122 mm
  • Kezdeti tömeg: 67 kg
  • Kiéget tömeg: 46 kg
  • Égés időtartalma: 1,7 s
  • Torkolati sebesség: 690 m/s
  • Légellenállási együtthatója: 0,2

És most nézzük a konkrét kódot.

Modellezés

Először is építenünk kell egy modellt, amit kezdjük az állandókkal:

from math import radians, pi

# állandok
#   kiéget tömeg kg
m = 46
#   átmérő, m
d = .122
# alak ellenállása (Drag Coefficient)
Cd = .2
# légürüsség kg/m^3
rho= 1.2
# gravitáció
g = 9.8
# radar pontossága
accuracy = 125
# menyi időnk van arra, hogy eldöntsük cselekedünk kell-e, másodperc
hatarido = 10
# mérési időpillanatok közt eltelt idő
delta_t = 0.1
# prior
#   feltételezett kezdő sebesség
v_b = 1000
#   feltételezett kilövési szög
theta_b = radians(45)
#   valós sebesség a kiégés pillanatában
v = 1000
#   valós kilövési szög
theta = radians(50)


# Számított állandók
#   kertszemecet területe
A = pi*(d**2)
#   terminal vericity
D = rho*Cd*A/2
#   menyi mérési pont van a határidőn belül, lefelé kerekítve
deadline = int(hatarido/delta_t)

Most pedig készítsük el a szimulált radarjelet.  Ehhez először modellezzük a pályát; ez lesz a „valós pálya”; majd zajt adunk az x és z tengelyekhez :

from math import sqrt
import numpy as np

def statespacemodel(X, dt):
    # pozíciók az adott időpillanatban
    x_t = []
    z_t = []
    v_x = []
    v_z = []    
    
    A = np.array( [ [1,0,delta_t,0], [0,1,0,delta_t], [0,0,1,0], [0,0,0,1]] )
    B = np.array( [ [ 0.5*delta_t**2, 0], [ 0, 0.5*delta_t**2], [delta_t, 0], [0, delta_t] ] )
    
    while True:
        U = np.array([
            [ -1.*(D/m) * X[2][0] * sqrt( X[2][0]**2 + X[3][0]**2 ) ], 
            [ -g-(D/m) * X[3][0] * sqrt( X[2][0]**2 + X[3][0]**2 ) ]]
        )
        X = A.dot(X) + B.dot(U)
    
        # mentjük a szimulált poziciót
        if X[1][0] >= 0:
            x_t.append( X[0][0] )
            z_t.append( X[1][0] )
            v_x.append( X[2][0] )
            v_z.append( X[3][0] )
        else:
            break
    
    return(x_t,z_t,v_x,v_z)

# modellezés
x_t,z_t,v_x,v_z = statespacemodel( np.array( [ [0], [1345], [ v * np.cos(theta) ], [ v * np.sin(theta) ]  ] ), delta_t )

# zaj hozzáadása
zaj = np.random.normal(0,125,len(x_t))
xt_z = np.add(x_t, zaj)
zaj = np.random.normal(0,125,len(z_t))
zt_z = np.add(z_t, zaj)

Az eredmény pedig valahogy így fog kinézni:

A megvalósítás

Először készítsük egy függvényt az ismert modellről. Ezt már korábban megtettem. Vegyük észre, hogy lényegében csak ezt a függvényt kell átírni, ha a modellünk változik:

# a modellhez használt függvény
def fx(X, dt):
    # Újraformázzuk a mátrixot mivel a saját state space megvalósításuk ezt a formátumot használja
    X = X.reshape(4,1)
    # a state space modell amit ismerünk
    #  az A-t és a B-t ki lehetne mozgatni mivel állandók, 
    #  de itt hagyom, hogy együtt legyen látható az egész modell
    A = np.array( [ [1,0,delta_t,0], [0,1,0,delta_t], [0,0,1,0], [0,0,0,1]] )
    B = np.array( [ [ 0.5*delta_t**2, 0], [ 0, 0.5*delta_t**2], [delta_t, 0], [0, delta_t] ] )
    U = np.array([
            [ -1.*(D/m) * X[2][0] * sqrt( X[2][0]**2 + X[3]**2 ) ], 
            [ -g-(D/m) * X[3][0] * sqrt( X[2][0]**2 + X[3]**2 ) ]]
                )
    X = A.dot(X) + B.dot(U)
    # vissza formázzuk az eredményt
    return X.flatten()

Még egy függvényt kell készítenünk, ami a megfigyelt változókat az előbb létrehozott modell State Mátrix megfelelő részére linkeli. Esetünkben ez egyszerű, mivel csak két megfigyelt változónk van: a magasság és a távolság.  Ezek pedig a modell 1 és 0 elemei. Szóval:

# mérési függvény, lényegében a mért értékeket az X mátrix megfelelő pontjára illeszti
def hx(x):
    return np.array([x[0], x[1]])

Ezután állítsuk be a mintavétel paramétereit. A MerweScaledSigmaPoints mintavételt fogjuk alkalmazni.3 Mint általában minden hiperparaméter beállítás esetén, ajánlott ezeket optimalizálni pl. “grid search” vagy valamilyen Bayesian optimalizációs eljárással.  Ettől most eltekintek, mert az α, β és κ értékeknek egy konzervatív4 beállítása is elég a szűrő hatékonyságának általános szemléltetéshez.  Szintén elég konzervatív a 4 elemes mintavétel ciklusonként. Ezt növelhetjük vagy csökkenthetjük az igényeiknek megfelelően.  Értelemszerűen, minél több mintát veszünk annál nagyobb a számítási igény, de annál pontosabbak is leszünk:

# ezek a legáltalánosabb eredmények Normál zajhoz
from filterpy.kalman import MerweScaledSigmaPoints
sigmas = MerweScaledSigmaPoints(4, alpha=.1, beta=2., kappa=0)

Most hozzuk létre a Kálmán szűrő objektumát. Ehhez a következő változókat kell beállítani:

  • dim_x — a modellünkben használt X nagysága.  Esetünkben ez a x és z pozíció és sebességek, szóval 4. 
  • dim_z — mennyi valós megfigyelés áll rendelkezésünkre, ami 2, mert a radarkép csak az x és z pozíciókról szolgáltat információt , a sebességről nem.
  • dt — a Δt, vagyis az egyes radarmegfigyelések közt eltelt idő
  • hx — a mért értékeket az X mátrix megfelelő pontjára illesztő függvény
  • fx — a state-space modellünk
  • points –a már említett mintavételi eljárás
# a Kálmán szűrő maga
from filterpy.kalman import UnscentedKalmanFilter
ukf = UnscentedKalmanFilter(dim_x=4, dim_z=2, dt=delta_t, hx=hx, fx=fx, points=sigmas)

Majd állítsunk be néhány egyéb tulajdonságát:  

  • Az R azt hívatott érzékeltetni, ennyire zajos a mintavétel.  Ehhez szerencsére a radar ismert zajszintje jó támpontot szolgáltat. 
  • A Q a rendszer stabilitását írja le, azt, mekkora az érték változása az egyes mérési időpontok között.
  • A P pedig a Covariance Estimate Mátrix.
# Covariance Estimate Mátrix
ukf.P *= 200
# mennyire zajosak a mérések
ukf.R *= np.array([accuracy,accuracy,sqrt(2*(accuracy**2)),sqrt(2*(accuracy**2))]).reshape(2,2)
# menyit változik a rendszer az egyes időpillanatok között
from filterpy.common import Q_discrete_white_noise
ukf.Q = Q_discrete_white_noise(2, dt=delta_t, var=0.03, block_size=2)

Ezek után lépkedjünk végig a megfigyelésünkön és végezzük el a számítást:


# a Kálmán szűrő által számított eredmény
#   pálya adatok:
zp = []
#   földet érés 
f = []

# 1 ponttól indulunk mivel megfigyelőként nem tudható a 0.
for i in range(1,len(xt_z)):
    # aktuális helyzet
    z = np.array([ xt_z[i], zt_z[i] ])
    ukf.predict()
    # prior 
    if i == 1:
        ukf.x = [ukf.x[0], ukf.x[1], 1000*cos(radians(45)), 1000*sin(radians(45))]
    ukf.update(z)
    zp.append(ukf.x)
    
    # földet élés előrejelezése csak az deadline előtt
    if i <= deadline:
        x_f, z_f, _ , _ = statespacemodel(ukf.x.reshape(4,1), delta_t) 
        if len(x_f) > 0:
            f.append(x_f[-1])
           
fs.append(f)
zps.append(zp)

Az eredmény

Először is nézzük meg, mi a helyzet a célpont előrejelzéssel.  Az alábbi ábrán látható, hogyan alakul a Kálmán szűrő előrejelzése a célpontra nézve az idő függvényében.  Négy független szimuláció eredményét is felrajzoltam, hogy jobban kirajzolódjon a trend.  Kb. 6 másodperc után olyan 200 m-es biztonsággal meg tudjuk határozni a várható becsapódás pontját és a szűrő nagyjából stabilizálódik.  Ez bőven a 10 másodperces határon belül van.

Ha lefuttatjuk több százszor a szimulációt akkor a következő Sűrűségfüggvényt kapjuk a célpontra nézve a 10. másodperc végén:

Az előrejelzések 90% nagyjából ±100 méteres távolságra van a tényleges becsapódási pontól.  Úgy gondolom, ennek elégnek kell lennie arra, hogy eldöntsük a rakéta lakott területet veszélyeztet e?

Ezután nézzük meg, hogy alakul a pálya előrejelzése csak egy Δt időben előre:

Ez elsőre jónak tűnik. Az idő előrehaladtával egyre kisebb a szűrő által számított és a valós pálya különbsége. Részletesebben: az előrejelzési hiba a z tengely mentén az idő múlásával így alakul:

Az x tengely mentén pedig:


50 másodperc után a medián és a valós pálya kb. egybe esik, és 90% az előrejelzéseknek ±25 m-en távolságon belüli a z tengelyen (zöld terület a fenti ábrán), és kb. ±15m-en belül az x tengely mentén. Ami jelentős előrelépésnek tűnik a piros színnel jelölt radar adatokhoz képest.

Összegzés

Mint fentebb látható a Kálmán szűrő bizonyos implementációival képesek vagyunk gyorsan, valós időben kezelni a a zajos megfigyelésekből eredő problémákat nem lineáris rendszerek esetén is. Viszont ehhez rendelkeznünk kell a környezet egy megbízható state-space modelljével.


Lábjegyzet

  1. Úgy tudom, a Vaskupola első lépésben a rakéta típusát próbálja meghatározni, és csak az így becsült kiégési idő után próbálja kalkulálni a pályaadatokat.  Ez lényegében leegyszerűsíti a pályaszámítást.
  2. Tudja valaki, hogy milyen eloszlást követ egy radar rendszer hibája?
  3. A filterpy könyvtár másik lehetősége a JulierSigmaPoints.  A MerweScaledSigmaPoints eljárás Van der Merwe tanulmánya alapján készült, míg a JulierSigmaPoints padig S. Julier cikke nyomán.
  4. Ha valakit jobban érdekel a téma lsd. Sebastian Bitzer: The UKF exposed
Hírdetés

Vélemény, hozzászólás?

Adatok megadása vagy bejelentkezés valamelyik ikonnal:

WordPress.com Logo

Hozzászólhat a WordPress.com felhasználói fiók használatával. Kilépés /  Módosítás )

Twitter kép

Hozzászólhat a Twitter felhasználói fiók használatával. Kilépés /  Módosítás )

Facebook kép

Hozzászólhat a Facebook felhasználói fiók használatával. Kilépés /  Módosítás )

Kapcsolódás: %s