見出し画像

インバージョンの基礎理論④非線形最小二乗法(数式あり)【私の備忘録】

みなさん,こんにちは.とある地学屋の日常ことゆうすけです.
今回は,インバージョンの基礎理論シリーズの最終回となります.前回は,非線形最小二乗法の流れについて説明しました.今回は,数式を用いて非線形最小二乗法について説明していきます.

4-1. どんなときに非線形最小二乗法を使うのか?

 前回,ヤコビアン行列の中に求めたいモデル(未知数)が含まれる場合は,②で説明した最小二乗法のような計算では最適なモデルを求めることができないため非線形最小二乗法を適用しないといけないということを述べました.非線形問題となる場合の例を(4-1)式に示します.

画像1

                           (応用地質(株)HP トモグラフィー解析の基礎 資料⑤より)
ここで,x1, x2, x3が求めたいモデル,Zは変数であり,下記のように表せます.

画像2

したがって,Yは(4-2)に示すように,方程式はm本出てくることになります.

画像3

ヤコビアン行列を求めると(4-3)式のようになります.

画像4

(4-3)式より,ヤコビアン行列の中に求めたいモデルXが入っていることがわかります.このような場合,非線形最小二乗法によりモデルを求めることになります.

4-2. 非線形最小二乗法によるインバージョンの計算

 上記に示した数理モデル(4-1)を例に,非線形最小二乗法によりモデルを求める方法を以下に説明していきます.
 初期モデルと理論値の残差を求めるまでの計算手順を示したものを図5に示します.また,最小二乗法によりモデルを修正し,新しいモデルで観測値との残差を求める計算手順を示したものを図6に示します.手順(ⅰ)~(ⅱ)は図5の赤線,(ⅲ)~(ⅴ)を図6の赤線に示します.

(ⅰ)初期モデルX0に対する理論値Y0を求めます.
(ⅱ)理論値Y0と観測値Yの残差を求めます.
 残差を(4-4)式に,二乗平均平方根(RMSE: Root Mean Square)を(4-5)式に示します.通常,残差の大小はRMSEをもとに判断します.

画像5

画像6

画像7

図5. 非線形最小二乗法によるインバージョンの計算手順(はじめの部分)

(ⅲ)最小二乗法によりモデル修正量を求めます.
 求めたいモデル修正量をX0とすると,残差二乗和Eは(2-6)式に倣って,
(4-6)式のように表せます.

画像8

(4-6)が最小になるとき,(4-7)式に示す正規方程式が得られます.

画像9

(4-7)をX0について解くと,(4-8)式のように表せます.

画像10

ここで,新しいモデルX1は,(4-9)式で求めることができます.

画像11

(ⅳ)修正したモデルから新しい理論値を求めます.
 X1から求まる新しい理論値Y1を求めます.

(ⅴ) 残差が十分小さくなるまで(ⅱ)~(ⅳ)を繰り返し計算(反復計算)します.
(ⅵ)残差が十分小さくなったら,計算を終了します.
 反復計算の回数は最初に決めておく必要があります.N回反復計算した場合,XNが最適なモデルとなります.

画像12

図6. 非線形最小二乗法によるインバージョン計算手順(反復計算部分)

4-3. 非線形最小二乗法の応用(より実践的な解析のために)

 非線形問題において最適なモデルを求めるためには,上述したような計算が基本になりますが,実務で使用されるソフトウェアでは計算をうまく収束させるためにいろいろな工夫がされています.例として,レーベンバーグ・マルカート法(Levenberg-Marquardt法),平滑化制約付非線形最小二乗法を紹介します.

1)レーベンバーグ・マルカート法(Levenberg-Marquardt法)

 レーベンバーグ・マルカート法(単にマルカート法)では,最小二乗法によりモデル修正量を求めるとき,(4-10)に示す汎関数Φが最小になるように計算します.計算結果が安定して収束するために,係数λを調整します.(4-10)は(4-6)の残差二乗和に収束条件(λの項)が加わった形となります.

画像13

(4-10)が最小になるように整理して得られた正規方程式(4-11)をΔXについて解くことでモデルを求めます.

画像14

2) 平滑化制約付非線形最小二乗法

 モデルパラメータ数が多い場合やデータに多くのノイズが含まれる場合に,詳細な地下構造を求めようとしてモデルの未知数を増やすと解を安定して求めることが難しいという問題が起きます.例として,地下構造モデルを沢山のブロックに区分させるとしましょう.隣接するブロック間では求めたい物性値が急激に変化しないという平滑化制約条件を設けます.ブロック間の求めたい物性値の差分をとる行列をCとしたとき,(4-12)に示す汎関数Uが最小になるように計算します.

画像15

(4-12)を最小になるように整理して得られた正規方程式(4-13)をΔXについて解くことでモデルを求めます.

画像16

 参考資料

・応用地質(株)HP トモグラフィー解析の基礎 資料⑤
https://www.oyo.co.jp/oyocms_hq/wp-content/uploads/2015/01/text_tomography.pdf

・産総研資料 比抵抗法電気探査の測定・解析と電極配置に関する研究
https://staff.aist.go.jp/takakura-s/electric/el-pri.pdf


※専門家の方へ,誤った記述がありましたらご指摘ください.

この記事が気に入ったらサポートをしてみませんか?