264 KiB
264 KiB
🎯 [OPS02] 다중 센서 경보 중복 항적 통합 시뮬레이터
Google Colab 실습 환경 (작전·전투 수행 / Level 3)
📌 1. 프로젝트 개요 및 배경
| 항목 | 내용 |
|---|---|
| 과정 분류 | 작전·전투 수행 (Level 3) |
| 학습 소요 시간 | 10개 미션 (5일 / 총 27.5시간) |
| 환경 사양 | Python 3.10+, numpy, pandas, matplotlib, scipy, scikit-learn |
📖 프로젝트 목표
미국-이란 간 무력 충돌 사례에서는 드론·미사일이 여러 방향에서 동시에 들어오면서 레이더, 광학(EO/IR), 음향 센서가 같은 비행체를 각자 경보해 상황판에 중복 항적이 쏟아지고, 이 때문에 위협 수를 잘못 세거나 대응이 늦어지는 문제가 주목받았다. 본 프로젝트에서는 가상 방어 지역에 배치한 가상 센서 4기와 가상 비행체 궤적으로 교육용 경보 데이터를 코드로 만들어, 연동 규격과 통합 스키마 설계, 시각 동기화와 ENU 좌표 변환 전처리, 최단거리 매칭 Baseline, 확장 칼만 필터(EKF) 항적 추적·통합, 융합 신뢰도 지수, 전술 상황판(COP), 조작관 HITL 검토, 통신 단절·재밍 스트레스 시험, AAR 재보정, 최종 통합까지 다중 센서 중복 항적 통합 시뮬레이터를 단계별로 구현한다. 최종 목표는 비행체 하나를 하나의 항적으로 정확히 보여 주어 지휘관의 인지 부하를 줄이고, 항적 통합·분리 같은 판단은 조작관이 근거와 함께 확정하도록 설계하는 것이다.
⚠️ [교육용 가상 데이터셋 원칙]
- 본 노트북의 모든 데이터(제원, 로그, 시나리오 등)는 보안 규정을 준수하여 코드 내 난수 기반 가상 데이터로 생성됩니다.
- 실제 부대 위치, 작전 인가, 실제 전장 정보, 군사 기밀 데이터는 일체 포함되지 않습니다.
In [ ]:
# ============================================================================
# [Google Colab 환경 설정] 필수 라이브러리 설치 및 임포트
# ============================================================================
!pip install -q numpy pandas matplotlib scipy scikit-learn
import os
import sys
import time
import math
import json
import re
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# 그래프 영문 폰트 기본 설정 (환경 독립성 보장)
plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
print(f"실행 환경 준비 완료 (Python {sys.version.split()[0]})")
In [ ]:
"""레이더·EO/IR·음향 센서 경보 메시지와 통합 항적(Fused Track)의 데이터 스키마를 정의하고, 가상 경보 샘플 20건을 검증하며 오류 주입 시험으로 검증기를 확인한다."""
import json
import re
import numpy as np
import pandas as pd
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 1] 다종 센서(레이더·EO/IR·음향) 연동 규격 분석 및 중복 항적 데이터셋 스키마 설계 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# 1) 스키마 정의: 센서 경보(원시), 전술 파라미터, 통합 항적(출력)
SCHEMA = {
"report": {
"msg_id": {"type": str, "pattern": r"^M\d{5}$"},
"sensor": {"type": str, "enum": list(SENSORS)},
"stype": {"type": str, "enum": ["radar", "eoir", "acoustic"]},
"t_sensor": {"type": (int, float), "min": -5, "max": 86400},
"range_km": {"type": (int, float), "min": 0, "max": 60},
"az_deg": {"type": (int, float), "min": 0, "max": 360},
"sr": {"type": (int, float), "min": 0.001, "max": 5},
"saz": {"type": (int, float), "min": 0.01, "max": 10},
},
"tactical": {
"weather": {"type": str, "enum": list(WEATHER)},
"position_error_m": {"type": (int, float), "min": 0, "max": 500},
"comm_delay_ms": {"type": (int, float), "min": 0, "max": 5000},
},
"fused_track": {
"track_id": {"type": str, "pattern": r"^FT\d{3}$"},
"state_enu": {"type": list, "length": 4},
"sensors": {"type": list, "min_len": 1},
"confidence": {"type": (int, float), "min": 0, "max": 1},
"status": {"type": str, "enum": ["TENTATIVE", "CONFIRMED", "HOLD", "SPLIT", "MERGED"]},
"source_msgs": {"type": list, "min_len": 1},
},
}
def check(name, value, rule):
# [TODO : 학생 실습 1] 다종 센서(레이더·EO/IR·음향) 연동 규격 분석 및 중복 항적 데이터셋 스키마 설계 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def validate(rec, group):
errs = []
for k, rule in SCHEMA[group].items():
errs += [f"{k}: 필수 필드 누락"] if k not in rec else check(k, rec[k], rule)
if group == "report" and not errs and SENSORS[rec["sensor"]]["type"] != rec["stype"]:
errs.append(f"stype: 센서 {rec['sensor']}의 유형과 불일치")
return errs
# 2) [교육용가상자료] 가상 시나리오에서 경보 샘플 20건 추출
scn = make_scenario(seed=1, n=4, T=10)
rng = np.random.default_rng(1)
samples = [{k: r[k] for k in SCHEMA["report"]} for r in rng.choice(scn["reports"], 20, replace=False)]
tactical = {"weather": "clear", "position_error_m": 35.0, "comm_delay_ms": 120}
track = {"track_id": "FT001", "state_enu": [3.2, 4.1, 0.05, -0.02], "sensors": ["RDR-A", "EOIR-C"],
"confidence": 0.87, "status": "CONFIRMED", "source_msgs": ["M00012", "M00015"]}
n_err = sum(len(validate(s, "report")) for s in samples) + len(validate(tactical, "tactical")) + len(validate(track, "fused_track"))
print(f"경보 샘플 {len(samples)}건 + 전술 파라미터 1건 + 통합 항적 1건 검증: 오류 {n_err}건")
# 3) 오류 주입 시험
bad = json.loads(json.dumps(samples[0]))
bad["msg_id"] = "MSG-1" # 형식 오류
bad["sensor"] = "RDR-Z" # 등록되지 않은 센서
bad["az_deg"] = 412.0 # 범위 오류
del bad["sr"] # 필수 필드 누락
bad_trk = dict(track, sensors=[], status="AUTO_FIRE") # 센서 없는 항적, 허용되지 않은 상태
print("\n[오류 주입 결과]")
for e in validate(bad, "report") + validate(bad_trk, "fused_track"):
print(" -", e)
# 4) 명세서 요약: 필드 정의 + 센서 등록부
rows = []
for g, fields in SCHEMA.items():
for k, r in fields.items():
t = r["type"].__name__ if isinstance(r["type"], type) else "number"
rows.append((g, k, t, ", ".join(f"{x}={y}" for x, y in r.items() if x != "type")))
spec = pd.DataFrame(rows, columns=["그룹", "필드", "자료형", "규칙"])
print(f"\n[데이터 명세 요약] 필드 {len(spec)}개 정의 (정의 완성도 {spec['규칙'].str.len().gt(0).mean():.0%})")
print(spec.to_string(index=False))
reg = pd.DataFrame(SENSORS).T[["type", "pos", "max_r", "sr", "saz", "pd"]]
reg.columns = ["유형", "위치(ENU km)", "최대거리(km)", "거리오차(km)", "방위오차(도)", "탐지확률"]
print("\n[가상 센서 등록부]")
print(reg.to_string())
# 5) 중복 경보 현황: 같은 시각 같은 비행체를 몇 개 센서가 동시에 보고하는가
df = pd.DataFrame(scn["reports"])
df["sec"] = df["t_sensor"].round().astype(int)
dup = df[df["truth"] != "FA"].groupby(["truth", "sec"])["sensor"].nunique()
print(f"\n비행체 1개당 같은 초에 평균 {dup.mean():.1f}개 센서가 경보(최대 {dup.max()}개) -> 통합 항적 스키마가 필요한 이유")
In [ ]:
"""이종 센서 경보의 정제(중복·결측·반사파 제거), 시각 동기화(Time-sync) 오차 보정, 센서 국소 극좌표 -> 공통 ENU 좌표 변환(가상 기준점의 위경도·ECEF 변환 포함)을 구현하고 처리 속도·정확도를 단위 시험으로 확인한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 2] 이종 센서 시각 동기화(Time-sync) 및 공간 좌표계(ENU/ECEF) 변환 전처리 구현 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 2] 이종 센서 시각 동기화(Time-sync) 및 공간 좌표계(ENU/ECEF) 변환 전처리 구현 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
# 1) 위경도(LLA) <-> ECEF <-> ENU 변환 (WGS84 타원체 공식)
A, E2 = 6378137.0, 6.69437999014e-3
ORIGIN = (10.0, 20.0, 0.0) # [교육용가상자료] 가상 기준점 위경도(실제 위치 아님)
def lla_to_ecef(lat, lon, h):
la, lo = np.radians(lat), np.radians(lon)
n = A / np.sqrt(1 - E2 * np.sin(la) ** 2)
return np.array([(n + h) * np.cos(la) * np.cos(lo), (n + h) * np.cos(la) * np.sin(lo), (n * (1 - E2) + h) * np.sin(la)])
def _rot(lat, lon):
la, lo = np.radians(lat), np.radians(lon)
return np.array([[-np.sin(lo), np.cos(lo), 0],
[-np.sin(la) * np.cos(lo), -np.sin(la) * np.sin(lo), np.cos(la)],
[np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)]])
def ecef_to_enu(p, origin=ORIGIN):
return _rot(*origin[:2]) @ (p - lla_to_ecef(*origin))
def enu_to_ecef(enu, origin=ORIGIN):
return _rot(*origin[:2]).T @ enu + lla_to_ecef(*origin)
def ecef_to_lla(p):
lon = np.arctan2(p[1], p[0]); r = np.hypot(p[0], p[1]); lat = np.arctan2(p[2], r * (1 - E2))
for _ in range(5):
n = A / np.sqrt(1 - E2 * np.sin(lat) ** 2)
h = r / np.cos(lat) - n
lat = np.arctan2(p[2], r * (1 - E2 * n / (n + h)))
return np.degrees(lat), np.degrees(lon), h
print("[센서 설치 지점(가상 기준점 기준, 실제 위치 아님): 가상 ENU -> 위경도 -> ECEF -> ENU 왕복 변환]")
rt = []
for sid, s in SENSORS.items():
enu = np.array([s["pos"][0] * 1000, s["pos"][1] * 1000, 30.0])
lla = ecef_to_lla(enu_to_ecef(enu))
back = ecef_to_enu(lla_to_ecef(*lla))
rt.append(np.linalg.norm(back - enu))
print(f"{sid:7s} 위도 {lla[0]:.5f} 경도 {lla[1]:.5f} 고도 {lla[2]:6.1f} m 왕복 오차 {rt[-1] * 1000:.3f} mm")
# 2) [교육용가상자료] 원시 경보 + 전술 잡음 주입(재전송 중복, 결측, 반사파 거리값, 오차값 누락)
scn = make_scenario(seed=2, n=8, T=60, special="fast")
rng = np.random.default_rng(2)
raw = [dict(r) for r in scn["reports"]]
n0 = len(raw)
raw += [dict(raw[i]) for i in rng.choice(n0, int(0.03 * n0), replace=False)] # 재전송 중복 3%
for i in rng.choice(n0, int(0.02 * n0), replace=False):
raw[i]["az_deg"] = None # 측정값 결측 2%
for i in rng.choice(n0, int(0.01 * n0), replace=False):
raw[i]["range_km"] = SENSORS[raw[i]["sensor"]]["max_r"] * rng.uniform(1.5, 3) # 반사파(다중경로) 1%
miss_sigma = rng.choice(n0, int(0.05 * n0), replace=False)
for i in miss_sigma:
raw[i]["sr"] = None # 오차값 누락 5% -> 등록부 값으로 보정
rng.shuffle(raw) # 도착 순서 뒤섞임
t0 = time.perf_counter()
filled = 0
for r in raw:
if r["sr"] is None:
r["sr"] = SENSORS[r["sensor"]]["sr"]; filled += 1
cleaned, why = clean(raw)
off = estimate_offsets(scn["heartbeats"])
enu = sorted((to_enu(r, off) for r in cleaned), key=lambda r: r["t"])
dt = time.perf_counter() - t0
print(f"\n원시 경보 {len(raw)}건 -> 정제 {len(enu)}건, 제거 사유 {why}, 오차값 보정 {filled}건(보정률 {filled / len(miss_sigma):.0%})")
print(f"처리 속도 {len(raw) / dt:,.0f} 건/초 (목표 1,000건/초 이상: {'달성' if len(raw) / dt >= 1000 else '미달'})")
# 3) 시각 동기화: 기준 신호로 추정한 센서 시계 오차
print("\n센서 실제 시계오차(s) 추정값(s) 추정 오차(s)")
for sid, s in SENSORS.items():
print(f"{sid:7s} {s['clock']:14.3f} {off[sid]:10.3f} {abs(off[sid] - s['clock']):11.3f}")
def along_bias(reps, use_t):
"""고속 표적 경보가 정답 위치보다 진행 방향으로 얼마나 앞서거나 뒤처지는지(m, 평균). 시계 오차 x 속도만큼 생긴다."""
e = []
for r in reps:
if scn["kinds"].get(r["truth"]) == "fast":
tr, ts = scn["targets"][r["truth"]], np.arange(scn["T"])
g = np.array([np.interp(r[use_t], ts, tr[:, j]) for j in range(4)])
v = g[2:] / np.linalg.norm(g[2:])
e.append(np.dot([r["x"] - g[0], r["y"] - g[1]], v))
return float(np.mean(e)) * 1000
print("\n[고속 표적 경보의 진행 방향 위치 편향: 시계 오차 x 속도만큼 어긋난다]")
print("센서 동기화 전(m) 동기화 후(m)")
bias = {}
for sid in SENSORS:
rs = [r for r in enu if r["sensor"] == sid]
bias[sid] = (along_bias(rs, "t_sensor"), along_bias(rs, "t"))
print(f"{sid:7s} {bias[sid][0]:+11.0f} {bias[sid][1]:+12.0f}")
# 4) 단위 시험 5종 (pytest 없이 assert로 수행)
tests = {
"ENU 왕복 변환 오차 1 mm 미만": max(rt) < 1e-3,
"방위 0도 = 정북(+y) 방향": np.allclose([to_enu({"sensor": "RDR-A", "t_sensor": 0, "range_km": 5, "az_deg": 0}, {})[k] for k in "xy"], [0, 5]),
"중복 메시지 전량 제거": len({r["msg_id"] for r in enu}) == len(enu),
"반사파(최대거리 초과) 전량 제거": all(r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1 for r in enu),
"시계오차 추정 오차 0.05초 미만": all(abs(off[s] - SENSORS[s]["clock"]) < 0.05 for s in SENSORS),
"출력 시각 오름차순 정렬": all(a["t"] <= b["t"] for a, b in zip(enu, enu[1:])),
}
for k, v in tests.items():
print(f"[{'PASS' if v else 'FAIL'}] {k}")
print(f"단위 시험 {sum(tests.values())}/{len(tests)} 통과, 전처리 오류율 {1 - sum(tests.values()) / len(tests):.0%}")
# 5) 그림: 센서별 진행 방향 위치 편향(동기화 전/후)
sids = list(SENSORS)
x = np.arange(len(sids))
plt.figure(figsize=(7, 4))
plt.bar(x - 0.2, [bias[s][0] for s in sids], 0.4, label="before time-sync", color="tab:gray")
plt.bar(x + 0.2, [bias[s][1] for s in sids], 0.4, label="after time-sync", color="tab:blue")
plt.axhline(0, color="k", lw=0.8)
plt.xticks(x, sids); plt.ylabel("along-track bias (m)"); plt.legend()
plt.title("Fast-target along-track bias by sensor (synthetic)")
plt.tight_layout(); plt.show()
In [ ]:
"""매 스캔(1초)마다 여러 센서 경보를 최단거리(유클리드/마할라노비스) 매칭으로 묶어 비행체 1개당 1개의 융합 경보로 만드는 1차 항적 융합(Fusion) Baseline을 구현하고, 통합 정확도와 한계를 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 3] 최단거리 매칭 기반 단일 항적-센서 경보 1차 융합(Fusion) 알고리즘 구현 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 3] 최단거리 매칭 기반 단일 항적-센서 경보 1차 융합(Fusion) 알고리즘 구현 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
def fuse_scan(reps, gate, metric):
"""정밀한 센서 경보부터 차례로, 같은 센서가 아직 없는 가장 가까운 묶음에 넣는다(게이트 밖이면 새 묶음).
묶음 위치는 측정 공분산의 역수(정보량) 가중 평균으로 갱신한다."""
clusters = []
for r in sorted(reps, key=lambda r: SENSORS[r["sensor"]]["sr"]):
p, C = np.array([r["x"], r["y"]]), meas_cov(r)
best, bd = None, gate
for c in clusters:
if r["sensor"] in c["sensors"]:
continue
d = p - c["pos"]
dist = np.hypot(*d) if metric == "euclid" else np.sqrt(d @ np.linalg.solve(C + c["cov"], d))
if dist < bd:
best, bd = c, dist
if best is None:
clusters.append({"pos": p, "cov": C, "sensors": {r["sensor"]}, "members": [r], "dist": []})
else:
info = np.linalg.inv(best["cov"]) + np.linalg.inv(C)
best["pos"] = np.linalg.solve(info, np.linalg.solve(best["cov"], best["pos"]) + np.linalg.solve(C, p))
best["cov"] = np.linalg.inv(info)
best["sensors"].add(r["sensor"]); best["members"].append(r); best["dist"].append(round(float(bd), 2))
return clusters
def evaluate_scans(reps, gate, metric):
"""통합 정확도 = (스캔, 실제 비행체) 중 정확히 1개의 순수한 묶음으로 통합된 비율."""
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
ok = tot = n_cl = n_multi_fa = 0
t0 = time.perf_counter()
for k, rs in by.items():
cl = fuse_scan(rs, gate, metric)
n_cl += len(cl)
n_multi_fa += sum(1 for c in cl if len(c["sensors"]) >= 2 and all(m["truth"] == "FA" for m in c["members"]))
for g in {r["truth"] for r in rs if r["truth"] != "FA"}:
hold = [c for c in cl if any(m["truth"] == g for m in c["members"])]
tot += 1
ok += len(hold) == 1 and all(m["truth"] == g for m in hold[0]["members"])
ms = (time.perf_counter() - t0) * 1000 / len(by)
return {"경보": len(reps), "융합 묶음": n_cl, "통합 정확도": round(ok / tot, 3), "다중센서 거짓묶음": n_multi_fa, "스캔당 ms": round(ms, 2)}
# [교육용가상자료] 기상·기동이 다른 가상 시나리오 4개
scns = [make_scenario(30 + i, n=6, weather=w, special=sp) for i, (w, sp) in enumerate([("clear", None), ("rain", None), ("clear", "turn"), ("fog", "formation")])]
reps = [prepare(s) for s in scns]
rows = []
for name, gate, metric in [("유클리드 1.0km", 1.0, "euclid"), ("마할라노비스 3.0", 3.0, "mahal")]:
for s, r in zip(scns, reps):
rows.append({"방식": name, "기상": s["weather"], **evaluate_scans(r, gate, metric)})
df = pd.DataFrame(rows)
print(df.to_string(index=False))
summ = df.groupby("방식")[["통합 정확도", "스캔당 ms"]].mean().round(3)
print("\n" + summ.to_string())
best = summ["통합 정확도"].max()
print(f"\n목표(통합 정확도 75% 이상) 달성 여부: {'달성' if best >= 0.75 else '미달'} (최고 {best:.1%})")
# 게이트 크기 민감도(맑음 시나리오)
print("\n[게이트 민감도] 작으면 같은 비행체가 쪼개지고, 크면 다른 비행체·오경보가 섞인다")
for metric, gates in [("euclid", [0.3, 0.6, 1.0, 1.5, 2.5]), ("mahal", [1.5, 2.0, 3.0, 4.0, 6.0])]:
print(metric, " ".join(f"{g}:{evaluate_scans(reps[0], g, metric)['통합 정확도']:.3f}" for g in gates))
# 시각 동기화 효과
e0 = evaluate_scans(prepare(scns[2], sync=False), 3.0, "mahal")["통합 정확도"]
print(f"\n시각 동기화 생략 시(기동 시나리오) 통합 정확도 {e0:.3f} -> 동기화 적용 {df.iloc[6]['통합 정확도']:.3f}")
# 판단 근거(투명성): 한 스캔의 융합 결과
k = 25
cl = fuse_scan([r for r in reps[0] if int(round(r["t"])) == k], 3.0, "mahal")
print(f"\n[스캔 {k}초 융합 결과 예시: 묶음 {len(cl)}개]")
for c in sorted(cl, key=lambda c: -len(c["sensors"]))[:5]:
print(f"위치({c['pos'][0]:6.2f},{c['pos'][1]:6.2f}) 센서 {sorted(c['sensors'])} 결합 거리 {c['dist']} 정답 {sorted({m['truth'] for m in c['members']})}")
print("한계: 스캔마다 독립적으로 묶어 속도 정보를 쓰지 않으므로, 가까운 편대·교차 표적과 단일 센서 오경보를 구분하기 어렵다(미션 4에서 칼만 필터로 보완).")
# 그림: 원시 경보 vs 융합 묶음
fig, ax = plt.subplots(1, 2, figsize=(11, 5), sharex=True, sharey=True)
colors = dict(zip(SENSORS, ["tab:blue", "tab:green", "tab:orange", "tab:red"]))
for r in [r for r in reps[0] if int(round(r["t"])) == k]:
ax[0].scatter(r["x"], r["y"], color=colors[r["sensor"]], s=18)
for sid, c in colors.items():
ax[0].scatter([], [], color=c, label=sid)
for c in cl:
ax[1].scatter(*c["pos"], s=40 + 40 * len(c["sensors"]), color="tab:purple" if len(c["sensors"]) > 1 else "lightgray", edgecolor="k")
for g, tr in scns[0]["targets"].items():
ax[1].plot(tr[k, 0], tr[k, 1], "k+", ms=10)
ax[0].set_title(f"Raw alerts at t={k}s ({sum(1 for r in reps[0] if int(round(r['t'])) == k)})"); ax[0].legend(fontsize=8)
ax[1].set_title(f"Fused clusters ({len(cl)}), + = truth, size = #sensors")
for a in ax:
a.set_xlabel("East (km)"); a.set_ylabel("North (km)")
plt.tight_layout(); plt.show()
In [ ]:
"""확장 칼만 필터(EKF)와 전역 최근접 할당(GNN)으로 여러 센서 경보를 하나의 항적으로 추적·통합하는 추적기를 고도화하고, 선회·교차·편대·고속·기만 항적·재밍·통신 단절 등 복합 위협 시나리오 10종에서 재현율·처리 지연·거짓 항적 패턴을 검증한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 4] 칼만 필터(Kalman Filter) 기반 다중 센서 경보 중복 항적 추적 및 궤적 통합 고도화 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 4] 칼만 필터(Kalman Filter) 기반 다중 센서 경보 중복 항적 추적 및 궤적 통합 고도화 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
# ---- 확장 칼만 필터(EKF) 다중 센서 항적 추적기 공통 함수 ----
# 상태 x = [동(km), 북(km), 동속(km/s), 북속(km/s)], 측정 z = 센서 기준 [거리(km), 방위(rad)]
def ekf_predict(x, P, dt, q):
F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], float)
G = np.array([[dt ** 2 / 2, 0], [0, dt ** 2 / 2], [dt, 0], [0, dt]])
return F @ x, F @ P @ F.T + q * G @ G.T
def h_polar(x, spos):
dx, dy = x[0] - spos[0], x[1] - spos[1]
r2 = max(dx * dx + dy * dy, 1e-9); r = np.sqrt(r2)
return np.array([r, np.arctan2(dx, dy)]), np.array([[dx / r, dy / r, 0, 0], [dy / r2, -dx / r2, 0, 0]])
def innovation(x, P, rep):
s = SENSORS[rep["sensor"]]
zp, H = h_polar(x, s["pos"])
y = np.array([rep["range_km"], np.radians(rep["az_deg"])]) - zp
y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi # 방위 차이는 -pi~pi로 감싼다
S = H @ P @ H.T + np.diag([rep["sr"] ** 2, np.radians(rep["saz"]) ** 2])
return y, H, S, float(y @ np.linalg.solve(S, y))
class Tracker:
"""센서 경보를 스캔(1초) 단위로 받아 EKF 항적을 갱신한다. 같은 스캔 안에서 센서별로
순차 갱신하므로, 여러 센서가 같은 비행체를 보면 하나의 항적으로 통합된다."""
def __init__(self, gate=9.21, q=4e-4, confirm=(3, 5), max_miss=4, merge_gate=6.0):
self.gate, self.q, self.confirm, self.max_miss, self.merge_gate = gate, q, confirm, max_miss, merge_gate
self.tracks, self.dead, self.nid, self.n_merged = [], [], 0, 0
def _new(self, rep, k):
C = meas_cov(rep)
P = np.zeros((4, 4)); P[:2, :2] = C; P[2, 2] = P[3, 3] = 0.25 ** 2
self.tracks.append({"id": f"FT{self.nid:03d}", "x": np.array([rep["x"], rep["y"], 0, 0]), "P": P, "k": k,
"scans": [k], "sensors": {rep["sensor"]: 1}, "truths": [rep["truth"]], "nis": [],
"hist": [(k, rep["x"], rep["y"])], "confirmed": False, "miss": 0, "born": k, "upd": True})
self.nid += 1
def _assign(self, tracks, rs, sid):
"""게이트(마할라노비스 거리) 안에서 전역 최근접 할당(GNN). 할당된 보고서 번호 집합을 돌려준다."""
if not tracks or not rs:
return set()
spos = SENSORS[sid]["pos"]
Z = np.array([[r["range_km"], np.radians(r["az_deg"])] for r in rs])
R = np.diag([rs[0]["sr"] ** 2, np.radians(rs[0]["saz"]) ** 2])
cost = np.full((len(tracks), len(rs)), 1e6)
for i, tr in enumerate(tracks):
zp, H = h_polar(tr["x"], spos)
Y = Z - zp
Y[:, 1] = (Y[:, 1] + np.pi) % (2 * np.pi) - np.pi
Si = np.linalg.inv(H @ tr["P"] @ H.T + R)
d2 = np.einsum("ij,jk,ik->i", Y, Si, Y)
cost[i] = np.where(d2 < self.gate, d2, 1e6)
used = set()
for i, j in zip(*linear_sum_assignment(cost)):
if cost[i, j] < self.gate:
tr, r = tracks[i], rs[j]
y, H, S, d2 = innovation(tr["x"], tr["P"], r)
K = tr["P"] @ H.T @ np.linalg.inv(S)
tr["x"] = tr["x"] + K @ y
tr["P"] = (np.eye(4) - K @ H) @ tr["P"]
tr["nis"].append(d2); tr["truths"].append(r["truth"])
tr["sensors"][sid] = tr["sensors"].get(sid, 0) + 1
tr["upd"] = True; used.add(j)
return used
def _merge(self):
"""같은 비행체를 따라가는 중복 항적(위치·속도가 통계적으로 같음)을 오래된 항적으로 통합한다."""
drop = set()
for a in range(len(self.tracks)):
for b in range(a + 1, len(self.tracks)):
A, B = self.tracks[a], self.tracks[b]
if A["id"] in drop or B["id"] in drop:
continue
d = A["x"] - B["x"]
d2 = float(d[:2] @ np.linalg.solve(A["P"][:2, :2] + B["P"][:2, :2], d[:2]))
if d2 < self.merge_gate and np.hypot(d[2], d[3]) < 0.05:
keep, lose = (A, B) if len(A["scans"]) >= len(B["scans"]) else (B, A)
for sid, n in lose["sensors"].items():
keep["sensors"][sid] = keep["sensors"].get(sid, 0) + n
keep["truths"] += lose["truths"]; keep["scans"] = sorted(set(keep["scans"]) | set(lose["scans"]))
keep["merged"] = keep.get("merged", 0) + 1
drop.add(lose["id"])
self.tracks = [t for t in self.tracks if t["id"] not in drop]
self.n_merged += len(drop)
def step(self, k, reps):
for tr in self.tracks: # 1) 예측
tr["x"], tr["P"] = ekf_predict(tr["x"], tr["P"], k - tr["k"], self.q)
tr["k"], tr["upd"] = k, False
for sid in sorted({r["sensor"] for r in reps}, key=lambda s: SENSORS[s]["sr"]): # 2) 정밀 센서부터 순차 갱신
rs = [r for r in reps if r["sensor"] == sid]
used = self._assign([t for t in self.tracks if t["confirmed"]], rs, sid) # 확정 항적 우선
rest = [r for j, r in enumerate(rs) if j not in used]
used2 = self._assign([t for t in self.tracks if not t["confirmed"]], rest, sid)
for j, r in enumerate(rest):
if j not in used2:
self._new(r, k)
if self.merge_gate > 0:
self._merge()
keep = []
for tr in self.tracks: # 3) 확정·삭제 관리
if tr["upd"]:
tr["miss"] = 0
if tr["scans"][-1] != k:
tr["scans"].append(k)
tr["hist"].append((k, tr["x"][0], tr["x"][1]))
else:
tr["miss"] += 1
m, n = self.confirm
if not tr["confirmed"] and sum(1 for s in tr["scans"] if s > k - n) >= m:
tr["confirmed"] = True
limit = self.max_miss if tr["confirmed"] else 1
(keep if tr["miss"] <= limit else self.dead).append(tr)
self.tracks = keep
def all_tracks(self):
return self.dead + self.tracks
def run_tracker(reps, **kw):
trk = Tracker(**kw)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
for k in range(min(by), max(by) + 1):
trk.step(k, by.get(k, []))
return trk.all_tracks()
def label_track(tr):
vals, cnt = np.unique(tr["truths"], return_counts=True)
return str(vals[np.argmax(cnt)]), float(cnt.max() / cnt.sum())
def evaluate(tracks, scn, min_reports=5):
"""정답 연결 평가: 재현율, 항적 분할(한 표적에 항적 2개 이상), 거짓 항적, 병합 오류, 위치 오차."""
conf = [t for t in tracks if t["confirmed"]]
for t in conf:
t["label"], t["purity"] = label_track(t)
cnt = pd.Series([r["truth"] for r in scn["reports"] if r["truth"].startswith("TGT")]).value_counts()
truths = [k for k, v in cnt.items() if v >= min_reports]
per = {g: [t for t in conf if t["label"] == g] for g in truths}
found = [g for g in truths if per[g]]
split = sum(max(0, len(v) - 1) for v in per.values())
false = [t for t in conf if not t["label"].startswith("TGT")]
merge = [t for t in conf if t["label"].startswith("TGT") and t["purity"] < 0.7]
err = []
for g in found:
for t in per[g]:
for k, x, y in t["hist"][3:]:
if 0 <= k < scn["T"]:
err.append(np.hypot(x - scn["targets"][g][k, 0], y - scn["targets"][g][k, 1]))
n = len(truths)
acc = len(found) / max(n + split + len(false) + len(merge), 1)
return {"정답표적": n, "탐지표적": len(found), "재현율": round(len(found) / max(n, 1), 3), "분할": split,
"거짓항적": len(false), "병합오류": len(merge), "위치RMSE(km)": round(float(np.sqrt(np.mean(np.square(err)))) if err else np.nan, 3),
"항적정확도": round(acc, 3)}
# [교육용가상자료] 복합 위협 시나리오 10종
SCEN = [("S01 기본(맑음)", dict(seed=401)), ("S02 선회 기동", dict(seed=402, special="turn")),
("S03 교차 표적", dict(seed=403, special="cross")), ("S04 근접 편대", dict(seed=404, special="formation")),
("S05 고속 표적", dict(seed=405, special="fast")), ("S06 기만 항적(EO/IR)", dict(seed=406, ghost=True)),
("S07 레이더A 재밍", dict(seed=407, noise_mult={"RDR-A": 4.0})), ("S08 레이더B 단절", dict(seed=408, dropout=("RDR-B", 15, 45))),
("S09 우천+오경보 증가", dict(seed=409, weather="rain", fa_mult=2.0)), ("S10 안개+선회+기만", dict(seed=410, weather="fog", special="turn", ghost=True))]
rows, lat_all, falses, trues = [], [], [], []
for name, kw in SCEN:
scn = make_scenario(**kw)
reps = prepare(scn)
t0 = time.perf_counter()
trk = run_tracker(reps)
ms = (time.perf_counter() - t0) * 1000 / scn["T"]
ev = evaluate(trk, scn)
lat_all.append(ms)
for t in trk:
if t["confirmed"] and not t["label"].startswith("TGT"):
falses.append({"시나리오": name.split()[0], "원인": t["label"], "센서수": len(t["sensors"]), "센서": "+".join(sorted(t["sensors"])),
"수명(s)": len(t["scans"]), "스캔당 보고": len(t["truths"]) / len(t["scans"])})
elif t["confirmed"]:
trues.append({"수명(s)": len(t["scans"]), "스캔당 보고": len(t["truths"]) / len(t["scans"]), "센서수": len(t["sensors"])})
rows.append({"시나리오": name, "표적": ev["정답표적"], "재현율": ev["재현율"], "분할": ev["분할"], "거짓항적": ev["거짓항적"],
"위치RMSE(km)": ev["위치RMSE(km)"], "ms/스캔": round(ms, 1), "통과": "PASS" if ev["재현율"] >= 0.85 else "FAIL"})
R = pd.DataFrame(rows)
print("[복합 위협 시나리오 10종 EKF 추적·통합 결과 - 교육용가상자료]")
print(R.to_string(index=False))
rec = R["재현율"].mean()
print(f"\n평균 재현율 {rec:.1%} (목표 85% 이상: {'달성' if rec >= 0.85 else '미달'}), 시나리오 통과 {(R['통과'] == 'PASS').sum()}/10")
print(f"스캔당 처리 지연 최대 {max(lat_all):.1f} ms (목표 200ms 이하: {'달성' if max(lat_all) <= 200 else '미달'})")
F = pd.DataFrame(falses)
print(f"\n[거짓 항적 패턴] 총 {len(F)}개")
if len(F):
print(F.groupby(["원인", "센서수"]).size().rename("개수").reset_index().to_string(index=False))
T_ = pd.DataFrame(trues)
print(f"실제 항적: 평균 수명 {T_['수명(s)'].mean():.0f}s, 스캔당 보고 {T_['스캔당 보고'].mean():.1f}건, 센서 {T_['센서수'].mean():.1f}종")
print(f"거짓 항적: 평균 수명 {F['수명(s)'].mean():.0f}s, 스캔당 보고 {F['스캔당 보고'].mean():.1f}건, 센서 {F['센서수'].mean():.1f}종")
print("오경보 유발 패턴: 우천·오경보 증가 시 클러터가 우연히 연속 연결된 짧은 항적, EO/IR 단독 기만 항적 -> 수명·스캔당 보고·센서 종류 수를 신뢰도 특징으로 추가(미션 5)")
# 기만 항적 시나리오에서 궤적 그림
scn = make_scenario(seed=410, weather="fog", special="turn", ghost=True)
trk = run_tracker(prepare(scn)); evaluate(trk, scn)
plt.figure(figsize=(6, 6))
for tid, tr in scn["targets"].items():
plt.plot(tr[:, 0], tr[:, 1], "k-", lw=0.8)
for t in trk:
if t["confirmed"]:
h = np.array(t["hist"])
plt.plot(h[:, 1], h[:, 2], "--", c="tab:red" if not t["label"].startswith("TGT") else "tab:blue", lw=1.2)
for sid, s in SENSORS.items():
plt.plot(*s["pos"], "^", c="tab:green"); plt.text(s["pos"][0] + 0.3, s["pos"][1], sid, fontsize=8)
plt.xlabel("East (km)"); plt.ylabel("North (km)"); plt.title("EKF fused tracks (blue) / false tracks (red) - synthetic")
plt.axis("equal"); plt.tight_layout(); plt.show()
In [ ]:
"""통합 항적마다 센서별 신뢰도 가중치(기상 반영), 항적 수명, 스캔당 보고 수, 센서 종류 수, 측정 정합도(NIS)를 특징으로 융합 신뢰도 지수와 High/Med/Low 등급을 산출해, 모든 확정 항적을 경보하는 방식 대비 오경보율 감축과 등급 일관성을 확인한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 5] 센서별 오탐 가중치 반영 및 항적 통합 신뢰도 지수 산출 엔진 개발 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 5] 센서별 오탐 가중치 반영 및 항적 통합 신뢰도 지수 산출 엔진 개발 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
# ---- 확장 칼만 필터(EKF) 다중 센서 항적 추적기 공통 함수 ----
# 상태 x = [동(km), 북(km), 동속(km/s), 북속(km/s)], 측정 z = 센서 기준 [거리(km), 방위(rad)]
def ekf_predict(x, P, dt, q):
F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], float)
G = np.array([[dt ** 2 / 2, 0], [0, dt ** 2 / 2], [dt, 0], [0, dt]])
return F @ x, F @ P @ F.T + q * G @ G.T
def h_polar(x, spos):
dx, dy = x[0] - spos[0], x[1] - spos[1]
r2 = max(dx * dx + dy * dy, 1e-9); r = np.sqrt(r2)
return np.array([r, np.arctan2(dx, dy)]), np.array([[dx / r, dy / r, 0, 0], [dy / r2, -dx / r2, 0, 0]])
def innovation(x, P, rep):
s = SENSORS[rep["sensor"]]
zp, H = h_polar(x, s["pos"])
y = np.array([rep["range_km"], np.radians(rep["az_deg"])]) - zp
y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi # 방위 차이는 -pi~pi로 감싼다
S = H @ P @ H.T + np.diag([rep["sr"] ** 2, np.radians(rep["saz"]) ** 2])
return y, H, S, float(y @ np.linalg.solve(S, y))
class Tracker:
"""센서 경보를 스캔(1초) 단위로 받아 EKF 항적을 갱신한다. 같은 스캔 안에서 센서별로
순차 갱신하므로, 여러 센서가 같은 비행체를 보면 하나의 항적으로 통합된다."""
def __init__(self, gate=9.21, q=4e-4, confirm=(3, 5), max_miss=4, merge_gate=6.0):
self.gate, self.q, self.confirm, self.max_miss, self.merge_gate = gate, q, confirm, max_miss, merge_gate
self.tracks, self.dead, self.nid, self.n_merged = [], [], 0, 0
def _new(self, rep, k):
C = meas_cov(rep)
P = np.zeros((4, 4)); P[:2, :2] = C; P[2, 2] = P[3, 3] = 0.25 ** 2
self.tracks.append({"id": f"FT{self.nid:03d}", "x": np.array([rep["x"], rep["y"], 0, 0]), "P": P, "k": k,
"scans": [k], "sensors": {rep["sensor"]: 1}, "truths": [rep["truth"]], "nis": [],
"hist": [(k, rep["x"], rep["y"])], "confirmed": False, "miss": 0, "born": k, "upd": True})
self.nid += 1
def _assign(self, tracks, rs, sid):
"""게이트(마할라노비스 거리) 안에서 전역 최근접 할당(GNN). 할당된 보고서 번호 집합을 돌려준다."""
if not tracks or not rs:
return set()
spos = SENSORS[sid]["pos"]
Z = np.array([[r["range_km"], np.radians(r["az_deg"])] for r in rs])
R = np.diag([rs[0]["sr"] ** 2, np.radians(rs[0]["saz"]) ** 2])
cost = np.full((len(tracks), len(rs)), 1e6)
for i, tr in enumerate(tracks):
zp, H = h_polar(tr["x"], spos)
Y = Z - zp
Y[:, 1] = (Y[:, 1] + np.pi) % (2 * np.pi) - np.pi
Si = np.linalg.inv(H @ tr["P"] @ H.T + R)
d2 = np.einsum("ij,jk,ik->i", Y, Si, Y)
cost[i] = np.where(d2 < self.gate, d2, 1e6)
used = set()
for i, j in zip(*linear_sum_assignment(cost)):
if cost[i, j] < self.gate:
tr, r = tracks[i], rs[j]
y, H, S, d2 = innovation(tr["x"], tr["P"], r)
K = tr["P"] @ H.T @ np.linalg.inv(S)
tr["x"] = tr["x"] + K @ y
tr["P"] = (np.eye(4) - K @ H) @ tr["P"]
tr["nis"].append(d2); tr["truths"].append(r["truth"])
tr["sensors"][sid] = tr["sensors"].get(sid, 0) + 1
tr["upd"] = True; used.add(j)
return used
def _merge(self):
"""같은 비행체를 따라가는 중복 항적(위치·속도가 통계적으로 같음)을 오래된 항적으로 통합한다."""
drop = set()
for a in range(len(self.tracks)):
for b in range(a + 1, len(self.tracks)):
A, B = self.tracks[a], self.tracks[b]
if A["id"] in drop or B["id"] in drop:
continue
d = A["x"] - B["x"]
d2 = float(d[:2] @ np.linalg.solve(A["P"][:2, :2] + B["P"][:2, :2], d[:2]))
if d2 < self.merge_gate and np.hypot(d[2], d[3]) < 0.05:
keep, lose = (A, B) if len(A["scans"]) >= len(B["scans"]) else (B, A)
for sid, n in lose["sensors"].items():
keep["sensors"][sid] = keep["sensors"].get(sid, 0) + n
keep["truths"] += lose["truths"]; keep["scans"] = sorted(set(keep["scans"]) | set(lose["scans"]))
keep["merged"] = keep.get("merged", 0) + 1
drop.add(lose["id"])
self.tracks = [t for t in self.tracks if t["id"] not in drop]
self.n_merged += len(drop)
def step(self, k, reps):
for tr in self.tracks: # 1) 예측
tr["x"], tr["P"] = ekf_predict(tr["x"], tr["P"], k - tr["k"], self.q)
tr["k"], tr["upd"] = k, False
for sid in sorted({r["sensor"] for r in reps}, key=lambda s: SENSORS[s]["sr"]): # 2) 정밀 센서부터 순차 갱신
rs = [r for r in reps if r["sensor"] == sid]
used = self._assign([t for t in self.tracks if t["confirmed"]], rs, sid) # 확정 항적 우선
rest = [r for j, r in enumerate(rs) if j not in used]
used2 = self._assign([t for t in self.tracks if not t["confirmed"]], rest, sid)
for j, r in enumerate(rest):
if j not in used2:
self._new(r, k)
if self.merge_gate > 0:
self._merge()
keep = []
for tr in self.tracks: # 3) 확정·삭제 관리
if tr["upd"]:
tr["miss"] = 0
if tr["scans"][-1] != k:
tr["scans"].append(k)
tr["hist"].append((k, tr["x"][0], tr["x"][1]))
else:
tr["miss"] += 1
m, n = self.confirm
if not tr["confirmed"] and sum(1 for s in tr["scans"] if s > k - n) >= m:
tr["confirmed"] = True
limit = self.max_miss if tr["confirmed"] else 1
(keep if tr["miss"] <= limit else self.dead).append(tr)
self.tracks = keep
def all_tracks(self):
return self.dead + self.tracks
def run_tracker(reps, **kw):
trk = Tracker(**kw)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
for k in range(min(by), max(by) + 1):
trk.step(k, by.get(k, []))
return trk.all_tracks()
def label_track(tr):
vals, cnt = np.unique(tr["truths"], return_counts=True)
return str(vals[np.argmax(cnt)]), float(cnt.max() / cnt.sum())
def evaluate(tracks, scn, min_reports=5):
"""정답 연결 평가: 재현율, 항적 분할(한 표적에 항적 2개 이상), 거짓 항적, 병합 오류, 위치 오차."""
conf = [t for t in tracks if t["confirmed"]]
for t in conf:
t["label"], t["purity"] = label_track(t)
cnt = pd.Series([r["truth"] for r in scn["reports"] if r["truth"].startswith("TGT")]).value_counts()
truths = [k for k, v in cnt.items() if v >= min_reports]
per = {g: [t for t in conf if t["label"] == g] for g in truths}
found = [g for g in truths if per[g]]
split = sum(max(0, len(v) - 1) for v in per.values())
false = [t for t in conf if not t["label"].startswith("TGT")]
merge = [t for t in conf if t["label"].startswith("TGT") and t["purity"] < 0.7]
err = []
for g in found:
for t in per[g]:
for k, x, y in t["hist"][3:]:
if 0 <= k < scn["T"]:
err.append(np.hypot(x - scn["targets"][g][k, 0], y - scn["targets"][g][k, 1]))
n = len(truths)
acc = len(found) / max(n + split + len(false) + len(merge), 1)
return {"정답표적": n, "탐지표적": len(found), "재현율": round(len(found) / max(n, 1), 3), "분할": split,
"거짓항적": len(false), "병합오류": len(merge), "위치RMSE(km)": round(float(np.sqrt(np.mean(np.square(err)))) if err else np.nan, 3),
"항적정확도": round(acc, 3)}
# 센서 신뢰도 가중치 매트릭스(기상별) [교육용가상자료]: 기본 신뢰도 x 기상 탐지 계수
BASE_W = {"radar": 1.0, "eoir": 0.8, "acoustic": 0.6}
WMAT = {w: {sid: round(BASE_W[s["type"]] * WEATHER[w][s["type"]] / WEATHER[w]["fa"], 2) for sid, s in SENSORS.items()} for w in WEATHER}
def track_features(trk, weather):
rows = []
for t in trk:
if not t["confirmed"]:
continue
n = len(t["scans"])
rows.append({"id": t["id"], "life": min(n, 10) / 10, "rps": len(t["truths"]) / n, "nsens": len(t["sensors"]),
"wsup": sum(WMAT[weather][s] * min(c / n, 1) for s, c in t["sensors"].items()),
"nis": float(np.median(t["nis"])) if t["nis"] else 9.0, "eoir_only": float(set(t["sensors"]) == {"EOIR-C"}),
"label": t["label"], "true": t["label"].startswith("TGT")})
return pd.DataFrame(rows)
FEAT = ["life", "rps", "nsens", "wsup", "nis", "eoir_only"]
def dataset(seeds, fa=2.0, **kw):
out = []
for i, sd in enumerate(seeds):
w = ["clear", "night", "fog", "rain"][i % 4]
scn = make_scenario(sd, weather=w, fa_mult=fa, ghost=(i % 3 == 0), special=[None, "turn", "cross", "formation"][i % 4], **kw)
reps = prepare(scn)
for ts in (8, 15, 30): # 실시간 판단 시점: 해당 시각까지의 보고만 사용
trk = run_tracker([r for r in reps if r["t"] <= ts]); evaluate(trk, scn)
df = track_features(trk, w); df["weather"] = w; df["scn"] = sd; df["ts"] = ts; out.append(df)
return pd.concat(out, ignore_index=True)
print("[센서 신뢰도 가중치 매트릭스]")
print(pd.DataFrame(WMAT).T.to_string())
tr = dataset(range(500, 508))
te = dataset(range(600, 606), fa=4.0, noise_mult={"RDR-B": 3.0}) # 평가: 오경보 4배 + 레이더B 재밍(학습보다 가혹)
eng = LogisticRegression(max_iter=1000, C=1.0).fit(tr[FEAT], tr["true"])
def grade(F):
s = eng.predict_proba(F[FEAT])[:, 1]
g = np.where(s >= 0.8, "High", np.where(s >= 0.5, "Med", "Low"))
return s, g
te["score"], te["grade"] = grade(te)
false_n = (~te["true"]).sum()
far_old = false_n / len(te) # 기존: 확정 항적 전부 경보
alarm = te["grade"] == "High"
far_new = (~te["true"] & alarm).sum() / max(alarm.sum(), 1)
red = 1 - far_new / far_old
print(f"\n평가 항적(시점별) {len(te)}개(실제 {te['true'].sum()} / 거짓 {false_n}), 시나리오 6개 x 판단 시점 3개(8·15·30초), 기상 4종, 오경보 4배, 레이더B 재밍, 기만 항적 포함")
print(f"오경보율(경보 중 거짓 항적 비율): 기존 {far_old:.1%} -> High 경보 {far_new:.1%}, 감축 {red:.0%} (목표 40% 이상: {'달성' if red >= 0.4 else '미달'})")
print(f"실제 항적 High 비율 {(te['true'] & alarm).sum() / te['true'].sum():.1%}, 실제 항적 Low 판정 {(te['true'] & (te['grade'] == 'Low')).sum()}건, ROC-AUC {roc_auc_score(te['true'], te['score']):.3f}")
print("\n등급별 항적 수:")
print(pd.crosstab(te["grade"], te["true"].map({True: "실제", False: "거짓"})).reindex(["High", "Med", "Low"]).fillna(0).astype(int).to_string())
coef = pd.Series(eng.coef_[0], index=FEAT).round(2)
print("\n특징 가중치(로지스틱):", {k: float(v) for k, v in coef.items()})
# 일관성: 특징에 ±5% 측정 흔들림을 10회 주어도 등급이 유지되는 비율
rng = np.random.default_rng(7)
same = []
for _ in range(10):
P = te.copy()
for c in ["life", "rps", "wsup", "nis"]:
P[c] = P[c] * rng.uniform(0.95, 1.05, len(P))
same.append((grade(P)[1] == te["grade"]).mean())
cons = float(np.mean(same))
print(f"\n신뢰도 등급 일관성(특징 ±5% 흔들림 10회): {cons:.1%} (목표 90% 이상: {'달성' if cons >= 0.9 else '미달'})")
print("\n[임계값 민감도] High 기준 경보 수 경보 중 거짓")
for th in [0.6, 0.7, 0.8, 0.9]:
a = te["score"] >= th
print(f" {th:.1f} {a.sum():4d} {(~te['true'] & a).sum() / max(a.sum(), 1):6.1%}")
print("등급 운용: High=통합 항적 경보 / Med=조작관 확인(미션 7) / Low=감시 유지(삭제하지 않음). 교전 여부는 교전수칙·지휘계통에 따른다.")
plt.figure(figsize=(6, 3.5))
plt.hist([te["score"][te["true"]], te["score"][~te["true"]]], bins=20, label=["true track", "false track"], color=["tab:blue", "tab:red"])
plt.axvline(0.8, c="k", ls="--"); plt.axvline(0.5, c="gray", ls=":"); plt.xlabel("fusion confidence"); plt.legend(); plt.title("Fusion confidence (synthetic)")
plt.tight_layout(); plt.show()
In [ ]:
"""스캔(1초)마다 센서별 원시 경보와 EKF 융합 항적을 나란히 갱신하는 야간 모드 전술 상황판(COP)을 구현하고, 원시 경보 대비 표시 항목 감소(인지 부하), 방어 지역 접근 긴급 강조, 화면 갱신 시간과 항적 선택 응답 시간을 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 6] 다중 센서 융합 항적 및 센서별 원시 경보 비교 전술 상황판(COP) 대시보드 구현 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 6] 다중 센서 융합 항적 및 센서별 원시 경보 비교 전술 상황판(COP) 대시보드 구현 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
# ---- 확장 칼만 필터(EKF) 다중 센서 항적 추적기 공통 함수 ----
# 상태 x = [동(km), 북(km), 동속(km/s), 북속(km/s)], 측정 z = 센서 기준 [거리(km), 방위(rad)]
def ekf_predict(x, P, dt, q):
F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], float)
G = np.array([[dt ** 2 / 2, 0], [0, dt ** 2 / 2], [dt, 0], [0, dt]])
return F @ x, F @ P @ F.T + q * G @ G.T
def h_polar(x, spos):
dx, dy = x[0] - spos[0], x[1] - spos[1]
r2 = max(dx * dx + dy * dy, 1e-9); r = np.sqrt(r2)
return np.array([r, np.arctan2(dx, dy)]), np.array([[dx / r, dy / r, 0, 0], [dy / r2, -dx / r2, 0, 0]])
def innovation(x, P, rep):
s = SENSORS[rep["sensor"]]
zp, H = h_polar(x, s["pos"])
y = np.array([rep["range_km"], np.radians(rep["az_deg"])]) - zp
y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi # 방위 차이는 -pi~pi로 감싼다
S = H @ P @ H.T + np.diag([rep["sr"] ** 2, np.radians(rep["saz"]) ** 2])
return y, H, S, float(y @ np.linalg.solve(S, y))
class Tracker:
"""센서 경보를 스캔(1초) 단위로 받아 EKF 항적을 갱신한다. 같은 스캔 안에서 센서별로
순차 갱신하므로, 여러 센서가 같은 비행체를 보면 하나의 항적으로 통합된다."""
def __init__(self, gate=9.21, q=4e-4, confirm=(3, 5), max_miss=4, merge_gate=6.0):
self.gate, self.q, self.confirm, self.max_miss, self.merge_gate = gate, q, confirm, max_miss, merge_gate
self.tracks, self.dead, self.nid, self.n_merged = [], [], 0, 0
def _new(self, rep, k):
C = meas_cov(rep)
P = np.zeros((4, 4)); P[:2, :2] = C; P[2, 2] = P[3, 3] = 0.25 ** 2
self.tracks.append({"id": f"FT{self.nid:03d}", "x": np.array([rep["x"], rep["y"], 0, 0]), "P": P, "k": k,
"scans": [k], "sensors": {rep["sensor"]: 1}, "truths": [rep["truth"]], "nis": [],
"hist": [(k, rep["x"], rep["y"])], "confirmed": False, "miss": 0, "born": k, "upd": True})
self.nid += 1
def _assign(self, tracks, rs, sid):
"""게이트(마할라노비스 거리) 안에서 전역 최근접 할당(GNN). 할당된 보고서 번호 집합을 돌려준다."""
if not tracks or not rs:
return set()
spos = SENSORS[sid]["pos"]
Z = np.array([[r["range_km"], np.radians(r["az_deg"])] for r in rs])
R = np.diag([rs[0]["sr"] ** 2, np.radians(rs[0]["saz"]) ** 2])
cost = np.full((len(tracks), len(rs)), 1e6)
for i, tr in enumerate(tracks):
zp, H = h_polar(tr["x"], spos)
Y = Z - zp
Y[:, 1] = (Y[:, 1] + np.pi) % (2 * np.pi) - np.pi
Si = np.linalg.inv(H @ tr["P"] @ H.T + R)
d2 = np.einsum("ij,jk,ik->i", Y, Si, Y)
cost[i] = np.where(d2 < self.gate, d2, 1e6)
used = set()
for i, j in zip(*linear_sum_assignment(cost)):
if cost[i, j] < self.gate:
tr, r = tracks[i], rs[j]
y, H, S, d2 = innovation(tr["x"], tr["P"], r)
K = tr["P"] @ H.T @ np.linalg.inv(S)
tr["x"] = tr["x"] + K @ y
tr["P"] = (np.eye(4) - K @ H) @ tr["P"]
tr["nis"].append(d2); tr["truths"].append(r["truth"])
tr["sensors"][sid] = tr["sensors"].get(sid, 0) + 1
tr["upd"] = True; used.add(j)
return used
def _merge(self):
"""같은 비행체를 따라가는 중복 항적(위치·속도가 통계적으로 같음)을 오래된 항적으로 통합한다."""
drop = set()
for a in range(len(self.tracks)):
for b in range(a + 1, len(self.tracks)):
A, B = self.tracks[a], self.tracks[b]
if A["id"] in drop or B["id"] in drop:
continue
d = A["x"] - B["x"]
d2 = float(d[:2] @ np.linalg.solve(A["P"][:2, :2] + B["P"][:2, :2], d[:2]))
if d2 < self.merge_gate and np.hypot(d[2], d[3]) < 0.05:
keep, lose = (A, B) if len(A["scans"]) >= len(B["scans"]) else (B, A)
for sid, n in lose["sensors"].items():
keep["sensors"][sid] = keep["sensors"].get(sid, 0) + n
keep["truths"] += lose["truths"]; keep["scans"] = sorted(set(keep["scans"]) | set(lose["scans"]))
keep["merged"] = keep.get("merged", 0) + 1
drop.add(lose["id"])
self.tracks = [t for t in self.tracks if t["id"] not in drop]
self.n_merged += len(drop)
def step(self, k, reps):
for tr in self.tracks: # 1) 예측
tr["x"], tr["P"] = ekf_predict(tr["x"], tr["P"], k - tr["k"], self.q)
tr["k"], tr["upd"] = k, False
for sid in sorted({r["sensor"] for r in reps}, key=lambda s: SENSORS[s]["sr"]): # 2) 정밀 센서부터 순차 갱신
rs = [r for r in reps if r["sensor"] == sid]
used = self._assign([t for t in self.tracks if t["confirmed"]], rs, sid) # 확정 항적 우선
rest = [r for j, r in enumerate(rs) if j not in used]
used2 = self._assign([t for t in self.tracks if not t["confirmed"]], rest, sid)
for j, r in enumerate(rest):
if j not in used2:
self._new(r, k)
if self.merge_gate > 0:
self._merge()
keep = []
for tr in self.tracks: # 3) 확정·삭제 관리
if tr["upd"]:
tr["miss"] = 0
if tr["scans"][-1] != k:
tr["scans"].append(k)
tr["hist"].append((k, tr["x"][0], tr["x"][1]))
else:
tr["miss"] += 1
m, n = self.confirm
if not tr["confirmed"] and sum(1 for s in tr["scans"] if s > k - n) >= m:
tr["confirmed"] = True
limit = self.max_miss if tr["confirmed"] else 1
(keep if tr["miss"] <= limit else self.dead).append(tr)
self.tracks = keep
def all_tracks(self):
return self.dead + self.tracks
def run_tracker(reps, **kw):
trk = Tracker(**kw)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
for k in range(min(by), max(by) + 1):
trk.step(k, by.get(k, []))
return trk.all_tracks()
def label_track(tr):
vals, cnt = np.unique(tr["truths"], return_counts=True)
return str(vals[np.argmax(cnt)]), float(cnt.max() / cnt.sum())
def evaluate(tracks, scn, min_reports=5):
"""정답 연결 평가: 재현율, 항적 분할(한 표적에 항적 2개 이상), 거짓 항적, 병합 오류, 위치 오차."""
conf = [t for t in tracks if t["confirmed"]]
for t in conf:
t["label"], t["purity"] = label_track(t)
cnt = pd.Series([r["truth"] for r in scn["reports"] if r["truth"].startswith("TGT")]).value_counts()
truths = [k for k, v in cnt.items() if v >= min_reports]
per = {g: [t for t in conf if t["label"] == g] for g in truths}
found = [g for g in truths if per[g]]
split = sum(max(0, len(v) - 1) for v in per.values())
false = [t for t in conf if not t["label"].startswith("TGT")]
merge = [t for t in conf if t["label"].startswith("TGT") and t["purity"] < 0.7]
err = []
for g in found:
for t in per[g]:
for k, x, y in t["hist"][3:]:
if 0 <= k < scn["T"]:
err.append(np.hypot(x - scn["targets"][g][k, 0], y - scn["targets"][g][k, 1]))
n = len(truths)
acc = len(found) / max(n + split + len(false) + len(merge), 1)
return {"정답표적": n, "탐지표적": len(found), "재현율": round(len(found) / max(n, 1), 3), "분할": split,
"거짓항적": len(false), "병합오류": len(merge), "위치RMSE(km)": round(float(np.sqrt(np.mean(np.square(err)))) if err else np.nan, 3),
"항적정확도": round(acc, 3)}
ALERT_R = 4.0 # 가상 방어 지역 중심 반경(km) 이내 접근 시 긴급 강조
COLORS = {"RDR-A": "deepskyblue", "RDR-B": "cyan", "EOIR-C": "yellow", "ACU-D": "violet"}
def grade(t):
n = len(t["scans"])
if len(t["sensors"]) >= 2 and len(t["truths"]) / n >= 1.5 and n >= 5:
return "High"
return "Med"
scn = make_scenario(seed=611, n=6, T=40, weather="night", ghost=True)
reps = prepare(scn)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
trk = Tracker()
plt.style.use("dark_background") # 야간 작전 모드
fig, ax = plt.subplots(2, 2, figsize=(11, 8))
upd, n_raw, n_fused, emerg_log = [], [], [], []
for k in range(scn["T"]):
t0 = time.perf_counter()
rk = by.get(k, [])
trk.step(k, rk)
conf = [t for t in trk.tracks if t["confirmed"]]
for t in conf:
t["grade"] = grade(t)
t["emerg"] = np.hypot(*(t["x"][:2] - CENTER)) < ALERT_R and t["grade"] == "High"
n_raw.append(len(rk)); n_fused.append(len(conf))
for a in ax.flat:
a.clear()
for sid in SENSORS: # 패널 1: 센서별 원시 경보
pts = np.array([[r["x"], r["y"]] for r in rk if r["sensor"] == sid]).reshape(-1, 2)
ax[0, 0].scatter(pts[:, 0], pts[:, 1], s=10, c=COLORS[sid], label=sid)
ax[0, 0].set_title(f"Raw sensor alerts t={k}s ({len(rk)})"); ax[0, 0].legend(fontsize=7, loc="upper right")
for t in conf: # 패널 2: 융합 항적
h = np.array(t["hist"])
c = "red" if t["emerg"] else ("lime" if t["grade"] == "High" else "orange")
ax[0, 1].plot(h[-10:, 1], h[-10:, 2], c=c, lw=2 if t["emerg"] else 1)
ax[0, 1].text(t["x"][0], t["x"][1], t["id"], fontsize=7, color=c)
ax[0, 1].add_patch(plt.Circle(CENTER, ALERT_R, fill=False, ls="--", color="red"))
ax[0, 1].set_title(f"Fused tracks ({len(conf)})")
for a in ax[0]:
a.set_xlim(-10, 26); a.set_ylim(-14, 18); a.set_aspect("equal")
ax[1, 0].plot(n_raw, label="raw alerts"); ax[1, 0].plot(n_fused, label="fused tracks"); ax[1, 0].legend(); ax[1, 0].set_title("Items on screen per scan")
ax[1, 1].axis("off"); ax[1, 1].set_title("Track list (operator review)")
for i, t in enumerate(sorted(conf, key=lambda t: (not t["emerg"], t["grade"] != "High"))[:8]):
ax[1, 1].text(0, 0.9 - i * 0.11, f"{t['id']} {t['grade']:4s} sensors {len(t['sensors'])} pos({t['x'][0]:5.1f},{t['x'][1]:5.1f})"
+ (" !! APPROACH" if t["emerg"] else ""), color="red" if t["emerg"] else "white", fontsize=8, family="monospace")
if any(t["emerg"] for t in conf):
emerg_log.append(k)
for sp in ax[0, 1].spines.values():
sp.set_color("red"); sp.set_linewidth(3)
fig.canvas.draw()
upd.append((time.perf_counter() - t0) * 1000)
# 항적 선택 응답성: 목록에서 항적을 고르면 기여 센서·최근 궤적·속도를 보여 주는 콜백
def on_select(tid):
t = next(t for t in trk.all_tracks() if t["id"] == tid)
sp = np.hypot(t["x"][2], t["x"][3]) * 1000
return {"id": tid, "grade": t.get("grade"), "sensors": dict(t["sensors"]), "speed_mps": round(sp), "last_pos": [round(float(v), 2) for v in t["x"][:2]]}
resp = []
for t in [t for t in trk.tracks if t["confirmed"]]:
t0 = time.perf_counter(); info = on_select(t["id"]); resp.append((time.perf_counter() - t0) * 1000)
ev = evaluate(trk.all_tracks(), scn)
red = 1 - np.mean(n_fused[5:]) / np.mean(n_raw[5:])
print("[전술 상황판(COP) - 교육용가상자료, 가상 방어 지역·야간]")
print(f"패널 4종(원시 경보, 융합 항적, 표시 항목 추이, 항적 목록) {len(upd)}회 갱신")
print(f"화면 갱신 시간 최대 {max(upd):.0f} ms, 평균 {np.mean(upd):.0f} ms (목표 1초 이내: {'달성' if max(upd) < 1000 else '미달'})")
print(f"항적 선택 응답 최대 {max(resp):.2f} ms (목표 100ms 이내: {'달성' if max(resp) < 100 else '미달'})")
print(f"스캔당 표시 항목: 원시 경보 평균 {np.mean(n_raw[5:]):.1f}개 -> 융합 항적 {np.mean(n_fused[5:]):.1f}개 (인지 부하 {red:.0%} 감소)")
fk = [t for t in trk.all_tracks() if t["confirmed"] and not t["label"].startswith("TGT")]
print(f"융합 결과: 정답 표적 {ev['정답표적']}개 중 {ev['탐지표적']}개 추적, 거짓 항적 {ev['거짓항적']}개 중 Med 표시 {sum(t.get('grade') == 'Med' for t in fk)}개")
print(f"방어 지역 접근 긴급 강조 스캔: {len(emerg_log)}회 (첫 발생 t={emerg_log[0] if emerg_log else '-'}s)")
print("선택 예시:", info)
print("상황판은 판단 보조이며, 항적 분리·통합 확정은 조작관 검토(미션 7), 교전 여부는 교전수칙·지휘계통에 따른다.")
plt.tight_layout(); plt.show()
In [ ]:
# [TODO : 학생 실습 7] 중복 항적 강제 분리 및 통합 확정을 위한 조작관 HITL 검토/승인 체계 구축
# 미션 지침 및 검증 기준을 확인하고 코드를 실행 및 튜닝하세요.
"""센서 간 상충 정보로 생긴 항적 충돌(중복 의심, 병합 의심, 단일 센서 항적)을 조작관 검토 대기열에 올리고, 강제 통합·분리·추적 보류·현행 유지를 조작관이 근거와 함께 직접 결정하도록 하는 HITL 검토/승인 체계와 해시 체인 감사 로그를 구현한다."""
import hashlib
import json
import numpy as np
import pandas as pd
ACTIONS = {"MERGE": "강제 통합", "SPLIT": "항적 분리", "HOLD": "추적 보류(감시 유지)", "KEEP": "현행 유지"}
SEVERITY = {"DUP": 0, "MIXED": 1, "SINGLE": 2} # 중복 의심 > 병합 의심 > 단일 센서
HOLD_LIMIT_S = 60
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심(km)
class WorkflowError(ValueError):
pass
class ReviewQueue:
def __init__(self, tracks, clock):
self.tracks, self.items, self.log, self.clock, self.feedback = tracks, {}, [], clock, []
def _record(self, body):
prev = self.log[-1]["hash"] if self.log else "GENESIS"
body = {**body, "ts": round(self.clock(), 1)}
h = hashlib.sha256((prev + json.dumps(body, ensure_ascii=False, sort_keys=True)).encode()).hexdigest()
self.log.append({**body, "prev": prev, "hash": h})
def verify_log(self):
prev = "GENESIS"
for n, rec in enumerate(self.log, start=1):
body = {k: v for k, v in rec.items() if k not in ("prev", "hash")}
if rec["prev"] != prev or rec["hash"] != hashlib.sha256((prev + json.dumps(body, ensure_ascii=False, sort_keys=True)).encode()).hexdigest():
return False, n
prev = rec["hash"]
return True, -1
def submit(self, cid, kind, tids, note):
"""시스템은 충돌을 '제안'만 하고 항적 상태는 바꾸지 않는다(PENDING)."""
d = min(np.hypot(*(np.array(self.tracks[t]["pos"]) - CENTER)) for t in tids)
self.items[cid] = {"id": cid, "kind": kind, "tids": tids, "note": note, "dist_km": round(float(d), 1), "status": "PENDING",
"submitted": self.clock(), "hold_since": None, "escalated": False}
self._record({"event": "SUBMIT", "id": cid, "kind": kind, "tids": tids, "by": "FUSION"})
def pending(self):
todo = [i for i in self.items.values() if i["status"] in ("PENDING", "HOLD")]
return sorted(todo, key=lambda i: (not i["escalated"], SEVERITY[i["kind"]], i["dist_km"], i["submitted"]))
def decide(self, cid, action, operator, reason):
if action not in ACTIONS:
raise WorkflowError(f"허용되지 않은 조치: {action}")
if not operator:
raise WorkflowError("조작관 ID 필요")
it = self.items.get(cid)
if it is None or it["status"] in ("DONE",):
raise WorkflowError("대기열에 없거나 이미 종결된 충돌")
if len(reason.strip()) < 5:
raise WorkflowError("판단 근거 5자 이상 필요")
T = [self.tracks[t] for t in it["tids"]]
if action == "MERGE":
if len(T) < 2:
raise WorkflowError("통합은 항적 2개 이상 필요")
gap = np.hypot(*(np.array(T[0]["pos"]) - np.array(T[1]["pos"])))
if gap > 1.0 and len(reason.strip()) < 15:
raise WorkflowError(f"간격 {gap:.1f}km 항적 통합은 근거 15자 이상 필요(오통합 방지)")
T[1]["status"] = f"MERGED->{T[0]['id']}"; T[0]["sensors"] = sorted(set(T[0]["sensors"]) | set(T[1]["sensors"]))
elif action == "SPLIT":
if len(T[0]["sensors"]) < 2:
raise WorkflowError("단일 센서 항적은 분리 대상이 아님")
new = f"{T[0]['id']}b"
self.tracks[new] = {"id": new, "pos": T[0]["pos"], "sensors": T[0]["sensors"][1:], "status": "SPLIT"}
T[0]["sensors"] = T[0]["sensors"][:1]
elif action == "HOLD":
for t in T:
t["status"] = "HOLD"
it["status"] = "HOLD" if action == "HOLD" else "DONE"
it["hold_since"] = self.clock() if action == "HOLD" else None
it["decided"], it["action"] = self.clock(), action
self._record({"event": "DECIDE", "id": cid, "action": action, "operator": operator, "reason": reason})
def check_holds(self):
out = []
for i in self.items.values():
if i["status"] == "HOLD" and not i["escalated"] and self.clock() - i["hold_since"] > HOLD_LIMIT_S:
i["escalated"] = True; out.append(i["id"])
self._record({"event": "ESCALATE", "id": i["id"], "to": "방공 상황장교", "by": "SYSTEM"})
return out
def commander_feedback(self, cid, text):
self.feedback.append({"id": cid, "text": text}); self._record({"event": "FEEDBACK", "id": cid, "text": text})
class FakeClock:
t = 0.0
def __call__(self):
return self.t
# [교육용가상자료] 융합 엔진이 넘긴 가상 항적과 충돌 6건
tracks = {t: {"id": t, "pos": p, "sensors": s, "status": "ACTIVE"} for t, p, s in [
("FT003", (9.1, 3.0), ["RDR-A", "RDR-B"]), ("FT007", (9.4, 3.2), ["EOIR-C", "ACU-D"]), ("FT010", (14.0, -3.0), ["RDR-A", "RDR-B", "EOIR-C"]),
("FT012", (3.5, 12.0), ["EOIR-C"]), ("FT015", (12.0, 6.0), ["RDR-A"]), ("FT016", (13.5, 6.5), ["RDR-B"]), ("FT020", (6.0, 0.5), ["RDR-A", "ACU-D"])]}
clock = FakeClock()
q = ReviewQueue(tracks, clock)
for cid, kind, tids, note in [("C01", "DUP", ["FT003", "FT007"], "0.4km 간격, 속도 유사"), ("C02", "MIXED", ["FT010"], "EO/IR 방위가 레이더와 2도 어긋남"),
("C03", "SINGLE", ["FT012"], "EO/IR 단독 30초 유지"), ("C04", "DUP", ["FT015", "FT016"], "1.6km 간격, 레이더별 단독"),
("C05", "SINGLE", ["FT020"], "음향 소실 후 레이더 단독"), ("C06", "MIXED", ["FT020"], "음향·레이더 위치 불일치")]:
clock.t += 3; q.submit(cid, kind, tids, note)
print("검토 순서(충돌 유형·방어 지역 거리 순):", [i["id"] for i in q.pending()])
print("자동 적용된 충돌 수:", sum(i["status"] != "PENDING" for i in q.items.values()))
blocked = 0
for cid, a, op, r in [("C01", "MERGE", "", "근거 충분"), ("C02", "AUTO", "op1", "자동 처리"), ("C03", "SPLIT", "op1", "단독 항적 분리 시도"),
("C04", "MERGE", "op1", "위치 유사"), ("C05", "KEEP", "op1", "확인")]:
try:
q.decide(cid, a, op, r)
except WorkflowError as e:
blocked += 1; print(f"차단({cid} {a}): {e}")
clock.t += 10; q.decide("C01", "MERGE", "op1", "속도·고도 추정 일치, 같은 비행체로 판단")
clock.t += 6; q.decide("C02", "SPLIT", "op2", "EO/IR 영상상 2기 근접 비행 확인")
clock.t += 5; q.decide("C04", "MERGE", "op1", "레이더 A·B 시계 보정 후 궤적 일치, 동일 표적으로 판단")
clock.t += 4; q.decide("C03", "HOLD", "op2", "기만 가능성, 레이더 재확인 요청")
clock.t += 3; q.decide("C05", "KEEP", "op2", "레이더 연속 추적 정상")
clock.t += 90
print("보류 기한 초과 상신:", q.check_holds(), "/ 상신 후 검토 순서:", [i["id"] for i in q.pending()])
q.commander_feedback("C03", "EO/IR 단독 항적은 30초 이상 지속 시 상황장교에게 즉시 보고")
df = pd.DataFrame(q.items.values())
df["대기(s)"] = (df["decided"] - df["submitted"]).round(1)
print("\n" + df[["id", "kind", "tids", "dist_km", "status", "action", "escalated", "대기(s)"]].fillna("-").to_string(index=False))
print("\n항적 상태:", {k: v["status"] for k, v in tracks.items()})
n_dec = sum(r["event"] == "DECIDE" for r in q.log)
n_auto = sum(r["event"] == "DECIDE" and not r.get("operator") for r in q.log)
n_exp = sum("decided" in i for i in q.items.values())
print(f"\n자동 적용 {n_auto}건 / 규칙 위반 차단 {blocked}건 / 결정 로그 {n_dec}건(결정 {n_exp}건, 로깅 성공률 {n_dec / n_exp:.0%})")
print(f"평균 검토 대기 {df['대기(s)'].mean():.1f}s, 감사 로그 {len(q.log)}건, 해시 체인 검증: {q.verify_log()}")
q.log[6]["action"] = "KEEP"
print("7번째 기록 변조 후 검증(False, 깨진 위치):", q.verify_log())
print("\n책임 분계: 융합 엔진=충돌 제안, 조작관=통합·분리·보류 결정, 상황장교=상신 건 판단. 항적 결정은 교전 결정이 아니며 교전은 교전수칙·지휘계통에 따른다.")
In [ ]:
"""레이더 통신 단절과 전자전 재밍 세기를 단계적으로 높여 항적 융합 성능의 한계를 찾고, 비행체·오경보 트래픽을 평시의 5배로 늘리며 손상 메시지를 섞은 부하 시험에서 시스템 다운·데이터 유실 없이 기능 축소 모드로 버티는지와 병목 함수를 확인한다."""
import cProfile
import pstats
import time
import tracemalloc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 8] 다수 센서 통신 단절 및 전자전 재밍 환경 모의 항적 융합 한계치 스트레스 테스트 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 8] 다수 센서 통신 단절 및 전자전 재밍 환경 모의 항적 융합 한계치 스트레스 테스트 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
# ---- 확장 칼만 필터(EKF) 다중 센서 항적 추적기 공통 함수 ----
# 상태 x = [동(km), 북(km), 동속(km/s), 북속(km/s)], 측정 z = 센서 기준 [거리(km), 방위(rad)]
def ekf_predict(x, P, dt, q):
F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], float)
G = np.array([[dt ** 2 / 2, 0], [0, dt ** 2 / 2], [dt, 0], [0, dt]])
return F @ x, F @ P @ F.T + q * G @ G.T
def h_polar(x, spos):
dx, dy = x[0] - spos[0], x[1] - spos[1]
r2 = max(dx * dx + dy * dy, 1e-9); r = np.sqrt(r2)
return np.array([r, np.arctan2(dx, dy)]), np.array([[dx / r, dy / r, 0, 0], [dy / r2, -dx / r2, 0, 0]])
def innovation(x, P, rep):
s = SENSORS[rep["sensor"]]
zp, H = h_polar(x, s["pos"])
y = np.array([rep["range_km"], np.radians(rep["az_deg"])]) - zp
y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi # 방위 차이는 -pi~pi로 감싼다
S = H @ P @ H.T + np.diag([rep["sr"] ** 2, np.radians(rep["saz"]) ** 2])
return y, H, S, float(y @ np.linalg.solve(S, y))
class Tracker:
"""센서 경보를 스캔(1초) 단위로 받아 EKF 항적을 갱신한다. 같은 스캔 안에서 센서별로
순차 갱신하므로, 여러 센서가 같은 비행체를 보면 하나의 항적으로 통합된다."""
def __init__(self, gate=9.21, q=4e-4, confirm=(3, 5), max_miss=4, merge_gate=6.0):
self.gate, self.q, self.confirm, self.max_miss, self.merge_gate = gate, q, confirm, max_miss, merge_gate
self.tracks, self.dead, self.nid, self.n_merged = [], [], 0, 0
def _new(self, rep, k):
C = meas_cov(rep)
P = np.zeros((4, 4)); P[:2, :2] = C; P[2, 2] = P[3, 3] = 0.25 ** 2
self.tracks.append({"id": f"FT{self.nid:03d}", "x": np.array([rep["x"], rep["y"], 0, 0]), "P": P, "k": k,
"scans": [k], "sensors": {rep["sensor"]: 1}, "truths": [rep["truth"]], "nis": [],
"hist": [(k, rep["x"], rep["y"])], "confirmed": False, "miss": 0, "born": k, "upd": True})
self.nid += 1
def _assign(self, tracks, rs, sid):
"""게이트(마할라노비스 거리) 안에서 전역 최근접 할당(GNN). 할당된 보고서 번호 집합을 돌려준다."""
if not tracks or not rs:
return set()
spos = SENSORS[sid]["pos"]
Z = np.array([[r["range_km"], np.radians(r["az_deg"])] for r in rs])
R = np.diag([rs[0]["sr"] ** 2, np.radians(rs[0]["saz"]) ** 2])
cost = np.full((len(tracks), len(rs)), 1e6)
for i, tr in enumerate(tracks):
zp, H = h_polar(tr["x"], spos)
Y = Z - zp
Y[:, 1] = (Y[:, 1] + np.pi) % (2 * np.pi) - np.pi
Si = np.linalg.inv(H @ tr["P"] @ H.T + R)
d2 = np.einsum("ij,jk,ik->i", Y, Si, Y)
cost[i] = np.where(d2 < self.gate, d2, 1e6)
used = set()
for i, j in zip(*linear_sum_assignment(cost)):
if cost[i, j] < self.gate:
tr, r = tracks[i], rs[j]
y, H, S, d2 = innovation(tr["x"], tr["P"], r)
K = tr["P"] @ H.T @ np.linalg.inv(S)
tr["x"] = tr["x"] + K @ y
tr["P"] = (np.eye(4) - K @ H) @ tr["P"]
tr["nis"].append(d2); tr["truths"].append(r["truth"])
tr["sensors"][sid] = tr["sensors"].get(sid, 0) + 1
tr["upd"] = True; used.add(j)
return used
def _merge(self):
"""같은 비행체를 따라가는 중복 항적(위치·속도가 통계적으로 같음)을 오래된 항적으로 통합한다."""
drop = set()
for a in range(len(self.tracks)):
for b in range(a + 1, len(self.tracks)):
A, B = self.tracks[a], self.tracks[b]
if A["id"] in drop or B["id"] in drop:
continue
d = A["x"] - B["x"]
d2 = float(d[:2] @ np.linalg.solve(A["P"][:2, :2] + B["P"][:2, :2], d[:2]))
if d2 < self.merge_gate and np.hypot(d[2], d[3]) < 0.05:
keep, lose = (A, B) if len(A["scans"]) >= len(B["scans"]) else (B, A)
for sid, n in lose["sensors"].items():
keep["sensors"][sid] = keep["sensors"].get(sid, 0) + n
keep["truths"] += lose["truths"]; keep["scans"] = sorted(set(keep["scans"]) | set(lose["scans"]))
keep["merged"] = keep.get("merged", 0) + 1
drop.add(lose["id"])
self.tracks = [t for t in self.tracks if t["id"] not in drop]
self.n_merged += len(drop)
def step(self, k, reps):
for tr in self.tracks: # 1) 예측
tr["x"], tr["P"] = ekf_predict(tr["x"], tr["P"], k - tr["k"], self.q)
tr["k"], tr["upd"] = k, False
for sid in sorted({r["sensor"] for r in reps}, key=lambda s: SENSORS[s]["sr"]): # 2) 정밀 센서부터 순차 갱신
rs = [r for r in reps if r["sensor"] == sid]
used = self._assign([t for t in self.tracks if t["confirmed"]], rs, sid) # 확정 항적 우선
rest = [r for j, r in enumerate(rs) if j not in used]
used2 = self._assign([t for t in self.tracks if not t["confirmed"]], rest, sid)
for j, r in enumerate(rest):
if j not in used2:
self._new(r, k)
if self.merge_gate > 0:
self._merge()
keep = []
for tr in self.tracks: # 3) 확정·삭제 관리
if tr["upd"]:
tr["miss"] = 0
if tr["scans"][-1] != k:
tr["scans"].append(k)
tr["hist"].append((k, tr["x"][0], tr["x"][1]))
else:
tr["miss"] += 1
m, n = self.confirm
if not tr["confirmed"] and sum(1 for s in tr["scans"] if s > k - n) >= m:
tr["confirmed"] = True
limit = self.max_miss if tr["confirmed"] else 1
(keep if tr["miss"] <= limit else self.dead).append(tr)
self.tracks = keep
def all_tracks(self):
return self.dead + self.tracks
def run_tracker(reps, **kw):
trk = Tracker(**kw)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
for k in range(min(by), max(by) + 1):
trk.step(k, by.get(k, []))
return trk.all_tracks()
def label_track(tr):
vals, cnt = np.unique(tr["truths"], return_counts=True)
return str(vals[np.argmax(cnt)]), float(cnt.max() / cnt.sum())
def evaluate(tracks, scn, min_reports=5):
"""정답 연결 평가: 재현율, 항적 분할(한 표적에 항적 2개 이상), 거짓 항적, 병합 오류, 위치 오차."""
conf = [t for t in tracks if t["confirmed"]]
for t in conf:
t["label"], t["purity"] = label_track(t)
cnt = pd.Series([r["truth"] for r in scn["reports"] if r["truth"].startswith("TGT")]).value_counts()
truths = [k for k, v in cnt.items() if v >= min_reports]
per = {g: [t for t in conf if t["label"] == g] for g in truths}
found = [g for g in truths if per[g]]
split = sum(max(0, len(v) - 1) for v in per.values())
false = [t for t in conf if not t["label"].startswith("TGT")]
merge = [t for t in conf if t["label"].startswith("TGT") and t["purity"] < 0.7]
err = []
for g in found:
for t in per[g]:
for k, x, y in t["hist"][3:]:
if 0 <= k < scn["T"]:
err.append(np.hypot(x - scn["targets"][g][k, 0], y - scn["targets"][g][k, 1]))
n = len(truths)
acc = len(found) / max(n + split + len(false) + len(merge), 1)
return {"정답표적": n, "탐지표적": len(found), "재현율": round(len(found) / max(n, 1), 3), "분할": split,
"거짓항적": len(false), "병합오류": len(merge), "위치RMSE(km)": round(float(np.sqrt(np.mean(np.square(err)))) if err else np.nan, 3),
"항적정확도": round(acc, 3)}
# 1) 결함 허용 한계: 센서 단절·재밍 단계별 융합 성능
print("[결함 허용 한계 - 교육용가상자료, 시나리오별 표적 6개 x 60초]")
cases = [("정상", {}), ("RDR-B 단절 30초", dict(dropout=("RDR-B", 15, 45))), ("RDR-A 단절 60초", dict(dropout=("RDR-A", 0, 60))),
("RDR-A 재밍 x3", dict(noise_mult={"RDR-A": 3.0})), ("레이더 2대 재밍 x3", dict(noise_mult={"RDR-A": 3.0, "RDR-B": 3.0})),
("레이더 2대 재밍 x6", dict(noise_mult={"RDR-A": 6.0, "RDR-B": 6.0})), ("RDR-A 단절+RDR-B 재밍 x6", dict(dropout=("RDR-A", 0, 60), noise_mult={"RDR-B": 6.0}))]
rows = []
for i, (nm, kw) in enumerate(cases):
evs = [evaluate(run_tracker(prepare(s)), s) for s in [make_scenario(800 + 10 * i + j, **kw) for j in range(2)]]
rows.append({"조건": nm, "재현율": np.mean([e["재현율"] for e in evs]), "분할": sum(e["분할"] for e in evs),
"거짓항적": sum(e["거짓항적"] for e in evs), "위치RMSE(km)": round(np.mean([e["위치RMSE(km)"] for e in evs]), 3),
"항적정확도": round(np.mean([e["항적정확도"] for e in evs]), 3)})
L = pd.DataFrame(rows)
print(L.to_string(index=False))
lim = L[L["항적정확도"] < 0.8]["조건"].tolist()
print("항적 정확도 80% 미만(한계):", lim if lim else "없음", "-> 해당 조건은 상황판에 '융합 신뢰 저하' 표시, 통합·분리 결정을 조작관 검토로 전환")
# 2) 부하 시험: 평시(표적 6, 오경보 1배) vs 500%(표적 30, 오경보 5배), 손상 메시지 2% 주입
def safe_prepare(scn, rng):
raw = list(scn["reports"])
for k in np.flatnonzero(rng.random(len(raw)) < 0.02): # 손상 메시지: 결측·미등록 센서·중복
bad = dict(raw[k])
kind = k % 3
if kind == 0:
bad["range_km"] = None
elif kind == 1:
bad["sensor"] = "UNKNOWN"
bad["msg_id"] = raw[k]["msg_id"] if kind == 2 else bad["msg_id"] + "x"
raw.append(bad)
reg = [r for r in raw if r.get("sensor") in SENSORS]
q_unknown = len(raw) - len(reg)
reps, why = clean(reg)
off = estimate_offsets(scn["heartbeats"])
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"]), len(raw), q_unknown + sum(why.values())
def run_load(n, fa, seed, budget_ms=1000.0, degrade_ms=150.0):
rng = np.random.default_rng(seed)
scn = make_scenario(seed, n=n, T=40, fa_mult=fa)
reps, total, quar = safe_prepare(scn, rng)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
trk, done, crash, degraded, worst = Tracker(), 0, 0, 0, 0.0
tracemalloc.start(); m0 = tracemalloc.get_traced_memory()[0]
last = 0.0
for k in range(min(by), max(by) + 1):
trk.merge_gate = 0.0 if last > degrade_ms else 6.0 # 기능 축소: 직전 스캔이 느리면 중복 통합 단계를 1스캔 생략
degraded += trk.merge_gate == 0.0
t0 = time.perf_counter()
try:
trk.step(k, by.get(k, [])); done += len(by.get(k, []))
except Exception:
crash += 1
last = (time.perf_counter() - t0) * 1000; worst = max(worst, last)
mem = (tracemalloc.get_traced_memory()[0] - m0) / 1e6; tracemalloc.stop()
ev = evaluate(trk.all_tracks(), scn)
return {"입력": total, "처리": done, "격리": quar, "유실": total - done - quar, "다운": crash, "기능축소 스캔": degraded,
"최대 ms/스캔": round(worst, 1), "메모리증가MB": round(mem, 2), "재현율": ev["재현율"]}
print("\n[부하 시험]")
res = {"평시": run_load(6, 1.0, 901), "500%": run_load(30, 5.0, 902)}
print(pd.DataFrame(list(res.values()), index=list(res)).to_string())
hi = res["500%"]
print(f"트래픽 {hi['입력'] / res['평시']['입력']:.1f}배: 시스템 다운 {hi['다운']}건, 데이터 유실 {hi['유실'] / hi['입력']:.1%} "
f"(목표 0건/0%: {'달성' if hi['다운'] == 0 and hi['유실'] == 0 else '미달'}), 스캔 최대 {hi['최대 ms/스캔']} ms(1초 주기 이내: {'예' if hi['최대 ms/스캔'] < 1000 else '아니오'})")
# 3) 병목 프로파일링(500% 조건)
scn = make_scenario(902, n=30, T=20, fa_mult=5.0)
reps = prepare(scn)
pr = cProfile.Profile(); pr.enable(); run_tracker(reps); pr.disable()
st = pstats.Stats(pr).stats
tot = {f"{fn[2]}": v[2] for fn, v in st.items()} # 함수별 자체 실행 시간
top = sorted(((k, v) for k, v in tot.items() if not k.startswith("<")), key=lambda x: -x[1])[:4]
cum = {fn[2]: v[3] for fn, v in st.items() if fn[2] in ("_assign", "_merge", "step", "ekf_predict")}
allt = cum.get("step", 1e-9)
print("\n[병목 프로파일] 단계별 누적 시간 비율:", {k: f"{v / allt:.0%}" for k, v in cum.items() if k != "step"})
print("자체 시간 상위 함수:", [f"{k}({v * 1000:.0f}ms)" for k, v in top])
bn = max((k for k in cum if k != "step"), key=cum.get)
print(f"병목: {bn}", "-> 개선 방향: 중복 통합(_merge)은 격자 분할로 비교 쌍을 줄이고, 할당(_assign)은 항적별 역행렬 반복을 벡터화한다. 기능 축소 시에도 경보는 버리지 않는다(Fail-safe).")
plt.figure(figsize=(7, 3.5))
plt.plot(L["조건"], L["재현율"], "o-", label="recall"); plt.plot(L["조건"], L["항적정확도"], "s--", label="track accuracy")
plt.axhline(0.8, c="k", ls=":"); plt.xticks(range(len(L)), [f"C{i}" for i in range(len(L))]); plt.legend(); plt.title("Fusion limit under dropout / jamming (synthetic)")
plt.tight_layout(); plt.show()
In [ ]:
"""초기 배치 파라미터로 운용한 모의 AAR 일지에서 항적 분할(Track Split)·거짓 항적·병합 오류 사례를 찾아 원인을 분석하고, 프로세스 잡음·게이트·확정 규칙·삭제 기준·중복 통합 게이트를 재탐색해 융합 임계값을 재보정한 뒤 새 시나리오 회귀 시험으로 정확도 향상과 재발 오류를 확인한다."""
import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 9] AAR 일지 기반 항적 왜곡/분할(Track Split) 오류 분석 및 융합 임계값 최적화 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 9] AAR 일지 기반 항적 왜곡/분할(Track Split) 오류 분석 및 융합 임계값 최적화 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
# ---- 확장 칼만 필터(EKF) 다중 센서 항적 추적기 공통 함수 ----
# 상태 x = [동(km), 북(km), 동속(km/s), 북속(km/s)], 측정 z = 센서 기준 [거리(km), 방위(rad)]
def ekf_predict(x, P, dt, q):
F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], float)
G = np.array([[dt ** 2 / 2, 0], [0, dt ** 2 / 2], [dt, 0], [0, dt]])
return F @ x, F @ P @ F.T + q * G @ G.T
def h_polar(x, spos):
dx, dy = x[0] - spos[0], x[1] - spos[1]
r2 = max(dx * dx + dy * dy, 1e-9); r = np.sqrt(r2)
return np.array([r, np.arctan2(dx, dy)]), np.array([[dx / r, dy / r, 0, 0], [dy / r2, -dx / r2, 0, 0]])
def innovation(x, P, rep):
s = SENSORS[rep["sensor"]]
zp, H = h_polar(x, s["pos"])
y = np.array([rep["range_km"], np.radians(rep["az_deg"])]) - zp
y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi # 방위 차이는 -pi~pi로 감싼다
S = H @ P @ H.T + np.diag([rep["sr"] ** 2, np.radians(rep["saz"]) ** 2])
return y, H, S, float(y @ np.linalg.solve(S, y))
class Tracker:
"""센서 경보를 스캔(1초) 단위로 받아 EKF 항적을 갱신한다. 같은 스캔 안에서 센서별로
순차 갱신하므로, 여러 센서가 같은 비행체를 보면 하나의 항적으로 통합된다."""
def __init__(self, gate=9.21, q=4e-4, confirm=(3, 5), max_miss=4, merge_gate=6.0):
self.gate, self.q, self.confirm, self.max_miss, self.merge_gate = gate, q, confirm, max_miss, merge_gate
self.tracks, self.dead, self.nid, self.n_merged = [], [], 0, 0
def _new(self, rep, k):
C = meas_cov(rep)
P = np.zeros((4, 4)); P[:2, :2] = C; P[2, 2] = P[3, 3] = 0.25 ** 2
self.tracks.append({"id": f"FT{self.nid:03d}", "x": np.array([rep["x"], rep["y"], 0, 0]), "P": P, "k": k,
"scans": [k], "sensors": {rep["sensor"]: 1}, "truths": [rep["truth"]], "nis": [],
"hist": [(k, rep["x"], rep["y"])], "confirmed": False, "miss": 0, "born": k, "upd": True})
self.nid += 1
def _assign(self, tracks, rs, sid):
"""게이트(마할라노비스 거리) 안에서 전역 최근접 할당(GNN). 할당된 보고서 번호 집합을 돌려준다."""
if not tracks or not rs:
return set()
spos = SENSORS[sid]["pos"]
Z = np.array([[r["range_km"], np.radians(r["az_deg"])] for r in rs])
R = np.diag([rs[0]["sr"] ** 2, np.radians(rs[0]["saz"]) ** 2])
cost = np.full((len(tracks), len(rs)), 1e6)
for i, tr in enumerate(tracks):
zp, H = h_polar(tr["x"], spos)
Y = Z - zp
Y[:, 1] = (Y[:, 1] + np.pi) % (2 * np.pi) - np.pi
Si = np.linalg.inv(H @ tr["P"] @ H.T + R)
d2 = np.einsum("ij,jk,ik->i", Y, Si, Y)
cost[i] = np.where(d2 < self.gate, d2, 1e6)
used = set()
for i, j in zip(*linear_sum_assignment(cost)):
if cost[i, j] < self.gate:
tr, r = tracks[i], rs[j]
y, H, S, d2 = innovation(tr["x"], tr["P"], r)
K = tr["P"] @ H.T @ np.linalg.inv(S)
tr["x"] = tr["x"] + K @ y
tr["P"] = (np.eye(4) - K @ H) @ tr["P"]
tr["nis"].append(d2); tr["truths"].append(r["truth"])
tr["sensors"][sid] = tr["sensors"].get(sid, 0) + 1
tr["upd"] = True; used.add(j)
return used
def _merge(self):
"""같은 비행체를 따라가는 중복 항적(위치·속도가 통계적으로 같음)을 오래된 항적으로 통합한다."""
drop = set()
for a in range(len(self.tracks)):
for b in range(a + 1, len(self.tracks)):
A, B = self.tracks[a], self.tracks[b]
if A["id"] in drop or B["id"] in drop:
continue
d = A["x"] - B["x"]
d2 = float(d[:2] @ np.linalg.solve(A["P"][:2, :2] + B["P"][:2, :2], d[:2]))
if d2 < self.merge_gate and np.hypot(d[2], d[3]) < 0.05:
keep, lose = (A, B) if len(A["scans"]) >= len(B["scans"]) else (B, A)
for sid, n in lose["sensors"].items():
keep["sensors"][sid] = keep["sensors"].get(sid, 0) + n
keep["truths"] += lose["truths"]; keep["scans"] = sorted(set(keep["scans"]) | set(lose["scans"]))
keep["merged"] = keep.get("merged", 0) + 1
drop.add(lose["id"])
self.tracks = [t for t in self.tracks if t["id"] not in drop]
self.n_merged += len(drop)
def step(self, k, reps):
for tr in self.tracks: # 1) 예측
tr["x"], tr["P"] = ekf_predict(tr["x"], tr["P"], k - tr["k"], self.q)
tr["k"], tr["upd"] = k, False
for sid in sorted({r["sensor"] for r in reps}, key=lambda s: SENSORS[s]["sr"]): # 2) 정밀 센서부터 순차 갱신
rs = [r for r in reps if r["sensor"] == sid]
used = self._assign([t for t in self.tracks if t["confirmed"]], rs, sid) # 확정 항적 우선
rest = [r for j, r in enumerate(rs) if j not in used]
used2 = self._assign([t for t in self.tracks if not t["confirmed"]], rest, sid)
for j, r in enumerate(rest):
if j not in used2:
self._new(r, k)
if self.merge_gate > 0:
self._merge()
keep = []
for tr in self.tracks: # 3) 확정·삭제 관리
if tr["upd"]:
tr["miss"] = 0
if tr["scans"][-1] != k:
tr["scans"].append(k)
tr["hist"].append((k, tr["x"][0], tr["x"][1]))
else:
tr["miss"] += 1
m, n = self.confirm
if not tr["confirmed"] and sum(1 for s in tr["scans"] if s > k - n) >= m:
tr["confirmed"] = True
limit = self.max_miss if tr["confirmed"] else 1
(keep if tr["miss"] <= limit else self.dead).append(tr)
self.tracks = keep
def all_tracks(self):
return self.dead + self.tracks
def run_tracker(reps, **kw):
trk = Tracker(**kw)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
for k in range(min(by), max(by) + 1):
trk.step(k, by.get(k, []))
return trk.all_tracks()
def label_track(tr):
vals, cnt = np.unique(tr["truths"], return_counts=True)
return str(vals[np.argmax(cnt)]), float(cnt.max() / cnt.sum())
def evaluate(tracks, scn, min_reports=5):
"""정답 연결 평가: 재현율, 항적 분할(한 표적에 항적 2개 이상), 거짓 항적, 병합 오류, 위치 오차."""
conf = [t for t in tracks if t["confirmed"]]
for t in conf:
t["label"], t["purity"] = label_track(t)
cnt = pd.Series([r["truth"] for r in scn["reports"] if r["truth"].startswith("TGT")]).value_counts()
truths = [k for k, v in cnt.items() if v >= min_reports]
per = {g: [t for t in conf if t["label"] == g] for g in truths}
found = [g for g in truths if per[g]]
split = sum(max(0, len(v) - 1) for v in per.values())
false = [t for t in conf if not t["label"].startswith("TGT")]
merge = [t for t in conf if t["label"].startswith("TGT") and t["purity"] < 0.7]
err = []
for g in found:
for t in per[g]:
for k, x, y in t["hist"][3:]:
if 0 <= k < scn["T"]:
err.append(np.hypot(x - scn["targets"][g][k, 0], y - scn["targets"][g][k, 1]))
n = len(truths)
acc = len(found) / max(n + split + len(false) + len(merge), 1)
return {"정답표적": n, "탐지표적": len(found), "재현율": round(len(found) / max(n, 1), 3), "분할": split,
"거짓항적": len(false), "병합오류": len(merge), "위치RMSE(km)": round(float(np.sqrt(np.mean(np.square(err)))) if err else np.nan, 3),
"항적정확도": round(acc, 3)}
INIT = dict(gate=5.99, q=5e-5, confirm=(2, 3), max_miss=1, merge_gate=0.0) # 1차 배치 파라미터(가상)
AAR_SCN = [dict(seed=1001, special="turn", weather="rain"), dict(seed=1002, special="cross", fa_mult=2.0),
dict(seed=1003, dropout=("RDR-A", 20, 35), weather="fog"), dict(seed=1004, special="fast", ghost=True)]
def run(kw, params, T=40):
scn = make_scenario(T=T, **kw)
trk = run_tracker(prepare(scn), **params)
return evaluate(trk, scn), trk
def errors(trk):
"""오류 사례 목록: 분할(같은 표적 항적 2개 이상), 거짓 항적, 병합 오류."""
conf = [t for t in trk if t["confirmed"]]
out, seen = [], {}
for t in sorted(conf, key=lambda t: t["born"]):
if not t["label"].startswith("TGT"):
out.append(("거짓항적", t["label"], len(t["sensors"]), len(t["scans"])))
elif t["purity"] < 0.7:
out.append(("병합오류", t["label"], len(t["sensors"]), len(t["scans"])))
elif t["label"] in seen:
out.append(("분할", t["label"], len(t["sensors"]), t["born"] - seen[t["label"]]))
seen.setdefault(t["label"], t["born"])
return out
print("[AAR 일지 분석 - 1차 배치 파라미터]", INIT)
log, acc0 = [], []
for i, kw in enumerate(AAR_SCN):
ev, trk = run(kw, INIT)
acc0.append(ev["항적정확도"])
for e in errors(trk):
log.append({"일지": f"AAR-{i + 1}", "조건": kw.get("special") or kw.get("weather") or "-", "오류": e[0], "표적": e[1], "센서수": e[2], "수명/간격(s)": e[3]})
A = pd.DataFrame(log)
print(f"시나리오 {len(AAR_SCN)}개 평균 항적 정확도 {np.mean(acc0):.1%}, 오류 {len(A)}건")
print(A.groupby(["일지", "조건", "오류"]).size().unstack(fill_value=0).to_string())
print("\n[원인 분석: 파라미터 하나씩만 바꿔 AAR 일지 재실행] 가설 -> 정확도 / 분할 / 거짓항적")
HYP = [("gate", 9.21, "게이트(5.99)가 좁아 기동·잡음 시 경보가 항적 밖으로 빠져 새 항적 생성"),
("q", 4e-4, "프로세스 잡음이 작아 선회·고속 기동 예측이 어긋남"),
("confirm", (4, 6), "확정 규칙(2/3)이 느슨해 클러터·기만 경보가 거짓 항적으로 확정"),
("max_miss", 4, "삭제 기준이 짧아 단절·미탐 1~2초에 항적이 끊겼다 재생성"),
("merge_gate", 6.0, "중복 통합 미사용으로 센서별 항적이 합쳐지지 않음")]
for k, v, why in HYP:
evs = [run(kw, {**INIT, k: v})[0] for kw in AAR_SCN]
acc = np.mean([e["항적정확도"] for e in evs])
print(f" {k:10s}={str(v):8s} {acc:6.1%} / {sum(e['분할'] for e in evs):2d} / {sum(e['거짓항적'] for e in evs):2d} ({acc - np.mean(acc0):+.1%}p) {why}")
print(" -> 정확도 개선 폭이 큰 파라미터가 주원인이다. 개선이 없거나 악화되는 가설은 기각하고 조합 탐색으로 최종값을 정한다.")
# 재보정: 별도 튜닝 시나리오에서 격자 탐색
TUNE = [dict(seed=1101, special="turn"), dict(seed=1102, fa_mult=2.0, weather="rain"), dict(seed=1103, dropout=("RDR-B", 10, 25), ghost=True)]
grid = list(itertools.product([9.21], [5e-5, 4e-4], [(2, 3), (3, 5), (4, 6)], [1, 4], [0.0, 6.0])) # 게이트는 카이제곱 99%(9.21)로 고정
best = None
for g, q, c, mm, mg in grid:
p = dict(gate=g, q=q, confirm=c, max_miss=mm, merge_gate=mg)
a = np.mean([run(kw, p, T=30)[0]["항적정확도"] for kw in TUNE])
if best is None or a > best[0]:
best = (a, p)
tuned = best[1]
print(f"\n격자 탐색 {len(grid)}개 조합 -> 튜닝 시나리오 정확도 {best[0]:.1%}")
yaml = "# configs/ops02_tuned_params.yaml [교육용가상자료]\ntracker:\n" + "\n".join(f" {k}: {list(v) if isinstance(v, tuple) else v}" for k, v in tuned.items())
yaml += "\n sync: heartbeat_median\n sensor_order: precise_first\nreview:\n route_to_operator: [split_suspect, single_sensor_over_30s]"
print(yaml)
# 회귀 시험: AAR 시나리오 재발 + 새 시나리오 일반화
recur, fa_b, fa_a = 0, 0, 0
for kw in AAR_SCN:
ev0, trk0 = run(kw, INIT)
ev, trk = run(kw, tuned)
old = {(e[0], e[1]) for e in errors(trk0) if e[0] != "거짓항적"}
new = {(e[0], e[1]) for e in errors(trk) if e[0] != "거짓항적"}
recur += len(old & new); fa_b += ev0["거짓항적"]; fa_a += ev["거짓항적"]
NEW = [dict(seed=1201 + i, **c) for i, c in enumerate([dict(special="turn"), dict(special="formation"), dict(weather="rain", fa_mult=2.0),
dict(dropout=("RDR-A", 10, 30)), dict(ghost=True, weather="night"), dict(special="fast")])]
b = [run(kw, INIT)[0]["항적정확도"] for kw in NEW]
a = [run(kw, tuned)[0]["항적정확도"] for kw in NEW]
print(f"\n[회귀 시험] 새 시나리오 {len(NEW)}개 항적 정확도: 재보정 전 {np.mean(b):.1%} -> 후 {np.mean(a):.1%} ({(np.mean(a) - np.mean(b)) * 100:+.1f}%p, 목표 +10%p: {'달성' if np.mean(a) - np.mean(b) >= 0.1 else '미달'})")
print(f"AAR 일지 오류 재발(같은 표적 분할·병합) {recur}건 (목표 0건: {'달성' if recur == 0 else '미달'}), 같은 일지 거짓 항적 {fa_b}건 -> {fa_a}건")
print("시나리오별 정확도(전->후):", [f"{x:.2f}->{y:.2f}" for x, y in zip(b, a)])
print("지속 학습 연계: 조작관 통합·분리 결정(미션 7 로그)을 정답 라벨로 누적 -> 분기별 재탐색 -> 본 회귀 시험 통과 시에만 파라미터 교체")
plt.figure(figsize=(6, 3.5))
x = np.arange(len(NEW))
plt.bar(x - 0.2, b, 0.4, label="initial"); plt.bar(x + 0.2, a, 0.4, label="retuned")
plt.xticks(x, ["turn", "formation", "rain+FA", "dropout", "ghost", "fast"]); plt.ylabel("track accuracy"); plt.legend(); plt.title("AAR recalibration (synthetic)")
plt.tight_layout(); plt.show()
In [ ]:
"""센서 메시지 검증·정제, 시각 동기화·ENU 변환, 재보정 파라미터를 적용한 EKF 항적 융합, 융합 신뢰도 등급, 조작관 검토 대기열·감사 로그까지 한 번에 실행하는 통합 시뮬레이터를 만들고, 자체 점검(함수 커버리지 포함)과 통합 방공망 브리핑 요약을 출력한다."""
import hashlib
import json
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# ---------------------------------------------------------------------------
# [교육용가상자료] 다중 센서(레이더·EO/IR·음향) 경보 생성기 (실제 장비 제원·배치와 무관)
# 가상 방어 지역(원점 기준 동/북 km)에 비행체 궤적을 만들고, 센서마다 '자기 위치 기준
# 거리·방위'와 '자기 시계 기준 시각'으로 경보를 낸다. 같은 비행체가 여러 센서에 동시에
# 잡히므로 경보가 중복되고, 센서마다 시계 오차·측정 오차·오경보가 다르다.
# ---------------------------------------------------------------------------
SENSORS = {
"RDR-A": {"type": "radar", "pos": (0.0, 0.0), "sr": 0.08, "saz": 0.4, "pd": 0.95, "fa": 1.0, "max_r": 40.0, "clock": 0.0},
"RDR-B": {"type": "radar", "pos": (16.0, -6.0), "sr": 0.10, "saz": 0.5, "pd": 0.93, "fa": 1.0, "max_r": 40.0, "clock": 0.42},
"EOIR-C": {"type": "eoir", "pos": (5.0, 9.0), "sr": 0.40, "saz": 0.15, "pd": 0.88, "fa": 0.5, "max_r": 18.0, "clock": -0.65},
"ACU-D": {"type": "acoustic", "pos": (10.0, 2.0), "sr": 0.60, "saz": 2.0, "pd": 0.80, "fa": 0.6, "max_r": 9.0, "clock": 0.93},
}
CENTER = np.array([8.0, 2.0]) # 가상 방어 지역 중심
KINDS = {"uav": (0.03, 0.06), "heli": (0.04, 0.07), "fast": (0.15, 0.20)} # 속도 범위(km/s, 추상값)
WEATHER = {"clear": {"eoir": 1.0, "acoustic": 1.0, "radar": 1.0, "fa": 1.0},
"night": {"eoir": 0.9, "acoustic": 1.0, "radar": 1.0, "fa": 1.1},
"fog": {"eoir": 0.55, "acoustic": 1.0, "radar": 1.0, "fa": 1.2},
"rain": {"eoir": 0.7, "acoustic": 0.6, "radar": 0.97, "fa": 1.6}}
def _targets(rng, n, T, special):
# [TODO : 학생 실습 10] 다중 센서 중복 항적 통합 시뮬레이터 최종 빌드 및 통합 방공망 브리핑 패키지 완성 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def make_scenario(seed, n=6, T=60, weather="clear", fa_mult=1.0, dropout=None, ghost=False, special=None,
noise_mult=None, drift=0.0):
"""시나리오 1개: 정답 궤적, 센서 경보(보고) 목록, 시계 동기화용 기준 신호(heartbeat)를 만든다.
dropout=(센서, 시작, 끝): 통신 단절, ghost: 한 센서에만 보이는 거짓 항적(기만·스푸핑),
noise_mult={센서: 배수}: 재밍으로 측정 잡음 증가, drift: 음향 센서 시계가 초당 drift초씩 흐트러짐."""
rng = np.random.default_rng(seed)
targets, kinds = _targets(rng, n, T, special)
w = WEATHER[weather]
reports, hb, mid = [], [], 0
gpos = CENTER + np.array([-9.0, 7.0])
for t in range(T):
for sid, s in SENSORS.items():
clock = s["clock"] + (drift * t if sid == "ACU-D" else 0.0)
if t % 10 == 0:
hb.append({"sensor": sid, "t_ref": float(t), "t_sensor": t + clock + rng.normal(0, 0.02)})
if dropout and sid == dropout[0] and dropout[1] <= t < dropout[2]:
continue
nm = (noise_mult or {}).get(sid, 1.0)
objs = [(tid, tr[t, :2]) for tid, tr in targets.items()]
if ghost and sid == "EOIR-C":
objs.append(("GHOST", gpos + 0.12 * t * np.array([0.6, -0.8])))
for tid, xy in objs:
d = xy - np.array(s["pos"])
r = float(np.hypot(*d))
if r > s["max_r"] or rng.random() > s["pd"] * w[s["type"]]:
continue
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(r + rng.normal(0, s["sr"] * nm), 3),
"az_deg": round(float((np.degrees(np.arctan2(d[0], d[1])) + rng.normal(0, s["saz"] * nm)) % 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": tid})
mid += 1
for _ in range(rng.poisson(s["fa"] * fa_mult * w["fa"] * (3 if nm > 1 else 1))): # 오경보(클러터·반사)
reports.append({"msg_id": f"M{mid:05d}", "sensor": sid, "stype": s["type"],
"t_sensor": round(t + clock + rng.normal(0, 0.02), 3),
"range_km": round(rng.uniform(0.5, s["max_r"]), 3), "az_deg": round(rng.uniform(0, 360), 3),
"sr": s["sr"] * nm, "saz": s["saz"] * nm, "truth": "FA"})
mid += 1
return {"targets": targets, "kinds": kinds, "reports": reports, "heartbeats": hb, "T": T, "weather": weather}
# ---------------------------------------------------------------------------
# ---- 시공간 전처리 공통 함수: 정제 -> 시각 동기화 -> 공통 좌표(ENU, km) 변환 ----
def estimate_offsets(heartbeats):
"""센서별 시계 오차 = 기준 신호(heartbeat) 대비 센서 시각 차이의 중앙값(잡음에 강함)."""
# [TODO : 학생 실습 10] 다중 센서 중복 항적 통합 시뮬레이터 최종 빌드 및 통합 방공망 브리핑 패키지 완성 관련 핵심 로직을 구현하세요.
# 힌트: 상단 주석 및 데이터 스키마/지침을 참고하여 작성합니다.
# <<< 여기에 코드를 작성하세요 >>>
pass # 작성 후 pass를 제거하세요.
def clean(reports):
"""중복 메시지·결측·범위 밖(반사파 등 최대 탐지거리 초과) 경보를 걸러낸다. (정제 목록, 사유별 건수)"""
seen, out, why = set(), [], {"중복": 0, "결측": 0, "범위 밖": 0}
for r in reports:
if r["msg_id"] in seen:
why["중복"] += 1; continue
seen.add(r["msg_id"])
if any(r.get(k) is None or (isinstance(r.get(k), float) and np.isnan(r[k])) for k in ("t_sensor", "range_km", "az_deg")):
why["결측"] += 1; continue
if not 0 <= r["range_km"] <= SENSORS[r["sensor"]]["max_r"] * 1.1:
why["범위 밖"] += 1; continue
out.append(r)
return out, why
def to_enu(r, offsets):
"""센서 국소 극좌표(거리·방위) -> 공통 ENU(동·북 km), 센서 시각 -> 기준 시각."""
s, a = SENSORS[r["sensor"]], np.radians(r["az_deg"])
return {**r, "t": r["t_sensor"] - offsets.get(r["sensor"], 0.0),
"x": s["pos"][0] + r["range_km"] * np.sin(a), "y": s["pos"][1] + r["range_km"] * np.cos(a)}
def meas_cov(r):
"""거리·방위 오차를 ENU 공분산(2x2)으로 변환한다(1차 근사)."""
a, rr = np.radians(r["az_deg"]), r["range_km"]
J = np.array([[np.sin(a), rr * np.cos(a)], [np.cos(a), -rr * np.sin(a)]])
return J @ np.diag([r["sr"] ** 2, np.radians(r["saz"]) ** 2]) @ J.T
def prepare(scn, sync=True):
reps, _ = clean(scn["reports"])
off = estimate_offsets(scn["heartbeats"]) if sync else {}
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"])
# ---- 확장 칼만 필터(EKF) 다중 센서 항적 추적기 공통 함수 ----
# 상태 x = [동(km), 북(km), 동속(km/s), 북속(km/s)], 측정 z = 센서 기준 [거리(km), 방위(rad)]
def ekf_predict(x, P, dt, q):
F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], float)
G = np.array([[dt ** 2 / 2, 0], [0, dt ** 2 / 2], [dt, 0], [0, dt]])
return F @ x, F @ P @ F.T + q * G @ G.T
def h_polar(x, spos):
dx, dy = x[0] - spos[0], x[1] - spos[1]
r2 = max(dx * dx + dy * dy, 1e-9); r = np.sqrt(r2)
return np.array([r, np.arctan2(dx, dy)]), np.array([[dx / r, dy / r, 0, 0], [dy / r2, -dx / r2, 0, 0]])
def innovation(x, P, rep):
s = SENSORS[rep["sensor"]]
zp, H = h_polar(x, s["pos"])
y = np.array([rep["range_km"], np.radians(rep["az_deg"])]) - zp
y[1] = (y[1] + np.pi) % (2 * np.pi) - np.pi # 방위 차이는 -pi~pi로 감싼다
S = H @ P @ H.T + np.diag([rep["sr"] ** 2, np.radians(rep["saz"]) ** 2])
return y, H, S, float(y @ np.linalg.solve(S, y))
class Tracker:
"""센서 경보를 스캔(1초) 단위로 받아 EKF 항적을 갱신한다. 같은 스캔 안에서 센서별로
순차 갱신하므로, 여러 센서가 같은 비행체를 보면 하나의 항적으로 통합된다."""
def __init__(self, gate=9.21, q=4e-4, confirm=(3, 5), max_miss=4, merge_gate=6.0):
self.gate, self.q, self.confirm, self.max_miss, self.merge_gate = gate, q, confirm, max_miss, merge_gate
self.tracks, self.dead, self.nid, self.n_merged = [], [], 0, 0
def _new(self, rep, k):
C = meas_cov(rep)
P = np.zeros((4, 4)); P[:2, :2] = C; P[2, 2] = P[3, 3] = 0.25 ** 2
self.tracks.append({"id": f"FT{self.nid:03d}", "x": np.array([rep["x"], rep["y"], 0, 0]), "P": P, "k": k,
"scans": [k], "sensors": {rep["sensor"]: 1}, "truths": [rep["truth"]], "nis": [],
"hist": [(k, rep["x"], rep["y"])], "confirmed": False, "miss": 0, "born": k, "upd": True})
self.nid += 1
def _assign(self, tracks, rs, sid):
"""게이트(마할라노비스 거리) 안에서 전역 최근접 할당(GNN). 할당된 보고서 번호 집합을 돌려준다."""
if not tracks or not rs:
return set()
spos = SENSORS[sid]["pos"]
Z = np.array([[r["range_km"], np.radians(r["az_deg"])] for r in rs])
R = np.diag([rs[0]["sr"] ** 2, np.radians(rs[0]["saz"]) ** 2])
cost = np.full((len(tracks), len(rs)), 1e6)
for i, tr in enumerate(tracks):
zp, H = h_polar(tr["x"], spos)
Y = Z - zp
Y[:, 1] = (Y[:, 1] + np.pi) % (2 * np.pi) - np.pi
Si = np.linalg.inv(H @ tr["P"] @ H.T + R)
d2 = np.einsum("ij,jk,ik->i", Y, Si, Y)
cost[i] = np.where(d2 < self.gate, d2, 1e6)
used = set()
for i, j in zip(*linear_sum_assignment(cost)):
if cost[i, j] < self.gate:
tr, r = tracks[i], rs[j]
y, H, S, d2 = innovation(tr["x"], tr["P"], r)
K = tr["P"] @ H.T @ np.linalg.inv(S)
tr["x"] = tr["x"] + K @ y
tr["P"] = (np.eye(4) - K @ H) @ tr["P"]
tr["nis"].append(d2); tr["truths"].append(r["truth"])
tr["sensors"][sid] = tr["sensors"].get(sid, 0) + 1
tr["upd"] = True; used.add(j)
return used
def _merge(self):
"""같은 비행체를 따라가는 중복 항적(위치·속도가 통계적으로 같음)을 오래된 항적으로 통합한다."""
drop = set()
for a in range(len(self.tracks)):
for b in range(a + 1, len(self.tracks)):
A, B = self.tracks[a], self.tracks[b]
if A["id"] in drop or B["id"] in drop:
continue
d = A["x"] - B["x"]
d2 = float(d[:2] @ np.linalg.solve(A["P"][:2, :2] + B["P"][:2, :2], d[:2]))
if d2 < self.merge_gate and np.hypot(d[2], d[3]) < 0.05:
keep, lose = (A, B) if len(A["scans"]) >= len(B["scans"]) else (B, A)
for sid, n in lose["sensors"].items():
keep["sensors"][sid] = keep["sensors"].get(sid, 0) + n
keep["truths"] += lose["truths"]; keep["scans"] = sorted(set(keep["scans"]) | set(lose["scans"]))
keep["merged"] = keep.get("merged", 0) + 1
drop.add(lose["id"])
self.tracks = [t for t in self.tracks if t["id"] not in drop]
self.n_merged += len(drop)
def step(self, k, reps):
for tr in self.tracks: # 1) 예측
tr["x"], tr["P"] = ekf_predict(tr["x"], tr["P"], k - tr["k"], self.q)
tr["k"], tr["upd"] = k, False
for sid in sorted({r["sensor"] for r in reps}, key=lambda s: SENSORS[s]["sr"]): # 2) 정밀 센서부터 순차 갱신
rs = [r for r in reps if r["sensor"] == sid]
used = self._assign([t for t in self.tracks if t["confirmed"]], rs, sid) # 확정 항적 우선
rest = [r for j, r in enumerate(rs) if j not in used]
used2 = self._assign([t for t in self.tracks if not t["confirmed"]], rest, sid)
for j, r in enumerate(rest):
if j not in used2:
self._new(r, k)
if self.merge_gate > 0:
self._merge()
keep = []
for tr in self.tracks: # 3) 확정·삭제 관리
if tr["upd"]:
tr["miss"] = 0
if tr["scans"][-1] != k:
tr["scans"].append(k)
tr["hist"].append((k, tr["x"][0], tr["x"][1]))
else:
tr["miss"] += 1
m, n = self.confirm
if not tr["confirmed"] and sum(1 for s in tr["scans"] if s > k - n) >= m:
tr["confirmed"] = True
limit = self.max_miss if tr["confirmed"] else 1
(keep if tr["miss"] <= limit else self.dead).append(tr)
self.tracks = keep
def all_tracks(self):
return self.dead + self.tracks
def run_tracker(reps, **kw):
trk = Tracker(**kw)
by = {}
for r in reps:
by.setdefault(int(round(r["t"])), []).append(r)
for k in range(min(by), max(by) + 1):
trk.step(k, by.get(k, []))
return trk.all_tracks()
def label_track(tr):
vals, cnt = np.unique(tr["truths"], return_counts=True)
return str(vals[np.argmax(cnt)]), float(cnt.max() / cnt.sum())
def evaluate(tracks, scn, min_reports=5):
"""정답 연결 평가: 재현율, 항적 분할(한 표적에 항적 2개 이상), 거짓 항적, 병합 오류, 위치 오차."""
conf = [t for t in tracks if t["confirmed"]]
for t in conf:
t["label"], t["purity"] = label_track(t)
cnt = pd.Series([r["truth"] for r in scn["reports"] if r["truth"].startswith("TGT")]).value_counts()
truths = [k for k, v in cnt.items() if v >= min_reports]
per = {g: [t for t in conf if t["label"] == g] for g in truths}
found = [g for g in truths if per[g]]
split = sum(max(0, len(v) - 1) for v in per.values())
false = [t for t in conf if not t["label"].startswith("TGT")]
merge = [t for t in conf if t["label"].startswith("TGT") and t["purity"] < 0.7]
err = []
for g in found:
for t in per[g]:
for k, x, y in t["hist"][3:]:
if 0 <= k < scn["T"]:
err.append(np.hypot(x - scn["targets"][g][k, 0], y - scn["targets"][g][k, 1]))
n = len(truths)
acc = len(found) / max(n + split + len(false) + len(merge), 1)
return {"정답표적": n, "탐지표적": len(found), "재현율": round(len(found) / max(n, 1), 3), "분할": split,
"거짓항적": len(false), "병합오류": len(merge), "위치RMSE(km)": round(float(np.sqrt(np.mean(np.square(err)))) if err else np.nan, 3),
"항적정확도": round(acc, 3)}
CONFIG = {"tracker": dict(gate=9.21, q=5e-5, confirm=(4, 6), max_miss=1, merge_gate=0.0), # 미션 9 재보정값
"high": {"min_sensors": 2, "min_rps": 1.5}}
CALLED = set()
def traced(fn):
def w(*a, **k):
CALLED.add(fn.__name__); return fn(*a, **k)
w.__name__ = fn.__name__
return w
@traced
def validate(r):
errs = []
if r.get("sensor") not in SENSORS:
errs.append("미등록 센서")
for k in ("t_sensor", "range_km", "az_deg"):
if not isinstance(r.get(k), (int, float)) or not np.isfinite(r[k]):
errs.append(f"{k} 결측·형식"); break
return errs
@traced
def ingest(scn):
ok, bad = [], {}
for r in scn["reports"]:
e = validate(r)
(bad.__setitem__(r["msg_id"], e) if e else ok.append(r))
reps, why = clean(ok)
off = estimate_offsets(scn["heartbeats"])
return sorted((to_enu(r, off) for r in reps), key=lambda r: r["t"]), bad, why, off
@traced
def fuse(reps):
return run_tracker(reps, **CONFIG["tracker"])
@traced
def grade(t):
n = len(t["scans"])
hi = len(t["sensors"]) >= CONFIG["high"]["min_sensors"] and len(t["truths"]) / n >= CONFIG["high"]["min_rps"]
return "High" if hi else "Med"
class Queue:
def __init__(self):
self.items, self.log = [], []
@traced
def submit(self, rec):
prev = self.log[-1]["hash"] if self.log else "GENESIS"
body = {"event": "SUBMIT", "id": rec["id"], "grade": rec["grade"], "reason": rec["reason"]}
self.log.append({**body, "prev": prev, "hash": hashlib.sha256((prev + json.dumps(body, ensure_ascii=False, sort_keys=True)).encode()).hexdigest()})
self.items.append({**rec, "status": "PENDING"}) # 자동 통합·분리·교전 연계 없음
@traced
def verify(self):
prev = "GENESIS"
for r in self.log:
body = {k: v for k, v in r.items() if k not in ("prev", "hash")}
if r["prev"] != prev or r["hash"] != hashlib.sha256((prev + json.dumps(body, ensure_ascii=False, sort_keys=True)).encode()).hexdigest():
return False
prev = r["hash"]
return True
@traced
def main(scn):
t0 = time.perf_counter()
reps, bad, why, off = ingest(scn)
trk = fuse(reps)
ev = evaluate(trk, scn)
q, table = Queue(), []
for t in [t for t in trk if t["confirmed"]]:
g = grade(t)
rec = {"id": t["id"], "grade": g, "sensors": "+".join(sorted(t["sensors"])), "scans": len(t["scans"]), "label": t["label"]}
table.append(rec)
if g == "Med":
q.submit({**rec, "reason": "단일 센서·보고 부족(기만·클러터 가능성)"})
ms = (time.perf_counter() - t0) * 1000 / scn["T"]
return pd.DataFrame(table), ev, bad, why, off, q, ms, len(reps)
T0 = time.perf_counter()
# [교육용가상자료] 야간, 표적 8개, 기만 항적, RDR-B 20초 단절, 손상 메시지 4건
scn = make_scenario(1501, n=8, T=60, weather="night", ghost=True, dropout=("RDR-B", 25, 45))
for k, (f, v) in enumerate([("range_km", None), ("sensor", "RDR-Z"), ("az_deg", float("nan")), ("t_sensor", "12:00")]):
bad = dict(scn["reports"][50 + k]); bad[f] = v; bad["msg_id"] = f"BAD{k}"; scn["reports"].append(bad)
D, ev, bad, why, off, q, ms, n_used = main(scn)
elapsed = time.perf_counter() - T0
true_hi = D[(D["grade"] == "High")]["label"].str.startswith("TGT")
chk = {
"손상 메시지 4건 차단": len(bad) == 4,
"메시지 건수 보존(유실 0)": n_used + len(bad) + sum(why.values()) == len(scn["reports"]),
"시계 오차 추정 4개 센서": len(off) == 4,
"표적 재현율 90% 이상": ev["재현율"] >= 0.9,
"High 항적 중 거짓 0건": bool(true_hi.all()),
"Med 항적 전량 검토 대기": len(q.items) == (D["grade"] == "Med").sum(),
"자동 적용 0건": all(i["status"] == "PENDING" for i in q.items),
"감사 로그 무결성": q.verify() and len(q.log) == len(q.items),
"스캔당 200ms 이내": ms < 200,
"항적 분할 0건": ev["분할"] == 0,
}
passed = sum(chk.values())
for k, v in chk.items():
print(f"[{'PASS' if v else 'FAIL'}] {k}")
funcs = ["validate", "ingest", "fuse", "grade", "submit", "verify", "main"]
cov = len(CALLED & set(funcs)) / len(funcs)
print(f"자체 점검 {passed}/{len(chk)} 통과, E2E 실행 성공률 {'100%' if passed == len(chk) else f'{passed / len(chk):.0%}'}, 함수 커버리지 {cov:.0%} (목표 85% 이상: {'달성' if cov >= 0.85 else '미달'})")
print("차단 메시지:", bad)
print("\n" + D.drop(columns="label").to_string(index=False))
n_raw = len(scn["reports"])
print(f"""
[통합 방공망 브리핑 요약 - 교육용가상자료]
1. 처리: 센서 4종 경보 {n_raw}건 수신, 손상 {len(bad)}건 차단, 정제 후 {n_used}건 융합, 전체 {elapsed:.1f}s, 스캔당 {ms:.1f} ms
2. 융합: 원시 경보 {n_raw}건 -> 통합 항적 {len(D)}개 (High {(D['grade'] == 'High').sum()} / Med {(D['grade'] == 'Med').sum()}), 표적 재현율 {ev['재현율']:.0%}, 분할 {ev['분할']}건, 위치 오차 {ev['위치RMSE(km)']} km
3. 검토: Med 항적 {len(q.items)}건 조작관 검토 대기(기만·클러터 가능성), 자동 통합·분리 0건, 감사 로그 무결성 {'정상' if q.verify() else '이상'}
4. 운용: 레이더B 단절 구간도 다른 센서로 추적을 유지했다. 항적 결정은 조작관이, 교전 여부는 교전수칙·지휘계통이 정한다.
5. 한계·확장: 가상 센서·가상 좌표 기준 결과로 실제 센서 제원·배치에 맞춘 재보정이 필요하다. 인접 체계 연동은 통합 항적(JSON) 표준 메시지 API로 확장한다.""")
plt.figure(figsize=(6, 6))
for tid, tr in scn["targets"].items():
plt.plot(tr[:, 0], tr[:, 1], "k-", lw=0.7)
for t in fuse(ingest(scn)[0]):
if t["confirmed"]:
h = np.array(t["hist"]); plt.plot(h[:, 1], h[:, 2], "--", c="tab:blue" if grade(t) == "High" else "tab:orange")
for sid, s in SENSORS.items():
plt.plot(*s["pos"], "^", c="tab:green")
plt.xlabel("East (km)"); plt.ylabel("North (km)"); plt.title("Integrated fusion: High (blue) / Med review (orange) - synthetic"); plt.axis("equal")
plt.tight_layout(); plt.show()