Arbitrary high order shaped stencils for time domain finite difference schemes in seismic wave propagation

技術要約：地震波伝播における時間領域有限差分スキームのための任意の高次形状スタ stencil（ステンシル）

問題提起
地震探査イメージングは、数値的な波の方程式の解に大きく依存しており、そこでは偽の数値波を抑制し、物理的な分散関係を正確に保持するために高忠実度の解が求められる。有限差分（FD）スキームは、その計算効率の高さから地震探しのイメージングにおける主要な選択肢であるが、固有の数値分散の問題を抱えている。この誤差は、位相速度が格子間隔や離散化の方向に依存することを引き起こし（異方性）、多くの場合、高周波成分の遅延として現れる。古典的な手法はテイラー展開に基づくものが主流であるが、これらはしばしば不十分な波分散の表現しか提供できない。本論文は、標準的なクロス型（十字型）の幾何学形状を超えて、最適化された任意の形状および次数のFDステンシルを用いることで、計算コストの削減と現実的な速度モデルのより優れた表現という実質的な利点が得られるかどうかの問いに取り組んでいる。

手法
著者らは、任意の形状のステンシルを用いて音響波方程式に対する高次時間領域FDスキップルを導出するための、一般的な数学的フレームワークを提案する。本研究は等方的な音響ケースに焦点を当て、Devitoという記号計算フレームワーク内で実装されている。これにより、高度に最適化されたステンシルベースのコード生成が可能となる。

手法は主に以下の3つの構成要素からなる。

1. ステンシル幾何学： 論文では、以下のような様々なステンシル形状を調査している。

   • クロス・ステンシル（Cross Stencils）： 軸に沿った古典的なライン・ステンシル。
   • 菱形およびクロス・菱形ステンシル（Rhombus and Cross-Rhombus Stencils）： 軸から外れた「内部」の点を追加したもの。
   • 正方形およびクロス・正方形ステンシル（Square and Cross-Square Stencils）： ステンシルを正方形に近い幾何学形状で満たしたもの。
     本研究では、格子点の数と重みの対称性を制御する幾何学的パラメータ（$M, N, P, Q$）を定義している。

2. 導出アプローチ： 著者らは、ステンシル重みを決定するために、いくつかの導出戦略を統合および拡張している。

   • SpatTE： 空間的テイラー展開（古典的な手法）。
   • SpecTE： 空間演算子のスペクトル・テイラー展開。
   • SpecLS： スペクトル・テイラー級数の最小二乗近似。
   • DispTE： 分散に基づくテイラー・トランケーション（打ち切り）。時間微分を空間微分と交換することで、分散誤差を低減させる。
   • DispLS： 全時間・空間スペクトル演算子の最小二乗近似を行い、分散誤差を最小化する。

3. 評価指標： 古典的な $L_1, L_2, L_\infty$ 誤差ノルムに加え、離散フーリエ変換（DFT）および連続ウェーブレット変換（CWT）解析を用いている。これらのツールは、周波数および時間・周波数領域における重み付き位相差と振幅偏差を測定することにより、分散効果を定量化する。

主な貢献

• 統一されたフレームワーク： テイラー、スペクトル、および最小二乗法のアプローチを統合し、任意の次数のFDスキームを任意の形状のステンシルで導出するための一般的な数学的定式化を提供する。
• 体系的な比較： さまざまなステンシル幾何学（Cross, Rhombus, Square）と、さまざまな導出手法（SpatTE, SpecTE, DispTE, SpecLS, DispLS）を組み合わせた包括的な数値比較を行う。
• Devitoへの実装： これらの手法をDevitoに実装し、高レベルの記号的定義がいかに大規模な音響地震インバージョン問題に適した効率的な低レベルコードへとコンパイルされるかを実証する。
• 非クロス型ステンシルの分析： 非クロス（内部）の点をステンシルに追加することが精度を向上させるという仮説を厳密に検証する。

結果
数値実験は、4つの速度モデル（均質、二層、2D SEG/EAGE Salt、およびMarmousi）に対して実施された。

• 均質および単純な不均質モデル： 一定または単純な二層速度モデルにおいて、最適化されたスキーム（特に DispTE-CrossRb および DispLS-CrossRb）は、高次（$M \geq 4$）において、古典的なSpatTEスキームよりも一般に分散抑制において優れている。
• 複雑な不均質モデル（Salt/Marmousi）： 大きな速度コントラストを持つモデルでは、最大速度に基づいてグリッド全体に単一のステンシルを適用する場合、速度依存型の最適化スキーム（DispTEやDispLSなど）の性能が低下する。
• 幾何学に関する知見：
  • クロス・ステンシル vs 複雑な形状： より多くの点が精度を高めるという直感に反して、結果はクロス・ステンシルが精度と計算コストの間の最適なトレードオフを提供することを一般的に示している。
  • 収穫逓減： 非クロス（内部）の点を追加すること（菱形や正方形の形状）は、必ずしも精度の向上や分散の低減をもたらさない。多くの場合、特に正方形ステンシルにおいては、重みの増加が線形システムの不良設定化、数値的不安定性、あるいは増加した計算コストに対する無視できるほどの改善を招く。
  • 空間のみの手法の堅牢性： 高度に不均質な媒体において、速度依存の重みに依存せずに最適化を行うスキーム（具体的には SpecLS-CrossRb および SpatTE-Cross）は、速度依存型のスキームよりも高い堅牢性と低い分散アーティファクトを示した。
• トレードオフ： 均質領域における高次スキーム（$M \geq 4$）では、最適化されたクロス・ステンシルが大幅な分散低減を実現する。しかし、複雑な媒体においては、最適化スキームの速度変化に対する敏感さが、より単純な速度独立型のスキームや適応的な戦略の方が好ましい可能性を示唆している。

意義および主張
本論文は、精度向上のために複雑なステンシル幾何学を探索する文献が多いものの、非クロス型のステンシル形状（菱形および正方形）は、増加した計算コストにもかかわらず、必ずしも精度の向上や分散の低減をもたらさないと主張している。

本研究の主な意義は以下の通りである：

1. クロス・ステンシルの検証： 最適化された分散低減型クロス・ステンシルが、古典的なアプローチやより複雑な形状と比較して、減少した計算コストで十分な精度を提供することを確認した。
2. 速度感受性の強調： 分散最適化スキームの利点が、速度モデルに強く依存することを実証した。複雑で不均質な媒体では、速度依存の重み計算に依存するスキームは、空間的に最適化されたスキームや古典的なスキームよりも性能が劣る可能性がある。
3. 実用的な実装： Devitoへの実装を提供することで、これらの手法をフルスケールの音響地震インバージョン（FWI）および逆時間マイグレーション（RTM）問題へ適用することを容易にする。
4. 将来の適応性： 本フレームワークは異方性モデル（VTI/TTI）にも適応可能であり、また（SpecLSのような）ステンシルの事前計算が反復的なインバージョン・ワークフローにおいて効率性の利点を提供することを述べている。

著者らは、複雑な媒体における実用的なアプリケーションにおいては、時間および空間のパラメータの選択、ならびに速度モデルが解の品質を決定的に左右すること、そして「Cross」の幾何学形状がテストされたシナリオにおいて最も効率的かつ堅牢な選択肢であることを結論付けている。