앙상블 칼만 필터

**앙상블 칼만 필터**(Ensemble Kalman filter, **EnKF**)는 지구 물리 모델의 편미분 방정식의 이산화와 같이 많은 변수가 있는 문제에 적합한 재귀 필터이다. EnKF는 대규모 문제에 대한 칼만 필터의 한 버전으로 시작되었으며(본질적으로 공분산 행렬이 표본 공분산으로 대체됨), 현재 앙상블 예보의 중요한 데이터 동화 구성 요소이다. EnKF는 파티클 필터와 관련이 있지만(이 맥락에서 파티클은 앙상블 구성원과 동일하다), EnKF는 관련된 모든 확률 분포가 정규 분포라는 가정을 한다. 적용 가능할 때, EnKF는 파티클 필터보다 훨씬 효율적이다.
# 서론
앙상블 칼만 필터(EnKF)는 베이즈 추론 문제에 대한 몬테카를로 구현이다. 모델링된 시스템의 상태에 대한 확률 밀도 함수(PDF) (사전 확률, 지구 과학에서는 종종 예측이라고 불림)와 데이터 가능도가 주어지면, 베이즈 정리를 사용하여 데이터 가능도를 고려한 후의 PDF (사후 확률, 종종 분석이라고 불림)를 얻는다. 이를 베이즈 업데이트라고 한다. 베이즈 업데이트는 모델을 시간적으로 진행시키는 것과 결합되어 주기적으로 새로운 데이터를 통합한다. 1960년에 소개된 원래의 칼만 필터는doi=10.1115/1.3662552 모든 PDF가 정규 분포라는 가정을 하며 (가우스 가정), 베이즈 업데이트에 의한 평균 및 공분산 행렬의 변화에 대한 대수적 공식과 시스템이 선형일 경우 평균 및 공분산을 시간적으로 진행시키는 공식을 제공한다. 그러나 고차원 시스템의 경우 공분산 행렬을 유지하는 것은 계산적으로 불가능하다. 이러한 이유로 EnKF가 개발되었다.hdl-access=freebibcode=1998MWRv..126..796H EnKF는 앙상블이라고 불리는 상태 벡터 컬렉션을 사용하여 시스템 상태의 분포를 나타내고, 공분산 행렬을 앙상블에서 계산된 표본 공분산으로 대체한다. 앙상블은 마치 무작위 표본인 것처럼 작동하지만, 앙상블 구성원들은 EnKF를 공유하기 때문에 실제로는 독립적이지 않다. 그들은 대략 독립적이라고 간주되며, 모든 계산은 실제로 독립적인 것처럼 진행된다. EnKF의 한 가지 장점은 PDF를 시간적으로 진행시키는 것이 앙상블의 각 구성원을 단순히 진행시키는 것으로 달성된다는 것이다.EnKF 및 관련 데이터 동화 기법에 대한 조사는 다음을 참조. isbn=978-3-540-38300-0
# 유도

## 칼만 필터
칼만 필터 \mathbf{x}를 모델의 n차원 상태 벡터라고 하고, 평균 \mathbf{\mu}와 공분산 Q를 갖는 가우스 확률 분포를 따른다고 가정하면, 즉 PDF는 다음과 같다. : p(\mathbf{x})\propto\exp\left( -\frac{1}{2}(\mathbf{x}-\mathbf{\mu })^{\mathrm{TQ^{-1}(\mathbf{x}-\mathbf{\mu})\right) . 여기서 그리고 아래에서 \propto는 비례함을 의미한다. PDF는 항상 전체 공간에 대한 적분이 1이 되도록 스케일링된다. 이 p(\mathbf{x})를 사전 확률이라고 부르며, 모델을 실행하여 시간적으로 진화했으며 이제 새로운 데이터를 설명하기 위해 업데이트되어야 한다. 데이터의 오차 분포가 알려져 있다고 가정하는 것이 자연스럽다. 데이터는 오차 추정치와 함께 제공되어야 하며, 그렇지 않으면 무의미하다. 여기서 데이터 \mathbf{d}는 공분산 R과 평균 H\mathbf{x}를 갖는 가우스 PDF를 갖는다고 가정한다. 여기서 H는 소위 관측 행렬이다. 공분산 행렬 R은 데이터 오차의 추정치를 설명한다. 데이터 벡터 \mathbf{d}의 엔트리에서 무작위 오차가 독립적이라면, R은 대각선이며 대각선 엔트리는 데이터 벡터 \mathbf{d}의 해당 엔트리 오차의 표준 편차("오차 크기")의 제곱이다. 값 H\mathbf{x}는 데이터 오차가 없을 때 상태 \mathbf{x}에 대한 데이터 값이 될 것이다. 그러면 시스템 상태 \mathbf{x}에 조건부인 데이터 \mathbf{d}의 확률 밀도 p(\mathbf{d}|\mathbf{x}) (이것은 데이터 가능도라고 불린다)는 다음과 같다. : p\left( \mathbf{d}|\mathbf{x}\right) \propto\exp\left( -\frac{1}{2}(\mathbf{d}-H\mathbf{x})^{\mathrm{TR^{-1}(\mathbf{d}-H\mathbf{x})\right) . 상태의 PDF와 데이터 가능도는 베이즈 정리에 의해 데이터 \mathbf{d}의 값에 조건부인 시스템 상태 \mathbf{x}의 새로운 확률 밀도 (이것은 사후)를 제공하기 위해 결합된다. : p\left( \mathbf{x}|\mathbf{d}\right) \propto p\left( \mathbf{d}|\mathbf{x}\right) p(\mathbf{x}). 데이터 \mathbf{d}는 일단 수신되면 고정되므로, 사후 상태를 \mathbf{x}|\mathbf{d} 대신 \mathbf{\hat{x로 나타내고 사후 PDF를 p\left( \mathbf{\hat{x\right) 로 나타낸다. 대수적 조작을 통해isbn=978-0-13-638122-8 사후 PDF도 가우스 분포임을 보일 수 있다. : p\left( \mathbf{\hat{x\right) \propto\exp\left( -\frac{1}{2}(\mathbf{\hat{x-\mathbf{\hat{\mu)^{\mathrm{T\hat{Q}^{-1}(\mathbf{\hat{x-\mathbf{\hat{\mu)\right) , 사후 평균 \mathbf{\hat{\mu와 공분산 \hat{Q}는 칼만 업데이트 공식에 의해 주어진다. : \mathbf{\hat{\mu=\mathbf{\mu}+K\left( \mathbf{d}-H\mathbf{\mu}\right) ,\quad\hat{Q}=\left( I-KH\right) Q, 여기서 : K=QH^{\mathrm{T\left( HQH^{\mathrm{T+R\right) ^{-1} 는 소위 칼만 이득 행렬이다.
## 앙상블 칼만 필터
EnKF는 상태 벡터 \mathbf{x}의 PDF의 공분산 행렬을 진화시키는 것을 피하는 칼만 필터의 몬테카를로 근사이다. 대신, PDF는 앙상블에 의해 표현된다. : X=\left[ \mathbf{x}_{1},\ldots,\mathbf{x}_{N}\right] =\left[ \mathbf{x}_{i}\right]. X는 열이 앙상블 구성원인 n\times N 행렬이며, 사전 앙상블이라고 불린다. 이상적으로, 앙상블 구성원들은 사전 분포로부터의 무작위 표본을 형성할 것이다. 그러나 앙상블 구성원들은 모든 EnKF 단계가 그들을 묶기 때문에 초기 앙상블을 제외하고는 일반적으로 독립적이지 않다. 그들은 대략 독립적이라고 간주되며, 모든 계산은 실제로 독립적인 것처럼 진행된다. 데이터 \mathbf{d}를 m\times N 행렬로 복제한다. : D=\left[ \mathbf{d}_{1},\ldots,\mathbf{d}_{N}\right] =\left[ \mathbf{d}_{i}\right], \quad \mathbf{d}_{i}=\mathbf{d}+\mathbf{\epsilon_{i, \quad \mathbf{\epsilon_{i \sim N(0,R), 따라서 각 열 \mathbf{d}_{i}는 데이터 벡터 \mathbf{d}에 m차원 정규 분포 N(0,R)에서 온 무작위 벡터를 더한 것으로 구성된다. 또한, X의 열이 사전 확률 분포에서 온 표본이라면, : \hat{X}=X+K(D-HX) 의 열은 사후 확률 분포에서 온 표본을 형성한다. 스칼라 경우에 H=1일 때 이를 살펴보자. x_i = \mu + \xi_i, \; \xi_i \sim N(0, \sigma_x^2), 그리고 d_i = d + \epsilon_i, \; \epsilon_i \sim N(0, \sigma_d^2)라고 하자. 그러면 :\hat{x}_i = \left(\frac{1/\sigma_x^2}{1/\sigma_x^2 + 1/\sigma_d^2} \mu + \frac{1/\sigma_d^2}{1/\sigma_x^2 + 1/\sigma_d^2} d \right)+ \left(\frac{1/\sigma_x^2}{1/\sigma_x^2 + 1/\sigma_d^2} \xi_i + \frac{1/\sigma_d^2}{1/\sigma_x^2 + 1/\sigma_d^2} \epsilon_i \right) . 첫 번째 합은 사후 평균이고, 두 번째 합은 독립성을 고려할 때 다음과 같은 분산을 갖는다. :\left(\frac{1/\sigma_x^2}{1/\sigma_x^2 + 1/\sigma_d^2}\right)^2 \sigma_x^2 + \left(\frac{1/\sigma_d^2}{1/\sigma_x^2 + 1/\sigma_d^2}\right)^2 \sigma_d^2 = \frac{1}{1/\sigma_x^2 + 1/\sigma_d^2}, 이것이 사후 분산이다. EnKF는 칼만 이득 행렬 K에서 상태 공분산 Q를 앙상블 구성원으로부터 계산된 표본 공분산 C(앙상블 공분산이라고 함)로 대체함으로써 얻어진다.citeseerx=10.1.1.67.4916 즉: K=CH^{\mathrm{T\left( HCH^{\mathrm{T+R\right) ^{-1}
# 구현

## 기본 공식
여기서는 다음을 따른다.bibcode=1998MWRv..126.1719Bbibcode=2003OcDyn..53..343E 앙상블 행렬 X와 데이터 행렬 D가 위와 같다고 가정한다. 앙상블 평균과 공분산은 다음과 같다. : E\left( X\right) =\frac{1}{N}\sum_{k=1}^{N}\mathbf{x}_{k},\quad C=\frac{AA^{T{N-1}, 여기서 : A=X-E\left( X\right) \mathbf{e}_{1\times N} =X-\frac{1}{N}\left( X\mathbf{e}_{N\times1}\right) \mathbf{e}_{1\times N}, 이고 \mathbf{e}는 지정된 크기의 모든 원소가 1인 행렬을 나타낸다. 사후 앙상블 X^{p}는 다음과 같이 주어진다. : X^{p}=X+CH^{T}\left( HCH^{T}+R\right) ^{-1}(D-HX), 여기서 교란된 데이터 행렬 D는 위와 같다. R은 공분산 행렬이므로 항상 양의 준정부호이고 보통 양의 정부호이므로 위 역행렬이 존재하며 이 공식은 숄레스키 분해를 사용하여 구현할 수 있다.date=June 2006에서는 R이 표본 공분산 \tilde{D} \tilde{D}^{T}/\left( N-1\right) 로 대체되며, 여기서 \tilde{D} = D - \frac{1}{N} d \, \mathbf{e}_{1\times N}이고 역행렬은 의사 역행렬로 대체되며, 특이값 분해(SVD)를 사용하여 계산된다. 이 공식들은 레벨 3 연산이 주가 되는 행렬 연산이므로,isbn=978-0-8018-3772-2 LAPACK(직렬 및 공유 메모리 컴퓨터) 및 ScaLAPACK(분산 메모리 컴퓨터)과 같은 소프트웨어 패키지를 사용하여 효율적인 구현에 적합하다. 행렬의 역행렬을 계산하고 곱하는 대신, 행렬의 숄레스키 분해를 계산하고 역행렬 곱셈을 많은 동시 우변을 가진 선형 시스템의 해로 처리하는 것이 훨씬 더 좋다 (몇 배 더 저렴하고 더 정확하다).
## 관측 행렬이 없는 구현
공분산 행렬을 앙상블 공분산으로 대체했으므로, 행렬 H를 명시적으로 지정하지 않고도 앙상블 관측치를 직접 사용하는 더 간단한 공식이 나온다. 더 구체적으로, 다음과 같은 형태의 함수 h(\mathbf{x})를 정의한다. : h(\mathbf{x})=H\mathbf{x}. 함수 h는 관측 함수 또는 역문제 맥락에서는 순방향 연산자라고 불린다. h(\mathbf{x})의 값은 측정값이 정확하다고 가정할 때 상태 \mathbf{x}에 대한 데이터 값이 될 것이다. 그러면 사후 앙상블은 다음과 같이 다시 쓸 수 있다. : X^{p}=X+\frac{1}{N-1}A\left( HA\right) ^{T}P^{-1}(D-HX) 여기서 : HA=HX-\frac{1}{N}\left( \left( HX\right) \mathbf{e}_{N\times1}\right) \mathbf{e}_{1\times N}, 그리고 : P=\frac{1}{N-1}HA\left( HA\right) ^{T}+R, 와 함께 :\left[ HA\right] _{i} =H\mathbf{x}_{i}-H\frac{1}{N}\sum_{j=1}^{N}\mathbf{x}_{j}\ =h\left( \mathbf{x}_{i}\right) -\frac{1}{N}\sum_{j=1}^{N}h\left( \mathbf{x}_{j}\right) . 결과적으로 앙상블 업데이트는 각 앙상블 구성원에 대해 관측 함수 h를 한 번 평가함으로써 계산될 수 있으며, 행렬 H를 명시적으로 알 필요가 없다. 이 공식은 고정된 오프셋 \mathbf{f}를 갖는 관측 함수 h(\mathbf{x})=H\mathbf{x+f}에 대해서도 성립하며, 이 또한 명시적으로 알 필요가 없다. 위 공식은 허리케인 소용돌이의 위치와 같은 비선형 관측 함수 h에 일반적으로 사용되어 왔다.doi-access=free 이 경우 관측 함수는 앙상블 구성원의 값으로부터 선형 함수로 본질적으로 근사된다.
## 대량의 데이터 포인트에 대한 구현
많은 수의 데이터 포인트 m의 경우 P^{-1}을 곱하는 것이 병목 현상이 된다. 다음 대안 공식은 데이터 포인트 수 m이 많고(예: 그리드 또는 픽셀 데이터 동화) 데이터 오류 공분산 행렬 R이 대각 행렬(데이터 오류가 상관되지 않는 경우)이거나 분해하기 저렴할 때(제한된 공분산 거리로 인해 밴드 행렬인 경우) 유리하다. 셔먼-모리슨-우드베리 공식을 사용하면doi=10.1137/1031049 : (R+UV^{T})^{-1}=R^{-1}-R^{-1}U(I+V^{T}R^{-1}U)^{-1}V^{T}R^{-1}, 와 함께 : U=\frac{1}{N-1}HA,\quad V=HA, 다음을 얻는다. :\begin{align} P^{-1} & =\left( R+\frac{1}{N-1}HA\left( HA\right) ^{T}\right) ^{-1}\ = \\\ & =R^{-1}\left[ I-\frac{1}{N-1}\left( HA\right) \left( I+\left( HA\right) ^{T}R^{-1}\frac{1}{N-1}\left( HA\right) \right) ^{-1}\left( HA\right) ^{T}R^{-1}\right] , \end{align} 이는 행렬 R을 포함하는 시스템(저렴하다고 가정)과 m개의 우변을 가진 크기 N 시스템의 해만 필요로 한다.
# 추가 확장
여기서 설명된 EnKF 버전은 데이터의 무작위화를 포함한다. 데이터의 무작위화가 없는 필터에 대해서는 다음을 참조하라.bibcode=2001MWRv..129.2884Abibcode=2004OcDyn..54..539Ebibcode=2003MWRv..131.1485T 앙상블 공분산은 랭크 부족하다 (상태 변수 수가 앙상블 멤버 수보다 훨씬 많다. 일반적으로 수백만 대 수십 미만). 이로 인해 공간적으로 멀리 떨어진 지점 쌍에 대한 큰 항을 갖는다. 실제로는 멀리 떨어진 위치의 물리적 필드 값이 그렇게 많이 상관되지 않으므로, 거리에 따라 공분산 행렬을 인위적으로 테이퍼링하여 지역화된 EnKF 알고리즘이 발생한다.bibcode=2003MWRv..131..634Abibcode=2004TellA..56..415O 이러한 방법은 계산에 사용되는 공분산 행렬을 수정하며, 결과적으로 사후 앙상블은 더 이상 사전 앙상블의 선형 조합으로만 이루어지지 않는다. 비선형 문제의 경우, EnKF는 비물리적인 상태를 가진 사후 앙상블을 생성할 수 있다. 이는 큰 공간 기울기를 가진 상태에 대한 벌점화와 같은 티호노프 정칙화에 의해 완화될 수 있다. 허리케인, 뇌우, 산불선, 스콜 라인, 비 전선과 같은 코히어런트 피쳐가 있는 문제의 경우, 상태 진폭을 가산적으로 보정하는 것뿐만 아니라 공간(그리드)에서 상태를 변형하여 수치 모델 상태를 조정할 필요가 있다. 2007년에 라벨라 외 다수는 앙상블을 사용하여 공동 위치-진폭 조정 모델을 도입하고, EnKF 및 기타 공식에 적용할 수 있는 순차적 근사치를 체계적으로 도출했다.bibcode=2007PhyD..230..127R 그들의 방법은 다른 사람들이 하는 것처럼 진폭 및 위치 오류가 독립적이거나 공동 가우스 분포를 따른다고 가정하지 않는다. 모핑 EnKF는 상태의 선형 조합 대신 영상 정합 및 모핑에서 빌려온 기술로 얻은 중간 상태를 사용한다.bibcode=2008TellA..60..131Bconference=11th Symposium on Integrated Observing and Assimilation Systems for the Atmosphere, Oceans, and Land Surface (IOAS-AOLS), CD-ROM, Paper 4.12, 87th American Meteorological Society Annual Meeting, San Antonio, TX, January 2007 형식적으로 EnKF는 가우스 가정을 따른다. 실제로는 가우스 가정이 충족되지 않을 수 있는 비선형 문제에도 사용될 수 있다. EnKF의 장점을 유지하면서 가우스 가정을 완화하려는 관련 필터에는 여러 가우스 커널로 상태 PDF를 피팅하는 필터,bibcode=1999MWRv..127.2741A 가우스 혼합으로 상태 PDF를 근사하는 필터,doi-access=free 밀도 추정에 의한 파티클 가중치 계산을 사용하는 파티클 필터의 변형, 및 파티클 필터 퇴화를 완화하기 위해 코시 분포를 갖는 두꺼운 꼬리 데이터 PDF를 사용하는 파티클 필터의 변형이 포함된다.bibcode=2003MWRv..131.2071V