気象情報解析特論　第8回　特異値分解を用いた気象データの解析（MCA）

担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import math
from matplotlib.colors import Normalize

下準備：SST・海面更正気圧のデータセットの読み込みと関数の定義

今回は，海面水温（SST）と海面更正気圧（SLP，以下単に気圧と呼ぶ）の両方のデータを用います。

SLPのデータ： http://web.is.ocha.ac.jp/~tsubasa/env_info/slpa_erai.npz

In [2]:
loadfile = 'detrended_ssta_OISST.npz' # デトレンド済みの海面水温偏差の入力ファイル名を定義
ssta_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
ssta = ssta_dataset['ssta'] # 海面水温偏差を変数sstaに保存
lon2_sst = ssta_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2_sst = ssta_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = ssta_dataset['y'] # 年(year)を変数yに保存

# 2つのデータセットに共通な1982年から2018年のみのデータを用いる
ssta = ssta[:, :, (1982 <= y)*(y <= 2018)]

loadfile = 'slpa_erai.npz' # 海面更正気圧偏差の入力ファイル名を定義
slpa_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
slpa = slpa_dataset['slpa'] # 海面更正気圧偏差を変数slpaに保存
lon2_slp = slpa_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2_slp = slpa_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = slpa_dataset['y'] # 年(year)を変数yに保存
m = slpa_dataset['m'] # 月(month)を変数mに保存
slpa = signal.detrend(slpa)# デトレンド（signal.detrendは一番後ろの次元についてdetrendしてくれる）

# 2つのデータセットに共通な1982年から2018年のみのデータを用いる
slpa = slpa[:, :, (1982 <= y)*(y <= 2018)]
m = m[(1982 <= y)*(y <= 2018)]
y = y[(1982 <= y)*(y <= 2018)] # y自身のサイズが変わってしまうので，yは一番最後に書き換えないとダメ

# 気象場を2つ描画する関数
def draw_two_fields(field, lon2, lat2, field_2, lon2_2, lat2_2, fig_title, vmin = -3, vmax = 3, vint = 1):        
    plt.figure()
    cm = plt.get_cmap('seismic')
    cs = plt.contourf(lon2, lat2, field,\
                  cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),\
                  levels=np.arange(vmin,vmax+vint,vint), extend='both')
    plt.colorbar(cs)
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    title = fig_title
    plt.title(title)
    plt.contour(lon2_2, lat2_2, field_2, colors=['black'])
    plt.xlim(100, 250)
    plt.ylim(10, 70)
    plt.show()
    
# 月別の時系列を2つ描画する関数
def plot_2_mon_time(time_series1, time_series2, lower = -3, upper = 3, init_year=1982, fin_year=2018):
    mon = np.arange(init_year, fin_year+1, 1/12)
    plt.plot(mon, time_series1)
    plt.plot(mon, time_series2, 'r')
    plt.plot(mon, 0*time_series1, 'k')
    plt.xlim(init_year, fin_year)
    plt.ylim(lower, upper)
    plt.show()
    
# 回帰図を計算する関数    
def reg_map(ref_time, var):
    [imt, jmt, tmt] = var.shape 

    # 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
    a_var = np.zeros((imt, jmt)) # 回帰係数a
    b_var = np.zeros((imt, jmt)) # 切片b
    
    # np.polyfitがエラーを吐かないようにするために，陸地の場所(nanが入っている)に一度ゼロを入れておきたい
    # （nanmeanのようなnanpolyfitという関数はない）
    is_land_grids_3D = (np.isnan(var)==True) 
    var[is_land_grids_3D]=0

    # 回帰図の計算
    # 経度方向にimt回，緯度方向にjmt回forループを回す
    for ii in range(0, imt):
        for jj in range(0, jmt):
            [a_var[ii, jj], b_var[ii, jj]] = np.polyfit(ref_time, var[ii, jj, :], 1)

        if (ii % 100 == 0):
            print(ii) # ちゃんと計算が進んでるかチェックするために，100回に1回iiを出力する
    
    # さっきゼロにしておいた陸地の場所にもう一度nanを戻す
    var[is_land_grids_3D]=np.nan
    is_land_grids_2D = np.squeeze(is_land_grids_3D[:, :, 0])
    a_var[is_land_grids_2D]=np.nan
    b_var[is_land_grids_2D]=np.nan
    
    return a_var

前回までのあらすじ

前回の最後で，特異値分解が気象データ解析の文脈で出てくるのは，主に3通りの場合であると言いました。

そのうち，3つ目に当たる「最大共分散解析」が今回のテーマです。

二つのデータ行列$X$と$Y$の共分散行列 $C_{XY} = XY^T/N$を特異値分解＝最大共分散分析(MCA)

すごく大雑把にいうと主成分分析の2変数バージョンです。日本の大気海洋研究コミュニティで「SVD解析を行った」という場合は，大抵この最大共分散分析（Maximum Covariance Analysis; MCA）を指します。

EOF解析は，データ行列$X$の分散を最大にするような基底を見つける分析方法でした。

それに対しMCAは，（サンプリング次元の長さが等しい）データ行列$X$と$Y$の共分散を最大にするような基底を見つける分析方法です。たとえば，海面水温と気圧のSVD1は，そのSVD1に射影して求まる時系列の共分散が最大となるような空間パターンです。

主成分分析は「1変数の卓越した変動」を探す解析でした。特異値分解は「2変数の関係が深い変動」を探す解析，といってもいいでしょう。

今回は，Wallace et al. 1992を参考に，SSTとSLP（海面更正気圧; sea level pressure）でMCAを行ってみましょう。

MCAの流れ

主成分分析の流れと比べながら見てみてください。わかりやすいように$X$は2箇所，$Y$は3箇所の観測地点があったことにして説明します。

実際にデータ解析をする際には，たとえば$X$がSST，$Y$がSLPのデータだと思っていただければ良いです（が，抽象的で分かりづらいと思うので，1回目はざっと目を通して実践例に進んでくれても良いです）。

I. データ行列の定義（XとYの大きさは，時間方向だけあっていれば良い）

$$X := \left( \begin{array}{c} \vec{x_1}\\ \vec{x_2} \end{array} \right) = \left( \begin{array}{cccc} x_{1, 1} & x_{1, 2} & \cdots & x_{1, N}\\ x_{2, 1} & x_{2, 2} & \cdots & x_{2, N} \end{array} \right)$$$$Y := \left( \begin{array}{c} \vec{y_1}\\ \vec{y_2}\\ \vec{y_3} \end{array} \right) = \left( \begin{array}{cccc} y_{1, 1} & y_{1, 2} & \cdots & y_{1, N}\\ y_{2, 1} & y_{2, 2} & \cdots & y_{2, N}\\ y_{3, 1} & y_{3, 2} & \cdots & y_{3, N} \end{array} \right)$$

II. 共分散行列の計算

$$C_{XY} := \left( \begin{array}{cc} \mathrm{cov}(\vec{x_1}, \vec{y_1}) & \mathrm{cov}(\vec{x_1}, \vec{y_2}) & \mathrm{cov}(\vec{x_1}, \vec{y_3})\\ \mathrm{cov}(\vec{x_2}, \vec{y_1}) & \mathrm{cov}(\vec{x_2}, \vec{y_2}) & \mathrm{cov}(\vec{x_1}, \vec{y_3}) \end{array} \right) = XY^T/N$$

III. 共分散行列を特異値分解する

$$C_{XY} = U \Sigma V^T$$

IV. $U$の特異ベクトルは$X$の空間構造を与え，$V$の特異ベクトルは$Y$の空間構造を与えるので，それぞれ特異値の大きい特異ベクトルから順にセットとして，SVD第1モード(SVD1)，SVD第2モード(SVD2)，...の空間構造の組とする（例：SVD1は {$\vec{u_1}$, $\vec{v_1}$ }の組）。

$$ U = \left( \begin{array}{cc} \vec{u_1} & \vec{u_2} \end{array} \right)$$$$ V = \left( \begin{array}{ccc} \vec{v_1} & \vec{v_2} & \vec{v_3} \end{array} \right)$$$$ \Sigma = \left( \begin{array}{ccc} \sigma_1 & 0 & 0\\ 0 & \sigma_2 & 0 \end{array} \right), \ \ \sigma_1 \geq \sigma_2 > 0$$

V. 特異ベクトルを並べた行列$U$に，元の行列$X$を射影すると，SVD1$_X$時系列，SVD2$_X$時系列，...を得る

$$ X^\ast = U^{T}X$$

ただし

$$X^\ast := \left( \begin{array}{c} \vec{SVD1}_X\\ \vec{SVD2}_X \end{array} \right)$$

VI. 特異ベクトルを並べた行列$V$に，元の行列$Y$を射影すると，SVD1$_Y$時系列，SVD2$_Y$時系列，...を得る

$$ Y^\ast = V^{T}Y$$

ただし

$$Y^\ast := \left( \begin{array}{c} \vec{SVD1}_Y\\ \vec{SVD2}_Y \end{array} \right)$$

VII. SVD1はSVD1$_X$とSVD1$_Y$の共分散が最大となるような空間構造の組であり，SVD2以降は「特異ベクトルが直交する」という制約の中で共分散の大きい順となる。特異値は各モードのSVD時系列同士の共分散を与える。それぞれのモードのSVD時系列は共分散を共有せず，MCA前後でデータ全体の共分散の和は保存される。

$$ X^\ast Y^{\ast T} = N \Sigma $$

実践例：海面水温（SST）と海面更正気圧（SLP）でMCAしてみる

さて，ここからは北太平洋のSSTと，北半球のSLPについてMCAを行っていきます。

データ行列の準備

まず，それぞれの緯度経度を抜き出します。

In [3]:
# 北太平洋(North Pacific; NP)のSSTは10N-70N, 100E-110Wで定義
ssta_NP = ssta[(100<=lon2_sst[:, 1])*(lon2_sst[:, 1]<=250),:, :]
ssta_NP = ssta_NP[:, (10<=lat2_sst[1, :])*(lat2_sst[1, :]<=70), :]
lon2_sst_NP = lon2_sst[(100<=lon2_sst[:, 1])*(lon2_sst[:, 1]<=250),:]
lon2_sst_NP = lon2_sst_NP[:, (10<=lat2_sst[1, :])*(lat2_sst[1, :]<=70)]
lat2_sst_NP = lat2_sst[(100<=lon2_sst[:, 1])*(lon2_sst[:, 1]<=250),:]
lat2_sst_NP = lat2_sst_NP[:, (10<=lat2_sst[1, :])*(lat2_sst[1, :]<=70)]
# 北半球(Northern Hemisphere; NH)のSLPは10N-90Nで定義
slpa_NH = slpa[:, (10<=lat2_slp[1, :]), :]
lon2_slp_NH = lon2_slp[:, (10<=lat2_slp[1, :])]
lat2_slp_NH = lat2_slp[:, (10<=lat2_slp[1, :])]

まず，EOF解析のときと同じように，latlon_to_Xを用いて，データ行列$X$と$Y$を作ります。

In [4]:
# latlon格子のデータをデータ行列Xに変換する関数
## var: 特異値分解したい3次元配列(lon, lat, time)
## lat2: 各格子の緯度が入った2次元行列
## mask: NaNではないところに1,NaNのところは0が入っている2次元行列
def latlon_to_X(var, lat2, mask):
    
    # 行列のサイズを得る
    imt, jmt, tmt = var.shape
    
    # 出力するデータ行列を初期化しておく（このとき，最大サイズ（すべて海）を仮定）
    X = np.zeros((imt*jmt, tmt))
    
    # 緯度のコサインの平方根をかけておく
    for tt in range(0, tmt):
        var[:, :, tt] = var[:, :, tt] * np.sqrt(np.cos(lat2*math.pi/180))
    
    # データ行列の作成
    cc = 0 # 何番目のNaNではない格子なのか数えるカウンター
    for ii in range(0, imt):
        for jj in range(0, jmt):
            if mask[ii, jj] == 1 :  # もしii, jj番地の格子がNaNではなければ
                X[cc, :] = np.squeeze(var[ii, jj, :]) # cc行目に，ii, jj番地の経度緯度のSST時系列を入れて
                cc = cc + 1 # カウンターをインクリ
                
    X = np.squeeze(X[0:cc, :]) # すべて海を仮定したが，実際は陸があるので結局cc-1行目までしかデータは入っていない
    
    return X

# NaNではないところに1,NaNのところは0が入っている2次元行列
mask_sst = np.isfinite(np.squeeze(ssta_NP[:, :, 0])) 
mask_slp = np.isfinite(np.squeeze(slpa_NH[:, :, 0])) # slpの場合はNaNはないので，単にnp.onesでも良い

# SST偏差のデータについて，データ行列Xの計算
X = latlon_to_X(ssta_NP, lat2_sst_NP, mask_sst)

# SLP偏差のデータについて，データ行列Yの計算
Y = latlon_to_X(slpa_NH, lat2_slp_NH, mask_slp)

出来上がった行列のサイズを見てみましょう。

In [5]:
X.shape
Out[5]:
(6075, 444)
In [6]:
Y.shape
Out[6]:
(3240, 444)

このように，時間方向（サンプリング次元）の方だけ揃った二つのデータ行列に，MCAは適用できます。

特異値分解

早速，$X$と$Y$の共分散行列$C_{XY} = XY^T/N$について，特異値分解していきましょう。

（大きな行列の特異値分解なので，少し計算に時間がかかります）

In [7]:
# 時間方向のサンプルサイズ
N = X.shape[1]
# 共分散行列の計算
C_XY = X@Y.T/N

# 特異値分解
# full_matricesをFalseにすると，min(M, N)番目の特異値・特異ベクトルまで
# しか計算しないので計算量を節約できる。
U, S, VT = np.linalg.svd(C_XY, full_matrices=False) 

# Vの転置
V = VT.T 

# Sigmaの計算
Sigma = np.zeros(C_XY.shape) # Xと同じ大きさのfloat64型行列Sigmaを用意
for ss in range(0, S.shape[0]):
    Sigma[ss, ss] = S[ss] #特異値を対角行列に並べる

# U, Vに元の行列を射影することでSVD時系列の計算
U_star = U.T@X
V_star = V.T@Y

SVDの描画

ここで，$U$と$V$の構造を見てみるためにlatlon格子に戻してみましょう。

In [8]:
# 直交行列U, Vをlatlon格子のデータに変換する関数
## U, V: SVDの空間構造の入った展開行列U(space, mode), V(space, mode)
## mask: NaNではないところに1,NaNのところは0が入っている2次元行列
def UV_to_latlon(U, mask):
        
    # 行列のサイズを得る
    smt, mmt = U.shape
    imt, jmt = mask.shape
    
    # 出力するデータ行列を初期化しておく（このとき，最大サイズ（すべて海）を仮定）
    SVDs = np.zeros((imt, jmt, mmt))
    
    # lat-lon行列の作成
    cc = 0 # 何番目のNaNではない格子なのか数えるカウンター
    for ii in range(0, imt):
        for jj in range(0, jmt):
            if mask[ii, jj] == 1 :  # もしii, jj番地の格子がNaNではなければ
                SVDs[ii, jj, :] = U[cc, :] # ii, jj番地に，cc列目の格子のデータを入れる
                cc = cc + 1 # カウンターをインクリ
            else:  # もしii, jj番地の格子がNaNであれば
                SVDs[ii, jj, :] = np.nan # EOFsの方のii, jj番地にもNaNを入れる
                 
    return SVDs

sst_SVDs = UV_to_latlon(U, mask_sst)
slp_SVDs = UV_to_latlon(V, mask_slp)

それでは，SVD第1モードを描画してみましょう！

In [9]:
svd1_sst_map = sst_SVDs[:, :, 0]
svd1_slp_map = slp_SVDs[:, :, 0]
draw_two_fields(svd1_sst_map, lon2_sst_NP, lat2_sst_NP,svd1_slp_map, lon2_slp_NH, lat2_slp_NH,'SVD1 of SST and SLP', -0.04, 0.04, 0.004)

ご覧の通り，暖かい海水と高気圧が少しずれて並んでいるパターンが見えてきました。

北半球では，大雑把には高気圧を右に見て等圧線に沿って風が吹くことが知られています（地衡風）。そのため，暖かい水があるところは南風，冷たい水があるところは北風になっていると考えることもできるでしょう。逆に，暖かい水によって偏西風ジェットが北に上がることも近年わかってきました(e.g., Ogawa et al. 2012)。

いずれにせよMCAによって，海面水温と気圧が最も密接に結びついている（=互いに最も共分散をよく説明する）パターンを取り出すことができました。

SVD1の時系列

次に，SVD1の時系列を見てみます。EOF解析の時と同じように，時系列は正規化しましょう。

In [10]:
svd1_sst_time = U_star[0, :]
svd1_sst_time = svd1_sst_time/np.std(svd1_sst_time) # 正規化
svd1_slp_time = V_star[0, :]
svd1_slp_time = svd1_slp_time/np.std(svd1_slp_time) # 正規化

plot_2_mon_time(svd1_sst_time, svd1_slp_time) # SVD1_SSTの時系列を青で，SVD1_SLPの時系列を赤で描画する

# SVD第一モードは共分散全体のうち何%を説明するか
print(Sigma[0, 0]/np.trace(Sigma)*100)
# 規格化後のSVD第一モードは，互いの分散を何%説明するか（相関係数の2乗で計算できる）
print((np.corrcoef(svd1_sst_time, svd1_slp_time)[0, 1])**2*100)

16.165042575120495
20.94533523012629

ちょっと見づらいので，二つの時系列にバタワースローパスフィルタをかけてみましょう。

In [11]:
# 引数は，時系列，サンプリング間隔（月別データだと1ヶ月），カットオフ周波数，フィルタオーダー
def butterworth_lowpass(time_series, delta_t, f_cut, order):
    b, a = signal.butter(order, f_cut, btype='low', fs = 1/delta_t) # フィルタオーダーは1
    time_series_lp = signal.filtfilt(b, a, time_series) # 与えられた時系列に大してフィルターを実行する
    return time_series_lp

svd1_sst_time_low = butterworth_lowpass(svd1_sst_time, 1, 0.12, 1)
svd1_slp_time_low = butterworth_lowpass(svd1_slp_time, 1, 0.12, 1)
plot_2_mon_time(svd1_sst_time_low, svd1_slp_time_low)

# 規格化後のSVD第一モードは，互いの分散を何%説明するか（相関係数の2乗で計算できる）
## 論文等では，普通はこちらの値を示すことが多い気がします
print((np.corrcoef(svd1_sst_time_low, svd1_slp_time_low)[0, 1])**2*100)

42.05564209992879

完璧とは言わないまでも，そこそこ良い相関に見えます。定量的には，（ローパスフィルター後の値で）互いの分散を42%説明するようなSVD1を取り出すことができました。

2つのデータから取り出しうる空間パターンの中で，この相関（厳密には，規格化前の時系列の共分散）が最大になるように選んだ空間パターンこそが，先ほど描画したSVD1なのです！

SVD時系列への回帰図（振幅の情報を持つSVD）

最後に，EOFのときと同じように，SVD時系列への回帰図を示します。これにより，SVD時系列が1標準偏差と等しくなったときにSSTやSLPがどの程度の振幅を持つのかがわかるようになります（研究では普通こちらの回帰図を論文に載せるため，実践上は先ほどのように特異ベクトルをわざわざlatlon格子に戻す必要はないです）。

ただし，MCAではSVD時系列が2つありますので，回帰図も2種類考えることができます。

同種回帰図(homogeneous regression map)： 「SSTのSVD1_SSTの時系列に対する回帰図」，「SLPのSVD1_SLPの時系列に対する回帰図」のように，同じ変数同士で回帰させる回帰図のこと。

異種回帰図(heterogeneous regression map)： 「SSTのSVD1_SLPの時系列に対する回帰図」，「SLPのSVD1_SSTの時系列に対する回帰図」のように，SVDした相手側の変数で回帰させる回帰図のこと。

同種回帰図は，特異ベクトルから得られた空間構造と同じパターンを示しますので，EOFのときと同様に振幅付SVDだと思えば良いです。このパターンにどれだけ意味があるかを調べるために，異種回帰図を作成します。MCAは，共分散が最大のパターンでさえあれば，どんなに小さなシグナルでも検出してしまいます。 異種回帰図は相手側の変数で回帰しますので，もし他の変動に埋もれてしまうような小さなシグナルだと，同種回帰図と同じようなパターンを得ることはできません。それゆえ，異種回帰図をつくることで「本当に注目すべき共通変動であるか」を確かめることができます。

つまり，同種と異種の2つの回帰図が似たようなパターンになれば，得られたモードはデータに内在する十分卓越した変動であることがわかり，MCAパターンが統計的にたまたま取り出されたノイズではなく，物理的な意味のあるものである可能性が高まります。

（それぞれのEOFで似たようなパターンが出るかどうか，というのも良い確認方法です）

では早速，2つの回帰図を描画してみましょう。

In [12]:
# 同種回帰図
print('Calculating homogeneous regression maps...')
a_sst_on_svd1_sst = reg_map(svd1_sst_time, ssta)
a_slp_on_svd1_slp = reg_map(svd1_slp_time, slpa)
# 異種回帰図
print('Calculating heterogeneous regression maps...')
a_sst_on_svd1_slp = reg_map(svd1_slp_time, ssta)
a_slp_on_svd1_sst = reg_map(svd1_sst_time, slpa)
# 描画
draw_two_fields(a_sst_on_svd1_sst, lon2_sst, lat2_sst, \
                a_slp_on_svd1_slp, lon2_slp, lat2_slp, \
                'Homogeneous Regression map of SVD1', vmin = -1.2, vmax = 1.2, vint = 0.1)
draw_two_fields(a_sst_on_svd1_slp, lon2_sst, lat2_sst, \
                a_slp_on_svd1_sst, lon2_slp, lat2_slp, \
                'Heterogeneous Regression map of SVD1', vmin = -1.2, vmax = 1.2, vint = 0.1)

Calculating homogeneous regression maps...
0
100
200
300
0
100
Calculating heterogeneous regression maps...
0
100
200
300
0
100

二つの回帰図が似たようなパターンになっているので，この共通変動パターンに注目する価値はありそうです。MCA大成功と言っていいでしょう！

課題8（締切：次回の授業開始前まで）

この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。

A問題. 最大共分散分析（MCA）とは何か，簡潔に説明してください。

B問題. SVD第2モードについて，SVD2時系列と同種回帰図を計算し，どのような現象を示しているか簡単に説明してください（難しく考えず，見たまま書いてくれればそれで良いです）。また，同種回帰図と異種回帰図は同じようなパターンになっているでしょうか？

C問題. 熱帯太平洋のSSTと熱帯全体のSLPの共分散が最大になるパターンを調べてみようと思ったとします。領域をそれぞれ自由に定めてMCAを行い，自由に考察してください。

D問題. 拙著 Kohyama et al. (2021) https://www.science.org/doi/10.1126/science.abh3295 （大学のネットワーク等でのみアクセス可）のFig. 2では，どのようにMCAが用いられていますか？

提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論8，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。