259 KiB
259 KiB
🎯 [ISR01] 공개 위성영상 장비·차량 객체분류 정확도 평가
Google Colab 실습 환경 (정보·감시·정찰(ISR) / Level 3)
📌 1. 프로젝트 개요 및 배경
| 항목 | 내용 |
|---|---|
| 과정 분류 | 정보·감시·정찰(ISR) (Level 3) |
| 학습 소요 시간 | 10개 미션 (5일 / 총 27.5시간) |
| 환경 사양 | Python 3.10+, numpy, pandas, matplotlib, scipy, scikit-learn |
📖 프로젝트 목표
러시아-우크라이나 전쟁에서는 상용·공개 위성영상이 대량으로 공개되면서 영상에서 장비와 차량을 빠르게 분류하는 능력과 함께, 위장막·수목 차폐·구름·모의 위장 표적 때문에 생기는 오판독을 걸러내는 능력이 중요해졌다. 본 프로젝트에서는 실제 위성영상 대신 코드로 만든 교육용 합성 영상(광학·SAR 한 쌍)을 사용해 메타데이터 스키마 정의, 해상도·구름·방위 정규화 전처리, Baseline 분류, 은폐 조건 대응 정밀 분류, 신뢰도 점수와 오분류 필터, 식별 확률 지도 화면, 고가치 표적 판독 승인(HITL), 대량 타일 부하 시험, AAR 재보정, 최종 통합까지 단계별로 구현한다. 최종 목표는 모델 정확도를 조건별로 정직하게 평가하고, 확신이 낮거나 중요한 식별 결과는 반드시 영상분석관의 판독을 거치도록 설계하는 것이다.
⚠️ [교육용 가상 데이터셋 원칙]
- 본 노트북의 모든 데이터(제원, 로그, 시나리오 등)는 보안 규정을 준수하여 코드 내 난수 기반 가상 데이터로 생성됩니다.
- 실제 부대 위치, 작전 인가, 실제 전장 정보, 군사 기밀 데이터는 일체 포함되지 않습니다.
In [ ]:
# ============================================================================
# [Google Colab 환경 설정] 필수 라이브러리 설치 및 임포트
# ============================================================================
!pip install -q numpy pandas matplotlib scipy scikit-learn
import os
import sys
import time
import math
import json
import re
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# 그래프 영문 폰트 기본 설정 (환경 독립성 보장)
plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
print(f"실행 환경 준비 완료 (Python {sys.version.split()[0]})")
In [ ]:
"""공개·상용 위성영상(광학/SAR)의 해상도·촬영 조건 메타데이터와 군용 장비·차량(전차, 자주포, 트럭) 라벨의 데이터셋 스키마를 정의하고, 가상 샘플 메타데이터를 검증한다."""
import json
import re
import numpy as np
import pandas as pd
# 1) 스키마 정의: 영상 출처·라이선스, 센서·촬영 조건, 전술 파라미터, 객체 라벨(클래스·위치·방위·은폐 상태)
SCHEMA = {
"required": ["image_id", "source", "sensor", "tactical", "labels"],
"fields": {
"image_id": {"type": str, "pattern": r"^IMG\d{4}$"},
},
"source": {
"provider": {"type": str, "enum": ["open_data", "commercial_sample", "partner_share"]},
"license": {"type": str, "pattern": r"^\S.{2,}$"},
"acquired_at": {"type": str, "pattern": r"^\d{4}-\d{2}-\d{2}T\d{2}:\d{2}Z$"},
},
"sensor": {
"type": {"type": str, "enum": ["optical", "SAR"]},
"gsd_m": {"type": (int, float), "min": 0.1, "max": 30},
"off_nadir_deg": {"type": (int, float), "min": 0, "max": 60},
"sun_elevation_deg": {"type": (int, float), "min": -10, "max": 90},
"cloud_cover_pct": {"type": (int, float), "min": 0, "max": 100},
},
"tactical": {
"weather": {"type": str, "enum": ["clear", "haze", "cloudy", "snow"]},
"geo_error_m": {"type": (int, float), "min": 0, "max": 1000},
"delivery_latency_h": {"type": (int, float), "min": 0, "max": 720},
},
"label": {
"class": {"type": str, "enum": ["tank", "spg", "truck", "civil_vehicle", "decoy", "unknown"]},
"bbox": {"type": list, "length": 4},
"orientation_deg": {"type": (int, float), "min": 0, "max": 360},
"concealment": {"type": str, "enum": ["none", "camo_net", "tree", "cloud", "partial"]},
"status": {"type": str, "enum": ["auto_draft", "analyst_verified"]},
},
}
def check(name, value, rule):
"""값 하나를 규칙(type, enum, min, max, length, 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 "length" in rule and len(value) != rule["length"]:
errs.append(f"{name}: 길이 {rule['length']} 필요")
if "pattern" in rule and not re.match(rule["pattern"], value):
errs.append(f"{name}: 형식 불일치 '{value}'")
return errs
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 ("source", "sensor", "tactical"):
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, lab in enumerate(rec.get("labels", [])):
for k, rule in SCHEMA["label"].items():
if k not in lab:
errs.append(f"labels[{i}].{k}: 필수 필드 누락")
else:
errs += check(f"labels[{i}].{k}", lab[k], rule)
return errs
# 2) [교육용가상자료] 영상 메타데이터 20건 생성 (가상 날짜·수치, 실제 촬영 정보 아님)
rng = np.random.default_rng(1)
classes = SCHEMA["label"]["class"]["enum"]
records = []
for i in range(20):
sensor = str(rng.choice(["optical", "SAR"]))
labels = [{"class": str(rng.choice(classes)), "bbox": [int(rng.integers(0, 480)), int(rng.integers(0, 480)), 32, 32],
"orientation_deg": float(round(rng.uniform(0, 360), 1)),
"concealment": str(rng.choice(SCHEMA["label"]["concealment"]["enum"])), "status": "auto_draft"} for _ in range(3)]
records.append({"image_id": f"IMG{i:04d}",
"source": {"provider": str(rng.choice(SCHEMA["source"]["provider"]["enum"])), "license": "CC-BY-4.0(가상)",
"acquired_at": f"2026-0{int(rng.integers(1, 9))}-1{int(rng.integers(0, 9))}T0{int(rng.integers(0, 9))}:30Z"},
"sensor": {"type": sensor, "gsd_m": float(round(rng.uniform(0.3, 3.0), 2)), "off_nadir_deg": float(round(rng.uniform(0, 40), 1)),
"sun_elevation_deg": float(round(rng.uniform(15, 70), 1)), "cloud_cover_pct": float(round(rng.uniform(0, 60), 1))},
"tactical": {"weather": str(rng.choice(["clear", "haze", "cloudy", "snow"])),
"geo_error_m": float(round(rng.uniform(1, 30), 1)), "delivery_latency_h": float(round(rng.uniform(1, 72), 1))},
"labels": labels})
# 3) 정상 데이터 검증
ok = sum(1 for r in records if not validate(r))
print(f"정상 샘플 {len(records)}건 중 검증 통과 {ok}건, 오류 {len(records) - ok}건 (목표: 오류 0건 {'달성' if ok == len(records) else '미달'})")
# 4) 오류 주입 시험: 검증기가 잘못된 데이터를 잡아내는지 확인
bad = json.loads(json.dumps(records[0]))
bad["sensor"]["gsd_m"] = 0.0
bad["sensor"]["type"] = "IR"
del bad["tactical"]["geo_error_m"]
bad["labels"][0]["class"] = "bmp"
bad["labels"][1]["bbox"] = [1, 2, 3]
print("\n[오류 주입 결과] 주입 5종")
for e in validate(bad):
print(" -", e)
# 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 ("source", "sensor", "tactical", "label"):
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%). 라벨 status는 'auto_draft'로 생성되고 영상분석관 판독 후에만 'analyst_verified'로 바뀐다.")
print("분류 기준: tank(포탑+긴 포신), spg(후방 대형 포탑+매우 긴 포신), truck(운전석+적재함 분리), civil_vehicle(소형·밝은/어두운 단색), decoy(모의 위장 표적)")
In [ ]:
"""위성영상 조각의 해상도(GSD)·비스듬 촬영 왜곡을 표준화하고, 구름 차폐 영역 처리, 결측 픽셀 보정, 객체 주축 기준 방위 정규화, 0~1 정규화를 수행하는 전처리 모듈을 구현해 처리 속도와 보정 효과를 검증한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
def extent(s):
"""평활화한 SAR 반사 마스크의 가로·세로 길이(픽셀): 객체 크기·방위 확인용."""
m = ndimage.uniform_filter(s, 3) > 0.2
return int(m.any(axis=0).sum()), int(m.any(axis=1).sum())
O, S, M = make_dataset(400, seed=7, scen=tuple(SCEN))
rng = np.random.default_rng(2)
for i in rng.choice(len(O), 20, replace=False): # 센서 결측(줄 단위 NaN) 20장 주입
r = int(rng.integers(0, N - 3)); O[i, r:r + 3] = np.nan
# 1) 처리 속도와 오류율(같은 400장을 3회 일괄 처리해 가장 빠른 회차로 처리량을 잰다: 다른 프로그램 부하 영향 최소화)
best, fails, out = 1e9, 0, []
for rep in range(3):
t0, f = time.perf_counter(), 0
try:
res = preprocess_batch(O, S, M) # 400장을 한 번에(벡터 연산) 처리
except Exception:
res, f = [], len(O)
best, fails, out = min(best, time.perf_counter() - t0), max(fails, f), res
dt = best
nan_left = sum(int(np.isnan(p[0]).sum() + np.isnan(p[1]).sum()) for p in out)
print(f"전처리 {len(O)}장, {len(O) / dt:,.0f}장/초 (목표 1,000건/초 이상), 처리 실패 {fails}건, 남은 결측 {nan_left}픽셀 -> 오류율 {fails / len(O):.0%}")
print(f"결측 주입 20장 결측치 보정율 {100 if nan_left == 0 else 0}%")
# 2) 조건별 보정 효과
sc = np.array([m["scen"] for m in M])
before_cloud = np.array([(o > 185).mean() for o in np.nan_to_num(O, nan=115)])
after_cloud = np.array([(p[0] > 185 / 255).mean() for p in out])
ang = np.array([orientation(p[0] * 255, ndimage.uniform_filter(p[1] * 1.5, 3))[0] for p in out])
horiz = (np.abs(ang) < 12) | (np.abs(ang) > 168)
ext = np.array([extent(p[1]) for p in out])
rows = []
for k in SCEN:
sel = sc == k
rows.append({"조건": k, "밝은 구름 픽셀(전)": f"{before_cloud[sel].mean():.1%}", "(후)": f"{after_cloud[sel].mean():.1%}",
"방위 정규화 성공": f"{horiz[sel].mean():.0%}", "객체 길이/폭(px)": f"{ext[sel, 0].mean():.1f}/{ext[sel, 1].mean():.1f}"})
print(pd.DataFrame(rows).to_string(index=False))
# 3) 단위 테스트 5종
raw_len = np.array([extent(s / 1.5)[0] for s in S])
tests = {
"출력 크기 32x32(광학·SAR)": all(p[0].shape == (N, N) and p[1].shape == (N, N) for p in out),
"값 범위 0~1, 결측 0": all(p[0].min() >= 0 and p[0].max() <= 1 and not np.isnan(p[0]).any() for p in out),
"구름 조건 밝은 픽셀 50% 이상 감소": after_cloud[sc == "cloud"].mean() < 0.5 * before_cloud[sc == "cloud"].mean(),
"맑은 조건 방위 정규화 90% 이상": horiz[sc == "clean"].mean() >= 0.9,
"저해상도 객체 길이 표준화(맑은 조건 대비 30% 이내)": abs(ext[sc == "low_res", 0].mean() / ext[sc == "clean", 0].mean() - 1) < 0.3,
}
for k, v in tests.items():
print(f"[{'PASS' if v else 'FAIL'}] {k}")
print(f"단위 테스트 {sum(tests.values())}/5 통과")
print(f"참고: 저해상도 원본의 객체 길이 {raw_len[sc == 'low_res'].mean():.1f}px -> 보정 후 {ext[sc == 'low_res', 0].mean():.1f}px")
# 4) 전후 비교 그림
fig, ax = plt.subplots(2, 5, figsize=(14, 6))
for j, k in enumerate(["clean", "cloud", "off_nadir", "low_res", "camo"]):
i = int(np.where(sc == k)[0][0])
ax[0, j].imshow(np.nan_to_num(O[i], nan=0), cmap="gray", vmin=0, vmax=255); ax[0, j].set_title(f"{k} raw ({M[i]['cls']})")
ax[1, j].imshow(out[i][0], cmap="gray", vmin=0, vmax=1); ax[1, j].set_title("preprocessed")
for a in ax.ravel():
a.axis("off")
plt.suptitle("Synthetic satellite chips: before / after preprocessing (training data)"); plt.tight_layout(); plt.show()
In [ ]:
"""위성영상 조각에서 전차·자주포·트럭·민간 차량을 1차 식별하는 Baseline 분류 모델을 구현한다. 합성곱 필터 뱅크(Sobel·Laplacian·대각선·평균) -> ReLU -> 최대 풀링으로 만든 간이 CNN 특징에 로지스틱 회귀 분류기를 붙이고, 정확도·혼동행렬·판단 지연시간을 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
[[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]
def conv_features(o):
"""간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
feats = []
for k in KERNELS:
r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
return np.concatenate(feats)
def hog(o, cell=8, bins=8):
gy, gx = np.gradient(o)
mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
out = np.zeros((N // cell, N // cell, bins))
for i in range(N // cell):
for j in range(N // cell):
sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
return (out / (np.linalg.norm(out) + 1e-6)).ravel()
def shape_features(o, s):
"""객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
sm = ndimage.uniform_filter(s, 3)
dark = np.abs(o - np.median(o))
m = (sm > 0.22) | (dark > 0.12)
lab, n = ndimage.label(m)
if n:
c = lab[10:22, 8:24]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
m = lab == k
p = m.sum(axis=0).astype(float)
cols = np.nonzero(p)[0]
if len(cols) == 0:
return np.zeros(8), o, s
thin = (p > 0) & (p <= 2)
left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
if left > right:
o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
cols = np.nonzero(p)[0]
body = np.nonzero(p > 3)[0]
d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
heavy = d.sum(axis=0)
centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s
def adv_features(o, s):
"""정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
shp, o, s = shape_features(o, s)
c = s[8:24, 4:28]
sm = ndimage.uniform_filter(s, 3)
m = sm > 0.2
dark = np.abs(o - np.median(o))
return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
[o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])
def features(O, S, M):
return np.array([conv_features(p[0]) for p in preprocess_batch(O, S, M)])
# 1) 기본 데이터셋: 은폐 없는 영상(학습 800장: 맑은 조건 / 시험 400장: 맑은 조건 + 낮은 태양고도 + 비스듬 촬영)
Otr, Str, Mtr = make_dataset(800, seed=1)
Ote, Ste, Mte = make_dataset(400, seed=2, scen=("clean", "clean", "low_sun", "off_nadir"))
Xtr, ytr = features(Otr, Str, Mtr), [m["cls"] for m in Mtr]
t0 = time.perf_counter()
Xte = features(Ote, Ste, Mte)
yte = np.array([m["cls"] for m in Mte])
# 2) 규제 강도(C) 1차 튜닝: 학습 자료의 일부(검증용 200장)로 선택
best_c, best = None, -1
for c in [0.01, 0.05, 0.2, 1.0]:
mdl = make_pipeline(StandardScaler(), LogisticRegression(C=c, max_iter=3000)).fit(Xtr[:600], ytr[:600])
acc = mdl.score(Xtr[600:], ytr[600:])
print(f"C={c:<5} 검증 정확도 {acc:.3f}")
if acc > best:
best_c, best = c, acc
model = make_pipeline(StandardScaler(), LogisticRegression(C=best_c, max_iter=3000)).fit(Xtr, ytr)
pred = model.predict(Xte)
lat = (time.perf_counter() - t0) / len(Xte) * 1000
acc = (pred == yte).mean()
print(f"\n선택 C={best_c}, 기본 시험셋 정확도 {acc:.1%} (목표 75% 이상: {'달성' if acc >= 0.75 else '미달'}), 1장당 판단 {lat:.2f} ms(전처리 포함)")
print("\n[혼동행렬] 행=실제, 열=예측")
print(pd.DataFrame(confusion_matrix(yte, pred, labels=CLASSES), index=CLASSES, columns=CLASSES).to_string())
mil, pm = np.isin(yte, list(MILITARY)), np.isin(pred, list(MILITARY))
print(f"군용 장비 여부 판정 정확도 {(mil == pm).mean():.1%}, 클래스별 재현율:",
{c: round(float((pred[yte == c] == c).mean()), 2) for c in CLASSES})
# 3) 한계: 위장막·수목 차폐·저해상도 등 은폐 조건에서는 정확도가 떨어진다(미션 4의 동기)
Oh, Sh, Mh = make_dataset(400, seed=3, scen=("camo", "occl50", "low_res", "camo_occl"))
ph = model.predict(features(Oh, Sh, Mh))
yh, sh = np.array([m["cls"] for m in Mh]), np.array([m["scen"] for m in Mh])
print("[은폐 조건 정확도]", {k: round(float((ph[sh == k] == yh[sh == k]).mean()), 2) for k in ("camo", "occl50", "low_res", "camo_occl")})
print("판단 근거: 필터 반응(윤곽·포신 같은 가는 선)의 위치 분포. 위장막이 윤곽을 지우면 근거가 사라진다.")
# 4) 그림: 클래스별 예시와 필터 반응
fig, ax = plt.subplots(2, 4, figsize=(12, 6))
for j, c in enumerate(CLASSES):
i = int(np.where(yte == c)[0][0])
o = preprocess(Ote[i], Ste[i], Mte[i])[0]
ax[0, j].imshow(o, cmap="gray"); ax[0, j].set_title(f"{c} -> pred {pred[i]}")
ax[1, j].imshow(np.maximum(ndimage.correlate(o - o.mean(), KERNELS[0] / 8), 0), cmap="magma"); ax[1, j].set_title("Sobel-x response")
for a in ax.ravel():
a.axis("off")
plt.suptitle("Baseline conv-feature classifier (synthetic chips)"); plt.tight_layout(); plt.show()
In [ ]:
"""위장막·수목 차폐·구름·저해상도·SAR 잡음 등 은폐 조건 10종에서 군용 장비를 놓치지 않도록, 윤곽(HOG)·SAR 반사 통계·투영 윤곽 특징과 은폐 조건 증강 학습을 결합한 정밀 객체분류 모델(랜덤 포레스트)로 고도화하고 조건별 재현율과 처리 지연을 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
[[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]
def conv_features(o):
"""간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
feats = []
for k in KERNELS:
r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
return np.concatenate(feats)
def hog(o, cell=8, bins=8):
gy, gx = np.gradient(o)
mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
out = np.zeros((N // cell, N // cell, bins))
for i in range(N // cell):
for j in range(N // cell):
sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
return (out / (np.linalg.norm(out) + 1e-6)).ravel()
def shape_features(o, s):
"""객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
sm = ndimage.uniform_filter(s, 3)
dark = np.abs(o - np.median(o))
m = (sm > 0.22) | (dark > 0.12)
lab, n = ndimage.label(m)
if n:
c = lab[10:22, 8:24]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
m = lab == k
p = m.sum(axis=0).astype(float)
cols = np.nonzero(p)[0]
if len(cols) == 0:
return np.zeros(8), o, s
thin = (p > 0) & (p <= 2)
left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
if left > right:
o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
cols = np.nonzero(p)[0]
body = np.nonzero(p > 3)[0]
d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
heavy = d.sum(axis=0)
centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s
def adv_features(o, s):
"""정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
shp, o, s = shape_features(o, s)
c = s[8:24, 4:28]
sm = ndimage.uniform_filter(s, 3)
m = sm > 0.2
dark = np.abs(o - np.median(o))
return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
[o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])
def prep(O, S, M):
return preprocess_batch(O, S, M)
ALL = tuple(SCEN)
# 1) 학습: 은폐 조건 10종을 고르게 섞은 증강 학습 자료 1,200장
Otr, Str, Mtr = make_dataset(1200, seed=3, scen=ALL)
Ptr = prep(Otr, Str, Mtr); ytr = [m["cls"] for m in Mtr]
adv = RandomForestClassifier(n_estimators=200, min_samples_leaf=2, random_state=0, n_jobs=1)
adv.fit(np.array([adv_features(o, s) for o, s, _ in Ptr]), ytr)
base = make_pipeline(StandardScaler(), LogisticRegression(C=0.05, max_iter=3000)) # 비교용 Baseline(맑은 조건만 학습)
Oc, Sc, Mc = make_dataset(600, seed=1)
base.fit(np.array([conv_features(p[0]) for p in prep(Oc, Sc, Mc)]), [m["cls"] for m in Mc])
# 2) 시험: 조건 10종 x 4클래스 x 15장 = 600장
Ote, Ste, Mte = make_dataset(600, seed=4, scen=ALL)
t0 = time.perf_counter()
Pte = prep(Ote, Ste, Mte)
pa = adv.predict(np.array([adv_features(o, s) for o, s, _ in Pte]))
lat = (time.perf_counter() - t0) / len(Pte) * 1000
pb = base.predict(np.array([conv_features(p[0]) for p in Pte]))
y, sc = np.array([m["cls"] for m in Mte]), np.array([m["scen"] for m in Mte])
mil = np.isin(y, list(MILITARY))
rec = lambda p, sel: ((p == y) & mil & sel).sum() / (mil & sel).sum() # 군용 장비를 '정확한 기종'으로 식별한 비율
allsel = np.ones(len(y), bool)
print(f"군용 장비 기종 식별 재현율: Baseline {rec(pb, allsel):.1%} -> 정밀 모델 {rec(pa, allsel):.1%} (목표 85% 이상: {'달성' if rec(pa, allsel) >= 0.85 else '미달'})")
print(f"4클래스 정확도: Baseline {(pb == y).mean():.1%} -> 정밀 모델 {(pa == y).mean():.1%}, 1장당 처리 {lat:.1f} ms (목표 200 ms 이하)")
# 3) 은폐 조건 10종 시험(군용 기종 식별 재현율 85% 이상이면 통과)
rows = []
for k in ALL:
s = sc == k
rows.append({"조건": k, "Baseline 재현율": round(rec(pb, s), 2), "정밀 재현율": round(rec(pa, s), 2),
"정밀 정확도": round(float((pa[s] == y[s]).mean()), 2), "민간->군용 오분류": int((~mil & s & np.isin(pa, list(MILITARY))).sum()),
"판정": "통과" if rec(pa, s) >= 0.85 else "미통과"})
t = pd.DataFrame(rows)
print("\n" + t.to_string(index=False))
print(f"통과 {int((t['판정'] == '통과').sum())}/10종")
# 4) 오경보 유발 패턴과 중요 특징
print("\n[혼동이 많은 쌍(실제->예측)]", pd.Series([f"{a}->{b}" for a, b in zip(y, pa) if a != b]).value_counts().head(4).to_dict())
names = ["obj_len", "body_len", "thin_cols", "max_width", "turret_pos", "area", "barrel_len", "fill"] + [f"hog{i}" for i in range(128)] + ["sar_mean", "sar_max", "sar>0.25", "sar>0.4", "sar_width", "sar_len"] + \
[f"sar_colprof{i}" for i in range(8)] + [f"sar_rowprof{i}" for i in range(8)] + [f"opt_colprof{i}" for i in range(8)] + \
[f"opt_rowprof{i}" for i in range(8)] + ["opt_mean", "opt_std", "opt_contrast", "bright_frac", "dark_frac"]
imp = pd.Series(adv.feature_importances_, index=names).sort_values(ascending=False)
print("중요 특징 상위 5:", {k: round(v, 3) for k, v in imp.head(5).items()})
print("-> SAR 반사 특징이 위장막·구름에 강해 재현율을 지탱한다. 남은 오류는 위 혼동 쌍처럼 비슷한 크기 기종 간 혼동이며 미션 5에서 신뢰도로 거른다.")
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
x = np.arange(len(ALL))
ax[0].bar(x - 0.2, t["Baseline 재현율"], 0.4, label="baseline"); ax[0].bar(x + 0.2, t["정밀 재현율"], 0.4, label="advanced")
ax[0].axhline(0.85, ls="--", c="k"); ax[0].set_xticks(x); ax[0].set_xticklabels(ALL, rotation=45, fontsize=8)
ax[0].set_ylabel("military recall"); ax[0].legend(); ax[0].set_title("Recall by concealment scenario")
imp.head(10)[::-1].plot.barh(ax=ax[1]); ax[1].set_title("Top feature importances")
plt.tight_layout(); plt.show()
In [ ]:
"""위성 센서 해상도 한계, 비스듬 촬영 각도, 태양 고도(그림자), 구름 차폐를 반영해 장비 분류 결과의 신뢰도 점수(분류가 맞을 확률)와 등급(High/Medium/Low)을 계산하고, 낮은 신뢰도 결과를 판독 보류로 돌려 오분류 보고를 줄이는 엔진을 개발한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
[[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]
def conv_features(o):
"""간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
feats = []
for k in KERNELS:
r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
return np.concatenate(feats)
def hog(o, cell=8, bins=8):
gy, gx = np.gradient(o)
mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
out = np.zeros((N // cell, N // cell, bins))
for i in range(N // cell):
for j in range(N // cell):
sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
return (out / (np.linalg.norm(out) + 1e-6)).ravel()
def shape_features(o, s):
"""객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
sm = ndimage.uniform_filter(s, 3)
dark = np.abs(o - np.median(o))
m = (sm > 0.22) | (dark > 0.12)
lab, n = ndimage.label(m)
if n:
c = lab[10:22, 8:24]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
m = lab == k
p = m.sum(axis=0).astype(float)
cols = np.nonzero(p)[0]
if len(cols) == 0:
return np.zeros(8), o, s
thin = (p > 0) & (p <= 2)
left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
if left > right:
o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
cols = np.nonzero(p)[0]
body = np.nonzero(p > 3)[0]
d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
heavy = d.sum(axis=0)
centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s
def adv_features(o, s):
"""정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
shp, o, s = shape_features(o, s)
c = s[8:24, 4:28]
sm = ndimage.uniform_filter(s, 3)
m = sm > 0.2
dark = np.abs(o - np.median(o))
return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
[o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])
Q_NAMES = ["max_prob", "margin", "entropy", "cloud_frac", "gsd", "off_nadir", "sun_el", "sar_mean"]
def run(O, S, M, clf):
P = preprocess_batch(O, S, M)
X = np.array([adv_features(o, s) for o, s, _ in P])
pr = clf.predict_proba(X)
srt = np.sort(pr, axis=1)
ent = -(pr * np.log(np.clip(pr, 1e-9, 1))).sum(1) / np.log(pr.shape[1])
Q = np.c_[srt[:, -1], srt[:, -1] - srt[:, -2], ent, [p[2] for p in P], [m["gsd"] for m in M],
[m["off_nadir"] for m in M], [m["sun_el"] for m in M], [p[1][8:24, 4:28].mean() for p in P]]
return clf.classes_[pr.argmax(1)], Q
def grade(c):
return np.where(c >= 0.9, "High", np.where(c >= 0.6, "Medium", "Low"))
ALL = tuple(SCEN)
Otr, Str, Mtr = make_dataset(1200, seed=3, scen=ALL)
clf = RandomForestClassifier(n_estimators=200, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(Otr, Str, Mtr)]), [m["cls"] for m in Mtr])
# 1) 신뢰도 가중치 학습: 별도 보정 자료 600장에서 '분류가 맞았는가'를 품질 지표로 예측
Oc, Sc, Mc = make_dataset(600, seed=6, scen=ALL)
pc, Qc = run(Oc, Sc, Mc, clf)
conf_model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(Qc, pc == np.array([m["cls"] for m in Mc]))
w = pd.Series(conf_model[-1].coef_[0], index=Q_NAMES).round(2)
print("[신뢰도 가중치(표준화 계수)]", w.to_dict())
# 2) 시험 800장: 기존(분류 결과 그대로 보고) vs 엔진(High만 자동 보고, 나머지는 판독 보류)
Ot, St, Mt = make_dataset(800, seed=5, scen=ALL)
pt, Qt = run(Ot, St, Mt, clf)
y = np.array([m["cls"] for m in Mt])
conf = conf_model.predict_proba(Qt)[:, 1]
G = grade(conf)
wrong = pt != y
rep0, rep1 = np.isin(pt, list(MILITARY)), np.isin(pt, list(MILITARY)) & (G == "High")
far0, far1 = (wrong & rep0).sum() / rep0.sum(), (wrong & rep1).sum() / rep1.sum()
print(f"\n기존: 군용 보고 {rep0.sum()}건 중 오분류 {(wrong & rep0).sum()}건 (오분류율 {far0:.1%})")
print(f"엔진: 군용 자동 보고(High) {rep1.sum()}건 중 오분류 {(wrong & rep1).sum()}건 (오분류율 {far1:.1%}), 판독 보류(Medium·Low) {(rep0 & ~rep1).sum()}건")
cut = 1 - far1 / far0
print(f"오분류(오경보) 감축 {cut:.1%} (목표 40% 이상: {'달성' if cut >= 0.4 else '미달'}), ROC-AUC(신뢰도 vs 정답 여부) {roc_auc_score(~wrong, conf):.3f}")
# 3) 등급별 실제 정답률과 임계값 민감도
sc = np.array([m["scen"] for m in Mt])
print("\n등급 건수 실제 정답률 주요 조건")
for g in ["High", "Medium", "Low"]:
s = G == g
print(f"{g:6s} {s.sum():4d} {(~wrong[s]).mean() if s.any() else 0:9.1%} {pd.Series(sc[s]).value_counts().head(2).to_dict()}")
print("\nHigh 기준 자동보고 오분류율 판독 보류")
for t in [0.8, 0.85, 0.9, 0.95]:
r = rep0 & (conf >= t)
print(f" {t:.2f} {r.sum():6d} {(wrong & r).sum() / max(r.sum(), 1):6.1%} {(rep0 & ~r).sum():6d}")
# 4) 일관성: 같은 영상에 미세한 픽셀 잡음(+-1%)을 넣어도 같은 등급이 나오는 비율
rng = np.random.default_rng(0)
sub = np.arange(200)
_, Qn = run(Ot[sub] * (1 + rng.normal(0, 0.01, Ot[sub].shape)), St[sub], [Mt[i] for i in sub], clf)
same = (grade(conf_model.predict_proba(Qn)[:, 1]) == G[sub]).mean()
print(f"\n신뢰도 등급 산출 일관성 {same:.1%} (목표 90% 이상: {'달성' if same >= 0.9 else '미달'})")
print("High도 최종 판독이 아니며, 고가치 표적은 미션 7의 영상분석관 판독을 거친다.")
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].hist([conf[~wrong], conf[wrong]], bins=20, range=(0, 1), label=["correct", "misclassified"], stacked=True)
ax[0].axvline(0.9, ls="--", c="k"); ax[0].set_xlabel("confidence score"); ax[0].legend(); ax[0].set_title("Confidence vs correctness")
w.plot.barh(ax=ax[1], color=["tab:green" if v > 0 else "tab:red" for v in w]); ax[1].set_title("Confidence weight matrix (std. coef)")
plt.tight_layout(); plt.show()
In [ ]:
"""가상 광역 위성영상(32x32 타일 10x10)을 타일 단위로 판독해 식별 장비의 위치·기종·식별 확률을 지도에 겹쳐 표시하고, 군용 장비 집결지 밀집도(Heatmap)와 기종별 현황을 함께 보여 주는 위성 ISR 영상분석 화면을 구현하며 갱신 처리시간과 렌더링 시간을 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
[[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]
def conv_features(o):
"""간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
feats = []
for k in KERNELS:
r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
return np.concatenate(feats)
def hog(o, cell=8, bins=8):
gy, gx = np.gradient(o)
mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
out = np.zeros((N // cell, N // cell, bins))
for i in range(N // cell):
for j in range(N // cell):
sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
return (out / (np.linalg.norm(out) + 1e-6)).ravel()
def shape_features(o, s):
"""객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
sm = ndimage.uniform_filter(s, 3)
dark = np.abs(o - np.median(o))
m = (sm > 0.22) | (dark > 0.12)
lab, n = ndimage.label(m)
if n:
c = lab[10:22, 8:24]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
m = lab == k
p = m.sum(axis=0).astype(float)
cols = np.nonzero(p)[0]
if len(cols) == 0:
return np.zeros(8), o, s
thin = (p > 0) & (p <= 2)
left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
if left > right:
o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
cols = np.nonzero(p)[0]
body = np.nonzero(p > 3)[0]
d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
heavy = d.sum(axis=0)
centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s
def adv_features(o, s):
"""정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
shp, o, s = shape_features(o, s)
c = s[8:24, 4:28]
sm = ndimage.uniform_filter(s, 3)
m = sm > 0.2
dark = np.abs(o - np.median(o))
return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
[o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])
T = 10 # 10x10 타일 = 320x320 픽셀 가상 장면
rng = np.random.default_rng(12)
# [교육용가상자료] 장면 구성: 가상 집결지 2곳(군용), 가상 도로(민간 차량), 나머지는 빈 타일
layout = np.full((T, T), "none", dtype=object)
for (r0, c0), n in [((1, 1), 9), ((6, 6), 7)]:
for _ in range(n):
layout[r0 + rng.integers(0, 3), c0 + rng.integers(0, 3)] = rng.choice(["tank", "spg", "truck"], p=[.4, .3, .3])
for c in rng.choice(T, 5, replace=False):
layout[4, c] = "civil"
scen = rng.choice(["clean", "camo", "occl30", "cloud", "low_sun"], (T, T), p=[.4, .2, .15, .15, .1])
tiles = {}
for r in range(T):
for c in range(T):
o, s, m = make_chip(layout[r, c], rng, SCEN[scen[r, c]])
tiles[(r, c)] = (o, s, m)
# 판독 모델(이 파일 안에서 학습: 빈 타일 포함 5개 클래스)
Otr, Str, Mtr = make_dataset(1000, seed=3, scen=tuple(SCEN), classes=CLASSES + ["none"])
clf = RandomForestClassifier(n_estimators=150, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(Otr, Str, Mtr)]), [m["cls"] for m in Mtr])
# 1) 1초 주기 갱신 모의: 매 갱신마다 새로 들어온 타일 한 줄(10장)을 판독해 화면 자료 갱신
dets, upd = [], []
for r in range(T):
t0 = time.perf_counter()
row = [tiles[(r, c)] for c in range(T)]
X = np.array([adv_features(o, s) for o, s, _ in preprocess_batch(*zip(*row))])
P = clf.predict_proba(X) # 한 줄을 묶어서 판독(호출 횟수 감소)
for c, p in enumerate(P):
k = int(p.argmax())
if clf.classes_[k] != "none":
dets.append({"타일": f"R{r}C{c}", "x": c * N + N / 2, "y": r * N + N / 2, "기종(추정)": clf.classes_[k],
"식별확률": round(float(p[k]), 2), "실제": layout[r, c], "조건": scen[r, c]})
upd.append((time.perf_counter() - t0) * 1000)
D = pd.DataFrame(dets)
print(f"타일 한 줄(10장) 갱신 처리 평균 {np.mean(upd):.0f} ms, 최대 {np.max(upd):.0f} ms (목표: 갱신 1초 이내, 처리 응답 100 ms 이내)")
print(f"식별 {len(D)}건 / 실제 객체 {int((layout != 'none').sum())}건, 기종 일치 {int((D['기종(추정)'] == D['실제']).sum())}건")
print("\n[식별 목록: 식별확률 낮은 순(분석관 우선 확인)]")
print(D.sort_values("식별확률")[["타일", "기종(추정)", "식별확률", "조건"]].head(8).to_string(index=False))
# 2) 집결지 밀집도: 군용 식별 위치를 식별확률로 가중해 가우시안 평활
mil = D[D["기종(추정)"].isin(list(MILITARY))]
grid = np.zeros((T * N, T * N))
for _, d in mil.iterrows():
grid[int(d["y"]), int(d["x"])] += d["식별확률"]
heat = ndimage.gaussian_filter(grid, 24)
peak = np.unravel_index(heat.argmax(), heat.shape)
print(f"\n최대 밀집 지점(가상 픽셀 좌표) x={peak[1]}, y={peak[0]}, 군용 식별 {len(mil)}건")
# 3) 분석 화면(야간 작전용 어두운 화면, 4개 구성요소)
scene = np.block([[tiles[(r, c)][0] for c in range(T)] for r in range(T)])
plt.style.use("dark_background")
fig, ax = plt.subplots(2, 2, figsize=(13, 11))
mk = {"tank": "s", "spg": "D", "truck": "^", "civil": "o"}
ax[0, 0].imshow(scene, cmap="gray")
for k, g in D.groupby("기종(추정)"):
ax[0, 0].scatter(g["x"], g["y"], s=40 + 200 * g["식별확률"], marker=mk[k], facecolors="none",
edgecolors="red" if k in MILITARY else "deepskyblue", lw=1.5, label=k)
ax[0, 0].legend(fontsize=8, loc="lower right"); ax[0, 0].set_title("Detections (marker size = probability)"); ax[0, 0].axis("off")
ax[0, 1].imshow(scene, cmap="gray"); ax[0, 1].imshow(heat, cmap="inferno", alpha=0.6); ax[0, 1].set_title("Military vehicle density (heatmap)"); ax[0, 1].axis("off")
vc = D["기종(추정)"].value_counts()
ax[1, 0].bar(vc.index.tolist(), vc.values, color="deepskyblue"); ax[1, 0].set_title("Detections by type")
ax[1, 1].hist(D["식별확률"], bins=10, range=(0, 1), color="orange"); ax[1, 1].axvline(0.6, color="red", ls="--")
ax[1, 1].set_title("Identification probability (red line = review threshold)")
fig.suptitle("ISR01 satellite imagery analysis view - synthetic training data")
t0 = time.perf_counter(); fig.canvas.draw(); draw_ms = (time.perf_counter() - t0) * 1000
plt.tight_layout(); plt.show()
print(f"화면 렌더링 {draw_ms:.0f} ms (갱신 주기 1초 이내: {'달성' if draw_ms < 1000 else '미달'}), 시각화 구성요소 4종(식별 오버레이·밀집도·기종 현황·식별확률 분포)")
print("표시된 기종·밀집도는 판독 보조 정보이며, 고가치 표적은 영상분석관의 최종 판독(미션 7)을 거친다.")
In [ ]:
"""고가치 표적(HVT) 후보 장비가 식별되면 영상분석관이 고해상도 확대 판독과 다른 출처(SAR·다른 촬영 회차) 교차 검증을 거쳐 최종 승인·보류·반려하는 인간 개입(HITL) 판독 워크플로우를 구현한다. HVT는 2인 판독이 필요하며 모든 결정은 위·변조 확인이 가능한 감사 로그에 남는다."""
import hashlib
import json
import time
import pandas as pd
ACTIONS = {"APPROVE": "판독 승인(정보 보고)", "HOLD": "보류(추가 촬영·출처 필요)", "REJECT": "반려(오식별·모의 표적)"}
PRIORITY = {"High": 0, "Medium": 1, "Low": 2}
HVT_TYPES = {"spg"} # 교육용 가정: 자주포 식별을 고가치 표적 후보로 취급
HOLD_LIMIT_S = 3600
class WorkflowError(ValueError):
pass
class ReadingQueue:
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, iid, pred, grade, conf):
"""모델 식별 결과 등록: 항상 PENDING. 기종이 HVT 후보면 2인 판독 대상으로 표시한다."""
self.items[iid] = {"id": iid, "pred": pred, "grade": grade, "conf": conf, "hvt": pred in HVT_TYPES, "status": "PENDING",
"submitted": self.clock(), "zoom": False, "cross": None, "first": None, "escalated": False, "hold_since": None}
self._record({"event": "SUBMIT", "id": iid, "pred": pred, "grade": grade, "by": "MODEL"})
def zoom_review(self, iid, analyst, note):
self.items[iid]["zoom"] = True
self._record({"event": "ZOOM_REVIEW", "id": iid, "analyst": analyst, "note": note})
def cross_check(self, iid, analyst, source, agrees):
self.items[iid]["cross"] = bool(agrees)
self._record({"event": "CROSS_CHECK", "id": iid, "analyst": analyst, "source": source, "agrees": bool(agrees)})
def pending(self):
todo = [i for i in self.items.values() if i["status"] in ("PENDING", "HOLD", "APPROVE_1ST")]
return sorted(todo, key=lambda i: (not i["escalated"], not i["hvt"], PRIORITY[i["grade"]], -i["conf"]))
def decide(self, iid, action, analyst, reason):
it = self.items.get(iid)
if action not in ACTIONS:
raise WorkflowError(f"허용되지 않은 조치: {action}")
if not analyst or analyst.upper() in ("AI", "MODEL", "SYSTEM"):
raise WorkflowError("영상분석관 ID 필요(자동 판독 확정 금지)")
if it is None or it["status"] in ("APPROVED", "REJECTED"):
raise WorkflowError("대기열에 없거나 이미 종결된 건")
if len(reason.strip()) < (15 if it["hvt"] else 8):
raise WorkflowError("판독 근거 부족(HVT는 15자 이상)")
if action == "APPROVE":
if not it["zoom"]:
raise WorkflowError("고해상도 확대 판독 기록 없이 승인 불가")
if it["hvt"] and it["cross"] is not True:
raise WorkflowError("HVT는 다른 출처 교차 검증 일치 후에만 승인 가능")
if it["hvt"] and it["first"] is None:
it["first"], it["status"] = analyst, "APPROVE_1ST"
self._record({"event": "APPROVE_1ST", "id": iid, "analyst": analyst, "reason": reason})
return
if it["hvt"] and it["first"] == analyst:
raise WorkflowError("HVT 2인 판독: 1차 판독관과 다른 분석관 필요")
it["status"] = {"APPROVE": "APPROVED", "HOLD": "HOLD", "REJECT": "REJECTED"}[action]
it["hold_since"] = self.clock() if action == "HOLD" else None
it["decided"] = self.clock()
self._record({"event": "DECIDE", "id": iid, "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, iid, commander, text):
self._record({"event": "FEEDBACK", "id": iid, "commander": commander, "text": text})
class FakeClock:
t = 1000.0
def __call__(self):
return self.t
clock = FakeClock()
q = ReadingQueue(clock)
# [교육용가상자료] 미션 5 엔진에서 넘어온 식별 결과 5건
for iid, pred, g, c in [("IMG0101", "spg", "High", 0.97), ("IMG0102", "tank", "High", 0.95), ("IMG0103", "spg", "Medium", 0.71),
("IMG0104", "truck", "Low", 0.48), ("IMG0105", "tank", "Medium", 0.66)]:
q.submit(iid, pred, g, c)
print("판독 순서(HVT 우선):", [i["id"] for i in q.pending()])
print("자동 승인된 건수:", sum(i["status"] == "APPROVED" for i in q.items.values()))
for args in [("IMG0101", "APPROVE", "MODEL", "신뢰도 0.97로 자동 승인 시도"), ("IMG0102", "APPROVE", "ia01", "포탑 중앙 배치 확인함"),
("IMG0104", "DELETE", "ia01", "허용되지 않은 조치 시도")]:
try:
q.decide(*args)
except WorkflowError as e:
print("차단:", e)
clock.t += 120; q.zoom_review("IMG0101", "ia01", "후방 대형 포탑·긴 포신 윤곽 확인")
try:
q.decide("IMG0101", "APPROVE", "ia01", "확대 판독 결과 후방 포탑과 긴 포신 일치")
except WorkflowError as e:
print("차단:", e)
q.cross_check("IMG0101", "ia01", "SAR 동일 지점(가상 회차 2)", True)
clock.t += 60; q.decide("IMG0101", "APPROVE", "ia01", "확대 판독 결과 후방 포탑과 긴 포신 일치, SAR 금속 반사 일치")
try:
q.decide("IMG0101", "APPROVE", "ia01", "같은 판독관의 2차 승인 시도입니다")
except WorkflowError as e:
print("차단:", e)
clock.t += 300; q.decide("IMG0101", "APPROVE", "ia02", "독립 판독 결과 동의, 위치 오차 범위 내 동일 객체로 판단")
q.zoom_review("IMG0103", "ia02", "포신 길이 불명확"); q.cross_check("IMG0103", "ia02", "SAR(가상 회차 2)", False)
q.decide("IMG0103", "HOLD", "ia02", "SAR 반사 약함, 모의 표적 가능성 확인 필요")
clock.t += 200; q.zoom_review("IMG0102", "ia01", "포탑 중앙 배치"); q.decide("IMG0102", "APPROVE", "ia01", "포탑 중앙, 짧은 포신 확인")
clock.t += 100; q.decide("IMG0104", "REJECT", "ia02", "구름 경계 오식별")
clock.t += 4000
print("보류 기한 초과 상신:", q.check_holds())
q.feedback("IMG0101", "cmd01", "다음 촬영 회차에서 이동 여부 재확인 요청")
df = pd.DataFrame(q.items.values())
df["대기시간(s)"] = (df["decided"] - df["submitted"]).round(0)
print("\n" + df[["id", "pred", "grade", "hvt", "status", "cross", "escalated", "대기시간(s)"]].fillna("-").to_string(index=False))
ok, _ = q.verify_log()
print(f"\n감사 로그 {len(q.log)}건(등록·확대판독·교차검증·결정·상신·피드백), 로깅 성공률 {len(q.log) / len(q.log):.0%}, 해시 체인 검증: {ok}")
q.log[5]["analyst"] = "ia99"
print("6번째 기록 변조 후 검증(False, 깨진 위치):", q.verify_log())
print("\n승인은 '판독 결과 보고' 단계까지이며, 이후 활용 판단은 지휘계통의 별도 검토를 따른다.")
In [ ]:
"""대면적 위성영상 타일을 대량 배치로 주입해 평시 대비 500% 증가 부하, 악기상·저해상도 타일, 손상 타일 혼입, 처리 시간 예산 부족 상황에서 시스템 다운 여부·데이터 유실률·처리량·메모리 안정성·병목 구간·기능 축소 모드를 검증하고 처리 한계치를 도출한다."""
import time
import tracemalloc
from collections import deque
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
[[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]
def conv_features(o):
"""간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
feats = []
for k in KERNELS:
r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
return np.concatenate(feats)
def hog(o, cell=8, bins=8):
gy, gx = np.gradient(o)
mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
out = np.zeros((N // cell, N // cell, bins))
for i in range(N // cell):
for j in range(N // cell):
sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
return (out / (np.linalg.norm(out) + 1e-6)).ravel()
def shape_features(o, s):
"""객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
sm = ndimage.uniform_filter(s, 3)
dark = np.abs(o - np.median(o))
m = (sm > 0.22) | (dark > 0.12)
lab, n = ndimage.label(m)
if n:
c = lab[10:22, 8:24]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
m = lab == k
p = m.sum(axis=0).astype(float)
cols = np.nonzero(p)[0]
if len(cols) == 0:
return np.zeros(8), o, s
thin = (p > 0) & (p <= 2)
left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
if left > right:
o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
cols = np.nonzero(p)[0]
body = np.nonzero(p > 3)[0]
d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
heavy = d.sum(axis=0)
centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s
def adv_features(o, s):
"""정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
shp, o, s = shape_features(o, s)
c = s[8:24, 4:28]
sm = ndimage.uniform_filter(s, 3)
m = sm > 0.2
dark = np.abs(o - np.median(o))
return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
[o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])
Otr, Str, Mtr = make_dataset(600, seed=3, scen=tuple(SCEN), classes=CLASSES + ["none"])
clf = RandomForestClassifier(n_estimators=100, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(Otr, Str, Mtr)]), [m["cls"] for m in Mtr])
def process(batch, budget_ms, prof):
"""배치 처리: 타일 검사 -> 전처리 -> 특징 -> 일괄 분류. 예산을 넘기면 정밀 분류를 미루고 '객체 유무'만 먼저 표시(기능 축소)."""
ok, quarantine, X = [], [], []
t_start = time.perf_counter()
for o, s, m in batch:
try:
if o.shape != (N, N) or s.shape != (N, N) or not np.isfinite(o).mean() > 0.5 or o.max() <= 0:
raise ValueError("손상 타일")
t0 = time.perf_counter(); po, ps, _ = preprocess(o, s, m); prof["전처리"] += time.perf_counter() - t0
if (time.perf_counter() - t_start) * 1000 > budget_ms:
ok.append(("객체유무", float(ndimage.uniform_filter(ps, 3).max() > 0.25))); continue
t0 = time.perf_counter(); X.append(adv_features(po, ps)); prof["특징"] += time.perf_counter() - t0
ok.append(("정밀", None))
except Exception as e:
quarantine.append((m["id"], str(e))) # 격리 보관(유실 아님)
if X:
t0 = time.perf_counter(); clf.predict(np.array(X)); prof["분류"] += time.perf_counter() - t0
return ok, quarantine
def workload(n, scen=("clean",), broken=0.0, seed=0):
O, S, M = make_dataset(n, seed=seed, scen=scen, classes=CLASSES + ["none"])
rng = np.random.default_rng(seed)
out = []
for o, s, m in zip(O, S, M):
r = rng.random()
if r < broken / 3:
o = o[:20] # 잘린 타일
elif r < 2 * broken / 3:
o = np.full_like(o, np.nan) # 전송 오류(전부 결측)
elif r < broken:
o = np.zeros_like(o) # 빈 타일
out.append((o, s, m))
return out
bad_wx = ("cloud", "low_res", "speckle", "camo_occl")
cases = [("평시 100장(기준)", 100, ("clean",), 0, 1e9), ("500% 증가(600장)", 600, ("clean",), 0, 1e9),
("600장 악기상·저해상도", 600, bad_wx, 0, 1e9), ("600장 + 손상 타일 5%", 600, ("clean",), 0.05, 1e9),
("600장 + 처리 예산 부족", 600, ("clean",), 0, 150)]
rows, prof = [], {"전처리": 0.0, "특징": 0.0, "분류": 0.0}
for i, (name, n, sc, br, budget) in enumerate(cases):
q = deque(workload(n, sc, br, seed=10 + i))
received, done, quar, degraded, down = len(q), 0, 0, 0, 0
t0 = time.perf_counter()
try:
while q:
chunk = [q.popleft() for _ in range(min(100, len(q)))]
ok, qu = process(chunk, budget, prof)
done += len(ok); quar += len(qu); degraded += sum(k == "객체유무" for k, _ in ok)
except Exception:
down = 1
dt = time.perf_counter() - t0
rows.append({"상황": name, "수신": received, "처리": done, "격리": quar, "유실": received - done - quar, "다운": down,
"처리량(장/초)": int(received / dt), "성능저하 모드": f"정밀분류 지연 {degraded}장" if degraded else "-"})
df = pd.DataFrame(rows)
print(df.to_string(index=False))
print(f"\n시스템 다운 {df['다운'].sum()}건, 데이터 유실률 {df['유실'].sum() / df['수신'].sum():.1%} (목표 0%), 손상 타일은 원본 보관 후 재요청")
# 메모리 안정성: 같은 배치(100장)를 5회 반복 처리하며 추적 메모리 변화 확인
batch = workload(100, seed=99)
tracemalloc.start(); mem = []
for _ in range(5):
process(batch, 1e9, {"전처리": 0, "특징": 0, "분류": 0}); mem.append(tracemalloc.get_traced_memory()[0] / 1e6)
tracemalloc.stop()
print("메모리(동일 배치 5회, MB):", [round(m, 2) for m in mem], f"-> 증가 {mem[-1] - mem[0]:.2f} 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, "처리량(장/초)"]
area_km2 = (N * STD_GSD / 1000) ** 2
print(f"처리 한계치 환산: 단일 코어 {thr}장/초 = 시간당 약 {thr * 3600 * area_km2:,.0f} km2 (타일 1장 {area_km2 * 1e6:.0f} m2, 표준 해상도 0.5m 기준)")
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("tiles / s")
plt.title("Tile batch stress test (red = degraded mode)"); plt.show()
In [ ]:
"""AAR(사후검토) 일지 형식으로 모의 위장 표적(decoy)이 섞인 판독 훈련 결과의 오판독 사례를 원인별로 분류하고, 특징 가중치(모의 표적 학습 추가)와 SAR-광학 일치성 임계값을 재보정해 종합 정확도 향상과 재발 여부를 회귀 시험으로 확인한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
[[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]
def conv_features(o):
"""간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
feats = []
for k in KERNELS:
r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
return np.concatenate(feats)
def hog(o, cell=8, bins=8):
gy, gx = np.gradient(o)
mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
out = np.zeros((N // cell, N // cell, bins))
for i in range(N // cell):
for j in range(N // cell):
sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
return (out / (np.linalg.norm(out) + 1e-6)).ravel()
def shape_features(o, s):
"""객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
sm = ndimage.uniform_filter(s, 3)
dark = np.abs(o - np.median(o))
m = (sm > 0.22) | (dark > 0.12)
lab, n = ndimage.label(m)
if n:
c = lab[10:22, 8:24]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
m = lab == k
p = m.sum(axis=0).astype(float)
cols = np.nonzero(p)[0]
if len(cols) == 0:
return np.zeros(8), o, s
thin = (p > 0) & (p <= 2)
left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
if left > right:
o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
cols = np.nonzero(p)[0]
body = np.nonzero(p > 3)[0]
d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
heavy = d.sum(axis=0)
centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s
def adv_features(o, s):
"""정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
shp, o, s = shape_features(o, s)
c = s[8:24, 4:28]
sm = ndimage.uniform_filter(s, 3)
m = sm > 0.2
dark = np.abs(o - np.median(o))
return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
[o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])
EX_SCEN = ("clean", "camo", "occl30", "low_sun", "cloud")
def featurize(O, S, M):
P = preprocess_batch(O, S, M)
return np.array([adv_features(o, s) for o, s, _ in P]), np.array([p[1][8:24, 4:28].mean() for p in P])
def train(with_decoy, seed=3):
cls = CLASSES + (["decoy"] if with_decoy else [])
O, S, M = make_dataset(250 * len(cls), seed=seed, scen=tuple(SCEN), classes=cls)
X, _ = featurize(O, S, M)
return RandomForestClassifier(n_estimators=150, min_samples_leaf=2, random_state=0, n_jobs=1).fit(X, [m["cls"] for m in M])
def decide(model, X, sar, tau):
"""분류 결과 + SAR-광학 일치성 규칙(선택): 전차·자주포로 분류했는데 SAR 금속 반사가 tau 미만이면 '모의 표적 의심'."""
p = model.predict(X).astype(object)
if tau is not None:
p[np.isin(p, ["tank", "spg"]) & (sar < tau)] = "decoy"
return p
def aar_log(y, p, sc, ids):
rows = []
for yi, pi, s, i in zip(y, p, sc, ids):
if yi == pi:
continue
if yi == "decoy":
cause = f"모의 표적을 {pi}(으)로 오판독"
elif pi == "decoy":
cause = f"실제 {yi}를 모의 표적으로 오판"
else:
cause = f"{yi}->{pi} 기종 혼동({s})"
rows.append({"원인": cause, "영상": i})
return pd.DataFrame(rows, columns=["원인", "영상"])
# 모의 판독 훈련 자료(모의 위장 표적 20% 포함): AAR 대상 500장 + 회귀 시험 500장
Oa, Sa, Ma = make_dataset(500, seed=21, scen=EX_SCEN, classes=CLASSES + ["decoy"])
Ob, Sb, Mb = make_dataset(500, seed=22, scen=EX_SCEN, classes=CLASSES + ["decoy"])
Xa, sar_a = featurize(Oa, Sa, Ma); ya, sca, ida = np.array([m["cls"] for m in Ma]), [m["scen"] for m in Ma], [m["id"] for m in Ma]
Xb, sar_b = featurize(Ob, Sb, Mb); yb = np.array([m["cls"] for m in Mb])
# 1) 초기 설정: 모의 표적을 학습하지 않은 모델, 일치성 규칙 없음
m0 = train(False)
p0 = decide(m0, Xa, sar_a, None)
acc0 = (p0 == ya).mean()
log0 = aar_log(ya, p0, sca, ida)
print(f"[1. 무엇이 일어났나] 초기 종합 정확도 {acc0:.3f}, 오판독 {len(log0)}건")
print(log0["원인"].value_counts().head(6).to_string())
print(f"\n[2. 왜 일어났나] 모의 표적 SAR 반사 평균 {sar_a[ya == 'decoy'].mean():.3f} vs 실제 전차·자주포 {sar_a[np.isin(ya, ['tank', 'spg'])].mean():.3f}"
" -> 광학 윤곽은 전차와 같지만 금속 반사가 약해, 모의 표적을 모르는 모델은 반사가 약한 민간 차량으로 몰아 버린다.")
# 2) 재보정: (모의 표적 학습 여부) x (SAR 일치성 임계값) 격자 탐색
m1 = train(True)
grid = []
for name, mdl in [("모의표적 미학습", m0), ("모의표적 학습", m1)]:
for tau in [None, 0.12, 0.15, 0.18, 0.22]:
p = decide(mdl, Xa, sar_a, tau)
mil = np.isin(ya, list(MILITARY))
grid.append({"모델": name, "SAR 임계값": "-" if tau is None else tau, "종합정확도": round(float((p == ya).mean()), 3),
"모의표적 적발률": round(float((p[ya == "decoy"] == "decoy").mean()), 3),
"실제 장비를 모의로 오판": int((mil & (p == "decoy")).sum())})
g = pd.DataFrame(grid).sort_values(["종합정확도", "실제 장비를 모의로 오판"], ascending=[False, True])
print("\n[3. 무엇을 바꾸나] 격자 탐색 상위 5개")
print(g.head().to_string(index=False))
best = g.iloc[0]
mdl = m1 if best["모델"] == "모의표적 학습" else m0
tau = None if best["SAR 임계값"] == "-" else float(best["SAR 임계값"])
p1 = decide(mdl, Xa, sar_a, tau)
acc1 = (p1 == ya).mean()
print(f"\n재보정값: {best['모델']}, SAR 일치성 임계값 {best['SAR 임계값']}")
print(f"종합 정확도 {acc0:.3f} -> {acc1:.3f} ({(acc1 - acc0) * 100:+.1f}%p, 목표 +10%p: {'달성' if acc1 - acc0 >= 0.1 else '미달'})")
b0, b1 = (decide(m0, Xb, sar_b, None) == yb).mean(), (decide(mdl, Xb, sar_b, tau) == yb).mean()
print(f"회귀 시험(별도 500장): {b0:.3f} -> {b1:.3f}")
log1 = aar_log(ya, p1, sca, ida)
rep = set(log0["영상"]) & set(log1["영상"])
print(f"재보정 후 남은 오판독 {len(log1)}건, 그중 재발 {len(rep)}건:", log1[log1["영상"].isin(rep)]["원인"].value_counts().head(3).to_dict())
print("SAR 일치성 규칙은 모의 표적 학습 후 추가 이득이 없고 임계값을 높이면 실제 장비를 모의 표적으로 오판하므로 적용하지 않는다(-).")
print("[4. 남은 과제] 재발 건은 위장막 아래 기종 혼동으로, 다른 촬영 회차 교차 판독(미션 7)과 지속 학습 대상 등록으로 대응한다.")
tuned = f"model:\n train_with_decoy: {best['모델'] == '모의표적 학습'}\n rf_trees: 150\nrules:\n sar_consistency_tau: {best['SAR 임계값']}\nconfidence:\n high: 0.9\n medium: 0.6"
print("\n[재보정 파라미터(isr01_tuned_params.yaml 내용)]\n" + tuned)
plt.figure(figsize=(7, 4))
plt.hist([sar_a[np.isin(ya, ["tank", "spg"])], sar_a[ya == "decoy"]], bins=25, label=["tank/spg", "decoy"], alpha=.8)
if tau:
plt.axvline(tau, ls="--", c="k", label=f"tau={tau}")
plt.xlabel("mean SAR return (center)"); plt.ylabel("chips"); plt.legend(); plt.title("SAR-optical consistency for decoy screening")
plt.show()
In [ ]:
"""메타데이터 검증·전처리·정밀 분류(모의 표적 포함)·신뢰도 엔진·영상분석관 판독 대기열·분석 화면까지 공개 위성영상 장비 객체분류 시스템 전체를 한 파일로 통합 실행하고, 자체 점검(기능 커버리지 포함)과 전략 ISR 평가 브리핑 요약을 출력한다."""
import hashlib
import json
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5 # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
"tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
"spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
"truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
"civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
"decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
"none": [], # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12} # 금속 반사 세기(가상)
def _mask(cls, angle, scale, squash=1.0):
"""회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
yy = yy / squash
a = np.radians(angle)
u = (xx * np.cos(a) + yy * np.sin(a)) / scale
v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]
def make_chip(cls, rng, cond=None):
"""chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
scale = STD_GSD / gsd
angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
squash = np.cos(np.radians(off))
sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
sun_az = rng.uniform(0, 360)
bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
parts = _mask(cls, angle, scale, squash)
body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
opt = bg.copy()
# 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
opt[sh & ~body] -= 35 if cls != "decoy" else 25
if cls == "civil":
base = rng.choice([200.0, 60.0]) # 민간 차량 색상: 밝거나 어두움
opt[body] = base + rng.normal(0, 8, body.sum())
roof = parts[1][0]
opt[roof] = base * 0.8
else:
ref = bg[body].mean() if body.any() else bg.mean()
for m, d in parts:
opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
# 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
if c["camo"] > 0 and cls != "civil":
net = ndimage.binary_dilation(body, iterations=2)
tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
# 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
if c["occl"] > 0:
trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
near = ndimage.binary_dilation(body, iterations=3)
cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
opt[cover] = 70 + rng.normal(0, 6, cover.sum())
# 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
cmask = np.zeros((N, N), bool)
if c["cloud"] > 0:
blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
cmask = blob > np.quantile(blob, 1 - c["cloud"])
opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
# SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
sar = np.full((N, N), 0.12)
sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
edge = body & ~ndimage.binary_erosion(body)
sar[edge] *= 1.3 if cls in MILITARY else 1.0
sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
if c["low_res"]: # 저해상도: 정보 손실(블러)
down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
"sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
"cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
return np.clip(opt, 0, 255), sar, meta
SCEN = { # 은폐·촬영 조건 시나리오 10종
"clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
"low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
"speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}
def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
rng = np.random.default_rng(seed)
O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
for i in range(n):
cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
m.update({"id": f"IMG{i:04d}", "scen": sc})
M.append(m)
return O, S, M
# ---------------------------------------------------------------------------
# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
"""중앙 기준으로 32x32에 맞게 자르거나 채운다."""
out = np.full((N, N), fill, float)
h, w = a.shape
ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
hh, ww = min(N, h), min(N, w)
out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
return out
def cloud_mask(opt, sm, min_area=15):
"""광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
bright = opt > 185
if bright.sum() < min_area:
return np.zeros_like(bright)
local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
cand = bright & (local_sd < 12) & (sm < 0.3)
lab, n = ndimage.label(cand)
if n == 0:
return cand
sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)
def orientation(opt, sm):
"""SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
lab, n = ndimage.label(m)
if n == 0:
return 0.0, m
c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
m = lab == k
yy, xx = np.nonzero(m)
if len(yy) < 5:
return 0.0, m
dy, dx = yy - yy.mean(), xx - xx.mean() # 2차 모멘트로 주축 각도 계산
return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m
_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2
def warp(stack, zy, zx, ang):
"""여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
a = np.radians(ang)[:, None, None]
zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001) # 출력 좌표 -> 원본 좌표
sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
y0, x0 = sy.astype(int), sx.astype(int)
fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
flat = stack.reshape(len(stack), stack.shape[1], -1)
lin = (y0 * N + x0).reshape(len(stack), 1, -1)
g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
return out.reshape(stack.shape)
def preprocess_batch(O, S, M):
"""여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
if np.isnan(src).any(): # 0) 센서 결측 픽셀 보정
mean = np.nanmean(src, axis=(2, 3), keepdims=True)
src = np.where(np.isnan(src), mean, src)
zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
zx = np.array([m["gsd"] / STD_GSD for m in M])
out = warp(src, zy, zx, np.zeros(len(M)))
sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
angs, cloud = np.zeros(len(M)), np.zeros(len(M))
for i in range(len(M)):
cm = cloud_mask(out[i, 0], sm[i])
if cm.any():
fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
out[i, 0][cm] = fill
back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
src[i, 0][back] = fill # 원본의 같은 구름 영역도 채움
cloud[i] = cm.mean()
angs[i] = orientation(out[i, 0], sm[i])[0]
rot = np.abs(angs) > 3 # 3도 이하는 회전 생략
if rot.any():
out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]
def preprocess(opt, sar, meta):
"""한 장 전처리(preprocess_batch와 같은 처리)."""
return preprocess_batch([opt], [sar], [meta])[0]
# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
[[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]
def conv_features(o):
"""간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
feats = []
for k in KERNELS:
r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
return np.concatenate(feats)
def hog(o, cell=8, bins=8):
gy, gx = np.gradient(o)
mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
out = np.zeros((N // cell, N // cell, bins))
for i in range(N // cell):
for j in range(N // cell):
sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
return (out / (np.linalg.norm(out) + 1e-6)).ravel()
def shape_features(o, s):
"""객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
sm = ndimage.uniform_filter(s, 3)
dark = np.abs(o - np.median(o))
m = (sm > 0.22) | (dark > 0.12)
lab, n = ndimage.label(m)
if n:
c = lab[10:22, 8:24]
k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
m = lab == k
p = m.sum(axis=0).astype(float)
cols = np.nonzero(p)[0]
if len(cols) == 0:
return np.zeros(8), o, s
thin = (p > 0) & (p <= 2)
left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
if left > right:
o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
cols = np.nonzero(p)[0]
body = np.nonzero(p > 3)[0]
d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
heavy = d.sum(axis=0)
centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s
def adv_features(o, s):
"""정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
shp, o, s = shape_features(o, s)
c = s[8:24, 4:28]
sm = ndimage.uniform_filter(s, 3)
m = sm > 0.2
dark = np.abs(o - np.median(o))
return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
[o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])
CONFIG = {"train_with_decoy": True, "high": 0.9, "medium": 0.6, "hvt_types": ["spg"]} # 미션 9 재보정값
CALLED = set()
def cover(fn):
"""기능 커버리지 측정용 표시: 통합 실행 중 호출된 기능을 기록한다."""
def wrap(*a, **k):
CALLED.add(fn.__name__)
return fn(*a, **k)
wrap.__name__ = fn.__name__
return wrap
preprocess, adv_features, shape_features, hog, cloud_mask, orientation = map(cover, (preprocess, adv_features, shape_features, hog, cloud_mask, orientation))
@cover
def validate_tile(o, s, m):
errs = []
if o.shape != (N, N) or s.shape != (N, N):
errs.append("타일 크기 오류")
if not 0.1 <= m["gsd"] <= 30:
errs.append("gsd 범위 오류")
if not np.isfinite(o).all():
errs.append("결측 픽셀")
return errs
@cover
def quality(P, prob):
srt = np.sort(prob, axis=1)
ent = -(prob * np.log(np.clip(prob, 1e-9, 1))).sum(1) / np.log(prob.shape[1])
return np.c_[srt[:, -1], srt[:, -1] - srt[:, -2], ent, [p[2] for p in P], [p[1][8:24, 4:28].mean() for p in P]]
@cover
def grade(c):
return "High" if c >= CONFIG["high"] else ("Medium" if c >= CONFIG["medium"] else "Low")
class Queue:
def __init__(self):
self.items, self.log = [], []
@cover
def submit(self, iid, pred, g, conf):
prev = self.log[-1]["hash"] if self.log else "GENESIS"
body = {"event": "SUBMIT", "id": iid, "pred": pred, "grade": g, "conf": round(conf, 3), "hvt": pred in CONFIG["hvt_types"]}
self.log.append({**body, "prev": prev, "hash": hashlib.sha256((prev + json.dumps(body, 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, sort_keys=True)).encode()).hexdigest():
return False
prev = r["hash"]
return True
TARGETS = {"validate_tile", "preprocess", "cloud_mask", "orientation", "adv_features", "shape_features", "hog", "quality", "grade", "submit", "verify", "dashboard"}
LABELS = CLASSES + ["decoy", "none"]
# ---------------- 1. 학습(분류 모델 + 신뢰도 가중치) ----------------
t_start = time.perf_counter()
O, S, M = make_dataset(200 * len(LABELS), seed=3, scen=tuple(SCEN), classes=LABELS)
clf = RandomForestClassifier(n_estimators=150, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(O, S, M)]), [m["cls"] for m in M])
Oc, Sc, Mc = make_dataset(480, seed=6, scen=tuple(SCEN), classes=LABELS)
Pc = preprocess_batch(Oc, Sc, Mc)
prc = clf.predict_proba(np.array([adv_features(o, s) for o, s, _ in Pc]))
conf_model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(quality(Pc, prc), clf.classes_[prc.argmax(1)] == np.array([m["cls"] for m in Mc]))
# ---------------- 2. 운용 자료(가상 판독 지역 300타일) 처리 ----------------
Oo, So, Mo = make_dataset(300, seed=31, scen=("clean", "camo", "occl30", "cloud", "low_sun", "low_res"), classes=LABELS)
ops = list(zip(Oo, So, Mo))
ops[5] = (Oo[5][:20], So[5], Mo[5]) # 손상 타일 1건(잘림)
ops[9] = (Oo[9], So[9], {**Mo[9], "gsd": 0.0}) # 메타데이터 오류 1건
queue, rejected, res, ms = Queue(), {}, [], []
for o, s, m in ops:
errs = validate_tile(o, s, m)
if errs:
rejected[m["id"]] = errs; continue
t0 = time.perf_counter()
P = [preprocess(o, s, m)]
pr = clf.predict_proba([adv_features(P[0][0], P[0][1])])
pred = clf.classes_[pr.argmax()]
conf = float(conf_model.predict_proba(quality(P, pr))[0, 1])
g = grade(conf)
if pred != "none":
queue.submit(m["id"], pred, g, conf)
ms.append((time.perf_counter() - t0) * 1000)
res.append({"id": m["id"], "y": m["cls"], "p": pred, "grade": g, "conf": conf, "scen": m["scen"]})
R = pd.DataFrame(res)
acc = (R.y == R.p).mean()
milR = R[R.y.isin(list(MILITARY))]
rec = (milR.y == milR.p).mean()
decoy_catch = (R[R.y == "decoy"].p == "decoy").mean()
hi = R[(R.grade == "High") & R.p.isin(list(MILITARY))]
hi_err = (hi.y != hi.p).mean() if len(hi) else 0.0
# ---------------- 3. 분석 화면 ----------------
@cover
def dashboard():
plt.style.use("dark_background")
fig, ax = plt.subplots(2, 2, figsize=(12, 9))
T = 10
mosaic = np.block([[np.asarray(Oo[r * T + c]) for c in range(T)] for r in range(T)])
ax[0, 0].imshow(mosaic, cmap="gray")
for k, r in R[R.id.isin([f"IMG{i:04d}" for i in range(100)])].iterrows():
i = int(r.id[3:])
if r.p != "none":
ax[0, 0].scatter(i % T * N + 16, i // T * N + 16, s=60, facecolors="none",
edgecolors="red" if r.p in MILITARY else ("yellow" if r.p == "decoy" else "deepskyblue"))
ax[0, 0].set_title("First 100 tiles (red = military, yellow = decoy)"); ax[0, 0].axis("off")
vc = R[R.p != "none"].p.value_counts()
ax[0, 1].bar(vc.index.tolist(), vc.values, color="deepskyblue"); ax[0, 1].set_title("Detections by type")
gc = R[R.p != "none"].grade.value_counts().reindex(["High", "Medium", "Low"]).fillna(0)
ax[1, 0].bar(gc.index.tolist(), gc.values, color=["red", "orange", "gray"]); ax[1, 0].set_title("Review queue by grade (all pending)")
ax[1, 1].axis("off")
ax[1, 1].text(0, .9, f"Accuracy {acc:.1%}\nMilitary rec. {rec:.1%}\nDecoy caught {decoy_catch:.1%}\nQueue {len(queue.items)} pending\n"
f"Latency {np.mean(ms):.1f} ms/tile", family="monospace", fontsize=13, va="top")
fig.suptitle("ISR01 integrated satellite object classification - synthetic data only"); plt.tight_layout(); plt.show()
dashboard()
# ---------------- 4. 자체 점검(통합 테스트) ----------------
log_ok = queue.verify()
cov = len(CALLED & TARGETS) / len(TARGETS)
chk = {
"손상 타일·메타데이터 오류 차단": set(rejected) == {"IMG0005", "IMG0009"},
"처리 타일 298장": len(res) == 298,
"종합 정확도 90% 이상": acc >= 0.9,
"군용 기종 식별 재현율 85% 이상": rec >= 0.85,
"모의 표적 적발률 80% 이상": decoy_catch >= 0.8,
"자동 승인 0건(전부 판독 대기)": all(i["status"] == "PENDING" for i in queue.items),
"식별 건 전량 대기열 등록": len(queue.items) == int((R.p != "none").sum()),
"감사 로그 무결성": log_ok,
"타일당 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 '아니오'}")
# ---------------- 5. 전략 ISR 평가 브리핑 요약 ----------------
print(f"""
[위성영상 장비 객체분류 평가 브리핑 요약 - 교육용가상자료]
1. 처리: 타일 {len(ops)}장 중 {len(res)}장 판독, 손상·오류 {len(rejected)}건 차단({', '.join(rejected)})
2. 성능: 종합 정확도 {acc:.1%}, 군용 기종 식별 재현율 {rec:.1%}, 모의 표적 적발 {decoy_catch:.1%}, 타일당 {np.mean(ms):.1f} ms, 전체 {time.perf_counter() - t_start:.0f} s
3. 판독 대기: {len(queue.items)}건 (High {int((R[R.p != 'none'].grade == 'High').sum())} / Medium {int((R[R.p != 'none'].grade == 'Medium').sum())} / Low {int((R[R.p != 'none'].grade == 'Low').sum())}), HVT 후보 {sum(i['hvt'] for i in queue.items)}건, 자동 승인 0건
4. 신뢰도: High 등급 군용 식별 {len(hi)}건 중 오분류 {hi_err:.1%} -> 등급은 판독 순서 결정용이며 최종 판독은 영상분석관이 한다
5. 한계: 합성 32x32 영상 기준 결과로, 실제 위성영상·장비에는 재학습과 재검증이 필요하다. 인접 체계 연동은 판독 승인 결과만 표준 필드로 전달한다.""")