Ex.1 演習ノート ―Volveフィールドデータで学ぶ―¶

実践 貯留層工学 データ駆動型シミュレータ編 第1回 演習(Ex.1)

本ノートは データ駆動型シミュレータ編_第1回_演習.pptx の構成(環境構築→データ読込→前処理→可視化→特徴量整理)に沿って、 調査部門/参考資料/Volve/ に格納したEquinor Volveフィールドデータ(演習用抜粋)を実際に読み込み・可視化した記録である。

使用データ: Volve_BHP_15_9-F-1C.csv / Volve_BHFP_data.csv / Volve_Production_Data_Processed.csv

Step1: 環境構築確認¶

In [1]:
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Yu Gothic"
plt.rcParams["axes.unicode_minus"] = False

print("pandas", pd.__version__, "/ numpy", np.__version__, "/ matplotlib", matplotlib.__version__)

DATA_DIR = r"C:\kawata\testAI\調査部門\参考資料\Volve"

pandas 3.0.3 / numpy 2.5.0 / matplotlib 3.11.1

Step2: データ読込・前処理・可視化(1坑井から始める)¶

まずVolve_BHP_15_9-F-1C.csv(坑井15/9-F-1C単独、746行)のような小さいファイルから読み込み、操作に慣れる。

In [2]:
df1 = pd.read_csv(f"{DATA_DIR}/Volve_BHP_15_9-F-1C.csv")
print(df1.shape)
df1.head()

(746, 5)
Out[2]:
      DATEPRD NPD_WELL_BORE_NAME  ON_STREAM_HRS  AVG_DOWNHOLE_PRESSURE  AVG_DOWNHOLE_TEMPERATURE
0  07/04/2014         15/9-F-1 C           0.00                  0.000                     0.000
1  08/04/2014         15/9-F-1 C           0.00                    NaN                       NaN
2  09/04/2014         15/9-F-1 C           0.00                    NaN                       NaN
3  10/04/2014         15/9-F-1 C           0.00                    NaN                       NaN
4  11/04/2014         15/9-F-1 C          24.00                233.xxx                    96.876
In [3]:
df1["DATEPRD"] = pd.to_datetime(df1["DATEPRD"], dayfirst=True)
print(df1["DATEPRD"].min(), "〜", df1["DATEPRD"].max())
print(df1.isna().sum())

2014-04-07 00:00:00 〜 2016-04-21 00:00:00
DATEPRD                     0
NPD_WELL_BORE_NAME          0
ON_STREAM_HRS               0
AVG_DOWNHOLE_PRESSURE       3
AVG_DOWNHOLE_TEMPERATURE    3
dtype: int64
In [4]:
fig, ax = plt.subplots(figsize=(10, 4.5))
ax.plot(df1["DATEPRD"], df1["AVG_DOWNHOLE_PRESSURE"], color="#2A78D6", linewidth=1.2)
ax.set_title("坑井 15/9-F-1C: 坑底流動圧の推移")
ax.set_xlabel("日付"); ax.set_ylabel("坑底流動圧 (AVG_DOWNHOLE_PRESSURE)")
ax.grid(alpha=0.3)
fig.tight_layout()
plt.show()
No description has been provided for this image

気づき: 2016年4月末の最終行付近で値が0まで急落している。物理的にありえない挙動であり、 坑井のシャットイン(停止)または測定不良によるものと考えられる。実運用ではこの種の異常値を 特徴量整理の前に除外・フラグ付けする必要がある。

Step3: 7坑井全体のデータで欠損・傾向を確認する¶

In [5]:
df_all = pd.read_csv(f"{DATA_DIR}/Volve_BHFP_data.csv")
df_all["DATEPRD"] = pd.to_datetime(df_all["DATEPRD"], dayfirst=True)
print(df_all.shape)
print(df_all["NPD_WELL_BORE_NAME"].unique().tolist())

(15634, 5)
['15/9-F-1 C', '15/9-F-11', '15/9-F-12', '15/9-F-14', '15/9-F-15 D', '15/9-F-4', '15/9-F-5']
In [6]:
missing_by_well = df_all.groupby("NPD_WELL_BORE_NAME")["AVG_DOWNHOLE_PRESSURE"].apply(lambda s: s.isna().mean())
n_by_well = df_all.groupby("NPD_WELL_BORE_NAME").size()
summary_missing = pd.DataFrame({"行数": n_by_well, "欠損率": missing_by_well}).sort_values("欠損率", ascending=False)
print(summary_missing)
print("\n全体欠損率:", df_all["AVG_DOWNHOLE_PRESSURE"].isna().mean())

                      行数       欠損率
NPD_WELL_BORE_NAME                
15/9-F-5            3306  1.000000
15/9-F-4            3327  1.000000
15/9-F-11           1165  0.005150
15/9-F-1 C           746  0.004021
15/9-F-12           3056  0.001963
15/9-F-14           3056  0.001963
15/9-F-15 D          978  0.000000

全体欠損率: 0.42561084815146477

気づき(重要): pptx資料では「AVG_DOWNHOLE_PRESSURE列は全体の約42.6%が欠損」とだけ記載していたが、 実データを確認すると欠損はランダムではなく、坑井 15/9-F-4・15/9-F-5 の2本(水圧入井)で100%欠損していることが判明した。 残り5本の生産井ではほぼ欠損なし(0〜0.5%程度)。つまり「42.6%欠損」の実態は、 「坑底圧力ゲージが付いていない坑井が存在する」という坑井属性の問題であり、単純な補間では対処できない。 特徴量整理の際はこの2坑井を圧力系の特徴量から除外するか、稼働時間等の別変数で代替する必要がある。

In [7]:
fig, ax = plt.subplots(figsize=(9, 4.5))
summary_missing["欠損率"].plot(kind="bar", ax=ax, color="#EB6834")
ax.set_title("坑井別: 坑底圧力データの欠損率")
ax.set_ylabel("欠損率"); ax.set_xlabel("坑井名")
for label in ax.get_xticklabels():
    label.set_rotation(30); label.set_ha("right")
fig.tight_layout()
plt.show()
No description has been provided for this image
In [8]:
fig, ax = plt.subplots(figsize=(10, 5))
colors = ["#2A78D6", "#EB6834", "#1BAF7A", "#EDA100", "#E87BA4", "#008300", "#4A3AA7"]
for i, (well, g) in enumerate(df_all.groupby("NPD_WELL_BORE_NAME")):
    g_valid = g.dropna(subset=["AVG_DOWNHOLE_PRESSURE"])
    ax.plot(g_valid["DATEPRD"], g_valid["AVG_DOWNHOLE_PRESSURE"], label=well,
            color=colors[i % len(colors)], linewidth=1.0, alpha=0.85)
ax.set_title("坑井別 坑底流動圧の推移(7坑井重ね描き)")
ax.set_xlabel("日付"); ax.set_ylabel("坑底流動圧")
ax.legend(fontsize=8, ncol=2); ax.grid(alpha=0.3)
fig.tight_layout()
plt.show()
No description has been provided for this image

気づき: dropna後も一部の坑井で圧力が文字通り0まで落ちる期間がある(例: 15/9-F-12は2010年後半以降、 15/9-F-14は複数回)。NaNではなく数値としての0が入っているため、isna()だけでは検出できない。 探索的データ分析(EDA)では「NaN」と「物理的にありえない0」の両方を欠損・異常値として扱う必要がある、 という実データならではの教訓が得られた。

Step4: 生産量データ(油田全体)の読込・可視化¶

前処理の落とし穴: Volve_Production_Data_Processed.csvの日付は DD/MM/YYYY 表記(例: 01/09/2007 = 2007年9月1日)。 dayfirst=True を指定せずに読み込むと、月次データが疑似的な日次データとして誤解釈され、 時系列プロットが破綻する(実際に最初はこの誤りをそのまま可視化し、崩れたグラフになった)。

In [9]:
df_prod = pd.read_csv(f"{DATA_DIR}/Volve_Production_Data_Processed.csv")
df_prod["Date"] = pd.to_datetime(df_prod["Date"], dayfirst=True)  # dayfirst=True が必須
print(df_prod.shape)
print("monotonic:", df_prod["Date"].is_monotonic_increasing)
df_prod.head()

(112, 8)
monotonic: True
Out[9]:
        Date  p (psia)  Np (STB)  Gp (SCF)  Wp (STB)  Gi (SCF)  Wi (STB)  Rp (SCF/STB)
0 2007-09-01   4780.59         0         0       0.0         0         0           0.0
1 2007-10-01   4780.59         0         0       0.0         0         0           0.0
2 2007-11-01   4780.59         0         0       0.0         0         0           0.0
3 2007-12-01   4780.59         0         0       0.0         0         0           0.0
4 2008-01-01   4780.59         0         0       0.0         0         0           0.0
In [10]:
fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
axes[0].plot(df_prod["Date"], df_prod["p (psia)"], color="#2A78D6")
axes[0].set_ylabel("貯留層圧力 p (psia)")
axes[0].set_title("Volve: 貯留層圧力・累積生産量の推移(月次)")
axes[0].grid(alpha=0.3)

axes[1].plot(df_prod["Date"], df_prod["Np (STB)"], color="#1BAF7A", label="Np (累積油量)")
ax2 = axes[1].twinx()
ax2.plot(df_prod["Date"], df_prod["Wp (STB)"], color="#E87BA4", label="Wp (累積水量)")
axes[1].set_ylabel("Np (STB)", color="#1BAF7A")
ax2.set_ylabel("Wp (STB)", color="#E87BA4")
axes[1].set_xlabel("日付")
axes[1].grid(alpha=0.3)
fig.tight_layout()
plt.show()
No description has been provided for this image

読み取れる傾向: 2007年後半の生産開始直後は圧力が急低下し(初期ドローダウン)、 2010年頃から水圧入(Wi)の効果で圧力が回復・安定していく典型的な水圧入油田の挙動が確認できる。 累積油量Np(緑)はS字カーブ的に増加し、2013年以降は累積水量Wp(桃)の増加が加速している (水攻進行=ウォーターブレイクスルーの兆候)。

Step5: 特徴量整理(坑井別サマリテーブル)¶

In [11]:
feat = df_all.groupby("NPD_WELL_BORE_NAME").agg(
    観測日数=("DATEPRD", "count"),
    開始日=("DATEPRD", "min"),
    終了日=("DATEPRD", "max"),
    平均稼働時間=("ON_STREAM_HRS", "mean"),
    平均坑底圧力=("AVG_DOWNHOLE_PRESSURE", "mean"),
    坑底圧力欠損率=("AVG_DOWNHOLE_PRESSURE", lambda s: s.isna().mean()),
)
feat[["平均稼働時間", "平均坑底圧力", "坑底圧力欠損率"]] = feat[["平均稼働時間", "平均坑底圧力", "坑底圧力欠損率"]].round(2)
feat
Out[11]:
                    観測日数        開始日        終了日  平均稼働時間  平均坑底圧力  坑底圧力欠損率
NPD_WELL_BORE_NAME
15/9-F-1 C           746 2014-04-07 2016-04-21   13.38  246.67     0.00
15/9-F-11           1165 2013-07-08 2016-09-17   22.32  233.96     0.01
15/9-F-12           3056 2008-02-12 2016-09-17   21.34   80.73     0.00
15/9-F-14           3056 2008-02-12 2016-09-17   20.54  233.07     0.00
15/9-F-15 D          978 2014-01-12 2016-09-17   18.23  226.03     0.00
15/9-F-4            3327 2007-09-01 2016-12-01   20.24     NaN     1.00
15/9-F-5            3306 2007-09-01 2016-09-18   19.17     NaN     1.00
In [12]:
feat.to_csv("well_feature_table.csv", encoding="utf-8-sig")
print("saved: well_feature_table.csv")

saved: well_feature_table.csv

まとめ・第2回への申し送り¶

  • 環境構築(pandas/numpy/matplotlib)からデータ読込・前処理・可視化・特徴量整理まで、pptx資料(E1〜E15)の流れに沿って一通り実行できた。
  • 実データ特有のつまずきが3点見つかった。第2回の教材化・受講者への注意喚起に反映すべき:
    1. Volve_Production_Data_Processed.csvの日付はDD/MM/YYYY表記であり、dayfirst=Trueを付け忘れると月次データが壊れて見える(実際に一度失敗して気づいた)。
    2. 坑井15/9-F-4・15/9-F-5(水圧入井)は坑底圧力データが100%欠損。「欠損率42.6%」は坑井属性による構造的な欠損であり、ランダム欠損ではない。
    3. 欠損はNaNだけでなく、物理的にありえない「0」としても現れる(坑井停止区間等)。isna()だけでは検出できないため、EDAで実際にプロットして目視確認する価値がある。
  • 坑井別特徴量テーブル(well_feature_table.csv)を作成した。第2回『代理モデル構築の実践』では、この観測日数・稼働時間・坑底圧力の分布を入力として扱う想定。