환경 통계에서 Kriging과 기후 시계열을 결합하는 방법

관측소·위성·재분석 데이터는 관측 방식과 해상도가 다르고, 결측도 피하기 어렵다. 이를 연속된 공간 격자와 시간 예측으로 연결하려면 공간 자기상관을 다루는 Kriging과 기후 시계열 분석을 함께 설계해야 한다.

공간의 빈 곳과 시간의 변화를 함께 모델링하기

Kriging은 공간 자기상관을 이용하는 최적 선형 불편 예측기(BLUE)다. 반변량함수(variogram)로 공분산 구조를 추정한 뒤, 값이 관측되지 않은 위치를 계산한다. 이 과정에서는 평균과 공분산이 변하지 않는다는 정지성 가정을 검토해야 하며, 등방성 또는 이방성도 함께 고려한다. Nugget–Sill–Range 매개변수는 각각 미시적 잡음, 총분산, 공간 영향반경을 나타낸다.

기후 시계열은 추세, 계절성, 주기, 구조적 변화가 겹친 비정상(non-stationary) 시계열이다. 월별·계절별 반복과 저주파 경향이 공존하므로 ARIMA/SARIMA, 상태공간(Kalman), STL 분해 등을 이용해 신호와 잡음을 분리하고 예측한다.

분석 전에 맞춰야 할 데이터 기준

입력은 지상관측소의 포인트 데이터, 위성·레이더의 그리드 데이터, 재분석(ERA5 등) 데이터로 구성될 수 있다. 분석 전에는 좌표계를 통일하고 시간대와 캘린더를 정규화해야 한다.

결측치를 식별하고 단위를 변환하며, 스파이크나 평탄화 같은 이상치를 제거한다. 관측 품질 지표를 태깅해두면 이후 가중치 적용에도 활용할 수 있다. 시간 데이터에서는 결측 플래그와 윤일·서머타임 정규화가 필요하고, 공간 데이터에서는 위경도를 투영미터 좌표로 옮겨 거리 단위를 일관되게 맞춘다.

반변량함수와 보간 모델의 선택

경험적 반변량함수를 계산한 뒤 구형·지수·가우시안 모델을 적합한다. 방향에 따라 상관 구조가 달라지는 경우에는 방향별 반변량 분석을 수행하고, 좌표 회전이나 축척 변환을 적용한다.

Nugget/Sill/Range 추정은 잔차 반변량의 백색성, Cressie-Hawkins 견고 추정, CV 지표인 RMSE와 MSPE로 검증한다. 보간 단계에서는 평균이 상수라는 가정에 맞는 Ordinary Kriging과 드리프트·공변량을 포함하는 Universal Kriging 중 하나를 선택한다. 지형고도, 토지피복, 거리기반 기상후보변량도 함께 병합할 수 있다.

LOSO와 공간 블록 CV는 보간 성능의 낙관적 편향을 줄이는 데 쓰인다. 산출물은 추정값 격자뿐 아니라 Kriging 분산을 담은 불확실성 지도까지 포함해야 한다.

공간 격자와 시간 예측이 만나는 흐름

⤢✕

시계열은 STL/LOESS로 추세·계절·잔차를 분해한 뒤 SARIMA·Prophet·State-Space 모델로 단기·중기 예측을 수행한다. 고정 계절과 이동 계절성이 섞여 있다면 혼합 접근이 필요하다. 롤링 윈도우 CV로 예측 안정성을 평가하고, 예측구간(PIs)을 제공해 의사결정 리스크를 관리한다.

공간 성능은 RMSE/MAE, CRPS, 평균 Kriging 분산과 오차의 정합성(coverage)으로 평가한다. 시간 성능은 sMAPE, RMSE, 칼만 필터 잔차의 백색성·정규성 테스트로 확인한다. 운영 단계에서는 계산비용과 지연(latency)도 지속적으로 모니터링한다.

보간 방식별 운영 특성

방법	성능(정확도)	확장성(대용량)	일관성(이론)	안정성(잡음 민감)	운영 편의
IDW	중간, 구조 단순	높음, O(n log n)	낮음, 경험적	잡음 민감	매우 용이
Spline	곡면 부드러움 우수	중간	중간	외삽 불안정	보통
Kriging(Ordinary)	높음, 불확실성 제공	중간, O(n^3)	높음, 확률론 기반	견고(variogram 적정 시)	보통
Kriging(Universal)	매우 높음(공변량 포함)	중간	높음	모델링 복잡	보통-어려움
랜덤포레스트+좌표	높음(비선형)	높음(분산처리)	중간	과적합 주의	용이-보통

Kriging은 공분산 구조가 적절히 추정될 때 강점을 보인다. 대용량 데이터에서는 Nearest Neighbor GP나 tile+overlap 같은 근사 기법이 필요하다.

환경 데이터에 적용하는 장면

PM2.5 시공간 격자 예보

도시권 관측소 PM2.5에 고도·도로밀도 공변량과 재분석 바람장을 결합할 수 있다. QC를 거친 뒤 EPSG:5186으로 재투영하고, 방향성 반변량을 추정해 Universal Kriging을 적용한다.

결과물은 1 km 격자 PM2.5 지도와 Kriging 분산, 24~72시간 SARIMA 예측이다. 고정 스테이션이 누락되면 최근 kNN+Kalman 보정으로 대체하고, 공변량 최신화는 캐싱하며, 타일링으로 지연을 5분 이내로 유지한다.

일강수 그리드와 홍수 위험도

레이더 강수의 편향을 보정한 다음 관측소 데이터를 Kriging with external drift로 융합한다. 누적 강수 임계치 초과 확률은 Kriging 분산과 시계열 예측구간을 결합해 산출한다.

기온 정상기후와 파생 지수

30년 기준(예: 1991–2020)의 월 평균을 시공간 모델로 계산해 정상기후(Climatology)를 구축한다. 이 기반 위에서 일교차, 열지수, 폭염일수 같은 파생 지수의 공간 분포를 제공할 수 있다.

재현 가능한 운영 체계로 연결하기

공간 모델링에서는 Cressie 기반의 반변량 견고 추정과 로즈다이어그램을 이용한 이방성 테스트를 적용한다. 모델은 짧은 상관에는 구형, 매끄러운 표면에는 가우시안, 장거리 꼬리에는 지수 모델을 검토한다.

보간 결과는 LOSO/블록 CV와 RMSE·CRPS·Coverage로 보고한다. 대규모 처리에는 NNGP와 SKI(KISS-GP)를 고려할 수 있다. 시간 모델은 STL 분해 후 SARIMA(p,d,q)(P,D,Q)s를 탐색하고 롤링 CV로 안정성을 확인한다. 결측과 가변 간격은 상태공간 모델로 처리한다.

모델과 데이터는 MLflow/DVC로 버전 관리하고, 도커 기반 재현성 환경을 유지한다. SLO로 RMSE와 지연을 모니터링하며, 불확실성 시각화를 정책 의사결정과 연결한다.

Python으로 확인하는 최소 예제

전제: Python 3.10+, numpy 1.26, pandas 2.x, geopandas 0.14, pykrige 1.7, statsmodels 0.14, scikit-learn 1.4, pyproj, shapely, rasterio(옵션)

# pip install numpy pandas geopandas pykrige statsmodels scikit-learn pyproj shapely
import numpy as np, pandas as pd
from pykrige.ok import OrdinaryKriging
from sklearn.model_selection import KFold
import statsmodels.api as sm

# 1) 예시 데이터(관측소 포인트)
np.random.seed(42)
n = 60
lon = np.random.uniform(126.5, 127.5, n)
lat = np.random.uniform(37.0, 38.0, n)
# 숨은 공간장 + 잡음
true_field = lambda x,y: 20 + 3*np.sin((x-126.5)*3) + 2*np.cos((y-37.0)*4)
z = true_field(lon, lat) + np.random.normal(0, 0.5, n)

# 2) Kriging(경도/위도 직접 사용은 근사, 실제는 투영 좌표 사용 권장)
OK = OrdinaryKriging(
    lon, lat, z,
    variogram_model="spherical",
    variogram_parameters=None,  # 자동 추정
    enable_plotting=False
)
grid_lon = np.linspace(126.5, 127.5, 50)
grid_lat = np.linspace(37.0, 38.0, 50)
zk, ss = OK.execute("grid", grid_lon, grid_lat)  # zk: 추정, ss: Kriging 분산

# 3) 간단 교차검증(LOSO 근사)
kf = KFold(n_splits=5, shuffle=True, random_state=0)
errs = []
for tr, te in kf.split(lon):
    ok = OrdinaryKriging(lon[tr], lat[tr], z[tr], variogram_model="spherical")
    zk_te, _ = ok.execute("points", lon[te], lat[te])
    errs.extend(z[te] - zk_te.data)
rmse = np.sqrt(np.mean(np.array(errs)**2))
print("Spatial RMSE:", round(rmse, 3))

# 4) 시계열: 월별 기온 이상치 예시
idx = pd.date_range("2015-01-01", periods=120, freq="MS")
signal = 0.02*np.arange(120) + 1.5*np.sin(2*np.pi*idx.month/12)
y = signal + np.random.normal(0, 0.4, 120)
# STL 분해
stl = sm.tsa.STL(y, period=12).fit()
# SARIMA(간단 모형)
model = sm.tsa.SARIMAX(y, order=(1,1,1), seasonal_order=(1,1,1,12)).fit(disp=False)
fc = model.get_forecast(12)
pred = fc.predicted_mean
pi = fc.conf_int()
print("Forecast next 12 months head:\n", pred.head())

실제 운영에서는 위경도를 EPSG:5186, 3857 등의 투영 좌표로 변환한 뒤 반변량과 Kriging을 수행하는 편이 적합하다. PyKrige의 자동 variogram 추정은 단순화 옵션이므로, 정밀 업무에서는 scikit-gstat 등으로 견고 적합한 후 매개변수를 명시한다. 대용량 격자는 grid tiling과 경계 오버랩을 사용해 봉합 아티팩트를 줄인다.

정확도와 해석 가능성 사이의 선택

Kriging은 IDW보다 정확도와 불확실성 제공 측면에서 우수하지만, 계산복잡도와 모델링 부담이 커진다. 관측이 희소하거나 잡음이 큰 경우 nugget을 확대하고 견고 추정을 적용해 안정성을 확보한다.

Universal Kriging과 랜덤포레스트·GBDT는 모두 공변량을 사용할 수 있다. 해석가능성과 불확실성 정량화는 Kriging이, 비선형성과 상호작용은 ML이 강점을 가진다. ML로 드리프트를 추정하고 잔차를 Kriging으로 공간 보정하는 결합 전략도 가능하다.

SARIMA는 계절성이 뚜렷하고 설명변수가 적을 때 효율적이다. 상태공간·칼만 접근은 결측, 불규칙 간격, 동적 회귀 처리에 강점이 있다. 어느 모형을 택하더라도 예측 신뢰구간의 캘리브레이션과 롤링 CV가 필요하다.

타일링·근사 GP를 도입하면 공간 RMSE는 IDW 대비 1030% 감소하고, coverage는 90%±5%를 달성할 수 있으며, 운영 지연 시간은 3060% 단축될 수 있다. 불확실성 기반 의사결정 체계는 모델 설명가능성을 높이고, 데이터 거버넌스와 재현성은 규제 대응 및 협업 효율을 강화한다.