AKARI Tech Blog

燈株式会社のエンジニア・開発メンバーによる技術ブログです

トップ > AI・機械学習 > オープンデータからガウス過程回帰で地盤の硬さを予測してみた

オープンデータからガウス過程回帰で地盤の硬さを予測してみた

AI・機械学習 Light up Hack!

こんにちは、DX Solution 事業本部 Devの岡内です。 今週のAKARI Tech Blogを担当いたします！

本記事では、ガウス過程回帰（GPR: Gaussian Process Regression） を用いて、既存の地盤調査データから未知の地点における地盤の硬さ（N値）を推定するシステムの構築について解説します。地盤情報は建設やインフラ整備において極めて重要ですが、調査にはコストと時間がかかります。このシステムは、機械学習モデルによって効率的な地盤情報の推定と可視化を目指すものです。

システムの核となる技術、具体的な実装、構築したシステムの概要と予測結果、そして今後の展望について順を追って説明します。

ガウス過程回帰（GPR）の理論的背景

本システムの核となる技術が、ガウス過程回帰（GPR） です。GPRは機械学習の回帰手法の一つであり、大きな特徴として予測値とその不確実性（予測の信頼区間）を同時に定量評価できる点が挙げられます。これはベイジアン的なアプローチに由来するものです。

端的に言えば、「観測データ点（入力と出力のペア）を滑らかに補間する関数を、その関数が生成されうる確率（事後確率）とともに推定する」手法です。

ガウス過程：関数の分布

一般的な統計モデルが「データの分布」を仮定するのに対し、ガウス過程では「関数の分布」を考えます。これは、どのような形状の関数がどの程度の確率で出現するか、という確率分布を直接モデル化するアプローチです。この「関数の分布」は、以下の2つの関数によって完全に定義されます。

  • 平均関数 m(x): 関数の大域的なトレンドを表します。事前知識がない場合は、通常 m(x)=0 と仮定されます。
  • 共分散関数（カーネル関数）k(x,x′): 2つの入力点 x と x′ における関数値の相関の強さを示します。この関数が、生成される関数の滑らかさや周期性といった特性を決定します。

予測プロセス：事前分布から事後分布へ

GPRによる予測は、ベイズの定理に基づき、観測データを用いて「関数の分布」を更新していくプロセスとして解釈できます。

  1. 事前分布（データ観測前） 観測データが存在しない初期状態では、関数の分布はガウス過程（事前分布）に従います。この時点では、想定される関数は無限に存在し、下図の上側のように様々な形状の関数がサンプリングされます。
  2. 事後分布（データ観測後） 次に、いくつかの観測データ点を与えると、ベイズの定理によって、これらのデータ点を通過する、あるいはその近傍を通過する可能性が高い関数群に分布が更新されます。これが事後分布です。下の図を見ると、観測データ点の周辺では、関数の取りうる範囲（灰色の帯、信頼区間）が狭くなっていることが確認できます。これは、データによって局所的な関数の形状が特定され、不確実性が減少したことを意味します。対照的に、データが存在しない領域では信頼区間は広がり、不確実性が大きいことを示唆します。

上図：事前分布からサンプリングされた関数。 下図：観測データ（黒点）を与えた後の事後分布。太い黒線が予測平均値、薄い灰色の帯が不確実性（標準偏差）を表す。 https://github.com/scikit-learn/scikit-learn/blob/99bf3d8e4/sklearn/gaussian_process/kernels.py#L1445

このように予測の信頼度を空間的に評価できる点が、GPRを実用上非常に有用なものにしています。

実装：GPyTorchによるモデル定義

理論的な背景に続いて、実際にこのシステムで用いたGPRモデルの実装の一部をご紹介します。モデルの構築には、PyTorchベースのGPRライブラリであるGPyTorchを利用しました。GPyTorchはGPUを活用した大規模データへのスケーラブルな近似推論をサポートしており、今回のタスクに適しています。

以下に示すのは、モデル定義と学習・予測の処理をまとめたクラスの概要です。

import gpytorch
import joblib
import numpy as np
import torch
from sklearn.preprocessing import StandardScaler
from torch.utils.data import DataLoader, TensorDataset

# デバイスとデータ型を設定
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
dtype = torch.double
torch.set_default_dtype(dtype)

# GPyTorchの近似GPモデル定義
class ApproximateGPModel(gpytorch.models.ApproximateGP):
    def __init__(self, inducing_points, ard_num_dims):
        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(
            num_inducing_points=inducing_points.size(0), dtype=inducing_points.dtype
        )
        variational_strategy = gpytorch.variational.VariationalStrategy(
            self, inducing_points, variational_distribution, learn_inducing_locations=True
        )
        super().__init__(variational_strategy)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel(ard_num_dims=ard_num_dims)
        )

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

# N値予測のワークフローを管理するクラス
class NValuePredictor:
    def __init__(self, num_inducing_points=1000, batch_size=1024, num_epochs=50):
        self.model = None
        self.likelihood = gpytorch.likelihoods.GaussianLikelihood()
        self.scaler = StandardScaler()
        self.num_inducing_points = num_inducing_points
        self.batch_size = batch_size
        self.num_epochs = num_epochs
        self.feature_columns = None

    def train(self, X: np.ndarray, y: np.ndarray, feature_columns: list):
        """モデルを学習させる (標準化 + ミニバッチ学習)"""
        X_scaled = self.scaler.fit_transform(X)
        self.feature_columns = feature_columns
        
        X_tensor = torch.from_numpy(X_scaled).to(device)
        y_tensor = torch.from_numpy(y).to(device)

        train_loader = DataLoader(TensorDataset(X_tensor, y_tensor), batch_size=self.batch_size, shuffle=True)

        inducing_points_idx = np.random.choice(len(X_tensor), self.num_inducing_points, replace=False)
        inducing_points = X_tensor[inducing_points_idx, :]
        
        self.model = ApproximateGPModel(inducing_points=inducing_points, ard_num_dims=X_tensor.shape[1]).to(device)
        self.likelihood.to(device)

        optimizer = torch.optim.Adam(list(self.model.parameters()) + list(self.likelihood.parameters()), lr=0.01)
        mll = gpytorch.mlls.VariationalELBO(self.likelihood, self.model, num_data=y_tensor.size(0))

        self.model.train()
        self.likelihood.train()
        for i in range(self.num_epochs):
            for x_batch, y_batch in train_loader:
                optimizer.zero_grad()
                output = self.model(x_batch)
                loss = -mll(output, y_batch)
                loss.backward()
                optimizer.step()

    def predict(self, X_new: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        """新しい地点のN値を予測する"""
        X_new_scaled = self.scaler.transform(X_new)
        X_new_tensor = torch.from_numpy(X_new_scaled).to(device)
        
        self.model.eval()
        self.likelihood.eval()
        with torch.no_grad(), gpytorch.settings.fast_pred_var():
            preds = self.likelihood(self.model(X_new_tensor))
        
        return preds.mean.cpu().numpy(), preds.stddev.cpu().numpy()

    def save(self, model_path, scaler_path, features_path):
        """モデル、スケーラー、特徴量リストを保存する"""
        torch.save(self.model.state_dict(), model_path)
        joblib.dump(self.scaler, scaler_path)
        joblib.dump(self.feature_columns, features_path)

    @classmethod
    def load(cls, model_path, scaler_path, features_path):
        """学習済みモデルからインスタンスを生成して返す"""
        predictor = cls()
        predictor.feature_columns = joblib.load(features_path)
        
        ard_num_dims = len(predictor.feature_columns)
        inducing_points = torch.zeros(predictor.num_inducing_points, ard_num_dims, device=device)
        
        predictor.model = ApproximateGPModel(inducing_points=inducing_points, ard_num_dims=ard_num_dims)
        predictor.model.load_state_dict(torch.load(model_path, map_location=device))
        predictor.model.to(device)
        
        predictor.scaler = joblib.load(scaler_path)
        return predictor

このコードでは、主に2つのクラスを定義しています。

  • ApproximateGPModel: GPyTorchの近似GPモデルを継承したクラスです。大規模なデータセットに対応するため、誘導点（inducing points）を用いた変分推論（Variational Inference）を実装しています。カーネルには、一般的な選択肢であるRBFカーネルを採用し、各特徴量の重要度を学習できるようARD（Automatic Relevance Determination）を有効にしています。
  • NValuePredictor: モデルの学習、予測、保存、読み込みといった一連のワークフローを管理するラッパークラスです。trainメソッド内では、scikit-learnのStandardScalerによるデータの前処理、DataLoaderを用いたミニバッチ学習が行われます。predictメソッドでは、学習済みモデルを用いて新しいデータポイントに対するN値の平均と標準偏差（不確実性）を算出します。

システムの概要と予測結果

このGPRの原理を応用し、地盤のN値を予測するシステムを構築しました。

処理フロー

システムは以下の手順で予測を実行します。

  1. 特徴量収集: 予測対象地点の周囲から、既存のボーリングデータを収集します。
  2. データクレンジング: 収集したデータを整形します（例：超硬質地盤など特定の値の補完）。
  3. GPRモデル作成: 収集・整形したデータを基に、その都度GPRモデルを学習します。
  4. N値予測: ユーザーが指定した地点や断面に沿って、連続的に地下のN値を予測します。
  5. 結果の可視化: 予測結果を統合し、断面図としてヒートマップで可視化します。

予測結果

本システムにより、任意の地点におけるN値の深度方向の分布を予測することができます。下図は、ある1地点を対象とした予測結果の一例です。実線がN値の予測平均値を、薄い青色の帯が予測の標準偏差（±1σ）の範囲（信頼区間）を示しており、深さごとの予測値とその不確実性を視覚的に捉えることができます。

さらに、指定したラインに沿って連続的に予測を行うことで、地盤のN値断面図をヒートマップとして可視化することも可能です。これにより、地層の構造をより直感的に把握できます。

予測精度

現状のモデル精度は以下の通りです。改善の余地はありますが、一定の相関関係が確認でき、実用化に向けた基礎的な検証としては有意義な結果が得られました。

  • R2スコア (決定係数): 0.5904
  • MAE (平均絶対誤差): 7.6154
  • RMSE (二乗平均平方根誤差): 10.2706

予測に使用したデータ

本モデルの学習と予測には、主に国が公開している以下のオープンデータを利用しました。

ボーリングデータとN値

ボーリング調査は、地面を削孔して地中の土壌や岩石サンプルを採取する地盤調査手法です。その結果は「ボーリング柱状図」としてまとめられます。 N値（標準貫入試験打撃回数）は、このボーリング調査中に行われる標準貫入試験から得られる値で、地盤の硬軟や締まり具合を示す重要な指標です。N値が大きいほど、地盤が硬いことを意味します。

ボーリング柱状図の例。（https://www.kunijiban.pwri.go.jp/viewer/column/?id=149912）

国土地盤情報検索サイト「Kunijiban」

Kunijibanは、国土交通省が管理・公開している日本全国のボーリングデータを集約したデータベースです。誰でも無料でアクセスでき、地図上からボーリング柱状図データをダウンロードすることが可能です。今回は、本サイトの利用規約に準拠し、検証を目的として同サイトで公開されているボーリング柱状図データを活用させていただきました。

なお、本記事で示した検証では、埼玉県内の荒川周辺を対象としております。この地域はボーリングデータが豊富に公開されており、モデルの検証に適していると判断したためです。

国土地盤情報検索サイト「Kunijiban」のインターフェース。（https://www.kunijiban.pwri.go.jp/viewer/）

まとめと今後の展望

本稿では、オープンデータを活用し、ガウス過程回帰を用いて任意の地点・断面における地盤の硬さ(N値)を、その予測信頼度とともに可視化するシステムについて述べました。現状の予測精度（R2スコア: 0.59）は発展途上ですが、GPRの利点である「不確実性」の評価機能は、効率的な地盤調査計画の立案支援等への応用可能性を示唆します。

しかし、このアプローチはまだ発展の途上にあります。今後の展望として、先行研究で示されているような高度な空間推定手法を取り入れ、モデルをさらに進化させることを目指します。

間接データを活用した空間分布推定への展開

現在のモデルは、各地点のボーリング調査データを基にN値を推定する単一回帰の考え方に基づいています。しかし、地盤調査はコストや工期の制約から、限定された地点でしか実施できないのが実情です。そこで、より精緻な空間分布を推定するために、先行研究で提案されているアプローチを参考にします。

津田ほか(2025)の研究では、推定対象の物性値である「直接データ」（論文中ではN値）と、それと相関を有する別の種類の計測データである「間接データ」（論文中では杭施工時のトルク値）をGPRの枠組みで統合する手法が提案されています（ https://www.jstage.jst.go.jp/article/jscejj/81/15/81_24-15018/_article/-char/ja/）。この研究では、間接データを活用することで、直接データがない領域でもN値の空間分布をより正確に推定できることが示されています。

この考え方を我々のシステムにも応用します。

  • 直接データ: これまで通り、Kunijibanから取得するボーリング調査のN値。これは高精度ですが、空間的に疎（スパース）なデータです。
  • 間接データ: 今後の拡張として、防災科学技術研究所のJ-SHISが提供する全国の平均S波速度（AVS30）の活用を計画しています。AVS30は地盤の硬さに関連する指標であり、N値と一定の相関が期待されます。何より、これは日本全国をカバーする面的（連続的）なデータです（https://www.j-shis.bosai.go.jp/api-list）。

https://www.j-shis.bosai.go.jp/map/

AVS30という密な間接データを組み込むことで、N値という疎な直接データだけでは捉えきれなかった地盤の広域的な傾向をモデルに与えることができます。これにより、以下のような改善が期待されます。

  1. 推定精度の向上: N値の観測点がない地点においても、周辺のAVS30の情報を活用することで、より実現象に即した推定が可能になります。
  2. 不確実性の低減: 間接データによって空間的な拘束が加わるため、予測の信頼区間が狭まり、より信頼性の高い推定が期待できます。
  3. 詳細な空間変動の再現: 論文で示されているように、支持層のような急激な変化や、細かい地層の変動を捉える能力が向上し、より詳細な3次元地盤モデルの構築に繋がります。

この実現には、直接データと間接データの相互相関を考慮したGPRの定式化が必要となります。また、N値とAVS30では空間的な相関の度合い（SoF: Scale of Fluctuation）が異なる可能性があるため、論文で示されているように、異なるSoFを持つデータ同士を柔軟に扱えるモデルの構築が鍵となります。

これからは、点と点で予測するだけじゃなく、面や空間として地盤を捉えられるようにシステムをパワーアップさせていきます。より実務で使える、精度の高い地盤予測システムを完成させるのが目標です！

本取り組みの背景：Light up Hack!

実は、今回ご紹介したシステム開発は、社内の技術コンペ「Light up Hack!」での取り組みがきっかけとなっています。これは、私たちDX Solution 事業本部 Devチームのエンジニアが個人やチームを組み、業務の枠にとらわれずに自由な発想でプロダクトを開発するイベントです。

「Light up Hack!」は、新しい技術への挑戦や、普段の業務では扱わない領域の課題解決を通じて、技術力を高め、互いに学び合う貴重な機会となっています。今回のプロジェクトも、そうしたカルチャーから生まれた成果の一つです。

We’re Hiring!

株式会社AKARIでは、AI・機械学習技術を用いて社会課題の解決に取り組むエンジニアを募集しています。 ご興味をお持ちいただけましたら、下記採用ページよりご連絡ください。

akariinc.co.jp