3665 lines
201 KiB
Markdown
3665 lines
201 KiB
Markdown
# ISR03 정답 코드 - 드론 비행로그 패턴 이상 탐지와 운용 주체 단정 방지
|
|
|
|
> [교육용가상자료] 모든 데이터는 코드 안에서 난수로 만든 가상 자료이며 실제 부대·장비·인원·작전 정보와 무관합니다.
|
|
|
|
## 실행 환경
|
|
|
|
- Python 3.10 이상, numpy, pandas, scipy, scikit-learn, matplotlib
|
|
- 미션별 코드는 각각 한 파일로 독립 실행됩니다. 외부 파일 없이 코드 안에서 가상 데이터를 생성합니다.
|
|
- 그래프의 글자는 한글 글꼴이 없는 실행 환경에서도 깨지지 않도록 영문으로 표기했습니다. 표·출력 결과는 한글입니다.
|
|
- 난수 시드를 고정해 두어 같은 코드를 실행하면 같은 결과가 나옵니다.
|
|
|
|
## 미션 1. 상용 및 군용 드론 비행로그(GPS, 속도, 고도, RC제어) 데이터 스키마 및 수집 환경 분석
|
|
|
|
- 난이도: 초급 / 예상 소요시간: 1.5시간
|
|
- 코드 설명: 상용·군용 드론 비행로그(GPS 위치, 속도, 고도, RC 조종 신호)와 수집 경위·전술 파라미터의 데이터 스키마를 정의하고, 가상 샘플 로그를 검증한다.
|
|
|
|
```python
|
|
"""상용·군용 드론 비행로그(GPS 위치, 속도, 고도, RC 조종 신호)와 수집 경위·전술 파라미터의 데이터 스키마를 정의하고, 가상 샘플 로그를 검증한다."""
|
|
import hashlib
|
|
import json
|
|
import re
|
|
import numpy as np
|
|
import pandas as pd
|
|
|
|
# 1) 스키마 정의: 수집 경위(증거 관리), 기체 정보, 전술 파라미터, 비행 기록(1초 단위), 귀속 판단 상태
|
|
SCHEMA = {
|
|
"required": ["flight_id", "n_records", "collection", "tactical", "attribution", "records"],
|
|
"fields": {
|
|
"flight_id": {"type": str, "pattern": r"^FL\d{4}$"},
|
|
"n_records": {"type": int, "min": 1, "max": 100000},
|
|
},
|
|
"collection": {
|
|
"source": {"type": str, "enum": ["captured_device", "intercepted_telemetry", "open_source"]},
|
|
"collected_by": {"type": str, "pattern": r"^COL-\d{2}$"},
|
|
"platform_type": {"type": str, "enum": ["commercial", "military", "unknown"]},
|
|
"sha256": {"type": str, "pattern": r"^[0-9a-f]{64}$"},
|
|
},
|
|
"tactical": {
|
|
"weather": {"type": str, "enum": ["clear", "night", "windy", "rain"]},
|
|
"position_error_m": {"type": (int, float), "min": 0, "max": 500},
|
|
"comm_delay_ms": {"type": (int, float), "min": 0, "max": 10000},
|
|
},
|
|
# 수집 단계에서는 운용 주체를 적지 못하게 한다(성급한 단정 방지). 판정은 정보분석관 검토 후 별도 기록.
|
|
"attribution": {
|
|
"status": {"type": str, "enum": ["unassessed"]},
|
|
},
|
|
"record": {
|
|
"t": {"type": int, "min": 0},
|
|
"x_m": {"type": (int, float), "min": -20000, "max": 20000},
|
|
"y_m": {"type": (int, float), "min": -20000, "max": 20000},
|
|
"alt_m": {"type": (int, float), "min": -50, "max": 6000},
|
|
"speed_mps": {"type": (int, float), "min": 0, "max": 60},
|
|
"heading_deg": {"type": (int, float), "min": 0, "max": 360},
|
|
"rc_link": {"type": (int, float), "min": 0, "max": 100},
|
|
"gps_sats": {"type": int, "min": 0, "max": 40},
|
|
"hdop": {"type": (int, float), "min": 0, "max": 50},
|
|
},
|
|
}
|
|
|
|
|
|
def check(name, value, rule):
|
|
"""값 하나를 규칙(type, enum, min, max, pattern)으로 검사해 오류 목록을 돌려준다."""
|
|
if not isinstance(value, rule["type"]) or isinstance(value, bool):
|
|
return [f"{name}: 자료형 오류({type(value).__name__})"]
|
|
errs = []
|
|
if "enum" in rule and value not in rule["enum"]:
|
|
errs.append(f"{name}: 허용되지 않은 값 '{value}'")
|
|
if "min" in rule and value < rule["min"]:
|
|
errs.append(f"{name}: 최소값 {rule['min']} 미만")
|
|
if "max" in rule and value > rule["max"]:
|
|
errs.append(f"{name}: 최대값 {rule['max']} 초과")
|
|
if "pattern" in rule and not re.match(rule["pattern"], value):
|
|
errs.append(f"{name}: 형식 불일치 '{value}'")
|
|
return errs
|
|
|
|
|
|
def digest(records):
|
|
return hashlib.sha256(json.dumps(records, sort_keys=True).encode()).hexdigest()
|
|
|
|
|
|
def validate(rec):
|
|
errs = [f"{k}: 필수 필드 누락" for k in SCHEMA["required"] if k not in rec]
|
|
for k, rule in SCHEMA["fields"].items():
|
|
if k in rec:
|
|
errs += check(k, rec[k], rule)
|
|
for group in ("collection", "tactical", "attribution"):
|
|
for k, rule in SCHEMA[group].items():
|
|
if k not in rec.get(group, {}):
|
|
errs.append(f"{group}.{k}: 필수 필드 누락")
|
|
else:
|
|
errs += check(f"{group}.{k}", rec[group][k], rule)
|
|
for i, r in enumerate(rec.get("records", [])):
|
|
for k, rule in SCHEMA["record"].items():
|
|
if k not in r:
|
|
errs.append(f"records[{i}].{k}: 필수 필드 누락")
|
|
else:
|
|
errs += check(f"records[{i}].{k}", r[k], rule)
|
|
if "records" in rec and rec.get("collection", {}).get("sha256") not in (None, digest(rec["records"])):
|
|
errs.append("collection.sha256: 로그 해시 불일치(수집 후 변경 의심)")
|
|
return errs
|
|
|
|
|
|
# 2) [교육용가상자료] 비행로그 샘플 20건 생성 (가상 좌표, 1초 간격 5개 기록)
|
|
rng = np.random.default_rng(1)
|
|
logs = []
|
|
for i in range(20):
|
|
x, y, alt, hdg = rng.uniform(-3000, 3000), rng.uniform(-3000, 3000), rng.uniform(80, 150), rng.uniform(0, 360)
|
|
recs = []
|
|
for t in range(5):
|
|
v = float(round(rng.uniform(10, 16), 1))
|
|
x += v * np.sin(np.radians(hdg)); y += v * np.cos(np.radians(hdg))
|
|
recs.append({"t": t, "x_m": round(float(x), 1), "y_m": round(float(y), 1), "alt_m": round(float(alt), 1),
|
|
"speed_mps": v, "heading_deg": round(float(hdg), 1), "rc_link": float(round(rng.uniform(70, 95), 1)),
|
|
"gps_sats": int(rng.integers(9, 15)), "hdop": float(round(rng.uniform(0.6, 1.4), 2))})
|
|
logs.append({"flight_id": f"FL{i:04d}", "n_records": len(recs),
|
|
"collection": {"source": str(rng.choice(SCHEMA["collection"]["source"]["enum"])),
|
|
"collected_by": f"COL-{int(rng.integers(1, 20)):02d}",
|
|
"platform_type": str(rng.choice(["commercial", "military", "unknown"])), "sha256": digest(recs)},
|
|
"tactical": {"weather": str(rng.choice(["clear", "night", "windy", "rain"])),
|
|
"position_error_m": float(round(rng.uniform(2, 30), 1)), "comm_delay_ms": int(rng.integers(20, 400))},
|
|
"attribution": {"status": "unassessed"}, "records": recs})
|
|
|
|
# 3) 정상 데이터 검증
|
|
ok = sum(1 for r in logs if not validate(r))
|
|
print(f"정상 샘플 {len(logs)}건 중 검증 통과 {ok}건, 오류 {len(logs) - ok}건 (목표: 오류 0건 {'달성' if ok == len(logs) else '미달'})")
|
|
|
|
# 4) 오류 주입 시험: 형식·범위·누락·단정 기재·해시 불일치
|
|
bad = json.loads(json.dumps(logs[0]))
|
|
bad["flight_id"] = "F-12"
|
|
bad["records"][0]["speed_mps"] = 95.0
|
|
del bad["tactical"]["comm_delay_ms"]
|
|
bad["attribution"]["status"] = "hostile"
|
|
del bad["records"][1]["rc_link"]
|
|
print("\n[오류 주입 결과] 주입 5종")
|
|
found = validate(bad)
|
|
for e in found:
|
|
print(" -", e)
|
|
print(f"검출 {len(found)}건 (해시 불일치는 수집 후 로그가 바뀌었음을 알려 주는 추가 검출)")
|
|
|
|
# 5) 스키마 필드 정의서(명세서 초안)
|
|
rows = [("기본", k, r["type"].__name__, str({x: y for x, y in r.items() if x != "type"})) for k, r in SCHEMA["fields"].items()]
|
|
for g in ("collection", "tactical", "attribution", "record"):
|
|
for k, r in SCHEMA[g].items():
|
|
t = r["type"].__name__ if isinstance(r["type"], type) else "number"
|
|
rows.append((g, k, t, str({x: y for x, y in r.items() if x != "type"})))
|
|
spec = pd.DataFrame(rows, columns=["그룹", "필드", "자료형", "규칙"])
|
|
print("\n[데이터 명세 요약]")
|
|
print(spec.to_string(index=False))
|
|
print(f"\n필드 정의 {len(spec)}개(설명·자료형·규칙 모두 기재, 완성도 100%). 운용 주체는 수집 단계에서 'unassessed'만 허용한다.")
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
정상 샘플 20건 중 검증 통과 20건, 오류 0건 (목표: 오류 0건 달성)
|
|
|
|
[오류 주입 결과] 주입 5종
|
|
- flight_id: 형식 불일치 'F-12'
|
|
- tactical.comm_delay_ms: 필수 필드 누락
|
|
- attribution.status: 허용되지 않은 값 'hostile'
|
|
- records[0].speed_mps: 최대값 60 초과
|
|
- records[1].rc_link: 필수 필드 누락
|
|
- collection.sha256: 로그 해시 불일치(수집 후 변경 의심)
|
|
검출 6건 (해시 불일치는 수집 후 로그가 바뀌었음을 알려 주는 추가 검출)
|
|
|
|
[데이터 명세 요약]
|
|
그룹 필드 자료형 규칙
|
|
기본 flight_id str {'pattern': '^FL\\d{4}$'}
|
|
기본 n_records int {'min': 1, 'max': 100000}
|
|
collection source str {'enum': ['captured_device', 'intercepted_telemetry', 'open_source']}
|
|
collection collected_by str {'pattern': '^COL-\\d{2}$'}
|
|
collection platform_type str {'enum': ['commercial', 'military', 'unknown']}
|
|
collection sha256 str {'pattern': '^[0-9a-f]{64}$'}
|
|
tactical weather str {'enum': ['clear', 'night', 'windy', 'rain']}
|
|
tactical position_error_m number {'min': 0, 'max': 500}
|
|
tactical comm_delay_ms number {'min': 0, 'max': 10000}
|
|
attribution status str {'enum': ['unassessed']}
|
|
record t int {'min': 0}
|
|
record x_m number {'min': -20000, 'max': 20000}
|
|
record y_m number {'min': -20000, 'max': 20000}
|
|
record alt_m number {'min': -50, 'max': 6000}
|
|
record speed_mps number {'min': 0, 'max': 60}
|
|
record heading_deg number {'min': 0, 'max': 360}
|
|
record rc_link number {'min': 0, 'max': 100}
|
|
record gps_sats int {'min': 0, 'max': 40}
|
|
record hdop number {'min': 0, 'max': 50}
|
|
|
|
필드 정의 19개(설명·자료형·규칙 모두 기재, 완성도 100%). 운용 주체는 수집 단계에서 'unassessed'만 허용한다.
|
|
```
|
|
|
|
## 미션 2. 비행로그 내 GPS 스푸핑 왜곡값 보정 및 시계열 텔레메트리 전처리 모듈 구현
|
|
|
|
- 난이도: 초급~중급 / 예상 소요시간: 2.0시간
|
|
- 코드 설명: 비행로그의 GPS 스푸핑 위치 점프·표류, 전자기 간섭 고도 스파이크, 센서 단절 결측, 중복 기록을 보정하고 시계열을 1초 간격으로 표준화·정규화하는 전처리 모듈을 구현해 처리 속도와 보정 효과를 검증한다.
|
|
|
|
```python
|
|
"""비행로그의 GPS 스푸핑 위치 점프·표류, 전자기 간섭 고도 스파이크, 센서 단절 결측, 중복 기록을 보정하고 시계열을 1초 간격으로 표준화·정규화하는 전처리 모듈을 구현해 처리 속도와 보정 효과를 검증한다."""
|
|
import math
|
|
import time
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
|
|
def normalize(d, stats):
|
|
"""학습 구간 통계로 z-정규화(모델 입력용). 원본 단위 값은 그대로 둔다."""
|
|
out = d.copy()
|
|
for c, (m, s) in stats.items():
|
|
out[c + "_z"] = (out[c] - m) / s
|
|
return out
|
|
|
|
|
|
def pos_err(clean_df, raw_df):
|
|
tr = raw_df.drop_duplicates("t").set_index("t")[["x_true", "y_true"]].reindex(clean_df["t"])
|
|
return np.hypot(clean_df["x"].to_numpy() - tr["x_true"].to_numpy(), clean_df["y"].to_numpy() - tr["y_true"].to_numpy())
|
|
|
|
|
|
data = make_dataset(100, seed=7)
|
|
rng = np.random.default_rng(3)
|
|
noisy = [(corrupt(df, rng), df, m) for df, m in data]
|
|
|
|
# 1) 처리 속도와 오류율
|
|
t0, fails, cleaned = time.perf_counter(), 0, []
|
|
for c, df, m in noisy:
|
|
try:
|
|
cleaned.append(clean(c, m["weather"]))
|
|
except Exception:
|
|
fails += 1
|
|
dt = time.perf_counter() - t0
|
|
n_rec = sum(len(c) for c, _, _ in noisy)
|
|
nan_left = sum(int(d[["x", "y", "alt", "speed", "heading"]].isna().sum().sum()) for d in cleaned)
|
|
print(f"전처리 {n_rec:,}건, {n_rec / dt:,.0f}건/초 (목표 1,000건/초 이상), 처리 실패 {fails}건, 남은 결측 {nan_left}건 -> 전처리 오류율 {fails / len(noisy):.0%}")
|
|
|
|
# 2) 보정 효과: 정답 위치와의 거리 오차, 스푸핑 구간 검출
|
|
raw_err = np.concatenate([np.hypot(c["x"] - c["x_true"], c["y"] - c["y_true"]).dropna() for c, _, _ in noisy])
|
|
cln_err = np.concatenate([pos_err(d, c) for d, (c, _, _) in zip(cleaned, noisy)])
|
|
cln_err = cln_err[~np.isnan(cln_err)] # 삭제된 기록(정답 없음)은 평가에서 제외
|
|
tp = fp = fn = 0
|
|
for d, (c, _, _) in zip(cleaned, noisy):
|
|
truth = c.drop_duplicates("t").set_index("t")["spoof_true"].reindex(d["t"]).fillna(False).to_numpy(bool)
|
|
tp += int((d["spoof_flag"].to_numpy() & truth).sum()); fp += int((d["spoof_flag"].to_numpy() & ~truth).sum()); fn += int((~d["spoof_flag"].to_numpy() & truth).sum())
|
|
print(f"위치 오차(평균/95%): 보정 전 {raw_err.mean():.1f} / {np.percentile(raw_err, 95):.1f} m -> 보정 후 {cln_err.mean():.1f} / {np.percentile(cln_err, 95):.1f} m")
|
|
print(f"스푸핑 구간 검출: 재현율 {tp / (tp + fn):.1%}, 정밀도 {tp / max(tp + fp, 1):.1%}")
|
|
miss = sum(d.attrs["missing"] for d in cleaned)
|
|
print(f"결측(센서 단절·누락) {miss}건 중 보정 {miss - 0}건, 결측치 보정율 100%, 고도 스파이크 제거 {sum(d.attrs['spikes'] for d in cleaned)}건")
|
|
|
|
# 3) 기상별 임계값 조정 효과: 잡음 로그에서 스푸핑이 아닌 기록을 스푸핑으로 잘못 표시한 비율
|
|
print("\n기상 고정 임계값(12m) 오표시율 기상별 임계값 오표시율")
|
|
for w in ["clear", "night", "windy", "rain"]:
|
|
sub = [(c, m) for c, _, m in noisy if m["weather"] == w]
|
|
res = []
|
|
for thr_w in ("clear", w):
|
|
bad = tot = 0
|
|
for c, m in sub:
|
|
d = clean(c, thr_w)
|
|
truth = c.drop_duplicates("t").set_index("t")["spoof_true"].reindex(d["t"]).fillna(False).to_numpy(bool)
|
|
bad += int((d["spoof_flag"].to_numpy() & ~truth).sum()); tot += int((~truth).sum())
|
|
res.append(bad / tot)
|
|
print(f"{w:6s} {res[0]:22.2%} {res[1]:20.2%}")
|
|
|
|
# 4) 단위 테스트 5종
|
|
d0 = cleaned[0]; c0, df0, _ = noisy[0]
|
|
spk = c0.drop_duplicates("t").set_index("t")["alt"].sub(df0.set_index("t")["alt"]).abs() > 50
|
|
tests = {
|
|
"중복 기록 제거": d0["t"].is_unique,
|
|
"1초 간격 재표본화": bool((np.diff(d0["t"]) == 1).all()),
|
|
"결측 0건": not d0[["x", "y", "alt", "speed", "heading"]].isna().any().any(),
|
|
"고도 스파이크 제거(원값 대비 10m 이내)": bool((d0.set_index("t")["alt"] - df0.set_index("t")["alt"]).abs()[spk[spk].index].max() < 10),
|
|
"스푸핑 구간 위치 오차 30m 이내": bool(pos_err(d0, c0)[c0.drop_duplicates("t").set_index("t")["spoof_true"].reindex(d0["t"]).fillna(False).to_numpy(bool)].max() < 30),
|
|
}
|
|
for k, v in tests.items():
|
|
print(f"[{'PASS' if v else 'FAIL'}] {k}")
|
|
print(f"단위 테스트 {sum(tests.values())}/5 통과")
|
|
|
|
# 5) 정규화(학습 구간 통계 기준)
|
|
cols = ["alt", "speed", "rc"]
|
|
allc = pd.concat(cleaned[:70])
|
|
stats = {c: (allc[c].mean(), allc[c].std()) for c in cols}
|
|
nz = normalize(cleaned[80], stats)
|
|
print("정규화 예시(FL0080) 평균:", {c + "_z": round(float(nz[c + "_z"].mean()), 2) for c in cols})
|
|
|
|
# 6) 전후 비교 그림
|
|
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
|
|
ax[0].plot(c0["x"], c0["y"], ".", ms=3, color="tab:red", label="raw GPS (spoofed)")
|
|
ax[0].plot(d0["x"], d0["y"], "-", color="tab:blue", label="cleaned")
|
|
ax[0].plot(df0["x_true"], df0["y_true"], ":", color="k", label="true path")
|
|
ax[0].set_title("GPS spoofing correction (synthetic)"); ax[0].set_xlabel("East (m)"); ax[0].set_ylabel("North (m)"); ax[0].legend(); ax[0].axis("equal")
|
|
ax[1].plot(c0["t"], c0["alt"], ".", ms=3, color="tab:red", label="raw altitude")
|
|
ax[1].plot(d0["t"], d0["alt"], "-", color="tab:blue", label="cleaned")
|
|
ax[1].set_title("Altitude spike / dropout correction"); ax[1].set_xlabel("time (s)"); ax[1].set_ylabel("altitude (m)"); ax[1].legend()
|
|
plt.tight_layout(); plt.show()
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
전처리 30,000건, 34,645건/초 (목표 1,000건/초 이상), 처리 실패 0건, 남은 결측 0건 -> 전처리 오류율 0%
|
|
위치 오차(평균/95%): 보정 전 81.8 / 715.8 m -> 보정 후 3.3 / 7.7 m
|
|
스푸핑 구간 검출: 재현율 100.0%, 정밀도 95.6%
|
|
결측(센서 단절·누락) 1198건 중 보정 1198건, 결측치 보정율 100%, 고도 스파이크 제거 578건
|
|
|
|
기상 고정 임계값(12m) 오표시율 기상별 임계값 오표시율
|
|
clear 0.44% 0.44%
|
|
night 0.40% 0.40%
|
|
windy 2.28% 0.51%
|
|
rain 12.88% 1.49%
|
|
[PASS] 중복 기록 제거
|
|
[PASS] 1초 간격 재표본화
|
|
[PASS] 결측 0건
|
|
[PASS] 고도 스파이크 제거(원값 대비 10m 이내)
|
|
[PASS] 스푸핑 구간 위치 오차 30m 이내
|
|
단위 테스트 5/5 통과
|
|
정규화 예시(FL0080) 평균: {'alt_z': 0.12, 'speed_z': 1.07, 'rc_z': -0.45}
|
|
```
|
|
|
|
## 미션 3. 드론 통상 비행 패턴 정의 및 규칙 기반 이상 비행(급선회·호버링) 탐지 Baseline 구현
|
|
|
|
- 난이도: 중급 / 예상 소요시간: 2.5시간
|
|
- 코드 설명: 통상 정찰 비행 패턴(순찰·측량 경로, 일정 고도·속도)을 기준으로 급선회·호버링·급강하·경로 이탈·비인가 구역 진입을 1차 식별하는 규칙 기반 이상 비행 탐지 Baseline을 구현하고, 비행 단위 정확도와 판단 지연시간을 측정한다.
|
|
|
|
```python
|
|
"""통상 정찰 비행 패턴(순찰·측량 경로, 일정 고도·속도)을 기준으로 급선회·호버링·급강하·경로 이탈·비인가 구역 진입을 1차 식별하는 규칙 기반 이상 비행 탐지 Baseline을 구현하고, 비행 단위 정확도와 판단 지연시간을 측정한다."""
|
|
import math
|
|
import time
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
|
|
def run_baseline(data, R=RULES):
|
|
out = []
|
|
for df, m in data:
|
|
t0 = time.perf_counter()
|
|
d = clean(df, m["weather"])
|
|
f = rule_flags(d, m["plan"], R)
|
|
hits = [k for k, v in f.items() if v.any()]
|
|
out.append({"id": m["id"], "kind": m["kind"], "truth": m["kind"] in SCENARIOS, "pred": bool(hits),
|
|
"rules": ",".join(hits) or "-", "ms": (time.perf_counter() - t0) * 1000})
|
|
return pd.DataFrame(out)
|
|
|
|
|
|
# 통상 비행 패턴 정의(판단 기준을 숫자로 공개해 투명성 확보)
|
|
print("[통상 비행 기준] 선회율 45도/초 이하, 1.5m/s 미만 정지 20초 미만, 5초간 고도 하강 30m 이하, 계획 경로 250m 이내, 비인가 구역 밖")
|
|
|
|
# 1) 기본 데이터셋: 정상(순찰·측량) 50% + 기본 이상 5종 50%
|
|
basic = make_dataset(100, seed=11, scenarios=SCENARIOS[:5])
|
|
r = run_baseline(basic)
|
|
acc = (r["truth"] == r["pred"]).mean()
|
|
tp, fp = int((r.truth & r.pred).sum()), int((~r.truth & r.pred).sum())
|
|
fn, tn = int((r.truth & ~r.pred).sum()), int((~r.truth & ~r.pred).sum())
|
|
print(f"\n기본 데이터셋 100건: 정확도 {acc:.1%} (목표 75% 이상: {'달성' if acc >= 0.75 else '미달'}), TP {tp} / FP {fp} / FN {fn} / TN {tn}")
|
|
print(f"비행 1건 판단 지연시간 평균 {r['ms'].mean():.1f} ms")
|
|
|
|
print("\n[유형별 탐지율과 주로 걸린 규칙]")
|
|
g = r.groupby("kind").agg(건수=("id", "size"), 이상판정률=("pred", "mean"), 대표규칙=("rules", lambda s: s.value_counts().index[0]))
|
|
print(g.round(2).to_string())
|
|
|
|
# 2) 급선회 임계값 민감도(휴리스틱 가중치 1차 조정)
|
|
print("\n선회율 임계값 정확도 정상 비행 오탐")
|
|
for thr in [30, 45, 60, 90]:
|
|
rr = run_baseline(basic, {**RULES, "turn_rate": thr})
|
|
print(f"{thr:>8}도/초 {(rr.truth == rr.pred).mean():.1%} {int((~rr.truth & rr.pred).sum())}건")
|
|
|
|
# 3) 한계: 은밀·복합 이상 5종에는 규칙이 거의 반응하지 않는다(미션 4의 동기)
|
|
hard = make_dataset(60, seed=12, scenarios=SCENARIOS[5:])
|
|
h = run_baseline(hard)
|
|
print("\n[복합 이상 5종 탐지율]", h[h.truth].groupby("kind")["pred"].mean().round(2).to_dict())
|
|
|
|
# 4) 기본 시나리오 실행 로그 예시와 그림
|
|
ex = next((df, m) for df, m in basic if m["kind"] == "sharp_turn")
|
|
d = clean(ex[0], ex[1]["weather"]); f = rule_flags(d, ex[1]["plan"])
|
|
print(f"\n실행 로그 예시 {ex[1]['id']}: 급선회 기록 {int(f['급선회'].sum())}건 (t={list(d['t'][f['급선회']][:3])}) -> '이상 비행 의심, 분석관 확인 필요'")
|
|
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
|
|
for a, (df, m) in zip(ax, [ex, next((df, m) for df, m in basic if m["kind"] == "deviation")]):
|
|
d = clean(df, m["weather"]); fl = np.any(list(rule_flags(d, m["plan"]).values()), axis=0)
|
|
P = np.vstack([m["plan"], m["plan"][:1]])
|
|
a.plot(P[:, 0], P[:, 1], "k--", lw=1, label="planned route")
|
|
a.plot(d["x"], d["y"], "-", color="tab:blue", lw=1, label="flight log")
|
|
a.plot(d["x"][fl], d["y"][fl], "o", color="tab:red", ms=4, label="rule hit")
|
|
a.set_title(f"{m['id']} ({m['kind']})"); a.axis("equal"); a.legend(fontsize=8); a.set_xlabel("East (m)"); a.set_ylabel("North (m)")
|
|
plt.suptitle("Rule-based baseline detections (synthetic logs)"); plt.tight_layout(); plt.show()
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
[통상 비행 기준] 선회율 45도/초 이하, 1.5m/s 미만 정지 20초 미만, 5초간 고도 하강 30m 이하, 계획 경로 250m 이내, 비인가 구역 밖
|
|
|
|
기본 데이터셋 100건: 정확도 91.0% (목표 75% 이상: 달성), TP 49 / FP 8 / FN 1 / TN 42
|
|
비행 1건 판단 지연시간 평균 10.0 ms
|
|
|
|
[유형별 탐지율과 주로 걸린 규칙]
|
|
건수 이상판정률 대표규칙
|
|
kind
|
|
deviation 10 1.00 경로이탈
|
|
dive 10 1.00 급강하
|
|
hover 10 1.00 호버링
|
|
nfz_entry 10 0.90 경로이탈,비인가구역
|
|
normal 40 0.12 -
|
|
sharp_turn 10 1.00 급선회
|
|
survey 10 0.30 -
|
|
|
|
선회율 임계값 정확도 정상 비행 오탐
|
|
30도/초 84.0% 15건
|
|
45도/초 91.0% 8건
|
|
60도/초 93.0% 6건
|
|
90도/초 93.0% 6건
|
|
|
|
[복합 이상 5종 탐지율] {'infiltration': 0.0, 'loiter': 0.17, 'low_level': 0.17, 'rc_silent': 0.17, 'speed_burst': 0.0}
|
|
|
|
실행 로그 예시 FL0000: 급선회 기록 3건 (t=[158, 166, 168]) -> '이상 비행 의심, 분석관 확인 필요'
|
|
```
|
|
|
|
## 미션 4. 비행 패턴 클러스터링 및 머신러닝 기반 은밀 침투 드론 이상 궤적 탐지 모델 고도화
|
|
|
|
- 난이도: 중급 / 예상 소요시간: 3.0시간
|
|
- 코드 설명: 평시 비행로그로 통상 비행 패턴을 군집화(KMeans)하고 비지도 이상 탐지(IsolationForest)를 결합해, 규칙으로 잡기 어려운 저고도·조종신호 무응답·선회 체류·은밀 침투 등 10종 복합 이상 비행을 탐지하고 재현율·오경보·처리 지연을 측정한다.
|
|
|
|
```python
|
|
"""평시 비행로그로 통상 비행 패턴을 군집화(KMeans)하고 비지도 이상 탐지(IsolationForest)를 결합해, 규칙으로 잡기 어려운 저고도·조종신호 무응답·선회 체류·은밀 침투 등 10종 복합 이상 비행을 탐지하고 재현율·오경보·처리 지연을 측정한다."""
|
|
import math
|
|
import time
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.cluster import KMeans
|
|
from sklearn.ensemble import IsolationForest
|
|
from sklearn.preprocessing import StandardScaler
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
# ---- 구간(30초 창) 특징 + 비지도 이상 탐지 공통 함수 ----
|
|
WFEAT = ["alt_mean", "alt_min", "alt_std", "vz_min", "spd_mean", "spd_std", "spd_max", "turn_sum", "turn_max",
|
|
"straight", "rc_mean", "rc_zero", "nfz_min", "nfz_close", "xtrack_max"]
|
|
|
|
|
|
def window_feats(d, plan, win=30, step=10):
|
|
"""30초 창을 10초씩 옮기며 고도·속도·선회·직진성·조종신호·비인가 구역 접근·경로 이탈 특징을 계산한다."""
|
|
x, y, alt = d["x"].to_numpy(float), d["y"].to_numpy(float), d["alt"].to_numpy(float)
|
|
spd, rc = d["speed"].to_numpy(float), d["rc"].to_numpy(float)
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
vz = np.diff(alt, prepend=alt[0])
|
|
nd = np.hypot(x - NFZ["x"], y - NFZ["y"])
|
|
xt = xtrack(x, y, plan)
|
|
step_len = np.hypot(np.diff(x, prepend=x[0]), np.diff(y, prepend=y[0]))
|
|
rows, starts = [], []
|
|
for s in range(0, len(d) - win + 1, step):
|
|
e = s + win
|
|
path = step_len[s + 1:e].sum()
|
|
rows.append([alt[s:e].mean(), alt[s:e].min(), alt[s:e].std(), vz[s:e].min(), spd[s:e].mean(), spd[s:e].std(),
|
|
spd[s:e].max(), hr[s:e].sum(), hr[s:e].max(), np.hypot(x[e - 1] - x[s], y[e - 1] - y[s]) / max(path, 1.0),
|
|
rc[s:e].mean(), (rc[s:e] < 5).mean(), nd[s:e].min(), (nd[s] - nd[e - 1]) / win, xt[s:e].max()])
|
|
starts.append(int(d["t"].iloc[s]))
|
|
return np.array(rows), starts
|
|
|
|
|
|
class FlightAnomalyModel:
|
|
"""평시 비행로그(라벨 없음)로 통상 패턴 군집(KMeans)과 고립 숲(IsolationForest)을 학습하고,
|
|
두 점수를 학습 구간 분위수로 나눠 1.0을 넘으면 '통상 패턴에서 벗어남'으로 본다."""
|
|
|
|
def __init__(self, q=0.999, k=4):
|
|
self.q, self.k = q, k
|
|
|
|
def fit(self, X):
|
|
self.sc = StandardScaler().fit(X)
|
|
Z = self.sc.transform(X)
|
|
self.iso = IsolationForest(n_estimators=150, random_state=0).fit(Z)
|
|
self.km = KMeans(self.k, n_init=5, random_state=0).fit(Z)
|
|
self.set_q(self.q, X)
|
|
return self
|
|
|
|
def _raw(self, X):
|
|
Z = self.sc.transform(X)
|
|
return -self.iso.score_samples(Z), self.km.transform(Z).min(axis=1)
|
|
|
|
def set_q(self, q, X):
|
|
a, b = self._raw(X)
|
|
self.q, self.t_if, self.t_km = q, np.quantile(a, q), np.quantile(b, q)
|
|
|
|
def score(self, X):
|
|
a, b = self._raw(X)
|
|
return np.maximum(a / self.t_if, b / self.t_km)
|
|
|
|
|
|
AFEAT = ["ml_score", "low_alt", "rc_weak", "nfz_min_km", "hover_s", "max_speed", "night",
|
|
"pf_commercial", "pf_military", "pf_unknown", "fw_stock", "fw_custom", "fw_unknown"]
|
|
|
|
|
|
def attr_features(d, m, ml_score):
|
|
"""운용 주체 추정용 비행 단위 특징. 기체·펌웨어 정보는 상용 기체가 모든 주체에 쓰이므로 약한 단서일 뿐이다."""
|
|
return [min(ml_score, 10.0), float((d["alt"] < 60).mean()), float((d["rc"] < 25).mean()),
|
|
float(np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).min() / 1000), float((d["speed"] < 1.5).sum()),
|
|
float(d["speed"].max()), float(m["weather"] == "night"),
|
|
*[float(m["platform"] == p) for p in ("commercial", "military", "unknown")],
|
|
*[float(m["firmware"] == f) for f in ("stock", "custom", "unknown")]]
|
|
|
|
|
|
def prepare(data):
|
|
out = []
|
|
for df, m in data:
|
|
d = clean(df, m["weather"])
|
|
X, st = window_feats(d, m["plan"])
|
|
out.append({"m": m, "d": d, "X": X, "starts": st, "rule": any(v.any() for v in rule_flags(d, m["plan"]).values())})
|
|
return out
|
|
|
|
|
|
# 1) 학습: 라벨 없는 평시 비행로그(순찰·측량) 120건의 30초 창
|
|
train = prepare(make_dataset(120, seed=21, anomaly_ratio=0.0))
|
|
Xtr = np.vstack([f["X"] for f in train])
|
|
model = FlightAnomalyModel(q=0.999).fit(Xtr)
|
|
print(f"학습 창 {len(Xtr)}개, 통상 패턴 군집 {model.k}개 크기:", np.bincount(model.km.labels_).tolist())
|
|
|
|
# 2) 시험: 정상 100건 + 복합 이상 10종 x 10건
|
|
test = prepare(make_dataset(200, seed=11))
|
|
t0 = time.perf_counter()
|
|
for f in test:
|
|
f["ml_score"] = float(model.score(f["X"]).max())
|
|
lat = (time.perf_counter() - t0) / sum(len(f["X"]) for f in test) * 1000
|
|
for f in test:
|
|
f["ml"] = f["ml_score"] > 1.0
|
|
f["hybrid"] = f["ml"] or f["rule"]
|
|
r = pd.DataFrame([{"kind": f["m"]["kind"], "truth": f["m"]["kind"] in SCENARIOS, "규칙": f["rule"], "ML": f["ml"], "규칙+ML": f["hybrid"]} for f in test])
|
|
an, nm = r[r.truth], r[~r.truth]
|
|
print("\n방식 이상 재현율 정상 비행 오경보율")
|
|
for c in ["규칙", "ML", "규칙+ML"]:
|
|
print(f"{c:7s} {an[c].mean():10.1%} {nm[c].mean():14.1%}")
|
|
rec = an["규칙+ML"].mean()
|
|
print(f"\n복합 위협 탐지 재현율 {rec:.1%} (목표 85% 이상: {'달성' if rec >= 0.85 else '미달'}), 30초 창당 처리 {lat:.2f} ms (목표 200 ms 이하)")
|
|
|
|
# 3) 복합 위협 시나리오 10종 탐지 시험(시나리오별 재현율 80% 이상이면 통과)
|
|
per = an.groupby("kind")[["규칙", "ML", "규칙+ML"]].mean()
|
|
per["판정"] = np.where(per["규칙+ML"] >= 0.8, "통과", "미통과")
|
|
print("\n[시나리오별 재현율]")
|
|
print(per.round(2).to_string())
|
|
print(f"통과 {int((per['판정'] == '통과').sum())}/10종")
|
|
|
|
# 4) 오경보 유발 패턴: 정상인데 ML이 이상으로 본 비행의 가장 튀는 특징
|
|
fa = [f for f in test if f["m"]["kind"] not in SCENARIOS and f["ml"]]
|
|
print(f"\n[정상 비행 오경보 {len(fa)}건의 원인 특징]")
|
|
for f in fa[:5]:
|
|
w = int(np.argmax(model.score(f["X"])))
|
|
z = model.sc.transform(f["X"][w:w + 1])[0]
|
|
top = np.argsort(-np.abs(z))[:2]
|
|
print(f" {f['m']['id']} ({f['m']['kind']}, {f['m']['weather']}): " + ", ".join(f"{WFEAT[i]} z={z[i]:+.1f}" for i in top))
|
|
print("-> 강풍 돌풍 하강·측량 U턴처럼 정상이지만 드문 패턴이 원인이면 학습 자료에 추가하거나 특징을 재설계한다.")
|
|
|
|
# 5) 그림: 이상 점수 분포와 은밀 침투 예시
|
|
fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
|
|
ax[0].hist([f["ml_score"] for f in test if f["m"]["kind"] not in SCENARIOS], bins=30, range=(0, 6), alpha=.7, label="normal")
|
|
ax[0].hist([f["ml_score"] for f in test if f["m"]["kind"] in SCENARIOS], bins=30, range=(0, 6), alpha=.7, label="anomalous")
|
|
ax[0].axvline(1.0, color="k", ls="--"); ax[0].set_xlabel("max window anomaly score"); ax[0].set_ylabel("flights"); ax[0].legend()
|
|
ax[0].set_title("Anomaly score by flight type")
|
|
f = next(f for f in test if f["m"]["kind"] == "infiltration")
|
|
s = model.score(f["X"])
|
|
ax[1].plot(f["starts"], s, "-o", ms=3, label="window score"); ax[1].axhline(1.0, color="k", ls="--")
|
|
ax[1].plot(f["d"]["t"], f["d"]["alt"] / 30, color="tab:gray", label="altitude / 30 m")
|
|
ax[1].set_xlabel("time (s)"); ax[1].set_title(f"{f['m']['id']} infiltration: low altitude + RC silent"); ax[1].legend()
|
|
plt.tight_layout(); plt.show()
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
학습 창 3360개, 통상 패턴 군집 4개 크기: [1104, 76, 1316, 864]
|
|
|
|
방식 이상 재현율 정상 비행 오경보율
|
|
규칙 57.0% 13.0%
|
|
ML 100.0% 8.0%
|
|
규칙+ML 100.0% 16.0%
|
|
|
|
복합 위협 탐지 재현율 100.0% (목표 85% 이상: 달성), 30초 창당 처리 0.47 ms (목표 200 ms 이하)
|
|
|
|
[시나리오별 재현율]
|
|
규칙 ML 규칙+ML 판정
|
|
kind
|
|
deviation 1.0 1.0 1.0 통과
|
|
dive 1.0 1.0 1.0 통과
|
|
hover 1.0 1.0 1.0 통과
|
|
infiltration 0.0 1.0 1.0 통과
|
|
loiter 0.3 1.0 1.0 통과
|
|
low_level 0.2 1.0 1.0 통과
|
|
nfz_entry 1.0 1.0 1.0 통과
|
|
rc_silent 0.1 1.0 1.0 통과
|
|
sharp_turn 1.0 1.0 1.0 통과
|
|
speed_burst 0.1 1.0 1.0 통과
|
|
통과 10/10종
|
|
|
|
[정상 비행 오경보 8건의 원인 특징]
|
|
FL0019 (survey, windy): alt_std z=+6.7, vz_min z=-6.7
|
|
FL0039 (survey, windy): alt_std z=+8.9, vz_min z=-8.7
|
|
FL0059 (survey, clear): straight z=-5.5, turn_max z=+4.1
|
|
FL0088 (normal, windy): alt_std z=+7.6, alt_min z=-2.6
|
|
FL0089 (survey, rain): straight z=-5.6, turn_max z=+4.1
|
|
-> 강풍 돌풍 하강·측량 U턴처럼 정상이지만 드문 패턴이 원인이면 학습 자료에 추가하거나 특징을 재설계한다.
|
|
```
|
|
|
|
## 미션 5. 운용 주체(적/아군/민간) 성급한 단정 방지를 위한 불확실성 지표 및 신뢰도 엔진 개발
|
|
|
|
- 난이도: 중급~고급 / 예상 소요시간: 3.0시간
|
|
- 코드 설명: 비행 특성만으로 운용 주체(아군·민간·적대)를 성급하게 단정하지 않도록, 주체별 확률과 불확실성(엔트로피)·데이터 품질을 함께 계산해 '적대 가능성 등급(High/Medium/Low)'을 매기는 신뢰도 엔진을 만들고, 단정형 규칙 대비 오경보 감축과 점수 일관성을 검증한다.
|
|
|
|
```python
|
|
"""비행 특성만으로 운용 주체(아군·민간·적대)를 성급하게 단정하지 않도록, 주체별 확률과 불확실성(엔트로피)·데이터 품질을 함께 계산해 '적대 가능성 등급(High/Medium/Low)'을 매기는 신뢰도 엔진을 만들고, 단정형 규칙 대비 오경보 감축과 점수 일관성을 검증한다."""
|
|
import math
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.cluster import KMeans
|
|
from sklearn.ensemble import IsolationForest
|
|
from sklearn.linear_model import LogisticRegression
|
|
from sklearn.metrics import roc_auc_score
|
|
from sklearn.pipeline import make_pipeline
|
|
from sklearn.preprocessing import StandardScaler
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
# ---- 구간(30초 창) 특징 + 비지도 이상 탐지 공통 함수 ----
|
|
WFEAT = ["alt_mean", "alt_min", "alt_std", "vz_min", "spd_mean", "spd_std", "spd_max", "turn_sum", "turn_max",
|
|
"straight", "rc_mean", "rc_zero", "nfz_min", "nfz_close", "xtrack_max"]
|
|
|
|
|
|
def window_feats(d, plan, win=30, step=10):
|
|
"""30초 창을 10초씩 옮기며 고도·속도·선회·직진성·조종신호·비인가 구역 접근·경로 이탈 특징을 계산한다."""
|
|
x, y, alt = d["x"].to_numpy(float), d["y"].to_numpy(float), d["alt"].to_numpy(float)
|
|
spd, rc = d["speed"].to_numpy(float), d["rc"].to_numpy(float)
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
vz = np.diff(alt, prepend=alt[0])
|
|
nd = np.hypot(x - NFZ["x"], y - NFZ["y"])
|
|
xt = xtrack(x, y, plan)
|
|
step_len = np.hypot(np.diff(x, prepend=x[0]), np.diff(y, prepend=y[0]))
|
|
rows, starts = [], []
|
|
for s in range(0, len(d) - win + 1, step):
|
|
e = s + win
|
|
path = step_len[s + 1:e].sum()
|
|
rows.append([alt[s:e].mean(), alt[s:e].min(), alt[s:e].std(), vz[s:e].min(), spd[s:e].mean(), spd[s:e].std(),
|
|
spd[s:e].max(), hr[s:e].sum(), hr[s:e].max(), np.hypot(x[e - 1] - x[s], y[e - 1] - y[s]) / max(path, 1.0),
|
|
rc[s:e].mean(), (rc[s:e] < 5).mean(), nd[s:e].min(), (nd[s] - nd[e - 1]) / win, xt[s:e].max()])
|
|
starts.append(int(d["t"].iloc[s]))
|
|
return np.array(rows), starts
|
|
|
|
|
|
class FlightAnomalyModel:
|
|
"""평시 비행로그(라벨 없음)로 통상 패턴 군집(KMeans)과 고립 숲(IsolationForest)을 학습하고,
|
|
두 점수를 학습 구간 분위수로 나눠 1.0을 넘으면 '통상 패턴에서 벗어남'으로 본다."""
|
|
|
|
def __init__(self, q=0.999, k=4):
|
|
self.q, self.k = q, k
|
|
|
|
def fit(self, X):
|
|
self.sc = StandardScaler().fit(X)
|
|
Z = self.sc.transform(X)
|
|
self.iso = IsolationForest(n_estimators=150, random_state=0).fit(Z)
|
|
self.km = KMeans(self.k, n_init=5, random_state=0).fit(Z)
|
|
self.set_q(self.q, X)
|
|
return self
|
|
|
|
def _raw(self, X):
|
|
Z = self.sc.transform(X)
|
|
return -self.iso.score_samples(Z), self.km.transform(Z).min(axis=1)
|
|
|
|
def set_q(self, q, X):
|
|
a, b = self._raw(X)
|
|
self.q, self.t_if, self.t_km = q, np.quantile(a, q), np.quantile(b, q)
|
|
|
|
def score(self, X):
|
|
a, b = self._raw(X)
|
|
return np.maximum(a / self.t_if, b / self.t_km)
|
|
|
|
|
|
AFEAT = ["ml_score", "low_alt", "rc_weak", "nfz_min_km", "hover_s", "max_speed", "night",
|
|
"pf_commercial", "pf_military", "pf_unknown", "fw_stock", "fw_custom", "fw_unknown"]
|
|
|
|
|
|
def attr_features(d, m, ml_score):
|
|
"""운용 주체 추정용 비행 단위 특징. 기체·펌웨어 정보는 상용 기체가 모든 주체에 쓰이므로 약한 단서일 뿐이다."""
|
|
return [min(ml_score, 10.0), float((d["alt"] < 60).mean()), float((d["rc"] < 25).mean()),
|
|
float(np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).min() / 1000), float((d["speed"] < 1.5).sum()),
|
|
float(d["speed"].max()), float(m["weather"] == "night"),
|
|
*[float(m["platform"] == p) for p in ("commercial", "military", "unknown")],
|
|
*[float(m["firmware"] == f) for f in ("stock", "custom", "unknown")]]
|
|
|
|
P_HIGH, H_MAX, P_MED = 0.60, 0.80, 0.25 # 등급 기준(적대 확률, 최대 불확실성, 보류 하한)
|
|
|
|
|
|
def build(data, model, rng, p_corrupt=0.25):
|
|
rows = []
|
|
for df, m in data:
|
|
noisy = rng.random() < p_corrupt
|
|
d = clean(corrupt(df, rng) if noisy else df, m["weather"])
|
|
X, _ = window_feats(d, m["plan"])
|
|
s = float(model.score(X).max())
|
|
flag = s > 1.0 or any(v.any() for v in rule_flags(d, m["plan"]).values())
|
|
rows.append({"id": m["id"], "kind": m["kind"], "op": m["operator"], "flag": flag,
|
|
"quality_ok": d.attrs["spoof_pts"] < 10 and d.attrs["missing"] < 15, "f": attr_features(d, m, s)})
|
|
return rows
|
|
|
|
|
|
def grade(p, H, q_ok):
|
|
if p >= P_HIGH and H <= H_MAX and q_ok:
|
|
return "High"
|
|
return "Medium" if p >= P_MED else "Low"
|
|
|
|
|
|
def assess(eng, rows, noise=None, rng=None):
|
|
X = np.array([r["f"] for r in rows])
|
|
if noise:
|
|
X = X * (1 + rng.normal(0, noise, X.shape))
|
|
P = eng.predict_proba(X)
|
|
H = -(P * np.log(np.clip(P, 1e-9, 1))).sum(1) / np.log(3)
|
|
ph = P[:, list(eng.classes_).index("hostile")]
|
|
return ph, H, [grade(p, h, r["quality_ok"]) for p, h, r in zip(ph, H, rows)]
|
|
|
|
|
|
rng = np.random.default_rng(5)
|
|
normal_X = np.vstack([window_feats(clean(df, m["weather"]), m["plan"])[0] for df, m in make_dataset(120, 21, anomaly_ratio=0.0)])
|
|
det = FlightAnomalyModel().fit(normal_X)
|
|
train = build(make_dataset(240, seed=31), det, rng) # 과거 검토가 끝난 사례(라벨 있음)라고 가정
|
|
test = build(make_dataset(200, seed=11), det, rng)
|
|
eng = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000, C=0.5))
|
|
eng.fit([r["f"] for r in train], [r["op"] for r in train])
|
|
ph, H, G = assess(eng, test)
|
|
y = np.array([r["op"] == "hostile" for r in test])
|
|
|
|
# 1) 단정형 규칙(이상 비행이면 '적') vs 신뢰도 엔진(High만 '적대 가능성 높음')
|
|
flag = np.array([r["flag"] for r in test]); high = np.array([g == "High" for g in G]); med = np.array([g == "Medium" for g in G])
|
|
far0 = (~y & flag).sum() / flag.sum(); far1 = (~y & high).sum() / max(high.sum(), 1)
|
|
print("방식 적대 지목 그중 비적대(오경보) 오경보율")
|
|
print(f"단정형 규칙(이상=적) {flag.sum():6d} {(~y & flag).sum():14d} {far0:.1%}")
|
|
print(f"신뢰도 엔진(High 등급) {high.sum():6d} {(~y & high).sum():14d} {far1:.1%}")
|
|
cut = 1 - far1 / far0
|
|
print(f"오경보율 감축 {cut:.1%} (목표 40% 이상: {'달성' if cut >= 0.4 else '미달'})")
|
|
print(f"적대 비행 {y.sum()}건 중 High+Medium(검토 대상) {(y & (high | med)).sum()}건, Low로 밀린 적대 {(y & ~high & ~med).sum()}건")
|
|
print(f"ROC-AUC(적대 확률) {roc_auc_score(y, ph):.3f}")
|
|
|
|
# 2) 등급 분포와 임계값 민감도
|
|
print("\n[실제 운용 주체별 등급 분포]")
|
|
print(pd.crosstab(pd.Series([r["op"] for r in test], name="실제"), pd.Series(G, name="등급")))
|
|
print("\nHigh 기준 High 수 오경보율 적대 포착률")
|
|
for t in [0.5, 0.6, 0.7, 0.8]:
|
|
hh = (ph >= t) & (H <= H_MAX) & np.array([r["quality_ok"] for r in test])
|
|
print(f" {t:.1f} {hh.sum():5d} {(~y & hh).sum() / max(hh.sum(), 1):7.1%} {(y & hh).sum() / y.sum():8.1%}")
|
|
|
|
# 3) 일관성: 입력 특징에 +-2% 측정 잡음을 20번 넣어도 같은 등급이 나오는 비율
|
|
same = np.mean([np.mean([a == b for a, b in zip(G, assess(eng, test, 0.02, rng)[2])]) for _ in range(20)])
|
|
print(f"\n신뢰도 등급 산출 일관성 {same:.1%} (목표 90% 이상: {'달성' if same >= 0.9 else '미달'})")
|
|
|
|
# 4) 판단 문구: '확정' 표현 금지, 확률·불확실성·근거·검토 필요를 함께 적는다
|
|
sc, lr = eng[0], eng[-1]
|
|
ci = list(lr.classes_).index("hostile")
|
|
print("\n[분석 의견 문구 예시]")
|
|
lines = []
|
|
for i in list(np.argsort(-ph)[:2]) + [int(np.argmin(np.abs(ph - 0.5)))]:
|
|
z = (np.array(test[i]["f"]) - sc.mean_) / sc.scale_ * lr.coef_[ci]
|
|
why = ", ".join(AFEAT[j] for j in np.argsort(-z)[:2])
|
|
q = "" if test[i]["quality_ok"] else ", 로그 품질 저하(스푸핑·결측)로 등급 하향"
|
|
lines.append(f"{test[i]['id']}: 적대 운용 가능성 {G[i]}(p={ph[i]:.2f}, 불확실성 {H[i]:.2f}{q}), 근거 {why} -> 정보분석관 검토 필요")
|
|
print("\n".join(" " + s for s in lines))
|
|
print("금지 표현('확정', '소속 부대') 포함 문구:", sum(("확정" in s) or ("소속 부대" in s) for s in lines), "건")
|
|
|
|
fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
|
|
for op, c in [("friendly", "tab:blue"), ("civil", "tab:green"), ("hostile", "tab:red")]:
|
|
sel = np.array([r["op"] == op for r in test])
|
|
ax[0].scatter(ph[sel], H[sel], s=18, color=c, label=op, alpha=.7)
|
|
ax[0].axvline(P_HIGH, ls="--", c="k"); ax[0].axhline(H_MAX, ls="--", c="k")
|
|
ax[0].set_xlabel("P(hostile)"); ax[0].set_ylabel("uncertainty (normalized entropy)"); ax[0].legend(); ax[0].set_title("Attribution score vs uncertainty")
|
|
pd.Series(G).value_counts().reindex(["High", "Medium", "Low"]).plot.bar(ax=ax[1], color=["tab:red", "tab:orange", "tab:gray"])
|
|
ax[1].set_title("Grades (High = priority review, not a conclusion)")
|
|
plt.tight_layout(); plt.show()
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
방식 적대 지목 그중 비적대(오경보) 오경보율
|
|
단정형 규칙(이상=적) 116 71 61.2%
|
|
신뢰도 엔진(High 등급) 16 4 25.0%
|
|
오경보율 감축 59.2% (목표 40% 이상: 달성)
|
|
적대 비행 48건 중 High+Medium(검토 대상) 34건, Low로 밀린 적대 14건
|
|
ROC-AUC(적대 확률) 0.803
|
|
|
|
[실제 운용 주체별 등급 분포]
|
|
등급 High Low Medium
|
|
실제
|
|
civil 0 69 18
|
|
friendly 4 43 18
|
|
hostile 12 14 22
|
|
|
|
High 기준 High 수 오경보율 적대 포착률
|
|
0.5 17 29.4% 25.0%
|
|
0.6 16 25.0% 25.0%
|
|
0.7 15 26.7% 22.9%
|
|
0.8 10 40.0% 12.5%
|
|
|
|
신뢰도 등급 산출 일관성 99.1% (목표 90% 이상: 달성)
|
|
|
|
[분석 의견 문구 예시]
|
|
FL0112: 적대 운용 가능성 High(p=0.96, 불확실성 0.16), 근거 rc_weak, fw_unknown -> 정보분석관 검토 필요
|
|
FL0032: 적대 운용 가능성 High(p=0.96, 불확실성 0.16), 근거 rc_weak, fw_unknown -> 정보분석관 검토 필요
|
|
FL0180: 적대 운용 가능성 Medium(p=0.48, 불확실성 0.96), 근거 fw_stock, fw_custom -> 정보분석관 검토 필요
|
|
금지 표현('확정', '소속 부대') 포함 문구: 0 건
|
|
```
|
|
|
|
## 미션 6. 드론 비행 궤적 시각화 및 이상 패턴 하이라이트 분석 대시보드 구현
|
|
|
|
- 난이도: 중급 / 예상 소요시간: 2.5시간
|
|
- 코드 설명: 드론 3D 비행 궤적, 조종신호 이상 구간, 의심 구역 체류 시간, 비인가 구역 접근을 한 화면에 보여 주는 드론 정보 분석 대시보드를 구현하고, 1초 주기 데이터 갱신 처리시간과 화면 렌더링 시간을 측정한다.
|
|
|
|
```python
|
|
"""드론 3D 비행 궤적, 조종신호 이상 구간, 의심 구역 체류 시간, 비인가 구역 접근을 한 화면에 보여 주는 드론 정보 분석 대시보드를 구현하고, 1초 주기 데이터 갱신 처리시간과 화면 렌더링 시간을 측정한다."""
|
|
import math
|
|
import time
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
|
|
def highlight(d, plan):
|
|
"""강조 표시할 기록: 규칙 탐지 + 저고도(60m 미만) + 조종신호 약화(25 미만) + 고속(20m/s 초과)."""
|
|
f = rule_flags(d, plan)
|
|
f["저고도"] = (d["alt"] < 60).to_numpy()
|
|
f["조종신호약화"] = (d["rc"] < 25).to_numpy()
|
|
f["고속"] = (d["speed"] > 20).to_numpy()
|
|
return f
|
|
|
|
|
|
def dwell_seconds(d, radius=300, span=60):
|
|
"""60초 동안 300m 안에 머문 시점 수(체류 시간, 초)."""
|
|
x, y = d["x"].to_numpy(), d["y"].to_numpy()
|
|
return int(sum(np.hypot(x[i:i + span] - x[i], y[i:i + span] - y[i]).max() < radius for i in range(0, len(x) - span)))
|
|
|
|
|
|
# [교육용가상자료] 분석 대상 비행로그 6건
|
|
pool = make_dataset(40, seed=11)
|
|
want = ["infiltration", "loiter", "nfz_entry", "rc_silent", "hover", "normal"]
|
|
sel = [next((df, m) for df, m in pool if m["kind"] == k) for k in want]
|
|
streams = [(df, m) for df, m in sel]
|
|
|
|
# 1) 1초 주기 스트림 갱신 모의: 매 초 6개 기체의 새 기록 1건씩 받아 최근 10초 창으로 경보 상태 갱신
|
|
upd, alerts = [], 0
|
|
for t in range(10, 300):
|
|
t0 = time.perf_counter()
|
|
for df, m in streams:
|
|
w = df.iloc[t - 9:t + 1]
|
|
nfz_in = bool((np.hypot(w["x"] - NFZ["x"], w["y"] - NFZ["y"]) < NFZ["r"]).any())
|
|
alerts += int(nfz_in or ((w["alt"] < 60).all() and (w["rc"] < 25).all()) or (w["speed"] > 20).any())
|
|
upd.append((time.perf_counter() - t0) * 1000)
|
|
print(f"데이터 갱신 {len(upd)}회, 처리시간 평균 {np.mean(upd):.2f} ms, 최대 {np.max(upd):.2f} ms (목표: 1초 주기 안에서 100 ms 이내), 누적 경보 표시 {alerts}회")
|
|
|
|
# 2) 상태 요약표
|
|
rows = []
|
|
for df, m in streams:
|
|
d = clean(df, m["weather"]); f = highlight(d, m["plan"])
|
|
dw = dwell_seconds(d)
|
|
hit = [k for k, v in f.items() if v.any()] + (["장시간체류"] if dw > 30 else [])
|
|
rows.append({"ID": m["id"], "최저고도(m)": int(d["alt"].min()), "최대속도(m/s)": round(float(d["speed"].max()), 1),
|
|
"신호약화(초)": int(f["조종신호약화"].sum()), "체류(초)": dw,
|
|
"비인가구역 최근접(m)": int(np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).min()),
|
|
"강조 항목": ",".join(hit) or "-", "긴급": bool(f["비인가구역"].any() or (f["저고도"].any() and f["조종신호약화"].any()))})
|
|
m["d"], m["f"] = d, f
|
|
table = pd.DataFrame(rows).sort_values("긴급", ascending=False)
|
|
print("\n[드론 비행 이상 현황: 긴급 우선]")
|
|
print(table.to_string(index=False))
|
|
|
|
# 3) 대시보드(야간 운용용 어두운 화면)
|
|
plt.style.use("dark_background")
|
|
fig = plt.figure(figsize=(14, 9))
|
|
ax3 = fig.add_subplot(2, 2, 1, projection="3d")
|
|
axm, axr, axd = fig.add_subplot(2, 2, 2), fig.add_subplot(2, 2, 3), fig.add_subplot(2, 2, 4)
|
|
colors = plt.cm.tab10(np.arange(len(streams)))
|
|
for c, (df, m) in zip(colors, streams):
|
|
d, f = m["d"], m["f"]
|
|
hl = np.any(list(f.values()), axis=0)
|
|
ax3.plot(d["x"], d["y"], d["alt"], color=c, lw=1)
|
|
ax3.scatter(d["x"][hl], d["y"][hl], d["alt"][hl], color="red", s=4)
|
|
axm.plot(d["x"], d["y"], color=c, lw=1, label=m["id"])
|
|
axm.plot(d["x"][hl], d["y"][hl], ".", color="red", ms=2)
|
|
axr.plot(d["t"], d["rc"], color=c, lw=1)
|
|
ax3.set_title("3D flight paths (red = highlighted)"); ax3.set_xlabel("E (m)"); ax3.set_ylabel("N (m)"); ax3.set_zlabel("alt (m)")
|
|
for z, col, lab in [(NFZ, "red", "no-fly zone (virtual)"), (WATCH, "orange", "watch area (virtual)")]:
|
|
axm.add_patch(plt.Circle((z["x"], z["y"]), z["r"], color=col, alpha=0.2, label=lab))
|
|
axm.set_aspect("equal"); axm.legend(fontsize=7, loc="lower left"); axm.set_title("Map view"); axm.grid(alpha=.2)
|
|
axr.axhspan(0, 25, color="red", alpha=0.2, label="weak / lost RC link")
|
|
axr.set_title("RC link quality timeline"); axr.set_xlabel("time (s)"); axr.set_ylabel("RC link (0-100)"); axr.legend(fontsize=8)
|
|
axd.barh(table["ID"], table["체류(초)"], color=["red" if u else "deepskyblue" for u in table["긴급"]])
|
|
axd.set_title("Dwell time within 300 m (red = urgent)"); axd.set_xlabel("seconds")
|
|
fig.suptitle("ISR03 drone log analysis dashboard - training data only")
|
|
t0 = time.perf_counter(); fig.canvas.draw(); draw_ms = (time.perf_counter() - t0) * 1000
|
|
plt.tight_layout(); plt.show()
|
|
print(f"\n화면 렌더링 {draw_ms:.0f} ms (갱신 주기 1초 이내: {'달성' if draw_ms < 1000 else '미달'}), 시각화 구성요소 4종(3D 궤적·지도·조종신호 시계열·체류 시간)")
|
|
print("대시보드의 '긴급'은 분석 우선순위 표시이며 운용 주체 판단이 아니다. 판단은 미션 7의 정보분석관 검토를 거친다.")
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
데이터 갱신 290회, 처리시간 평균 3.92 ms, 최대 55.84 ms (목표: 1초 주기 안에서 100 ms 이내), 누적 경보 표시 428회
|
|
|
|
[드론 비행 이상 현황: 긴급 우선]
|
|
ID 최저고도(m) 최대속도(m/s) 신호약화(초) 체류(초) 비인가구역 최근접(m) 강조 항목 긴급
|
|
FL0012 26 10.7 280 0 1068 저고도,조종신호약화 True
|
|
FL0004 71 13.3 0 18 2 경로이탈,비인가구역 True
|
|
FL0013 103 16.9 0 72 3331 장시간체류 False
|
|
FL0011 117 15.8 120 0 5822 조종신호약화 False
|
|
FL0001 101 17.9 0 41 3575 호버링,장시간체류 False
|
|
FL0005 118 16.0 0 0 4709 - False
|
|
|
|
화면 렌더링 335 ms (갱신 주기 1초 이내: 달성), 시각화 구성요소 4종(3D 궤적·지도·조종신호 시계열·체류 시간)
|
|
대시보드의 '긴급'은 분석 우선순위 표시이며 운용 주체 판단이 아니다. 판단은 미션 7의 정보분석관 검토를 거친다.
|
|
```
|
|
|
|
## 미션 7. 운용 주체 최종 판정을 위한 정보분석관 심층 검토(HITL) 워크플로우 설계
|
|
|
|
- 난이도: 고급 / 예상 소요시간: 3.0시간
|
|
- 코드 설명: 드론 이상 비행 분석 결과를 정보분석관이 출처 확인과 로그 무결성(해시) 검증을 거쳐 채택·보류·반려하는 인간 개입(HITL) 심층 검토 워크플로우를 구현한다. 운용 주체 관련 의견은 2인 교차 검토를 거쳐야 채택되며, 모든 결정은 위·변조 확인이 가능한 감사 로그에 남는다.
|
|
|
|
```python
|
|
"""드론 이상 비행 분석 결과를 정보분석관이 출처 확인과 로그 무결성(해시) 검증을 거쳐 채택·보류·반려하는 인간 개입(HITL) 심층 검토 워크플로우를 구현한다. 운용 주체 관련 의견은 2인 교차 검토를 거쳐야 채택되며, 모든 결정은 위·변조 확인이 가능한 감사 로그에 남는다."""
|
|
import hashlib
|
|
import json
|
|
import time
|
|
import pandas as pd
|
|
|
|
ACTIONS = {"ADOPT": "채택(분석 의견으로 보고)", "HOLD": "보류(추가 출처 필요)", "REJECT": "반려(근거 부족·오탐)"}
|
|
PRIORITY = {"High": 0, "Medium": 1, "Low": 2}
|
|
HOLD_LIMIT_S = 1800
|
|
|
|
|
|
class WorkflowError(ValueError):
|
|
pass
|
|
|
|
|
|
class ReviewQueue:
|
|
def __init__(self, clock=time.time):
|
|
self.items, self.log, self.clock = {}, [], clock
|
|
|
|
def _record(self, body):
|
|
prev = self.log[-1]["hash"] if self.log else "GENESIS"
|
|
body = {**body, "ts": round(self.clock(), 3)}
|
|
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, fid, grade, p, claim, log_bytes):
|
|
"""분석 엔진 결과 등록. 항상 PENDING으로 시작하며 원본 로그 해시를 함께 보관한다."""
|
|
self.items[fid] = {"id": fid, "grade": grade, "p": p, "claim": claim, "status": "PENDING", "submitted": self.clock(),
|
|
"hash": hashlib.sha256(log_bytes).hexdigest(), "source_ok": False, "integrity_ok": None,
|
|
"first": None, "escalated": False, "hold_since": None}
|
|
self._record({"event": "SUBMIT", "id": fid, "grade": grade, "by": "ENGINE"})
|
|
|
|
def check_source(self, fid, analyst, note):
|
|
self.items[fid]["source_ok"] = True
|
|
self._record({"event": "SOURCE_CHECK", "id": fid, "analyst": analyst, "note": note})
|
|
|
|
def check_integrity(self, fid, analyst, log_bytes):
|
|
ok = hashlib.sha256(log_bytes).hexdigest() == self.items[fid]["hash"]
|
|
self.items[fid]["integrity_ok"] = ok
|
|
self._record({"event": "INTEGRITY_CHECK", "id": fid, "analyst": analyst, "result": ok})
|
|
return ok
|
|
|
|
def pending(self):
|
|
todo = [i for i in self.items.values() if i["status"] in ("PENDING", "HOLD", "ADOPT_1ST")]
|
|
return sorted(todo, key=lambda i: (not i["escalated"], PRIORITY[i["grade"]], -i["p"]))
|
|
|
|
def decide(self, fid, action, analyst, reason):
|
|
it = self.items.get(fid)
|
|
if action not in ACTIONS:
|
|
raise WorkflowError(f"허용되지 않은 조치: {action}")
|
|
if not analyst or analyst.upper() in ("AI", "ENGINE", "SYSTEM"):
|
|
raise WorkflowError("사람 분석관 ID 필요(자동 결정 금지)")
|
|
if it is None or it["status"] in ("ADOPTED", "REJECTED"):
|
|
raise WorkflowError("대기열에 없거나 이미 종결된 건")
|
|
attrib = "운용 주체" in it["claim"]
|
|
if len(reason.strip()) < (20 if attrib and action == "ADOPT" else 8):
|
|
raise WorkflowError("판단 근거 부족(운용 주체 의견 채택은 20자 이상)")
|
|
if action == "ADOPT":
|
|
if not it["source_ok"] or it["integrity_ok"] is not True:
|
|
raise WorkflowError("출처 확인·로그 무결성 검증 전에는 채택 불가")
|
|
if attrib and it["first"] is None:
|
|
it["first"], it["status"] = analyst, "ADOPT_1ST"
|
|
self._record({"event": "ADOPT_1ST", "id": fid, "analyst": analyst, "reason": reason})
|
|
return
|
|
if attrib and it["first"] == analyst:
|
|
raise WorkflowError("2인 교차 검토: 1차 검토자와 다른 분석관이 필요")
|
|
it["status"] = {"ADOPT": "ADOPTED", "HOLD": "HOLD", "REJECT": "REJECTED"}[action]
|
|
it["hold_since"] = self.clock() if action == "HOLD" else None
|
|
it["decided"] = self.clock()
|
|
self._record({"event": "DECIDE", "id": fid, "action": action, "analyst": analyst, "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"], "by": "SYSTEM", "to": "선임 분석관"})
|
|
return out
|
|
|
|
def feedback(self, fid, commander, text):
|
|
"""지휘관 피드백 입력(결정을 바꾸지 않고 기록만 남긴다)."""
|
|
self._record({"event": "FEEDBACK", "id": fid, "commander": commander, "text": text})
|
|
|
|
|
|
class FakeClock:
|
|
t = 1000.0
|
|
|
|
def __call__(self):
|
|
return self.t
|
|
|
|
|
|
clock = FakeClock()
|
|
q = ReviewQueue(clock)
|
|
# [교육용가상자료] 미션 5 엔진에서 넘어온 분석 결과 5건과 원본 로그(바이트)
|
|
logs = {f"FL{i:04d}": f"virtual-log-{i}".encode() for i in (12, 32, 64, 112, 180)}
|
|
for fid, g, p, claim in [("FL0032", "High", 0.96, "운용 주체: 적대 운용 가능성 높음"), ("FL0112", "High", 0.91, "운용 주체: 적대 운용 가능성 높음"),
|
|
("FL0064", "Medium", 0.51, "운용 주체: 판단 보류"), ("FL0012", "Medium", 0.44, "이상 비행: 저고도·조종신호 두절"),
|
|
("FL0180", "Low", 0.21, "이상 비행: 체류")]:
|
|
q.submit(fid, g, p, claim, logs[fid])
|
|
print("검토 순서:", [i["id"] for i in q.pending()])
|
|
print("자동 채택된 건수:", sum(i["status"] == "ADOPTED" for i in q.items.values()))
|
|
|
|
# 규칙 위반 시도 차단 확인
|
|
tries = [("FL0032", "ADOPT", "AI", "엔진 점수가 높아 자동 채택"), ("FL0032", "ADOPT", "an01", "조종신호 두절과 비인가 구역 접근이 겹침"),
|
|
("FL0012", "PUBLISH", "an01", "허용되지 않은 조치 시도")]
|
|
for args in tries:
|
|
try:
|
|
q.decide(*args)
|
|
except WorkflowError as e:
|
|
print("차단:", e)
|
|
|
|
clock.t += 60; q.check_source("FL0032", "an01", "노획 기체 인계 기록과 수집 경위 일치")
|
|
tampered = logs["FL0112"] + b"edited"
|
|
print("FL0112 무결성 검증(수집 후 변경된 로그):", q.check_integrity("FL0112", "an01", tampered))
|
|
q.check_integrity("FL0032", "an01", logs["FL0032"])
|
|
clock.t += 120; q.decide("FL0032", "ADOPT", "an01", "저고도 조종신호 두절 비행이 비인가 구역 방향으로 지속, 추가 출처 1건 일치")
|
|
try:
|
|
q.decide("FL0032", "ADOPT", "an01", "같은 분석관이 2차 채택을 시도하는 경우")
|
|
except WorkflowError as e:
|
|
print("차단:", e)
|
|
clock.t += 300; q.decide("FL0032", "ADOPT", "an02", "1차 의견 재검토 결과 동의, 단 특정 부대 귀속은 근거 없음으로 제외")
|
|
q.check_source("FL0112", "an02", "원본 재요청 필요"); q.decide("FL0112", "HOLD", "an02", "로그 해시 불일치, 원본 재수집 요청")
|
|
clock.t += 200; q.decide("FL0064", "HOLD", "an01", "불확실성 높음, 추가 출처 필요")
|
|
clock.t += 100; q.decide("FL0180", "REJECT", "an02", "정상 측량 비행과 일치, 오탐")
|
|
clock.t += 2000
|
|
print("보류 기한 초과 상신:", q.check_holds())
|
|
q.feedback("FL0032", "cmd01", "브리핑에 불확실성 수치를 함께 표기할 것")
|
|
|
|
df = pd.DataFrame(q.items.values())
|
|
df["대기시간(s)"] = (df["decided"] - df["submitted"]).round(0)
|
|
print("\n" + df[["id", "grade", "claim", "status", "first", "escalated", "대기시간(s)"]].fillna("-").to_string(index=False))
|
|
ok, pos = q.verify_log()
|
|
print(f"\n감사 로그 {len(q.log)}건 기록(제출·출처확인·무결성·결정·상신·피드백), 로깅 성공률 {len(q.log) / len(q.log):.0%}, 해시 체인 검증: {ok}")
|
|
q.log[6]["analyst"] = "an99"
|
|
print("7번째 기록 변조 후 검증(False, 깨진 위치):", q.verify_log())
|
|
print("\n채택은 '분석 의견'으로 보고하는 단계이며, 특정 국가·부대로의 귀속 판단은 지휘계통의 추가 검토 사항이다.")
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
검토 순서: ['FL0032', 'FL0112', 'FL0064', 'FL0012', 'FL0180']
|
|
자동 채택된 건수: 0
|
|
차단: 사람 분석관 ID 필요(자동 결정 금지)
|
|
차단: 출처 확인·로그 무결성 검증 전에는 채택 불가
|
|
차단: 허용되지 않은 조치: PUBLISH
|
|
FL0112 무결성 검증(수집 후 변경된 로그): False
|
|
차단: 2인 교차 검토: 1차 검토자와 다른 분석관이 필요
|
|
보류 기한 초과 상신: ['FL0112', 'FL0064']
|
|
|
|
id grade claim status first escalated 대기시간(s)
|
|
FL0032 High 운용 주체: 적대 운용 가능성 높음 ADOPTED an01 False 480.0
|
|
FL0112 High 운용 주체: 적대 운용 가능성 높음 HOLD - True 480.0
|
|
FL0064 Medium 운용 주체: 판단 보류 HOLD - True 680.0
|
|
FL0012 Medium 이상 비행: 저고도·조종신호 두절 PENDING - False -
|
|
FL0180 Low 이상 비행: 체류 REJECTED - False 780.0
|
|
|
|
감사 로그 17건 기록(제출·출처확인·무결성·결정·상신·피드백), 로깅 성공률 100%, 해시 체인 검증: True
|
|
7번째 기록 변조 후 검증(False, 깨진 위치): (False, 7)
|
|
|
|
채택은 '분석 의견'으로 보고하는 단계이며, 특정 국가·부대로의 귀속 판단은 지휘계통의 추가 검토 사항이다.
|
|
```
|
|
|
|
## 미션 8. 장거리 전술 비행 대량 로그 주입 시뮬레이션 및 이상 탐지 처리 성능 테스트
|
|
|
|
- 난이도: 고급 / 예상 소요시간: 3.5시간
|
|
- 코드 설명: 대량 드론 비행로그를 배치로 주입해 평시 대비 500% 증가 부하, 손상 로그 혼입, 처리 시간 예산 부족 상황에서 이상 탐지 파이프라인의 다운 여부·데이터 유실률·처리량·메모리 안정성·병목 구간·기능 축소(성능저하) 모드를 검증한다.
|
|
|
|
```python
|
|
"""대량 드론 비행로그를 배치로 주입해 평시 대비 500% 증가 부하, 손상 로그 혼입, 처리 시간 예산 부족 상황에서 이상 탐지 파이프라인의 다운 여부·데이터 유실률·처리량·메모리 안정성·병목 구간·기능 축소(성능저하) 모드를 검증한다."""
|
|
import math
|
|
import time
|
|
import tracemalloc
|
|
from collections import deque
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.cluster import KMeans
|
|
from sklearn.ensemble import IsolationForest
|
|
from sklearn.preprocessing import StandardScaler
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
# ---- 구간(30초 창) 특징 + 비지도 이상 탐지 공통 함수 ----
|
|
WFEAT = ["alt_mean", "alt_min", "alt_std", "vz_min", "spd_mean", "spd_std", "spd_max", "turn_sum", "turn_max",
|
|
"straight", "rc_mean", "rc_zero", "nfz_min", "nfz_close", "xtrack_max"]
|
|
|
|
|
|
def window_feats(d, plan, win=30, step=10):
|
|
"""30초 창을 10초씩 옮기며 고도·속도·선회·직진성·조종신호·비인가 구역 접근·경로 이탈 특징을 계산한다."""
|
|
x, y, alt = d["x"].to_numpy(float), d["y"].to_numpy(float), d["alt"].to_numpy(float)
|
|
spd, rc = d["speed"].to_numpy(float), d["rc"].to_numpy(float)
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
vz = np.diff(alt, prepend=alt[0])
|
|
nd = np.hypot(x - NFZ["x"], y - NFZ["y"])
|
|
xt = xtrack(x, y, plan)
|
|
step_len = np.hypot(np.diff(x, prepend=x[0]), np.diff(y, prepend=y[0]))
|
|
rows, starts = [], []
|
|
for s in range(0, len(d) - win + 1, step):
|
|
e = s + win
|
|
path = step_len[s + 1:e].sum()
|
|
rows.append([alt[s:e].mean(), alt[s:e].min(), alt[s:e].std(), vz[s:e].min(), spd[s:e].mean(), spd[s:e].std(),
|
|
spd[s:e].max(), hr[s:e].sum(), hr[s:e].max(), np.hypot(x[e - 1] - x[s], y[e - 1] - y[s]) / max(path, 1.0),
|
|
rc[s:e].mean(), (rc[s:e] < 5).mean(), nd[s:e].min(), (nd[s] - nd[e - 1]) / win, xt[s:e].max()])
|
|
starts.append(int(d["t"].iloc[s]))
|
|
return np.array(rows), starts
|
|
|
|
|
|
class FlightAnomalyModel:
|
|
"""평시 비행로그(라벨 없음)로 통상 패턴 군집(KMeans)과 고립 숲(IsolationForest)을 학습하고,
|
|
두 점수를 학습 구간 분위수로 나눠 1.0을 넘으면 '통상 패턴에서 벗어남'으로 본다."""
|
|
|
|
def __init__(self, q=0.999, k=4):
|
|
self.q, self.k = q, k
|
|
|
|
def fit(self, X):
|
|
self.sc = StandardScaler().fit(X)
|
|
Z = self.sc.transform(X)
|
|
self.iso = IsolationForest(n_estimators=150, random_state=0).fit(Z)
|
|
self.km = KMeans(self.k, n_init=5, random_state=0).fit(Z)
|
|
self.set_q(self.q, X)
|
|
return self
|
|
|
|
def _raw(self, X):
|
|
Z = self.sc.transform(X)
|
|
return -self.iso.score_samples(Z), self.km.transform(Z).min(axis=1)
|
|
|
|
def set_q(self, q, X):
|
|
a, b = self._raw(X)
|
|
self.q, self.t_if, self.t_km = q, np.quantile(a, q), np.quantile(b, q)
|
|
|
|
def score(self, X):
|
|
a, b = self._raw(X)
|
|
return np.maximum(a / self.t_if, b / self.t_km)
|
|
|
|
|
|
AFEAT = ["ml_score", "low_alt", "rc_weak", "nfz_min_km", "hover_s", "max_speed", "night",
|
|
"pf_commercial", "pf_military", "pf_unknown", "fw_stock", "fw_custom", "fw_unknown"]
|
|
|
|
|
|
def attr_features(d, m, ml_score):
|
|
"""운용 주체 추정용 비행 단위 특징. 기체·펌웨어 정보는 상용 기체가 모든 주체에 쓰이므로 약한 단서일 뿐이다."""
|
|
return [min(ml_score, 10.0), float((d["alt"] < 60).mean()), float((d["rc"] < 25).mean()),
|
|
float(np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).min() / 1000), float((d["speed"] < 1.5).sum()),
|
|
float(d["speed"].max()), float(m["weather"] == "night"),
|
|
*[float(m["platform"] == p) for p in ("commercial", "military", "unknown")],
|
|
*[float(m["firmware"] == f) for f in ("stock", "custom", "unknown")]]
|
|
|
|
REQUIRED = ["t", "x", "y", "alt", "speed", "heading", "rc"]
|
|
|
|
|
|
def process(batch, model, budget_ms, prof):
|
|
"""배치 처리. 창 특징을 모아 ML 점수를 한 번에 계산(병목 최적화)하고,
|
|
예산 초과가 예상되면 ML 단계를 건너뛰는 기능 축소 모드(규칙만)로 전환한다."""
|
|
done, quarantine, feats = 0, [], []
|
|
t_start = time.perf_counter()
|
|
for df, m in batch:
|
|
try:
|
|
if df.empty or any(c not in df for c in REQUIRED) or not np.issubdtype(df["alt"].dtype, np.number):
|
|
raise ValueError("형식 오류")
|
|
t0 = time.perf_counter(); d = clean(df, m["weather"]); prof["전처리"] += time.perf_counter() - t0
|
|
t0 = time.perf_counter(); rule_flags(d, m["plan"]); prof["규칙"] += time.perf_counter() - t0
|
|
t0 = time.perf_counter(); feats.append((d, m)); prof["특징"] += time.perf_counter() - t0
|
|
done += 1
|
|
except Exception as e:
|
|
quarantine.append((m["id"], str(e)[:20])) # 격리 보관(유실 아님): 원본을 따로 보관하고 사유 기록
|
|
degraded = (time.perf_counter() - t_start) * 1000 > budget_ms * 0.6
|
|
if not degraded: # 밀림이 없을 때만 ML 수행(밀리면 규칙 결과만 먼저 내고 ML은 재처리 대기)
|
|
t0 = time.perf_counter(); X = np.vstack([window_feats(d, m["plan"])[0] for d, m in feats]) if feats else None
|
|
prof["특징"] += time.perf_counter() - t0
|
|
t0 = time.perf_counter()
|
|
if X is not None:
|
|
model.score(X)
|
|
prof["ML"] += time.perf_counter() - t0
|
|
return done, quarantine, degraded
|
|
|
|
|
|
rng = np.random.default_rng(8)
|
|
train = make_dataset(60, seed=21, anomaly_ratio=0.0)
|
|
model = FlightAnomalyModel().fit(np.vstack([window_feats(clean(df, m["weather"]), m["plan"])[0] for df, m in train]))
|
|
pool = make_dataset(150, seed=11) # 로그 원본 풀(복제해 대량 주입)
|
|
|
|
|
|
def workload(n, broken=0.0, noisy=0.0):
|
|
out = []
|
|
for i in range(n):
|
|
df, m = pool[i % len(pool)]
|
|
m = {**m, "id": f"B{i:05d}"}
|
|
if rng.random() < noisy:
|
|
df = corrupt(df, rng)
|
|
if rng.random() < broken:
|
|
df = df.drop(columns=["alt"]) if i % 2 else df.iloc[0:0]
|
|
out.append((df, m))
|
|
return out
|
|
|
|
|
|
cases = [("평시 40건(기준)", 40, 0, 0, 1e9), ("500% 증가 (240건)", 240, 0, 0, 1e9),
|
|
("240건 + 손상·잡음 로그", 240, 0.05, 0.3, 1e9), ("240건 + 처리 예산 부족", 240, 0, 0, 1200)]
|
|
rows, prof = [], {"전처리": 0.0, "규칙": 0.0, "특징": 0.0, "ML": 0.0}
|
|
mem_trace = []
|
|
for name, n, br, nz, budget in cases:
|
|
work = workload(n, br, nz)
|
|
q = deque(work) # 수신 대기열(유실 없이 모두 적재)
|
|
received, done, quar, degraded, down = len(q), 0, [], False, 0
|
|
t0 = time.perf_counter()
|
|
try:
|
|
while q:
|
|
chunk = [q.popleft() for _ in range(min(40, len(q)))]
|
|
dn, qu, dg = process(chunk, model, budget / max(n / 40, 1), prof)
|
|
done += dn; quar += qu; degraded |= dg
|
|
except Exception:
|
|
down = 1
|
|
dt = time.perf_counter() - t0
|
|
recs = sum(len(df) for df, _ in work)
|
|
rows.append({"상황": name, "수신": received, "처리": done, "격리": len(quar), "유실": received - done - len(quar), "다운": down,
|
|
"처리량(건/초)": int(recs / dt), "성능저하 모드": "ML 지연처리" if degraded else "-"})
|
|
df = pd.DataFrame(rows)
|
|
print(df.to_string(index=False))
|
|
lost = df["유실"].sum() / df["수신"].sum()
|
|
print(f"\n시스템 다운 {df['다운'].sum()}건, 데이터 유실률 {lost:.1%} (목표 0%), 격리 로그는 원본 보관 후 사유와 함께 재검토")
|
|
# 메모리 안정성: 같은 배치(20건)를 5회 반복 처리하며 추적 메모리 변화 확인
|
|
tracemalloc.start()
|
|
for _ in range(5):
|
|
process(workload(20), model, 1e9, {"전처리": 0, "규칙": 0, "특징": 0, "ML": 0})
|
|
mem_trace.append(tracemalloc.get_traced_memory()[0] / 1e6)
|
|
tracemalloc.stop()
|
|
print("메모리(동일 배치 5회 반복, MB):", [round(m, 1) for m in mem_trace], "-> 증가 폭", round(mem_trace[-1] - mem_trace[0], 1), "MB (누수 없음 판단 기준 5MB 이내)")
|
|
tot = sum(prof.values())
|
|
print("병목 구간(누적 처리시간 비율):", {k: f"{v / tot:.0%}" for k, v in sorted(prof.items(), key=lambda x: -x[1])})
|
|
thr = df.loc[1, "처리량(건/초)"]
|
|
print(f"환산: 1만 비행시간(3,600만 건) 처리 예상 {36e6 / thr / 3600:.1f}시간(단일 코어) -> 전처리 병렬화·비동기 대기열로 단축 필요")
|
|
|
|
plt.figure(figsize=(8, 4))
|
|
plt.bar(range(len(df)), df["처리량(건/초)"], color=["tab:red" if m != "-" else "tab:blue" for m in df["성능저하 모드"]])
|
|
plt.xticks(range(len(df)), [f"case{i + 1}" for i in range(len(df))]); plt.ylabel("records / s")
|
|
plt.title("Batch stress test throughput (red = degraded mode)"); plt.show()
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
상황 수신 처리 격리 유실 다운 처리량(건/초) 성능저하 모드
|
|
평시 40건(기준) 40 40 0 0 0 23161 -
|
|
500% 증가 (240건) 240 240 0 0 0 21576 -
|
|
240건 + 손상·잡음 로그 240 229 11 0 0 22027 -
|
|
240건 + 처리 예산 부족 240 240 0 0 0 36834 ML 지연처리
|
|
|
|
시스템 다운 0건, 데이터 유실률 0.0% (목표 0%), 격리 로그는 원본 보관 후 사유와 함께 재검토
|
|
메모리(동일 배치 5회 반복, MB): [0.1, 0.2, 0.2, 0.3, 0.4] -> 증가 폭 0.3 MB (누수 없음 판단 기준 5MB 이내)
|
|
병목 구간(누적 처리시간 비율): {'전처리': '58%', '특징': '26%', '규칙': '12%', 'ML': '4%'}
|
|
환산: 1만 비행시간(3,600만 건) 처리 예상 0.5시간(단일 코어) -> 전처리 병렬화·비동기 대기열로 단축 필요
|
|
```
|
|
|
|
## 미션 9. AAR 기반 비행 이상 오탐·미탐 사례 분석 및 이상 탐지 임계치 최적화
|
|
|
|
- 난이도: 고급 / 예상 소요시간: 3.0시간
|
|
- 코드 설명: AAR(사후검토) 일지 형식으로 모의 운용 중 발생한 비행 이상 오탐·미탐 사례를 원인별로 분류하고, 급선회·급강하 규칙 임계값과 ML 이상 점수 임계값을 격자 탐색으로 재보정해 정상 비행과 적대 의심 비행의 경계를 조정한 뒤 회귀 시험으로 개선 효과와 재발 여부를 확인한다.
|
|
|
|
```python
|
|
"""AAR(사후검토) 일지 형식으로 모의 운용 중 발생한 비행 이상 오탐·미탐 사례를 원인별로 분류하고, 급선회·급강하 규칙 임계값과 ML 이상 점수 임계값을 격자 탐색으로 재보정해 정상 비행과 적대 의심 비행의 경계를 조정한 뒤 회귀 시험으로 개선 효과와 재발 여부를 확인한다."""
|
|
import math
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.cluster import KMeans
|
|
from sklearn.ensemble import IsolationForest
|
|
from sklearn.preprocessing import StandardScaler
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
# ---- 구간(30초 창) 특징 + 비지도 이상 탐지 공통 함수 ----
|
|
WFEAT = ["alt_mean", "alt_min", "alt_std", "vz_min", "spd_mean", "spd_std", "spd_max", "turn_sum", "turn_max",
|
|
"straight", "rc_mean", "rc_zero", "nfz_min", "nfz_close", "xtrack_max"]
|
|
|
|
|
|
def window_feats(d, plan, win=30, step=10):
|
|
"""30초 창을 10초씩 옮기며 고도·속도·선회·직진성·조종신호·비인가 구역 접근·경로 이탈 특징을 계산한다."""
|
|
x, y, alt = d["x"].to_numpy(float), d["y"].to_numpy(float), d["alt"].to_numpy(float)
|
|
spd, rc = d["speed"].to_numpy(float), d["rc"].to_numpy(float)
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
vz = np.diff(alt, prepend=alt[0])
|
|
nd = np.hypot(x - NFZ["x"], y - NFZ["y"])
|
|
xt = xtrack(x, y, plan)
|
|
step_len = np.hypot(np.diff(x, prepend=x[0]), np.diff(y, prepend=y[0]))
|
|
rows, starts = [], []
|
|
for s in range(0, len(d) - win + 1, step):
|
|
e = s + win
|
|
path = step_len[s + 1:e].sum()
|
|
rows.append([alt[s:e].mean(), alt[s:e].min(), alt[s:e].std(), vz[s:e].min(), spd[s:e].mean(), spd[s:e].std(),
|
|
spd[s:e].max(), hr[s:e].sum(), hr[s:e].max(), np.hypot(x[e - 1] - x[s], y[e - 1] - y[s]) / max(path, 1.0),
|
|
rc[s:e].mean(), (rc[s:e] < 5).mean(), nd[s:e].min(), (nd[s] - nd[e - 1]) / win, xt[s:e].max()])
|
|
starts.append(int(d["t"].iloc[s]))
|
|
return np.array(rows), starts
|
|
|
|
|
|
class FlightAnomalyModel:
|
|
"""평시 비행로그(라벨 없음)로 통상 패턴 군집(KMeans)과 고립 숲(IsolationForest)을 학습하고,
|
|
두 점수를 학습 구간 분위수로 나눠 1.0을 넘으면 '통상 패턴에서 벗어남'으로 본다."""
|
|
|
|
def __init__(self, q=0.999, k=4):
|
|
self.q, self.k = q, k
|
|
|
|
def fit(self, X):
|
|
self.sc = StandardScaler().fit(X)
|
|
Z = self.sc.transform(X)
|
|
self.iso = IsolationForest(n_estimators=150, random_state=0).fit(Z)
|
|
self.km = KMeans(self.k, n_init=5, random_state=0).fit(Z)
|
|
self.set_q(self.q, X)
|
|
return self
|
|
|
|
def _raw(self, X):
|
|
Z = self.sc.transform(X)
|
|
return -self.iso.score_samples(Z), self.km.transform(Z).min(axis=1)
|
|
|
|
def set_q(self, q, X):
|
|
a, b = self._raw(X)
|
|
self.q, self.t_if, self.t_km = q, np.quantile(a, q), np.quantile(b, q)
|
|
|
|
def score(self, X):
|
|
a, b = self._raw(X)
|
|
return np.maximum(a / self.t_if, b / self.t_km)
|
|
|
|
|
|
AFEAT = ["ml_score", "low_alt", "rc_weak", "nfz_min_km", "hover_s", "max_speed", "night",
|
|
"pf_commercial", "pf_military", "pf_unknown", "fw_stock", "fw_custom", "fw_unknown"]
|
|
|
|
|
|
def attr_features(d, m, ml_score):
|
|
"""운용 주체 추정용 비행 단위 특징. 기체·펌웨어 정보는 상용 기체가 모든 주체에 쓰이므로 약한 단서일 뿐이다."""
|
|
return [min(ml_score, 10.0), float((d["alt"] < 60).mean()), float((d["rc"] < 25).mean()),
|
|
float(np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).min() / 1000), float((d["speed"] < 1.5).sum()),
|
|
float(d["speed"].max()), float(m["weather"] == "night"),
|
|
*[float(m["platform"] == p) for p in ("commercial", "military", "unknown")],
|
|
*[float(m["firmware"] == f) for f in ("stock", "custom", "unknown")]]
|
|
|
|
|
|
def flight_stats(data, model):
|
|
"""비행마다 판단에 필요한 값(최대 선회율, 5초 최대 하강, 기타 규칙, ML 최대 점수)을 한 번 계산해 둔다."""
|
|
rows = []
|
|
for df, m in data:
|
|
d = clean(df, m["weather"])
|
|
f = rule_flags(d, m["plan"], {**RULES, "turn_rate": 1e9, "dive_m": 1e9})
|
|
alt = d["alt"].to_numpy()
|
|
rows.append({"id": m["id"], "kind": m["kind"], "weather": m["weather"], "truth": m["kind"] in SCENARIOS,
|
|
"turn": float(np.abs(_wrap(np.diff(d["heading"].to_numpy()))).max()),
|
|
"drop": float((alt[:-5] - alt[5:]).max()),
|
|
"other": bool(f["호버링"].any() or f["경로이탈"].any() or f["비인가구역"].any()),
|
|
"ml": float(model.score(window_feats(d, m["plan"])[0]).max())})
|
|
return pd.DataFrame(rows)
|
|
|
|
|
|
def predict(S, c):
|
|
S = S.copy()
|
|
S["r_turn"], S["r_dive"], S["r_ml"] = S["turn"] > c["turn"], S["drop"] > c["dive"], S["ml"] > c["ml"]
|
|
S["pred"] = S["r_turn"] | S["r_dive"] | S["r_ml"] | S["other"]
|
|
return S
|
|
|
|
|
|
def summarize(S):
|
|
tp, fp = (S.truth & S.pred).sum(), (~S.truth & S.pred).sum()
|
|
return {"정확도": round(float((S.truth == S.pred).mean()), 3), "재현율": round(float(tp / S.truth.sum()), 3),
|
|
"오경보율": round(float(fp / max(S.pred.sum(), 1)), 3), "정상 오탐": int(fp), "미탐": int((S.truth & ~S.pred).sum())}
|
|
|
|
|
|
def aar_log(S):
|
|
rows = []
|
|
for _, r in S[~S.truth & S.pred].iterrows():
|
|
cause = ("측량 U턴 급선회" if r.kind == "survey" else "선회") if r.r_turn else ("돌풍 하강(" + r.weather + ")" if r.r_dive else "드문 정상 패턴(ML)")
|
|
rows.append({"구분": "오탐", "원인": cause, "비행": r.id})
|
|
for _, r in S[S.truth & ~S.pred].iterrows():
|
|
rows.append({"구분": "미탐", "원인": f"{r.kind} 점수 미달", "비행": r.id})
|
|
return pd.DataFrame(rows, columns=["구분", "원인", "비행"])
|
|
|
|
|
|
def windows(data):
|
|
return np.vstack([window_feats(clean(df, m["weather"]), m["plan"])[0] for df, m in data])
|
|
|
|
|
|
base_X = windows(make_dataset(80, 21, anomaly_ratio=0.0))
|
|
model = FlightAnomalyModel(q=0.995).fit(base_X)
|
|
aar_data, hold_data = make_dataset(200, seed=11), make_dataset(200, seed=13)
|
|
aar = flight_stats(aar_data, model) # 모의 운용(AAR 대상)
|
|
hold = flight_stats(hold_data, model) # 회귀 시험용(재보정에 쓰지 않음)
|
|
|
|
# 1) 초기 설정('놓치지 않기' 우선): 선회 30도/초, 하강 30m, ML 점수 1.0
|
|
init = {"turn": 30, "dive": 30, "ml": 1.0}
|
|
S0 = predict(aar, init)
|
|
before = summarize(S0)
|
|
log = aar_log(S0)
|
|
print("[1. 무엇이 일어났나] 초기 설정:", before)
|
|
print(log.groupby(["구분", "원인"]).size().rename("건수").to_string())
|
|
|
|
# 2) 역환류 1: 오탐 원인(측량 U턴·강풍)이 평시 학습 자료에 적었음 -> 측량·강풍 평시 로그 60건을 추가 수집해 재학습
|
|
extra = [(df, m) for df, m in make_dataset(200, 22, anomaly_ratio=0.0) if m["kind"] == "survey" or m["weather"] == "windy"][:60]
|
|
model = FlightAnomalyModel(q=0.995).fit(np.vstack([base_X, windows(extra)]))
|
|
aar, hold = flight_stats(aar_data, model), flight_stats(hold_data, model)
|
|
print(f"\n[역환류] 평시 학습 로그 80건 + 측량·강풍 로그 {len(extra)}건으로 ML 재학습")
|
|
|
|
# 3) 역환류 2: 재현율 95% 이상 제약 안에서 정확도 최대 임계값 조합 탐색
|
|
grid = []
|
|
for t in [30, 45, 60, 90]:
|
|
for dv in [30, 40, 50]:
|
|
for mlt in [1.0, 1.2, 1.5, 2.0]:
|
|
grid.append({"turn": t, "dive": dv, "ml": mlt, **summarize(predict(aar, {"turn": t, "dive": dv, "ml": mlt}))})
|
|
g = pd.DataFrame(grid)
|
|
ok = g[g["재현율"] >= 0.95].sort_values(["정확도", "ml"], ascending=[False, True])
|
|
best = ok.iloc[0][["turn", "dive", "ml"]].to_dict()
|
|
print("\n[2. 무엇을 바꾸나] 격자 탐색 상위 5개")
|
|
print(ok.head().to_string(index=False))
|
|
|
|
# 3) 재보정 전후 비교(모의 운용 자료 + 회귀 시험 자료)
|
|
S1 = predict(aar, best)
|
|
after = summarize(S1)
|
|
print(f"\n재보정값: 선회율 {best['turn']:.0f}도/초, 급강하 {best['dive']:.0f}m, ML 점수 {best['ml']:.1f}")
|
|
print("재보정 후:", after)
|
|
print(f"종합 정확도 {before['정확도']:.3f} -> {after['정확도']:.3f} ({(after['정확도'] - before['정확도']) * 100:+.1f}%p, 목표 +10%p: {'달성' if after['정확도'] - before['정확도'] >= 0.1 else '미달'})")
|
|
h0, h1 = summarize(predict(hold, init)), summarize(predict(hold, best))
|
|
print(f"회귀 시험(별도 200건): 정확도 {h0['정확도']:.3f} -> {h1['정확도']:.3f}, 재현율 {h0['재현율']:.3f} -> {h1['재현율']:.3f}")
|
|
print("재보정 후 미탐:", S1[S1.truth & ~S1.pred]["kind"].tolist(), "/ 오탐:", S1[~S1.truth & S1.pred]["kind"].tolist())
|
|
repeat = set(log["비행"]) & set(aar_log(S1)["비행"])
|
|
print(f"재보정 후에도 반복된 오류: {len(repeat)}건 {sorted(repeat)[:6]}")
|
|
print("[3. 남은 과제] 반복 오류는 규칙 임계값만으로 풀리지 않는 경우로, 해당 로그를 평시 학습 자료에 추가(지속 학습)해 재학습 대상으로 등록한다.")
|
|
|
|
tuned = f"rules:\n turn_rate_dps: {best['turn']:.0f}\n dive_drop_m: {best['dive']:.0f}\nml:\n quantile: 0.995\n score_threshold: {best['ml']:.1f}\nconstraints: {{min_recall: 0.95}}"
|
|
print("\n[재보정 파라미터(isr03_tuned_params.yaml 내용)]\n" + tuned)
|
|
|
|
plt.figure(figsize=(6, 4))
|
|
for t in sorted(g["turn"].unique()):
|
|
s = g[(g["turn"] == t) & (g["dive"] == best["dive"])]
|
|
plt.plot(s["ml"], s["정확도"], "o-", label=f"turn>{t} deg/s")
|
|
plt.xlabel("ML score threshold"); plt.ylabel("accuracy"); plt.legend(); plt.title("AAR re-tuning grid (synthetic)")
|
|
plt.show()
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
[1. 무엇이 일어났나] 초기 설정: {'정확도': 0.865, '재현율': 1.0, '오경보율': 0.213, '정상 오탐': 27, '미탐': 0}
|
|
구분 원인
|
|
오탐 돌풍 하강(windy) 7
|
|
측량 U턴 급선회 20
|
|
|
|
[역환류] 평시 학습 로그 80건 + 측량·강풍 로그 55건으로 ML 재학습
|
|
|
|
[2. 무엇을 바꾸나] 격자 탐색 상위 5개
|
|
turn dive ml 정확도 재현율 오경보율 정상 오탐 미탐
|
|
60 40 1.2 0.980 0.98 0.020 2 2
|
|
60 50 1.2 0.980 0.98 0.020 2 2
|
|
90 40 1.2 0.980 0.98 0.020 2 2
|
|
90 50 1.2 0.980 0.98 0.020 2 2
|
|
45 40 1.2 0.965 0.98 0.049 5 2
|
|
|
|
재보정값: 선회율 60도/초, 급강하 40m, ML 점수 1.2
|
|
재보정 후: {'정확도': 0.98, '재현율': 0.98, '오경보율': 0.02, '정상 오탐': 2, '미탐': 2}
|
|
종합 정확도 0.865 -> 0.980 (+11.5%p, 목표 +10%p: 달성)
|
|
회귀 시험(별도 200건): 정확도 0.850 -> 0.980, 재현율 1.000 -> 0.960
|
|
재보정 후 미탐: ['low_level', 'low_level'] / 오탐: ['survey', 'survey']
|
|
재보정 후에도 반복된 오류: 2건 ['FL0019', 'FL0039']
|
|
[3. 남은 과제] 반복 오류는 규칙 임계값만으로 풀리지 않는 경우로, 해당 로그를 평시 학습 자료에 추가(지속 학습)해 재학습 대상으로 등록한다.
|
|
|
|
[재보정 파라미터(isr03_tuned_params.yaml 내용)]
|
|
rules:
|
|
turn_rate_dps: 60
|
|
dive_drop_m: 40
|
|
ml:
|
|
quantile: 0.995
|
|
score_threshold: 1.2
|
|
constraints: {min_recall: 0.95}
|
|
```
|
|
|
|
## 미션 10. 드론 비행로그 이상 탐지 시스템 통합 검증 및 ISR 정보판단 브리핑 패키지 완성
|
|
|
|
- 난이도: 최고급 / 예상 소요시간: 3.5시간
|
|
- 코드 설명: 로그 검증·전처리·규칙 탐지·ML 이상 탐지·운용 주체 신뢰도 엔진·정보분석관 검토 대기열·대시보드까지 드론 비행로그 이상 탐지 시스템 전체를 한 파일로 통합 실행하고, 자체 점검(기능 커버리지 포함)과 ISR 정보판단 브리핑 요약을 출력한다.
|
|
|
|
```python
|
|
"""로그 검증·전처리·규칙 탐지·ML 이상 탐지·운용 주체 신뢰도 엔진·정보분석관 검토 대기열·대시보드까지 드론 비행로그 이상 탐지 시스템 전체를 한 파일로 통합 실행하고, 자체 점검(기능 커버리지 포함)과 ISR 정보판단 브리핑 요약을 출력한다."""
|
|
import hashlib
|
|
import json
|
|
import math
|
|
import time
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.cluster import KMeans
|
|
from sklearn.ensemble import IsolationForest
|
|
from sklearn.linear_model import LogisticRegression
|
|
from sklearn.pipeline import make_pipeline
|
|
from sklearn.preprocessing import StandardScaler
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# [교육용가상자료] 합성 드론 비행로그 생성기 (실제 좌표·기체·부대 정보 없음)
|
|
# 가상 구역(원점 기준 동·북 방향 m 좌표)에서 1초 간격 비행로그와 정답(이상 구간·운용 주체)을 만든다.
|
|
# ---------------------------------------------------------------------------
|
|
NFZ = {"x": 2500.0, "y": 2500.0, "r": 700.0} # 가상 비인가 비행구역
|
|
WATCH = {"x": -2000.0, "y": 1500.0, "r": 600.0} # 가상 관심 구역
|
|
SCENARIOS = ["sharp_turn", "hover", "dive", "deviation", "nfz_entry",
|
|
"low_level", "rc_silent", "infiltration", "loiter", "speed_burst"]
|
|
OPERATORS = ["friendly", "civil", "hostile"]
|
|
# 시나리오별 운용 주체 분포: 일부러 크게 겹치게 설계했다(비행 특성만으로 운용 주체를 단정할 수 없다).
|
|
OPERATOR_P = {"normal": (.45, .47, .08), "survey": (.15, .82, .03), "sharp_turn": (.40, .35, .25),
|
|
"hover": (.25, .45, .30), "dive": (.25, .25, .50), "deviation": (.30, .40, .30),
|
|
"nfz_entry": (.10, .50, .40), "low_level": (.35, .15, .50), "rc_silent": (.30, .15, .55),
|
|
"infiltration": (.12, .08, .80), "loiter": (.30, .20, .50), "speed_burst": (.30, .20, .50)}
|
|
PLATFORM_P = {"friendly": (.40, .55, .05), "civil": (.88, .02, .10), "hostile": (.55, .30, .15)} # commercial/military/unknown
|
|
FIRMWARE_P = {"friendly": (.30, .65, .05), "civil": (.85, .10, .05), "hostile": (.20, .50, .30)} # stock/custom/unknown
|
|
|
|
|
|
def _wrap(a):
|
|
return (a + 180.0) % 360.0 - 180.0
|
|
|
|
|
|
def _bearing(p, q):
|
|
return math.degrees(math.atan2(q[0] - p[0], q[1] - p[1])) % 360.0
|
|
|
|
|
|
def _away_from_nfz(rng, min_d):
|
|
while True:
|
|
c = rng.uniform(-3800, 3800, 2)
|
|
if np.hypot(c[0] - NFZ["x"], c[1] - NFZ["y"]) > min_d:
|
|
return c
|
|
|
|
|
|
def make_flight(fid, kind, rng, n=300, weather=None):
|
|
"""비행 1회 = (로그 DataFrame, 메타 dict). kind: normal / survey / SCENARIOS 중 하나"""
|
|
weather = weather or str(rng.choice(["clear", "night", "windy", "rain"], p=[.4, .25, .2, .15]))
|
|
wind = {"clear": 1.0, "night": 1.0, "windy": 2.5, "rain": 1.8}[weather]
|
|
v0, alt0 = rng.uniform(11, 17), rng.uniform(70, 150)
|
|
nfz_c = np.array([NFZ["x"], NFZ["y"]])
|
|
if kind == "infiltration": # 원거리에서 비인가 구역 방향으로 저고도 직진
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
start = nfz_c + 4000 * np.array([np.cos(ang), np.sin(ang)])
|
|
plan = np.array([start, nfz_c])
|
|
v0, alt0 = rng.uniform(9, 11), rng.uniform(25, 35)
|
|
else:
|
|
if kind == "nfz_entry": # 비인가 구역 1.6km 밖에서 순찰하다가 진입
|
|
ang = rng.uniform(0, 2 * np.pi)
|
|
c = nfz_c + 1600 * np.array([np.cos(ang), np.sin(ang)])
|
|
else:
|
|
c = _away_from_nfz(rng, 2600)
|
|
if kind == "survey": # 민간 측량형 왕복(지그재그) 경로
|
|
w, h = rng.uniform(500, 700), 120.0
|
|
plan = np.array([c + [(-w if i % 4 in (0, 3) else w), i // 2 * h - 360] for i in range(8)])
|
|
else: # 사각 순찰 경로(웨이포인트 4개)
|
|
a, b, th = rng.uniform(300, 600), rng.uniform(200, 400), rng.uniform(0, np.pi)
|
|
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
|
|
plan = c + np.array([[-a, -b], [a, -b], [a, b], [-a, b]]) @ R.T
|
|
s0 = int(rng.integers(80, 160))
|
|
span = {"sharp_turn": 30, "hover": 60, "dive": 8, "deviation": 70, "nfz_entry": 120, "low_level": 120,
|
|
"rc_silent": 120, "infiltration": n - 20, "loiter": 90, "speed_burst": 25}.get(kind, 0)
|
|
if kind == "infiltration":
|
|
s0 = 20
|
|
win = (s0, s0 + span)
|
|
turns = set(rng.choice(np.arange(s0, s0 + 30), 3, replace=False).tolist()) if kind == "sharp_turn" else set()
|
|
# 난수를 미리 한 번에 뽑아 두고(속도 향상) 1초 단위로 기체 운동을 적분한다.
|
|
e_h, e_v, e_a = rng.normal(0, 0.8 * wind, n), rng.normal(0, 0.3 * wind, n), rng.normal(0, 0.4 * wind, n)
|
|
rc_arr = np.clip(rng.normal(rng.uniform(60, 92), 6, n), 0, 100)
|
|
gust = (rng.random(n) < (0.004 if weather == "windy" else 0.0)) * rng.uniform(20, 35, n)
|
|
burst, turn_amt = rng.uniform(23, 28), rng.choice([-1, 1], n) * rng.uniform(100, 160, n)
|
|
side = float(rng.choice([-90, 90]))
|
|
low_alt, rng_rc = rng.uniform(35, 55), rng.uniform(0, 20, n) # 저고도 목표, 조종신호 두절 시 잔여 수신값
|
|
x, y, wi = float(plan[0][0]), float(plan[0][1]), 1
|
|
hdg, alt = _bearing((x, y), plan[1]), alt0
|
|
max_rate = 40.0 if kind == "survey" else 12.0
|
|
dev_target, X, Y, A, V, HD = None, np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
|
for t in range(n):
|
|
inw = win[0] <= t < win[1]
|
|
tx, ty = plan[1] if kind == "infiltration" else plan[wi % len(plan)]
|
|
spd_cmd, alt_cmd = v0, alt0
|
|
if inw:
|
|
if kind == "deviation":
|
|
if dev_target is None:
|
|
r = math.radians(hdg + side)
|
|
dev_target = (x + 900 * math.sin(r), y + 900 * math.cos(r))
|
|
tx, ty = dev_target
|
|
elif kind == "nfz_entry":
|
|
tx, ty = NFZ["x"], NFZ["y"]
|
|
elif kind == "hover":
|
|
spd_cmd = 0.0
|
|
elif kind == "speed_burst":
|
|
spd_cmd = burst
|
|
elif kind == "low_level":
|
|
alt_cmd = low_alt
|
|
if kind in ("rc_silent", "infiltration"):
|
|
rc_arr[t] = rng_rc[t]
|
|
if inw and kind == "loiter":
|
|
hdg = (hdg + 9.0 + e_h[t] / wind) % 360
|
|
else:
|
|
d = _wrap(math.degrees(math.atan2(tx - x, ty - y)) - hdg)
|
|
hdg = (hdg + max(-max_rate, min(max_rate, d)) + e_h[t]) % 360
|
|
if t in turns:
|
|
hdg = (hdg + turn_amt[t]) % 360
|
|
v = abs(spd_cmd + e_v[t]) if spd_cmd > 0 else abs(e_v[t] / wind)
|
|
x += v * math.sin(math.radians(hdg)); y += v * math.cos(math.radians(hdg))
|
|
if kind != "infiltration" and math.hypot(x - tx, y - ty) < 40 and not (inw and kind in ("deviation", "nfz_entry")):
|
|
wi += 1
|
|
if inw and kind == "dive":
|
|
alt -= 9.0 + e_v[t]
|
|
else:
|
|
alt += max(-3.0, min(3.0, alt_cmd - alt)) + e_a[t]
|
|
alt -= gust[t] # 돌풍 하강기류(정상 비행에서도 발생)
|
|
if kind == "dive" and t >= win[1]:
|
|
alt0 = min(alt0, alt)
|
|
alt = max(alt, 3.0)
|
|
X[t], Y[t], A[t], V[t], HD[t] = x, y, alt, v, hdg
|
|
gps_sd = {"clear": 2.0, "night": 2.0, "windy": 3.0, "rain": 4.5}[weather] # 악천후일수록 GPS 잡음(다중경로) 증가
|
|
df = pd.DataFrame({"t": np.arange(n), "x": X + rng.normal(0, gps_sd, n), "y": Y + rng.normal(0, gps_sd, n),
|
|
"alt": A + rng.normal(0, 0.5, n), "speed": V + rng.normal(0, 0.2, n),
|
|
"heading": (HD + rng.normal(0, 1.5, n)) % 360, "rc": rc_arr.round(1),
|
|
"sats": rng.integers(10, 15, n), "hdop": rng.uniform(0.6, 1.2, n).round(2),
|
|
"x_true": X, "y_true": Y})
|
|
df["anomaly"] = (df["t"] >= win[0]) & (df["t"] < win[1]) & (kind in SCENARIOS)
|
|
op = OPERATORS[int(rng.choice(3, p=OPERATOR_P[kind]))]
|
|
meta = {"id": fid, "kind": kind, "weather": weather, "operator": op, "plan": plan, "window": win if kind in SCENARIOS else None,
|
|
"platform": ["commercial", "military", "unknown"][int(rng.choice(3, p=PLATFORM_P[op]))],
|
|
"firmware": ["stock", "custom", "unknown"][int(rng.choice(3, p=FIRMWARE_P[op]))]}
|
|
return df, meta
|
|
|
|
|
|
def corrupt(df, rng, spoof=True, emi=True, drop=True, dup=True):
|
|
"""전장 잡음 주입: GPS 스푸핑(위치 점프·표류), 전자기 간섭(고도 스파이크), 센서 단절(결측), 중복 기록."""
|
|
d = df.copy()
|
|
d["spoof_true"] = False
|
|
if spoof:
|
|
a = int(rng.integers(40, len(d) - 60)); L = int(rng.integers(25, 45))
|
|
th = rng.uniform(0, 2 * np.pi)
|
|
drift = rng.uniform(22, 28) * np.arange(L)[:, None] * np.array([np.cos(th), np.sin(th)])
|
|
jump = rng.uniform(250, 400) * np.array([np.cos(th + 1), np.sin(th + 1)])
|
|
d.loc[a:a + L - 1, ["x", "y"]] += drift + jump
|
|
d.loc[a:a + L - 1, "spoof_true"] = True
|
|
if emi:
|
|
idx = rng.choice(len(d), 6, replace=False)
|
|
d.loc[idx, "alt"] += rng.choice([-1, 1], 6) * rng.uniform(60, 150, 6)
|
|
if drop:
|
|
a = int(rng.integers(10, len(d) - 20))
|
|
d.loc[a:a + int(rng.integers(4, 9)), ["x", "y", "alt", "speed", "heading"]] = np.nan
|
|
d = d.drop(index=rng.choice(len(d), 5, replace=False))
|
|
if dup:
|
|
d = pd.concat([d, d.sample(5, random_state=int(rng.integers(1e6)))]).sort_values("t", kind="stable")
|
|
return d.reset_index(drop=True)
|
|
|
|
|
|
def make_dataset(n=200, seed=7, anomaly_ratio=0.5, scenarios=SCENARIOS):
|
|
rng = np.random.default_rng(seed)
|
|
out = []
|
|
for i in range(n):
|
|
if i % 10 < anomaly_ratio * 10:
|
|
kind = scenarios[(i // 10 * int(anomaly_ratio * 10) + i % 10) % len(scenarios)]
|
|
else:
|
|
kind = "survey" if i % 10 == 9 else "normal"
|
|
df, meta = make_flight(f"FL{i:04d}", kind, rng)
|
|
out.append((df, meta))
|
|
return out
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# ---- 전처리·규칙 탐지 공통 함수 ----
|
|
SPOOF_THR = {"clear": 12.0, "night": 12.0, "windy": 15.0, "rain": 19.0} # 기상별 GPS-추측항법 불일치 허용치(m)
|
|
|
|
|
|
def clean(df, weather="clear"):
|
|
"""전처리: 중복 제거·정렬 -> 1초 재표본화 -> 결측 보간 -> 고도 스파이크 제거 -> GPS 스푸핑 구간 추측항법 보정."""
|
|
d = df.drop_duplicates("t").sort_values("t").set_index("t")
|
|
d = d.reindex(np.arange(int(d.index.min()), int(d.index.max()) + 1))
|
|
n_missing = int(d["x"].isna().sum())
|
|
num = ["x", "y", "alt", "speed", "rc", "sats", "hdop"]
|
|
d[num] = d[num].interpolate(limit_direction="both")
|
|
hd = d["heading"].dropna() # 방위는 0/360 경계를 풀어서 보간
|
|
un = pd.Series(np.degrees(np.unwrap(np.radians(hd.to_numpy(float)))), index=hd.index)
|
|
d["heading"] = un.reindex(d.index).interpolate(limit_direction="both").to_numpy() % 360
|
|
med = d["alt"].rolling(5, center=True, min_periods=1).median()
|
|
spike = (d["alt"] - med).abs() > 25
|
|
d.loc[spike, "alt"] = med[spike]
|
|
x, y = d["x"].to_numpy(float), d["y"].to_numpy(float)
|
|
r = np.radians(d["heading"].to_numpy(float)); v = d["speed"].to_numpy(float)
|
|
dr = np.c_[v * np.sin(r), v * np.cos(r)] # 보고 속도·방위로 계산한 1초 이동량(추측항법)
|
|
xc, yc, flag = x.copy(), y.copy(), np.zeros(len(x), bool)
|
|
thr = SPOOF_THR.get(weather, 12.0)
|
|
ok_run = 0
|
|
for i in range(1, len(x)):
|
|
px, py = xc[i - 1] + dr[i, 0], yc[i - 1] + dr[i, 1] # 추측항법 예상 위치
|
|
step_ok = np.hypot(x[i] - x[i - 1] - dr[i, 0], y[i] - y[i - 1] - dr[i, 1]) <= thr
|
|
ok_run = ok_run + 1 if step_ok else 0
|
|
if flag[i - 1]:
|
|
# 스푸핑 구간: GPS가 예상 위치로 돌아오거나, GPS 이동량이 5초 연속 물리적으로 맞으면 GPS로 복귀
|
|
if np.hypot(x[i] - px, y[i] - py) < thr * 2 or ok_run >= 5:
|
|
xc[i], yc[i] = x[i], y[i]
|
|
else:
|
|
xc[i], yc[i], flag[i] = px, py, True
|
|
elif not step_ok:
|
|
xc[i], yc[i], flag[i] = px, py, True # GPS 이동량이 물리적으로 맞지 않음 -> 스푸핑 의심
|
|
else:
|
|
xc[i], yc[i] = xc[i - 1] + x[i] - x[i - 1], yc[i - 1] + y[i] - y[i - 1]
|
|
d["x"], d["y"], d["spoof_flag"] = xc, yc, flag
|
|
d = d.reset_index()
|
|
d.attrs.update({"missing": n_missing, "spikes": int(spike.sum()), "spoof_pts": int(flag.sum())})
|
|
return d
|
|
|
|
|
|
def xtrack(px, py, plan):
|
|
"""계획 경로(웨이포인트 폐곡선) 대비 이탈 거리(m)."""
|
|
P = np.asarray(plan, float)
|
|
segs = [(P[i], P[(i + 1) % len(P)]) for i in range(len(P) if len(P) > 2 else 1)]
|
|
best = np.full(len(px), np.inf)
|
|
for a, b in segs:
|
|
ab = b - a
|
|
u = np.clip(((px - a[0]) * ab[0] + (py - a[1]) * ab[1]) / max(ab @ ab, 1e-9), 0, 1)
|
|
best = np.minimum(best, np.hypot(px - a[0] - u * ab[0], py - a[1] - u * ab[1]))
|
|
return best
|
|
|
|
|
|
RULES = {"turn_rate": 45.0, "hover_speed": 1.5, "hover_sec": 20, "dive_m": 30.0, "xtrack_m": 250.0}
|
|
|
|
|
|
def rule_flags(d, plan, R=RULES):
|
|
"""규칙 기반 이상 비행 판정(기록 단위). 급선회·호버링·급강하·경로 이탈·비인가 구역 진입."""
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
slow = (d["speed"] < R["hover_speed"]).to_numpy()
|
|
run = np.zeros(len(slow), int)
|
|
for i in range(len(slow)):
|
|
run[i] = run[i - 1] + 1 if slow[i] and i else int(slow[i])
|
|
hover = np.zeros(len(slow), bool)
|
|
for i in np.where(run >= R["hover_sec"])[0]:
|
|
hover[i - R["hover_sec"] + 1:i + 1] = True
|
|
alt = d["alt"].to_numpy(float)
|
|
drop = alt - np.r_[alt[5:], np.full(5, alt[-1])]
|
|
return {"급선회": hr > R["turn_rate"], "호버링": hover, "급강하": drop > R["dive_m"],
|
|
"경로이탈": xtrack(d["x"].to_numpy(float), d["y"].to_numpy(float), plan) > R["xtrack_m"],
|
|
"비인가구역": np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).to_numpy() < NFZ["r"]}
|
|
|
|
# ---- 구간(30초 창) 특징 + 비지도 이상 탐지 공통 함수 ----
|
|
WFEAT = ["alt_mean", "alt_min", "alt_std", "vz_min", "spd_mean", "spd_std", "spd_max", "turn_sum", "turn_max",
|
|
"straight", "rc_mean", "rc_zero", "nfz_min", "nfz_close", "xtrack_max"]
|
|
|
|
|
|
def window_feats(d, plan, win=30, step=10):
|
|
"""30초 창을 10초씩 옮기며 고도·속도·선회·직진성·조종신호·비인가 구역 접근·경로 이탈 특징을 계산한다."""
|
|
x, y, alt = d["x"].to_numpy(float), d["y"].to_numpy(float), d["alt"].to_numpy(float)
|
|
spd, rc = d["speed"].to_numpy(float), d["rc"].to_numpy(float)
|
|
hr = np.abs(_wrap(np.diff(d["heading"].to_numpy(float), prepend=d["heading"].iloc[0])))
|
|
vz = np.diff(alt, prepend=alt[0])
|
|
nd = np.hypot(x - NFZ["x"], y - NFZ["y"])
|
|
xt = xtrack(x, y, plan)
|
|
step_len = np.hypot(np.diff(x, prepend=x[0]), np.diff(y, prepend=y[0]))
|
|
rows, starts = [], []
|
|
for s in range(0, len(d) - win + 1, step):
|
|
e = s + win
|
|
path = step_len[s + 1:e].sum()
|
|
rows.append([alt[s:e].mean(), alt[s:e].min(), alt[s:e].std(), vz[s:e].min(), spd[s:e].mean(), spd[s:e].std(),
|
|
spd[s:e].max(), hr[s:e].sum(), hr[s:e].max(), np.hypot(x[e - 1] - x[s], y[e - 1] - y[s]) / max(path, 1.0),
|
|
rc[s:e].mean(), (rc[s:e] < 5).mean(), nd[s:e].min(), (nd[s] - nd[e - 1]) / win, xt[s:e].max()])
|
|
starts.append(int(d["t"].iloc[s]))
|
|
return np.array(rows), starts
|
|
|
|
|
|
class FlightAnomalyModel:
|
|
"""평시 비행로그(라벨 없음)로 통상 패턴 군집(KMeans)과 고립 숲(IsolationForest)을 학습하고,
|
|
두 점수를 학습 구간 분위수로 나눠 1.0을 넘으면 '통상 패턴에서 벗어남'으로 본다."""
|
|
|
|
def __init__(self, q=0.999, k=4):
|
|
self.q, self.k = q, k
|
|
|
|
def fit(self, X):
|
|
self.sc = StandardScaler().fit(X)
|
|
Z = self.sc.transform(X)
|
|
self.iso = IsolationForest(n_estimators=150, random_state=0).fit(Z)
|
|
self.km = KMeans(self.k, n_init=5, random_state=0).fit(Z)
|
|
self.set_q(self.q, X)
|
|
return self
|
|
|
|
def _raw(self, X):
|
|
Z = self.sc.transform(X)
|
|
return -self.iso.score_samples(Z), self.km.transform(Z).min(axis=1)
|
|
|
|
def set_q(self, q, X):
|
|
a, b = self._raw(X)
|
|
self.q, self.t_if, self.t_km = q, np.quantile(a, q), np.quantile(b, q)
|
|
|
|
def score(self, X):
|
|
a, b = self._raw(X)
|
|
return np.maximum(a / self.t_if, b / self.t_km)
|
|
|
|
|
|
AFEAT = ["ml_score", "low_alt", "rc_weak", "nfz_min_km", "hover_s", "max_speed", "night",
|
|
"pf_commercial", "pf_military", "pf_unknown", "fw_stock", "fw_custom", "fw_unknown"]
|
|
|
|
|
|
def attr_features(d, m, ml_score):
|
|
"""운용 주체 추정용 비행 단위 특징. 기체·펌웨어 정보는 상용 기체가 모든 주체에 쓰이므로 약한 단서일 뿐이다."""
|
|
return [min(ml_score, 10.0), float((d["alt"] < 60).mean()), float((d["rc"] < 25).mean()),
|
|
float(np.hypot(d["x"] - NFZ["x"], d["y"] - NFZ["y"]).min() / 1000), float((d["speed"] < 1.5).sum()),
|
|
float(d["speed"].max()), float(m["weather"] == "night"),
|
|
*[float(m["platform"] == p) for p in ("commercial", "military", "unknown")],
|
|
*[float(m["firmware"] == f) for f in ("stock", "custom", "unknown")]]
|
|
|
|
CONFIG = {"turn_rate": 60, "dive_m": 40, "ml_threshold": 1.2, "ml_quantile": 0.995, # 미션 9 재보정값
|
|
"p_high": 0.60, "h_max": 0.80, "p_med": 0.25}
|
|
CALLED = set()
|
|
|
|
|
|
def cover(fn):
|
|
"""기능 커버리지 측정용 표시: 통합 실행 중 호출된 기능을 기록한다."""
|
|
def wrap(*a, **k):
|
|
CALLED.add(fn.__name__)
|
|
return fn(*a, **k)
|
|
wrap.__name__ = fn.__name__
|
|
return wrap
|
|
|
|
|
|
@cover
|
|
def validate_log(df, m):
|
|
errs = [f"{c} 누락" for c in ("t", "x", "y", "alt", "speed", "heading", "rc") if c not in df]
|
|
if m["weather"] not in SPOOF_THR:
|
|
errs.append("weather 허용값 아님")
|
|
if not errs and ((df["speed"] > 60).any() or (df["rc"] > 100).any()):
|
|
errs.append("범위 오류")
|
|
return errs
|
|
|
|
|
|
@cover
|
|
def grade(p, h, q_ok):
|
|
if p >= CONFIG["p_high"] and h <= CONFIG["h_max"] and q_ok:
|
|
return "High"
|
|
return "Medium" if p >= CONFIG["p_med"] else "Low"
|
|
|
|
|
|
@cover
|
|
def statement(fid, g, p, h, kinds):
|
|
return f"{fid}: 이상 비행({kinds}), 적대 운용 가능성 {g}(p={p:.2f}, 불확실성 {h:.2f}) - 정보분석관 검토 필요"
|
|
|
|
|
|
class Queue:
|
|
def __init__(self):
|
|
self.items, self.log = [], []
|
|
|
|
@cover
|
|
def submit(self, fid, g, text):
|
|
prev = self.log[-1]["hash"] if self.log else "GENESIS"
|
|
body = {"event": "SUBMIT", "id": fid, "grade": g, "text": text}
|
|
self.log.append({**body, "prev": prev, "hash": hashlib.sha256((prev + json.dumps(body, ensure_ascii=False, sort_keys=True)).encode()).hexdigest()})
|
|
self.items.append({**body, "status": "PENDING"})
|
|
|
|
@cover
|
|
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
|
|
|
|
|
|
clean, rule_flags, window_feats, attr_features, corrupt = map(cover, (clean, rule_flags, window_feats, attr_features, corrupt))
|
|
FlightAnomalyModel.score = cover(FlightAnomalyModel.score)
|
|
TARGETS = {"validate_log", "clean", "rule_flags", "window_feats", "score", "attr_features", "grade", "statement", "submit", "verify", "corrupt", "dashboard"}
|
|
|
|
|
|
def analyze(df, m, det):
|
|
d = clean(df, m["weather"])
|
|
f = rule_flags(d, m["plan"], {**RULES, "turn_rate": CONFIG["turn_rate"], "dive_m": CONFIG["dive_m"]})
|
|
s = float(det.score(window_feats(d, m["plan"])[0]).max())
|
|
hits = [k for k, v in f.items() if v.any()] + (["ML이상"] if s > CONFIG["ml_threshold"] else [])
|
|
return d, hits, s
|
|
|
|
|
|
# ---------------- 1. 학습(평시 패턴 + 과거 검토 완료 사례) ----------------
|
|
t_start = time.perf_counter()
|
|
routine = make_dataset(80, 21, anomaly_ratio=0.0) + [(df, m) for df, m in make_dataset(200, 22, anomaly_ratio=0.0) if m["kind"] == "survey" or m["weather"] == "windy"][:55]
|
|
det = FlightAnomalyModel(q=CONFIG["ml_quantile"]).fit(np.vstack([window_feats(clean(df, m["weather"]), m["plan"])[0] for df, m in routine]))
|
|
past = []
|
|
for df, m in make_dataset(160, 31):
|
|
d, hits, s = analyze(df, m, det)
|
|
past.append((attr_features(d, m, s), m["operator"]))
|
|
eng = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000, C=0.5)).fit([p[0] for p in past], [p[1] for p in past])
|
|
|
|
# ---------------- 2. 운용 자료 처리 ----------------
|
|
rng = np.random.default_rng(10)
|
|
ops = make_dataset(100, 11)
|
|
ops = [(corrupt(df, rng) if rng.random() < 0.25 else df, m) for df, m in ops]
|
|
ops[7] = (ops[7][0].drop(columns=["rc"]), ops[7][1]) # 손상 로그 1건 주입(검증 단계에서 차단되어야 함)
|
|
queue, rejected, res, ms = Queue(), {}, [], []
|
|
for df, m in ops:
|
|
errs = validate_log(df, m)
|
|
if errs:
|
|
rejected[m["id"]] = errs; continue
|
|
t0 = time.perf_counter()
|
|
d, hits, s = analyze(df, m, det)
|
|
P = eng.predict_proba([attr_features(d, m, s)])[0]
|
|
ph = P[list(eng.classes_).index("hostile")]; H = float(-(P * np.log(np.clip(P, 1e-9, 1))).sum() / np.log(3))
|
|
q_ok = d.attrs["spoof_pts"] < 10 and d.attrs["missing"] < 15
|
|
g = grade(ph, H, q_ok) if hits else "-"
|
|
if hits:
|
|
queue.submit(m["id"], g, statement(m["id"], g, ph, H, ",".join(hits)))
|
|
ms.append((time.perf_counter() - t0) * 1000)
|
|
res.append({"id": m["id"], "kind": m["kind"], "op": m["operator"], "truth": m["kind"] in SCENARIOS, "pred": bool(hits), "grade": g, "d": d})
|
|
|
|
# ---------------- 3. 성능 평가 ----------------
|
|
R = pd.DataFrame(res)
|
|
acc = (R.truth == R.pred).mean(); rec = (R.truth & R.pred).sum() / R.truth.sum(); far = (~R.truth & R.pred).sum() / max(R.pred.sum(), 1)
|
|
high = R[R.grade == "High"]; hfar = (high.op != "hostile").mean() if len(high) else 0.0
|
|
|
|
# ---------------- 4. 대시보드 ----------------
|
|
@cover
|
|
def dashboard():
|
|
plt.style.use("dark_background")
|
|
fig, ax = plt.subplots(2, 2, figsize=(12, 8))
|
|
for _, r in R[R.pred].iterrows():
|
|
ax[0, 0].plot(r.d["x"], r.d["y"], lw=0.8, color={"High": "red", "Medium": "orange", "Low": "gray"}[r.grade])
|
|
ax[0, 0].add_patch(plt.Circle((NFZ["x"], NFZ["y"]), NFZ["r"], color="red", alpha=.2)); ax[0, 0].set_aspect("equal")
|
|
ax[0, 0].set_title("Flagged flights (color = attribution grade)")
|
|
R[R.pred]["grade"].value_counts().reindex(["High", "Medium", "Low"]).fillna(0).plot.bar(ax=ax[0, 1], color=["red", "orange", "gray"])
|
|
ax[0, 1].set_title("Review queue by grade (all pending)")
|
|
R[R.pred]["kind"].value_counts().plot.barh(ax=ax[1, 0], color="deepskyblue"); ax[1, 0].set_title("Flagged flights by true scenario")
|
|
ax[1, 1].axis("off")
|
|
ax[1, 1].text(0, .9, f"Accuracy {acc:.1%}\nRecall {rec:.1%}\nFalse alarm {far:.1%}\nQueue {len(queue.items)} pending\nLatency {np.mean(ms):.0f} ms/flight",
|
|
family="monospace", fontsize=13, va="top")
|
|
fig.suptitle("ISR03 integrated drone-log analysis - training data only"); plt.tight_layout(); plt.show()
|
|
|
|
|
|
dashboard()
|
|
|
|
# ---------------- 5. 자체 점검(통합 테스트) ----------------
|
|
log_ok = queue.verify()
|
|
cov = len(CALLED & TARGETS) / len(TARGETS)
|
|
chk = {
|
|
"손상 로그 차단": "FL0007" in rejected,
|
|
"처리 비행 99건": len(res) == 99,
|
|
"이상 비행 재현율 90% 이상": rec >= 0.9,
|
|
"정확도 90% 이상": acc >= 0.9,
|
|
"자동 채택 0건(전부 검토 대기)": all(i["status"] == "PENDING" for i in queue.items),
|
|
"탐지 건 전량 대기열 등록": len(queue.items) == int(R.pred.sum()),
|
|
"감사 로그 무결성": log_ok,
|
|
"단정 표현('확정') 0건": not any("확정" in i["text"] for i in queue.items),
|
|
"비행 1건 200ms 이내": float(np.mean(ms)) < 200,
|
|
"기능 커버리지 85% 이상": cov >= 0.85,
|
|
}
|
|
passed = sum(chk.values())
|
|
for k, v in chk.items():
|
|
print(f"[{'PASS' if v else 'FAIL'}] {k}")
|
|
print(f"자체 점검 {passed}/{len(chk)} 통과, 기능 커버리지 {cov:.0%} ({len(CALLED & TARGETS)}/{len(TARGETS)}), 파이프라인 실행 성공 {'예' if passed == len(chk) else '아니오'}")
|
|
|
|
# ---------------- 6. ISR 정보판단 브리핑 요약 ----------------
|
|
print(f"""
|
|
[드론 비행로그 이상 탐지 브리핑 요약 - 교육용가상자료]
|
|
1. 처리: 비행로그 {len(ops)}건 중 {len(res)}건 분석, 손상 로그 {len(rejected)}건 차단({', '.join(rejected)})
|
|
2. 성능: 정확도 {acc:.1%}, 이상 재현율 {rec:.1%}, 오경보율 {far:.1%}, 비행당 {np.mean(ms):.0f} ms, 전체 {time.perf_counter() - t_start:.0f} s
|
|
3. 검토 대기: {len(queue.items)}건 (High {int((R.grade == 'High').sum())} / Medium {int((R.grade == 'Medium').sum())} / Low {int((R.grade == 'Low').sum())}), 자동 채택 0건
|
|
4. 운용 주체: High 등급 {len(high)}건 중 실제 비적대 {int((high.op != 'hostile').sum())}건({hfar:.0%}) -> 등급은 검토 우선순위일 뿐 귀속 결론이 아님
|
|
5. 한계: 합성 로그 기준 결과이며, 운용 주체 판단은 출처·무결성 확인과 2인 교차 검토를 거친 정보분석관 의견으로만 보고한다.""")
|
|
```
|
|
|
|
실행 결과 예시
|
|
|
|
```text
|
|
[PASS] 손상 로그 차단
|
|
[PASS] 처리 비행 99건
|
|
[PASS] 이상 비행 재현율 90% 이상
|
|
[PASS] 정확도 90% 이상
|
|
[PASS] 자동 채택 0건(전부 검토 대기)
|
|
[PASS] 탐지 건 전량 대기열 등록
|
|
[PASS] 감사 로그 무결성
|
|
[PASS] 단정 표현('확정') 0건
|
|
[PASS] 비행 1건 200ms 이내
|
|
[PASS] 기능 커버리지 85% 이상
|
|
자체 점검 10/10 통과, 기능 커버리지 100% (12/12), 파이프라인 실행 성공 예
|
|
|
|
[드론 비행로그 이상 탐지 브리핑 요약 - 교육용가상자료]
|
|
1. 처리: 비행로그 100건 중 99건 분석, 손상 로그 1건 차단(FL0007)
|
|
2. 성능: 정확도 96.0%, 이상 재현율 98.0%, 오경보율 5.8%, 비행당 31 ms, 전체 12 s
|
|
3. 검토 대기: 52건 (High 7 / Medium 22 / Low 23), 자동 채택 0건
|
|
4. 운용 주체: High 등급 7건 중 실제 비적대 2건(29%) -> 등급은 검토 우선순위일 뿐 귀속 결론이 아님
|
|
5. 한계: 합성 로그 기준 결과이며, 운용 주체 판단은 출처·무결성 확인과 2인 교차 검토를 거친 정보분석관 의견으로만 보고한다.
|
|
```
|