다운스케일링과 기후 모델 에뮬레이션

다운스케일링

왜 필요한가. 저해상도 격자 칸의 기온은 칸 안의 산꼭대기와 골짜기를 평균한 값이다. 해발 1500 m 봉우리와 200 m 골짜기는 기온 감률(고도 1 km마다 약 6.5 °C)만으로도 8 °C 넘게 다르다. 칸 평균을 그대로 쓰면 봉우리의 눈과 골짜기의 서리를 모두 놓친다.

정의. 다운스케일링은 저해상도 기후 정보 $\mathbf{x}$(그리고 지형 같은 고해상도 고정 정보 $\mathbf{s}$)에서 고해상도 기후 정보 $\mathbf{y}$를 얻는 일이다. 두 갈래가 있다.

• 동역학적 다운스케일링(dynamical downscaling): 전 지구 모델의 출력을 경계 조건으로 지역 기후 모델(regional climate model)을 촘촘한 격자로 돌린다. 물리적으로 일관되지만 비싸고, 지역 모델 자체의 편향이 더해진다.
• 통계적 다운스케일링(statistical downscaling): 저해상도와 고해상도 사이의 통계적 관계를 자료에서 학습한다. 싸고 빠르지만, 학습한 관계가 미래에도 성립한다는 가정에 기댄다.

통계적 방법은 다시 무엇으로 학습하느냐에 따라 나뉜다. 관측 기반의 대규모 변수(재분석)와 지역 관측 사이의 관계를 배워 기후 모델 출력에 적용하는 방식을 perfect prognosis라 하고, 기후 모델 출력과 관측 사이의 관계를 직접 배우는 방식을 MOS(model output statistics)라 한다. 전자는 기후 모델의 대규모 변수가 정확하다고 가정하고, 후자는 모델의 편향까지 함께 보정한다. 편향 보정(bias correction)은 MOS의 가장 단순한 형태다.

초해상도와의 차이. 영상 초해상도(super-resolution)와 문제의 모양은 같지만 두 가지가 다르다. 첫째, 저해상도 칸 하나에는 서로 다른 무수한 고해상도 장이 대응한다. 지형이 정하는 부분은 예측할 수 있지만, 날마다 달라지는 국지 변동은 저해상도 값으로 정해지지 않는다. 그래서 목표는 한 장의 영상이 아니라 조건부 분포 $p(\mathbf{y} \mid \mathbf{x}, \mathbf{s})$다. 둘째, 기후 모델의 저해상도 출력은 관측을 성기게 만든 것이 아니라 편향이 있는 별개의 자료다.

회귀로 하는 다운스케일링. 산악 단면 하나로 시작한다. 2 km 간격 64칸의 고해상도 기온은 기온 감률, 대규모 기온 편차(고지대일수록 조금 크게 반응), 이웃 칸끼리 상관된 국지 변동의 합이다. 저해상도는 고해상도 8칸의 평균이다. 고해상도 칸마다 "자기가 속한 저해상도 값"에 대한 1차 회귀를 학습한다. 칸마다 다른 절편이 그 칸의 고도 효과를 배운다.

import numpy as np

rng = np.random.default_rng(0)
n_day, n_fine, blk = 4000, 64, 8                      # 일 수, 고해상도 칸 수, 저해상도 칸 하나에 든 고해상도 칸 수
x_km = np.arange(n_fine) * 2.0                        # (64,) 2 km 간격 단면
elev_m = 200 + 1400 * np.abs(np.sin(x_km / 9)) ** 3 + 600 * (x_km > 64)   # (64,) 봉우리와 골짜기, 동쪽 고원
large = np.zeros(n_day)                               # 대규모 기온 편차: 하루 간격 AR(1)
for t in range(1, n_day):
    large[t] = 0.8 * large[t - 1] + rng.normal(0, 2.0)
eps = rng.normal(0, 1, (n_day, n_fine))
local = np.empty_like(eps); local[:, 0] = eps[:, 0]
for j in range(1, n_fine):                            # 국지 변동: 이웃 칸끼리 상관(0.8)된 잡음
    local[:, j] = 0.8 * local[:, j - 1] + 0.6 * eps[:, j]
fine_c = 15 - 0.0065 * elev_m + large[:, None] * (1 + elev_m / 3000) + 1.5 * local   # (n_day, 64)
coarse_c = fine_c.reshape(n_day, -1, blk).mean(-1)   # (n_day, 8) 저해상도 = 칸 평균

tr, te = slice(0, 3000), slice(3000, None)            # 앞 3000일로 학습, 나머지 1000일로 시험
X = np.repeat(coarse_c, blk, axis=1)                  # (n_day, 64) 각 고해상도 칸이 속한 저해상도 값
coef = [np.polyfit(X[tr, j], fine_c[tr, j], 1) for j in range(n_fine)]      # 칸마다 1차 회귀
pred = np.stack([np.polyval(coef[j], X[te, j]) for j in range(n_fine)], 1)  # (1000, 64)
rmse = lambda a: np.sqrt(((a - fine_c[te]) ** 2).mean())
print(f"저해상도 값 그대로  RMSE {rmse(X[te]):.2f} °C")
print(f"칸별 회귀           RMSE {rmse(pred):.2f} °C")

저해상도 값 그대로  RMSE 3.02 °C
칸별 회귀           RMSE 0.95 °C

저해상도 값을 그대로 쓰면 고도 차이를 모두 놓쳐 오차가 3 °C다. 칸별 회귀는 고도 효과를 배워 오차를 3분의 1로 줄인다. 남은 0.95 °C는 대부분 저해상도 값으로 정해지지 않는 국지 변동이다. 이 숫자만 보면 회귀가 좋은 다운스케일링으로 보인다. 다음 절에서 이 결론이 왜 불완전한지 본다.

편향 보정과 quantile mapping. 기후 모델 출력은 관측과 분포가 다르다. 강수에서 흔한 편향은 약한 비가 너무 자주 오고(drizzle) 강한 비는 너무 약한 것이다. Quantile mapping(QM)은 모델 값이 과거 기간 모델 분포에서 차지하는 분위를 구하고, 같은 분위의 관측 값으로 바꾼다.

$\hat F_\text{mod}$와 $\hat F_\text{obs}$는 과거(보정) 기간의 모델과 관측의 경험적 누적분포함수다. 과거 기간에서는 정의상 보정한 모델의 분포가 관측과 같아진다. 문제는 이 변환을 미래에 그대로 쓸 때다. 강수 강도가 30% 늘어나는 미래 모델 자료에 적용해 본다.

def daily_precip(g, n, p_wet, shape, scale_mm):       # 강수일 확률 × 감마 분포 강수량
    return np.where(g.random(n) < p_wet, g.gamma(shape, scale_mm, n), 0.0)

g = np.random.default_rng(2)
n = 30 * 365                                          # 30년치 일강수 [mm/day]
obs_h = daily_precip(g, n, 0.30, 0.7, 9.0)            # 관측(과거)
mod_h = daily_precip(g, n, 0.50, 1.0, 3.5)            # 모델(과거): 약한 비가 너무 잦고 강한 비가 약하다
mod_f = daily_precip(g, n, 0.50, 1.0, 3.5 * 1.3)      # 모델(미래): 강수 강도 30% 증가

qs = np.linspace(0, 1, 1001)
mq, oq = np.quantile(mod_h, qs), np.quantile(obs_h, qs)    # (1001,) 과거 기간의 분위수
qm = lambda x: np.interp(x, mq, oq)                   # 모델 값 → 같은 분위의 관측 값(범위 밖은 끝값 고정)

stat = lambda x: f"강수일 {np.mean(x > 0.1):.2f}  평균 {x.mean():5.2f}  99%분위 {np.percentile(x, 99):6.2f}"
print("과거 관측   ", stat(obs_h))
print("과거 모델   ", stat(mod_h))
print("과거 모델+QM", stat(qm(mod_h)))
chg = lambda stat_fn, f, h: 100 * (stat_fn(f) / stat_fn(h) - 1)     # 미래/과거 변화율 [%]
for label, s in [("평균", np.mean), ("99.9%분위", lambda x: np.percentile(x, 99.9))]:
    print(f"미래 변화 {label}: 모델 {chg(s, mod_f, mod_h):+.0f}%  QM 후 {chg(s, qm(mod_f), qm(mod_h)):+.0f}%")
print(f"미래 값 중 과거 모델 최댓값을 넘는 비율 {np.mean(mod_f > mq[-1]):.4f}")

과거 관측    강수일 0.29  평균  1.93  99%분위  24.33
과거 모델    강수일 0.49  평균  1.78  99%분위  13.98
과거 모델+QM 강수일 0.29  평균  1.94  99%분위  24.33
미래 변화 평균: 모델 +31%  QM 후 +56%
미래 변화 99.9%분위: 모델 +33%  QM 후 +74%
미래 값 중 과거 모델 최댓값을 넘는 비율 0.0004

과거 기간에서 QM은 완벽해 보인다. 강수일 비율 0.49가 관측의 0.29로, 99% 분위수 13.98 mm가 24.33 mm로 맞춰졌다. 그러나 미래의 변화 신호가 바뀐다. 모델은 평균과 극한이 모두 약 30% 늘어난다고 말하는데, QM을 거치면 평균 +56%, 99.9% 분위수 +74%가 된다. 모델의 분위 변화가 관측의 더 긴 꼬리 위로 옮겨지면서 늘어난 것이다. 보정 기간 밖의 값(최댓값 초과)은 끝값에 고정되어 잘린다. 이 예에서는 0.04%뿐이지만, 극한이 관심사라면 바로 그 0.04%가 중요하다.

QM이 극한의 상대 변화를 부풀릴 수 있다는 것은 체계적으로 보고되었고, 분위별 변화량을 보존하도록 고친 변형(quantile delta mapping 등)이 제안되었다(Cannon 외, 2015). 편향 보정이 공간·시간 구조를 고치지 못하고 모델의 근본적 오류를 덮을 수 있다는 비판적 검토도 있다(Maraun, 2016). 어떤 보정이든 과거 기간에서 맞는다는 것은 미래 기간에서의 근거가 아니다.

생성모델 기반 다운스케일링

왜 필요한가. 회귀는 제곱 오차를 최소화한다. 제곱 오차를 최소화하는 예측은 조건부 평균 $\mathbb{E}[\mathbf{y} \mid \mathbf{x}]$다(12.3장). 저해상도 값으로 정해지지 않는 국지 변동은 평균을 내면 사라진다. 그래서 회귀의 출력은 실제 고해상도 장보다 매끄럽고, 변동성과 극한이 작다.

세부 구조의 확률적 생성. 해법은 조건부 평균 대신 조건부 분포에서 표본을 뽑는 것이다. 가장 단순한 형태는 회귀 예측에 잔차 분포의 표본을 더하는 것이고, 더 일반적으로는 조건부 생성모델 $p_\theta(\mathbf{y} \mid \mathbf{x}, \mathbf{s})$를 학습한다(19.1장). 조건부 GAN과 diffusion 모델(19.4장, 19.5장)이 고해상도 강수·바람 장을 만드는 데 쓰인다. 같은 저해상도 입력에서 여러 표본을 뽑으면 국지 구조의 불확실성까지 앙상블로 얻는다.

생성의 질은 무엇으로 재는가. 회귀, 독립 잡음을 더한 회귀, 잔차의 공간 상관까지 흉내 낸 잡음을 더한 회귀를 네 가지 기준으로 비교한다. 점별 RMSE, 국지 변동(저해상도 값을 뺀 부분)의 표준편차 비, 99% 분위수의 편향, 국지 변동의 이웃 칸 상관(공간 구조)이다.

res = fine_c[tr] - np.stack([np.polyval(coef[j], X[tr, j]) for j in range(n_fine)], 1)  # (3000, 64) 학습 잔차
sd = res.std(0)                                       # (64,) 칸별 잔차 표준편차
phi = np.corrcoef(res[:, :-1].ravel(), res[:, 1:].ravel())[0, 1]    # 잔차의 이웃 칸 상관
z = np.random.default_rng(1).normal(0, 1, pred.shape)
cor = np.empty_like(z); cor[:, 0] = z[:, 0]
for j in range(1, n_fine):                            # 이웃 상관 phi를 가진 잡음
    cor[:, j] = phi * cor[:, j - 1] + np.sqrt(1 - phi**2) * z[:, j]

obs = fine_c[te]                                      # (1000, 64)
def detail(f):                                        # 저해상도 값을 뺀 국지 편차에서 날마다 달라지는 부분
    d = f - X[te]
    return d - d.mean(0)

print(f"학습 잔차의 이웃 칸 상관 phi = {phi:.2f}")
print(" RMSE  국지변동 표준편차비  99%분위 편향  국지변동 이웃상관  방법")
for name, f in [("정답", obs), ("회귀", pred), ("회귀+독립잡음", pred + sd * z), ("회귀+상관잡음", pred + sd * cor)]:
    d, q99 = detail(f), (np.percentile(f, 99, axis=0) - np.percentile(obs, 99, axis=0)).mean()
    lag1 = np.corrcoef(d[:, :-1].ravel(), d[:, 1:].ravel())[0, 1]
    print(f"{rmse(f):5.2f}   {(d.std(0) / detail(obs).std(0)).mean():10.2f}"
          f"   {q99:+11.2f}   {lag1:14.2f}     {name}")

학습 잔차의 이웃 칸 상관 phi = 0.44
 RMSE  국지변동 표준편차비  99%분위 편향  국지변동 이웃상관  방법
 0.00         1.00         +0.00             0.51     정답
 0.95         0.37         -0.38             0.78     회귀
 1.34         1.00         -0.02             0.16     회귀+독립잡음
 1.34         0.99         +0.01             0.51     회귀+상관잡음

표는 평가 지표 하나로는 다운스케일링을 판단할 수 없음을 보여 준다.

• 회귀는 RMSE가 가장 좋고 기후는 가장 틀리다. 국지 변동의 표준편차가 정답의 37%뿐이고, 99% 분위수를 0.38 °C 과소평가하며, 공간적으로 지나치게 매끄럽다(이웃 상관 0.78 대 0.51).
• 독립 잡음은 분포를 고치고 공간 구조를 망친다. 변동성과 극한은 맞지만 이웃 칸이 서로 무관하게 흔들려 이웃 상관이 0.16으로 떨어진다. 격자마다 따로 흔들리는 "모래알" 같은 장이다.
• 공간 상관까지 흉내 낸 표본만 네 기준을 모두 맞춘다. 그 대가로 RMSE는 0.95에서 1.34로 커진다. 독립적인 두 표본의 차이는 분산이 두 배이므로 $\sqrt{2}$배 가까이 커지는 것이 정상이다.

연구에서의 의미. 생성모델 기반 다운스케일링을 RMSE로만 평가하면 회귀에 질 수밖에 없다. 평가표에 분포(분위수, 변동성), 극한, 공간 구조(이웃 상관, 공간 스펙트럼 — 27.4장)를 함께 넣는다. 표본이 여러 개라면 CRPS처럼 분포 전체를 재는 점수(21.2장)를 쓴다. 반대로 생성모델의 표본이 "그럴듯해 보인다"는 시각적 인상도 근거가 될 수 없다. 19.7장의 생성모델 평가와 같은 원칙이다.

기후 모델 에뮬레이션

왜 필요한가. 온실가스·에어로졸 배출 경로마다 기후 모델을 돌리는 데는 큰 계산 자원이 든다. 정책 분석은 수백에서 수천 개의 경로를 비교하고 싶어 한다. 에뮬레이터는 기후 모델의 입력–출력 관계를 흉내 내는 값싼 대리모델(surrogate, 24.5장·26.2장)이다.

정의. 에뮬레이터는 시나리오(시간에 따른 강제력 또는 배출량) $\mathbf{u}_{1:t}$에서 기후 반응(전구 평균 기온, 지역별 기온·강수 장) $\mathbf{y}_t$로의 사상

를 기후 모델의 모의 결과로 학습한다. 가장 오래된 형태는 전구 평균 기온 변화만 에너지 수지 모형으로 흉내 내는 간단한 기후 모형과, 지역 변화가 전구 평균 기온 변화에 비례한다고 보는 패턴 스케일링(pattern scaling)이다. 학습 기반 에뮬레이터는 Gaussian process, 회귀, 신경망으로 지역별 장까지 직접 예측한다.

문제의 성질. 기후 모델의 모의 결과에는 같은 시나리오라도 초기 상태에 따라 달라지는 내부 변동성(internal variability)이 섞여 있다. 에뮬레이터가 맞혀야 하는 것은 강제된 반응(forced response)이고, 내부 변동성은 잡음이거나 따로 분포로 모델링할 대상이다. 그래서 학습 자료의 모의가 여러 초기 상태(앙상블)로 반복되어 있는지가 중요하다.

외삽 위험. 에뮬레이터를 가장 쓰고 싶은 곳은 학습에 쓴 시나리오보다 강한 강제력, 즉 학습 범위 밖이다. 강제력 0–4.5 W/m²의 시나리오 200개로 세 에뮬레이터(선형 회귀, 2차 다항식, k-최근접 이웃 평균)를 학습한다. 참 반응은 되먹임이 강제력에 따라 조금씩 세지는 약한 비선형 함수다.

true_resp = lambda f: 0.75 * f + 0.04 * f**2          # '기후 모델'의 반응 [K], 되먹임이 강제력에 따라 커진다
g = np.random.default_rng(3)
f_tr = g.uniform(0, 4.5, 200)                         # (200,) 학습 시나리오의 강제력 [W/m^2]
t_tr = true_resp(f_tr) + g.normal(0, 0.15, 200)       # 내부 변동성이 섞인 모의 결과 [K]

def knn(f, k=10):                                     # 최근접 시나리오 k개의 평균
    idx = np.argsort(np.abs(f_tr[None] - np.atleast_1d(f)[:, None]), 1)[:, :k]
    return t_tr[idx].mean(1)

emus = {"선형": np.poly1d(np.polyfit(f_tr, t_tr, 1)),
        "2차": np.poly1d(np.polyfit(f_tr, t_tr, 2)),
        "k-NN": knn}
f_val = g.uniform(0, 4.5, 100)                        # 학습 범위 안의 검증 시나리오
t_val = true_resp(f_val) + g.normal(0, 0.15, 100)
f_te = np.array([2.0, 4.0, 6.0, 8.5])                 # 4.5를 넘는 값은 학습 범위 밖
print("범위 안 검증 RMSE |  강제력 " + "  ".join(f"{f:5.1f}" for f in f_te) + "  W/m^2")
print("                  |  참 반응 " + "  ".join(f"{v:5.2f}" for v in true_resp(f_te)) + "  K")
for name, emu in emus.items():
    r = np.sqrt(np.mean((emu(f_val) - t_val) ** 2))
    print(f"{r:10.3f}        |  {name:6s}" + "  ".join(f"{v:5.2f}" for v in emu(f_te)))

범위 안 검증 RMSE |  강제력   2.0    4.0    6.0    8.5  W/m^2
                  |  참 반응  1.66   3.64   5.94   9.27  K
     0.157        |  선형     1.72   3.57   5.42   7.73
     0.149        |  2차     1.66   3.64   5.99   9.44
     0.161        |  k-NN   1.63   3.59   4.02   4.02

범위 안의 검증 RMSE는 0.149–0.161 K로, 내부 변동성(0.15 K) 수준에서 구별되지 않는다. 범위 안의 예측(2, 4 W/m²)도 모두 참값에 가깝다. 범위 밖에서는 다르다. 8.5 W/m²에서 선형은 1.5 K 과소, k-최근접 이웃은 학습 범위 끝의 값(4.02 K)에 머물러 참값의 절반도 안 된다. 2차 다항식이 맞은 것은 참 반응이 마침 2차 함수였기 때문이고, 실제 기후 모델의 반응 형태는 모른다. 범위 안의 검증으로는 이 셋 중 무엇이 외삽에 맞는지 알 수 없다.

이것이 7.2장에서 본 보간과 외삽의 구분이 기후 문제에서 나타나는 모습이다. 처방은 세 가지다. 학습 시나리오가 쓰려는 범위를 덮게 한다(가장 강한 시나리오를 학습에 넣고 중간 시나리오로 검증한다). 에너지 수지처럼 외삽해도 성립하는 물리 구조로 함수 형태를 제약한다(24.5장). 그래도 범위 밖이라면, 에뮬레이터들 사이의 불일치를 불확실성으로 보고한다.

평가

왜 필요한가. 다운스케일링과 에뮬레이션의 산출물은 날씨 예보가 아니라 기후, 즉 통계다. 기후 모델의 2050년 7월 1일은 실제 그날과 대응하지 않으므로, 날짜별 오차가 아니라 분포를 비교해야 한다. 그리고 산출물이 쓰일 곳은 학습 기간이 아닌 미래다.

평가표에는 네 축을 넣는다.

미래 기후로의 외삽 검증. 미래의 정답은 없지만 미래를 흉내 낸 검증은 할 수 있다.

• 가상 현실 검증(perfect model test, pseudo-reality): 고해상도 기후 모델의 출력을 "정답"으로 보고, 그것을 성기게 만든 자료로 다운스케일링을 학습한다. 같은 고해상도 모델의 미래 모의에도 적용해, 미래에서도 관계가 성립하는지 확인한다. 미래 "정답"이 있으므로 변화 신호의 보존을 직접 잴 수 있다.
• 기후가 다른 기간으로 나누기: 관측 기간 안에서도 따뜻한 해들과 차가운 해들, 습한 해들과 건조한 해들로 나눠 학습·시험하면, 기후 조건이 바뀔 때의 성능 저하를 미리 볼 수 있다(15.4장의 분할을 기후 조건으로 한 것).
• 변화 신호의 보존: 다운스케일링·보정 전후로 미래–과거 변화량(평균과 극한)을 비교한다. 위 QM 예처럼 보정이 변화 신호를 바꾸면, 그 바뀐 신호가 물리적으로 정당한지 설명해야 한다.

연구에서의 의미. 다운스케일링 논문의 표가 과거 기간의 RMSE와 상관만 담고 있다면, 그 결과는 미래 기후 분석에 쓰기 위한 근거가 되지 못한다. 위의 네 축과 외삽 검증 중 하나 이상을 포함하는지 확인한다.

참고 문헌

• Maraun (2016). Bias correcting climate change simulations – a critical review. Current Climate Change Reports, 2, 211–220. — 편향 보정의 가정과 한계에 대한 검토
• Cannon, Sobie, Murdock (2015). Bias correction of GCM precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes? Journal of Climate, 28(17), 6938–6959. — quantile mapping이 극한의 변화 신호를 부풀리는 현상과 quantile delta mapping