【SIGNATE】ひろしまQuest2022：河川の水位予測

今回は、SIGNATEのデータ分析コンペ「ひろしまQuest2022：河川の水位予測」のデータを使って河川水予測モデルを作成しました。
手法としてはNearest-Neighbor法（k近傍法）を使用し、当日の24時間分の降水量や河川水位のデータから翌日の河川水位を予測する簡単なモデルを作成してみました。

コンペについて

【プラットフォーム】SIGNATE
【コンペ名】ひろしまQuest2022：河川の水位予測　広島県にある河川の水位の予測にチャレンジしよう！
【主催者】ひろしまサンドボックス推進協議会事務局
【期間】2022年11月15日(火)～2023年1月26日(木)
【情報公開ポリシー】モデル・分析結果ともにコンペ終了後に非営利目的での公開可
【概要】広島県内の河川の１時間ごとの観測データ（河川水位、雨量、潮位）から翌日（24時間分）の河川水位を予測
【評価指標】RMSE

データ

データとしては、①河川水位、②雨量、③潮位の3種類の時系列データ・観測地点情報に加え、④ダム情報のマートがあります。

①～③のデータは１日分（１時間ごとに24カラムのデータ）が横持ちで並んでいたので驚きました（すぐに縦持ちへ変換）。
①～③はそれぞれ上記のような観測値を記載したファイルと観測地点の情報（水系名や河川名、住所、経緯度情報など）が記載されたマートの２つのCSVからなります。

同じ河川でも観測地点は①～③でバラバラで、それぞれ全く別の場所で観測しています。
コンペのチュートリアルに記載の方法で地図上に可視化したところ、類似した地点名（たとえば、AA(国)とAA）でも別の観測場所でした。
潮位に関しては、河口付近の海岸のデータです。

県境をまたぐような河川については、広島県内の情報しかありません。
例えば江の川（ごうのかわ）は源流から広島県内区間のデータはありますが、県境を超えてから河口（島根県江津市）までのデータはありません。
ただし、中国山地の分水嶺は県境よりも広島県側にあるため、予測にはほとんど問題はありません。
上流の河川水位を予測するために下流のデータは不要です（逆は困る）。

③の潮位データに関しては、河口付近の水位変動は上流からの河川流量よりも海の潮の干潮の影響の方が大きいようです。
雨が降っていない日には、ほとんど潮の満ち引きだけで河川水位が変化するようです。
潮位データが影響するのは一部の河川の下流域（今回のデータでたまたま河口付近の観測地点がある場合）に限られます。
ただし、広島平野を流れる太田川では河口から数km上流でも潮位の影響があるようで、その点は注意が必要です。
河口付近以外では参照する意味がないデータなので、特定の観測地点の場合のみ条件分岐等で扱うデータになります。

④のダム情報は、わずか12地点のダムのデータで、ダムの情報（水系名や河川名、住所、経緯度情報、貯水容量など）のデータです。
広島県内でも有数の規模の太田川水系のダムは一箇所も記載がありません。
一方、河川水位や雨量のデータには太田川水系の観測地点も多数含まれます。
念のため調べてみましたが、当然、太田川水系にもダムが複数存在するため、かなり不完全なマートです。
また、ダムの大きさがわかっても貯水量の時系列データが無い限り河川水位予測には役に立たないため、何のために配布されているのかよくわからないデータです。

分析の流れ

ここからは分析の流れについて紹介します。

はじめにデータの基礎俯瞰を行った上で、先行研究の文献調査を行いました。
文献調査の結果を元に手法選択を行い、水系ごとにモデルを作成しました。
水系ごとに予測誤差の評価（RMSEの算出）を行い、モデルの精度評価を行いました。
ここで時間切れになってしまい、モデルの精度向上までは行えませんでした。

各ステップの詳細について、以下で説明します。

データの基礎俯瞰

基礎俯瞰では、先述の「データ」の項目に記載しているようなデータの特徴について確認した上で、水系ごとにデータの可視化を行いました。

河川水位や雨量の時系列データを可視化したところ、上流で河川水位が上昇したあとに、時間差で下流の河川水位が上昇しています。
ただし、上流と下流でほとんど差がない場合もありました。

一方、地理的に近い地点の雨量と河川水位のピークを比較したところ、雨量のピークが来たあとに遅れて河川水位のピークがやってきます。
河川水位は増水時に急に上昇し、減水時はゆるやかに低下していきます。
雨量はもっと極端に変動し、河川水位の増水時以上に急激に増加し、時間による上下動も激しいです。

このあたりの特徴は、事前知識が無くてもおおむね妥当だとわかるような事象ですが、データを把握する意味で文献調査の前に確認しました。

文献調査

次に、手法選択のために文献調査を行いました。
河川水位予測は学術的に研究されているので、やみくもに手法を考えるよりも過去の研究結果を前提にした方が分析の進め方として適切です。

理想的には英語の学術論文を中心に日本語・英語両方の文献を参照するべきです。
しかし、今回は
①そこまで時間をかけられない
②日本の河川は距離が短く勾配が急であり、他国の事例をそのまま適用できるとは限らない
という事情からから、日本語の文献を中心に参照しました。

文献調査の結果は、以下の記事にまとめています。

雨量から河川水位を予測する手法

はじめに、雨量から河川水位を予測する手法について調べました。
雨量から河川水位を予測する手法は大きく二つに分けらます。

一つは、流出モデルを用いた手法です。
降水量やダムの貯水量、河川流量などをモデル化して定量的に表現し、降水から海に到達するまでの水の流れを全て把握するような方法です。
この手法は学術的にはメジャーな手法ですが、河川の水深やダムの貯水量の時系列データなど、今回のデータセットには含まれていない多数の情報を必要とします。
そのため、今回は見送りです。

２つ目は、雨量から直接河川水位を予測する手法です。
降水が河川水位に反映されるまでには様々な過程を経るのですが、それらを飛ばして単純に過去の雨量と河川水位から未来の河川水位を予測する手法です。
算出根拠がブラックボックスになるため１つ目のアプローチよりも解釈可能性が落ちますが、今回のデータセットでも実装可能なアプローチです。

具体的な手法としては、Nearest-Neighbor法（k近傍法）やニューラルネットワークが使用されています。
Nearest-Neighbor法は2000年代を中心に文献が見つかるのに対し、ニューラルネットワークは2010年代中心に文献が見つかるため、k近傍法→ニューラルネットワークという流れで研究が進んでいるようです。

以上をふまえ、ひとまずNearest-Neighbor法で実装しました。
今回のように先行研究が豊富な事例に対しては、先行研究をふまえた手法を第一選択とした方がよいと考えています。
文献調査を後回しにすると何かあった時の手戻りが大変なので、早い段階で着手しておくと後からが楽になります。
Google Scholarにもあるように「巨人の肩の上に立つ」と遠くまで見渡せます。

感潮区間について

次に河口付近の水位の変動に潮位の影響があることから、感潮区間について調べました。
感潮区間とは、河口付近に位置するため海の潮の満ち引きの影響をうけて河川水位が変動するような区間のことです。
調査内容は次のページにまとめています。

感潮区間の潮位の影響をふまえる際も選択すべき手法は変わらず、Nearest-Neighbor法やニューラルネットワークを使うようです。
そのため、感潮区間の河川水位の予測の際には潮位の情報を加えるだけでよさそうです。

予測実施

今回使用するNearest-Neighbor法はk近傍法のk=1の場合の結果を利用した方法です。
まず、学習データとして使うために、過去データについて連続する24時間単位でデータを取得し、24時間単位の河川水位/雨量の推移パターンを学習させます。
一日単位で予測するのではなく、連続する24個（24時間分）のデータを読み枠として、１時間ずつ読み枠をずらしながら学習させていきます（スライド窓k近傍法）。

実際に予測したい一日の連続する24時間のデータに最も近い推移パターンを学習データの中から探します。
最も近い推移パターンを見つけたら、その推移パターンの次の24時間のデータを予測結果として使用します。

予測に使用した観測地点

予測自体は予測したい地点の河川水位データさえあれば可能です。
しかし、ある地点の水位は、過去にその河川の上流で降った雨量が影響します。
そのため、その地点の水位だけではなく、その水系全体の水位・雨量データを全て使って推移パターンを学習させました。
潮位に関しては、河口付近のみに影響するため、一旦無視して進めます。

同じ水系内であっても、その地点より下流の観測地点の水位や雨量は無関係です。
しかし、今回与えられた観測地点データは住所や経緯度情報があるものの、水系内での上流/下流の位置関係に関する情報はありません。

アイデアとしては、広島県内の河川はおおむね北から南に流れているので、緯度情報でソートして上流/下流を判定する方法などが考えられます。
しかしこの方法には問題もあり、①同じ水系でも合流前の支流上流部や②河川が蛇行してい部分、③日本海側に流れる江の川など単純に適用できない場合があります。

そのため、一旦上流/下流の関係性を無視して水系内の全ての水位・雨量データを予測に使用しました。

学習と予測

データは水系ごとに集約した上で、水系内の全観測地点の水位・雨量データを24時間単位で学習させました。
学習には全体のデータ（2191日分）の最初の８割を使用し、残り２割をテストデータとしました。

以下の学習モデル作成と予測のコードを記載します。

# 関数定義：連続するwidth個のデータのパターンを取得
def embed(df, width):
    cnt_cols = df.shape[1]
    dim = width*cnt_cols
    emb = np.empty((0, dim), float)
    for i in range(df.shape[0] - width +1):
        tmp = np.array(df[i:i+width])[::-1].reshape((1, -1))
        emb = np.append(emb, tmp, axis=0)
    return emb

# 関数定義：水系ごとにNN法にて予測実施
def peform_pred(suikei_name, width, nk):
    # 前処理済（水系別雨量＋水位、縦持ち変換済）データの読み込み
    file_path = 'data_' + suikei_name + '.csv'
    df_raw = pd.read_csv(file_path, encoding='cp932')

    # 訓練データとテストデータに分割
    train_ratio = 0.8
    train_test_boder = round(df_raw.shape[0]*train_ratio)
    
    df_train = df_raw.iloc[:train_test_boder]
    df_test = df_raw.iloc[train_test_boder:]

    df_train = df_train.set_index(['date', 'time'])
    df_test = df_test.set_index(['date', 'time'])

    # 使用する学習・テストデータ
    df_train_use = df_train.iloc[:-24]
    df_test_use = df_test.iloc[:-24]

    # 窓幅を使ってベクトルの集合を作成
    # 訓練データ
    train = embed(df_train_use, width=24)
    # テストデータ
    test = embed(df_test_use, width=24)

    # k近傍法実行
    neigh = NearestNeighbors(n_neighbors=nk)
    neigh.fit(train)
    # 学習データの中からテストデータの並びに最も近いパターンを取得
    # 距離を計算
    neigh_d, neigh_ind  = neigh.kneighbors(test)
    # 結果をdfに格納
    df_neigh = pd.DataFrame(neigh_d, columns=['distance'])
    df_neigh['neighbor_index'] = neigh_ind
    df_neigh['predict_index'] = df_neigh['neighbor_index'] + 24

    # knnの予測結果への結合用に加工
    # 結合用前処理
    #df_test_for_merge = df_test_use.reset_index().reset_index()
    df_test_for_merge = df_test.reset_index().reset_index()
    df_test_for_merge['index_tomorrow'] = df_test_for_merge['index'] + 24

    df_neigh_for_merge = df_neigh.reset_index()

    # knnの予測結果へテストデータの値を結合
    key_cols = 'index'
    df_neigh_merged = pd.merge(df_neigh_for_merge, df_test_for_merge, how='left', on=key_cols)

    # 学習データを結合用に加工
    df_train_reset = df_train.reset_index().reset_index()
    # カラムから水位データのみ取得
    obs_cols = df_train_use.columns.tolist()
    obs_cols = [col for col in obs_cols if (col.endswith('水位'))]

    # 水位について24時間分を横持ちで取得
    for obs_col in obs_cols:
        # 取得したい観測データ
        df_neigh_24 = df_neigh_merged.copy()

        for i in range(24):
            filled_i = str(i).zfill(2)
            # 0-23時までの水位データ取得
            # テストデータ
            df_neigh_24[filled_i + 'h_test'] = df_neigh_24['index_tomorrow'].apply(lambda x: df_test_for_merge[obs_col].loc[x+i])
            # テストデータに最近傍の学習データ
            df_neigh_24[filled_i + 'h_pred'] = df_neigh_24['predict_index'].apply(lambda x: df_train_reset[obs_col].loc[x+i])
        
        # 24時間で1行なので、時刻が0時のレコードのみ抽出
        df_neigh_daily = df_neigh_24[df_neigh_24['time']=='00:00:00']

        # 出力
        out_path = 'prediction_' + suikei_name + '_' + obs_col + '.csv'
        df_neigh_daily.to_csv(out_path, index=False, encoding='cp932')
    
    # 学習したモデルを返す
    return(train)

# 水系別にNN法にて予測実行
# 水系
suikei = '太田川'
# 窓幅
width = 24
# K近傍法のk
nk = 1

# 予測を実行して学習モデルを返す
trained_model = peform_pred(suikei_name=suikei, width=width, nk=nk)
# 学習モデルを水系ごとに出力
out_path = 'model_' + suikei + '.pickle'
with open(out_path, mode='wb') as f:  # with構文でファイルパスとバイナリ書き込みモードを設定
    pickle.dump(trained_model, f)

精度評価

最後にRMSEで精度評価を行いました。
全166の観測地点の水位の予測結果についてRMSEを求めたところ、RMSEの平均値は0.175でした。
RMSEの値の分布は下表のとおりです。

RMSEが大きい（=予測誤差が大きい）観測地点上位は、次のようになりました。

上位4地点はいずれも太田川河口のデルタ（広島平野）であり、いずれも上流の影響よりも潮の満ち引きの影響が大きい場所です。
今回のモデルは潮位データを使用していないので当然の結果です。
芦田川の大谷池（おおだにいけ）は名前のとおり、河川というよりも湖沼にあたる上、上流に位置するので下流域の水位や雨量データが役に立たなかったためと考えられます。
これらの地点はRMSEが大きく、全体の精度の足を引っ張っています。
今回のRMSEが大きい観測地点のうち河口付近の地点は潮位も考慮するようにモデルを変更することで、精度改善が見込まれます。

まとめ

今回はNN法を使って潮位を無視した簡単なモデルを作るところで終わりましたが、それでも一定の予測精度が出ました。
潮位の影響をうける感潮区間など誤差の原因がはっきりしている箇所の手直しをするだけで、更に精度改善が見込まれます。

加えて、今回は各観測地点の相対的な位置関係（上流側、下流側、別河川等）はデータが無いので無視しました。
それらを考慮して予測したい観測地点の上流側のデータだけを予測に使うといったことも改善点として挙げられます。
誤差が大きい河川上位には、太田川や江の川のような観測地点が多い河川が来ているため、予測に使用する観測地点の数を減らすことで精度向上が見込まれます。

参考文献

ひろしまQuest2022：河川の水位予測　広島県にある河川の水位の予測にチャレンジしよう！　SIGNATE　2023/2/4閲覧
【Python/Scikit-learn】k近傍法で時系列データの異常部位検出　2023/2/4閲覧