【画像認識】3次元地下画像の断層分布

地震探査による3次元地下画像内の断層分布の論文を見つけ、自前のPCで実装してみたので内容を以下にまとめます。なお、ここで載せている内容は、論文を参考にして休日に実装したものです。テスト用データもオープンソースのものを使用しています。

Jupyter Notebook
https://github.com/Jun-Tam/3D-Seismic-Image-Fault-Segmentation/blob/master/Fault_Detection.ipynb

地震探査とは?

油田・ガス田の探査に使われる技術で人工地震を発生させて、地下で反射した波を収録・処理して地下の画像を得る方法です。陸上ではダイナマイトやバイブロサイズと呼ばれる車が使われ、会場では巨大なエアガンで気泡を生じさせて、振動を起こします。

以下の図は、Googleの画像検索で適当に拾ってきたものですが、地下にある山のような構造の中に油が貯まる訳です。しかし、断層があると断層に沿って油やガスが地表付近へ抜ける事もあるので、断層分布の把握が重要となります。

従来は目で見ながら断層を解釈していましたが、近年は画像処理技術でよく使われるEdge Detectionに似た処理で断層を抽出する手法も一般的になっています。しかし、地層境界を示す反射面もエッジとして認識されるため、きれいな断層の分布を抽出するのはなかなか難しいです。

そこで、深層学習の適用が最近では研究トピックとなっているのですが、一番の問題は地下に本当にあるかどうか分からない断層に対して教師あり学習をしなければならないということです。手作業で断層のlabelingをするのは現実的ではなく、解釈者によっても異なるため、あまり良い方法ではありません。そこで、人工的に断層の入った地下の３次元画像を作り、断層の識別を学習させるというのが今回の試みです。

人工断層画像の作成

学習用データを人工的に作るのですが、実際の地層にできる限り近づけるため、以下のワークフローに従います。
ちなみに下の図は縦に切った断面図ですが、実際には３次元の画像です。

１．平行に広がった地層を仮定。地層の位置と硬さ（の差）はランダム。
２．２次元ガウシアン変形を与える。
３．平面変形を与える。(planar displacement)
４．断層を入れる。
５．ウェーブレットで畳み込む。
６．ノイズを加える。
７．中心部だけ抜き取る。
８．ゲインを上げる。

以上のワークフローをランダムなパラメータを与えながら必要な学習データ分作成します。
同時に断層を入れる際にラベルとして、断層の有無を格納したボリュームも作っていきます。

２，３だけ、CuPyを使って高速化しました。以下は２のCuPy使用箇所です。

def func_gauss2d(prm,refl,a,b,c,d,sigma):
    ''' Apply 2D Gaussian deformation '''
    xy_cp = cp.asarray(prm.xy)
    refl_cp = cp.asarray(refl)    
    a_cp = cp.asarray(a.astype('float64'))
    b_cp = cp.asarray(b.astype('float64'))
    c_cp = cp.asarray(c.astype('float64'))
    d_cp = cp.asarray(d.astype('float64'))
    sigma_cp = cp.asarray(sigma.astype('float64'))    
    z_cp = cp.asarray(prm.z)

    # Parallelize computation on GPU using cupy
    func_gauss2d = cp.ElementwiseKernel(
            in_params='T x, T y, T b, T c, T d, T sigma',
            out_params='T z',
            operation=''' z = b*expf(-(powf(x-c,2) + powf(y-d,2))/(2*powf(sigma,2))); ''',
            name='func_gauss2d'
            )

    gauss_2d_cp = cp.zeros_like(xy_cp[:,0])            
    for i in range(len(b)):
        gauss_2d_cp += func_gauss2d(xy_cp[:,0],xy_cp[:,1],b_cp[i],c_cp[i],d_cp[i],sigma_cp[i])
    s1_cp = a_cp +(1.5/z_cp)*cp.outer(cp.transpose(gauss_2d_cp),z_cp)

    for i in range(prm.nxy_tr):
        s = s1_cp[i,:] + z_cp
        mat = cp.tile(z_cp,(len(s),1)) - cp.tile(cp.expand_dims(s,1),(1,len(z_cp)))
        refl_cp[i,:] = cp.dot(refl_cp[i,:], cp.sinc(mat))

    refl = np.reshape(cp.asnumpy(refl_cp), [prm.nxy_tr, prm.nz_tr])    
    return refl

上の図では、断層が１つしか入っていませんが、学習を効率化させるため、実際には一つの画像に断層を５つ以上入れています。以下の図は出来上がった画像の例です。

このような３次元画像をトレーニング用200個、バリデーション用20個作成しました。

モデルの学習

モデルは層数を減らしたU-Netをkerasで実装し、GeForce GTX 1080 Tiを使って学習させました。
メモリ容量は11GBありますが、3次元のデータとモデルを扱うためメモリーが不足し、バッチ数は３までしか上げられませんでした。パッチサイズは 128 x 128 x 128、Data augmentationのため、学習データをランダムで回転（0, 90, 180, 270 deg)させています。

デモ

学習させたモデルを実データに適用し、スライスを表示したデモを作成しました。テスト用データは公開されているオランダの海上にあるF3という鉱区のデータです。

テストデータ
Netherlands Offshore F3 Block (dGB Open Seismic Repository)
https://terranubis.com/datainfo/Netherlands-Offshore-F3-Block-Complete

GPUのメモリを上げてバッチ数を増やせれば、もう少しはっきりとした分布が見えるかも知れません。

参照論文

Using synthetic data sets to train an end-to-end convolutional neural network for 3D seismic fault segmentation, Xinming Wu et al., Geophysics, 2019