(57)【要約】【課題】本発明は、リモートセンシングや水中音響情報処理の分野に関し、水中音場を高速にシミュレートすることを課題とする。【解決手段】音波伝播モデルに応じて処理を分割し、並列実行する。すなわち、ノーマルモード法については、多数の固有値・固有ベクトルを並列計算する。ＦＤＴＤ法については、領域を分割して適宜境界データを交換しながら並列計算する。【効果】本発明によれば、音場計算を高速化でき、反復シミュレーションや実時間シミュレーションの実現が容易になる。【選択図】  図１
2005024354.tif 000002
【特許請求の範囲】
【請求項１】
ＦＤＴＤ法の音場を並列計算する音場計算装置であって、
初期化処理部と分割処理部と配布処理部と部分計算処理部と統合処理部とからなり、
初期化処理部は媒質特性の場と格子化条件とを入力とし差分スキームおよび初期音場の格子データを出力し、
分割処理部は領域分割条件と上記差分スキームと初期音場とを入力とし
部分領域用に分割された差分スキームと初期音場と格子化条件とを出力し、
配布処理部は音源条件と領域分割条件と分割された差分スキームと初期音場と格子化条件とを入力とし
各部分計算処理および回収統合処理を起動して音源および上記入力データを該当するプロセスに出力し、
部分計算処理部はこれらのデータをもとに
所定の時間ステップごとに音場の更新および隣接する領域との必要な境界データの交換をおこない、
また音場を必要に応じて時間・空間方向にリサンプリングのうえ出力し、
統合処理部は各部分計算処理部からの出力を回収しつなぎ合わせてひとつの音場データを作成し・出力する
ことを特徴とする音場計算装置。
【請求項２】
請求項１の音場計算装置において、
各部分計算処理部が境界データ交換用のバッファを備え、
複数の時間ステップごとにまとめて境界データの交換をおこなう
ことを特徴とする音場計算装置。
【請求項３】
請求項１の音場計算装置において、
音場計算の中断要求を待ち受け、
中断要求があった場合各部分計算処理部を時刻をそろえたうえで止め、
その時点での各部分計算処理部の音場を回収・統合・保存して終了し、
また次回起動時に初期化処理部が上記で保存しておいた音場を取得し、
これにより音場計算の随時中断・再開を可能とする
ことを特徴とする音場計算装置。
【請求項４】
ノーマルモード法の音場計算を並列的に実施する音場計算装置であって、
初期化処理部と複数の固有ベクトル算出部と回収統合部とからなり、
初期化処理部は音源、海洋、格子化の各条件を入力として
差分スキームを算出・出力し、
配布部は上記差分スキームを入力とし
そのコピーをそれぞれの固有ベクトル算出処理部に出力し、
固有ベクトル算出処理部は差分スキームを入力とし
これに対応する対称行列の所定の固有ベクトルを算出・出力し、
回収統合部は上記の固有ベクトルを入力とし
固有ベクトルを所定の方法で重ね合わせることにより音場を算出・出力する
ことを特徴とする音場計算装置。
【請求項５】
請求項１〜４の音場計算装置において、
分散オブジェクトとの通信インタフェースを備え、
音源位置や観測位置などの条件を入力として随時受け付け、
その都度音場計算を実施または過去の音場計算結果を検索して
入力データに対応する音圧や伝播時間などの音場データを出力する
ことを特徴とする音場計算装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明はリモートセンシングや水中音響情報処理の分野に関する。
【０００２】
【従来の技術】
各種の音響モデルを説明するため、まず水中音響で用いられる波動の方程式について説明する。座標系を図５のとおりとし、記号はｐを音圧、ρを密度、ｃを音速、ｔを時刻、ｋを波数とする。次の式は通常の波動方程式である。
【０００３】
【数１】
2005024354.tif 000003
水中音響の分野では、これに音源の定常性と等方性を仮定して極座標表示した次の方程式（Ｂｅｓｓｅｌ  の方程式）がよく用いられる。
【０００４】
【数２】
2005024354.tif 000004
次に各種の音響モデルの概要を説明する。これらについては「海洋音響学会誌  Ｖｏｌ．２８  Ｎｏ．４」（平成１３年１０月）により詳しい説明がある。
【０００５】
ノーマルモード法とは、Ｂｅｓｓｅｌ  の方程式において
【０００６】
【数３】
2005024354.tif 000005
のように深度と距離とを変数分離して
【０００７】
【数４】
2005024354.tif 000006
の形の固有値問題に帰着する方法である。数値的には、上記の固有値問題を差分化した行列の固有値・固有ベクトルを多数計算し、上記の変数分離式のように重ね合わせることで音場ｐ＝ｐ（ｒ，ｚ）を算出する。
【０００８】
ＰＥ（Ｐａｒａｂｏｌｉｃ  Ｅｑｕａｔｉｏｎ；  放物型方程式）法とは、Ｂｅｓｓｅｌ  の方程式を
【０００９】
【数５】
2005024354.tif 000007
のように  Ｈａｎｋｅｌ  関数Ｈ０と摂動項Ψとに変数分離して変形し（ｋ０  は典型的な波数をあらわす定数）、
【００１０】
【数６】
2005024354.tif 000008
の方程式に帰着するものである。Ｘはｚのみについての２階の微分作用素である。ＰＥ法の名前は、上記方程式がｒを時間に見立てると熱方程式（放物型方程式）に近い形であることに由来する。Ｈａｎｋｅｌ  関数とは、Ｂｅｓｓｅｌ  の方程式の解としてよく知られたものである。これについては日本数学会編「数学事典」岩波書店に説明がある。数値的には、上記方程式を差分化した３重対角行列係数の１次方程式を解き、ｒが小さいほうから大きいほうへ摂動Ψ＝Ψ（ｒ，ｚ）を順次算出する。
【００１１】
ＦＤＴＤ（Ｆｉｎｉｔｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｔｉｍｅ  Ｄｏｍａｉｎ；  有限差分時間領域）法は、波動方程式を導くもととなる媒質の状態方程式
【００１２】
【数７】
2005024354.tif 000009
と媒質の運動方程式
【００１３】
【数８】
2005024354.tif 000010
とを連立させて解くものである。前者は音圧が変位に比例する  Ｈｏｏｋ  の法則をあらわす状態方程式、後者は媒質の加速度が音圧勾配に比例し粘性抵抗の分だけ減るすることを現す運動方程式である。ｐ＝ｐ（ｒ，ｚ，ｔ）は音圧、ν＝ν（ｒ，ｚ，ｔ）は媒質の速度、Ｋとηは別途定まる係数である。媒質速度νが決まると状態方程式により次の時刻の音圧ｐが決まり、音圧ｐ決まると運動方程式により次の時刻の媒質速度νが決まる構造になっている。数値的には、これらを差分化した  ｌｅａｐ  ｆｌｏｇ（蛙飛び）アルゴリズムと呼ばれるｐ，νの交互計算により音場ｐを時々刻々算出する。ノーマルモード法およびＰＥ法では時間非依存の音場ｐ＝ｐ（ｒ，ｚ）を算出するのに対し、ＦＤＴＤ法では時間依存の音場ｐ＝ｐ（ｒ，ｚ，ｔ）を算出する。
【００１４】
従来、電波伝播を対象にしたＦＤＴＤ法を並列化した例が打矢・柏「並列型スーパーコンピュータを用いたＦＤＴＤ並列計算」ＩＥＩＣＥ  Ｔｒａｎｓ．  ｖｏｌ．Ｊ８４．Ｃ  Ｎｏ．１１  ２００１  に記載されている。しかし基本方程式の異なる音波伝播を対象にしたＦＤＴＤ法を並列化した例は見当たらない。またノーマルモード法やＰＥ法を並列化した例は見当たらない。
【００１５】
【非特許文献１】
「海洋音響学会誌  Ｖｏｌ．２８  Ｎｏ．４」（平成１３年１０月）
【非特許文献２】
日本数学会編「数学事典」岩波書店
【非特許文献３】
打矢・柏「並列型スーパーコンピュータを用いたＦＤＴＤ並列計算」ＩＥＩＣＥ  Ｔｒａｎｓ．  ｖｏｌ．Ｊ８４．Ｃ  Ｎｏ．１１  ２００１
【００１６】
【発明が解決しようとする課題】
本発明は水中音響シミュレーションを高速におこなうことを課題とする。
【００１７】
【課題を解決するための手段】
ノーマルモード法については、多数の固有値・固有ベクトルを並列計算する。ＦＤＴＤ法については、領域を分割して適宜境界データを交換しながら並列計算する。
【００１８】
【発明の実施の形態】
第１の例として、ＦＤＴＤ法を並列化する装置を説明する。図１がその構成およびデータの流れをあらわしている。音源データ１１１は音源の周波数や位置などからなる。海洋１１２は音速プロファイルや吸収率プロファイルからなる。格子化の条件１１３は、解析の空間的および時間的範囲とその刻み幅、出力の空間的および時間的範囲とその刻み幅（通常は後者のほうがかなり粗い）などからなる。初期化処理部１０１では、海洋１１２と格子１１３の条件から定まる差分スキーム１１５と、格子から決まる初期音場ｐ・初期媒質速度場νの格子データ１１６とを出力する。ｐ、νの初期値は通常至る所０である。分割処理部１０２は、解析および出力の領域を分割し、また分割した解析領域を計算機に割り振る。ここで差分スキーム１１５は分割された領域に対応する分割された差分スキーム１１７に、同じく初期音場・初期媒質速度場１１６は分割された初期音場・媒質速度場１１８に、同じく解析用・出力用の格子１１６は分割された格子１１９に変換され、また分割された領域で音場計算をするためのＦＤＴＤサーバ１０６を起動する各計算機やそれらサーバ間の接続方法を記述するスクリプト１２０が生成・実行される。配布処理部では、分割された差分スキーム１１７、初期データ１１８、格子１１９を対応するＦＤＴＤサーバ１０６に配布する。部分音場計算処理部１０４内では、各ＦＤＴＤサーバ１０６が時々刻々の音場を計算し、隣接する部分領域から境界のデータを交換しながら全体領域の音場をもれなく計算する。回収統合処理部１０５では、上記各ＦＤＴＤサーバ１０６の音場計算結果を定期的に回収し、それらをつなげて１個の音場データに統合する。算出された音場１２１は必要に応じて加工・保存される。
【００１９】
境界データの交換の仕方を、図２をもちいて説明する。ある時刻ｔ（ｋ）において、各部分計算領域で正しい音場ｐ・媒質速度場νを保持しているとする。次の境界データ交換時刻ｔ（ｋ＋１）＝ｔ（ｋ）＋Δｔまでの処理にうつる前に、各部分計算領域は境界の内側の正しいデータを隣接領域（辺を接するもの）に提供し、同じく境界の外側の正しいデータを隣接領域から取得・保持する。音速をｃとすると、次のデータ交換時刻まで境界の内側で正しいデータを保持するには、境界からｃΔｔの厚さに相当する部分のデータを交換すればよい。ここで、まず上下方向のデータ交換をおこない、それが完了してから左右方向のデータ交換をおこなうようにする。こうすることで、図２の濃い色の部分のような、直接辺を接する領域が正しい値を保持していない部分についても、対角にある領域から正しい値を取得できる。なお境界データ交換の対象の領域の厚さ（符号２３０）〜ｃΔｔを大きくとるとデータ交換の回数が少なくて済むので、通信のオーバーヘッドが大きい場合に好適である。
【００２０】
第２の例として、ノーマルモード法を並列化する方法を図４をもちいて説明する。初期化処理４０１で差分スキームを決定したあとは、従来の技術の節でも説明したように、ノーマルモード法の計算の中心は３重対角行列の固有値・固有ベクトルを（小さいほうから）多数算出する部分４０２である。この部分を並列に実行し、結果を回収統合処理部４０３で集めて重ね合わせ、最終的に音場４１５を得る。したがってこの固有値・固有ベクトル４０２を互いに独立に算出できれば十分であり、以下ｋ番目の固有値・固有ベクトルを算出する方法の例を説明する。
【００２１】
対称行列Ａのλ’以下の固有値の数をＳ（Ａ，λ’）という記号であらわす。３重対角行列の場合、Ｓ（Ａ，λ’）の値は  Ｓｔｕｒｍ  の定理を用いて容易に計算できる。これについては  名取「数値解析とその応用」  １９９０  コロナ社  に説明があるとおりである。明らかにＳ（Ａ，−｜Ａ｜−ε）＝０、Ｓ（Ａ，｜Ａ｜＋ε）＝ｎ（次元）であり、λａ  ＝  −｜Ａ｜−ε、λｂ  ＝  ｜Ａ｜＋εはｋ番目の固有値の下界および上界になっている。そこでつぎの二分法の処理を繰り返す。
【００２２】
【数９】
2005024354.tif 000011
λａとλｂ  とが他の固有値に比べｋ番目の固有値に十分近づいたら、つぎの逆反復法の処理を収束まで繰り返す。
【００２３】
【数１０】
2005024354.tif 000012
これにより、ｋ番目の固有値が縮退していない場合はｋ番目の固有値・固有ベクトルが得られる。なお固有値が縮退している場合は、すでに求まっている固有ベクトルで所定の初期値を  Ｓｃｈｍｉｄｔ  直交化したものを逆反復法の初期値として用いれば、並列化の効果は損なわれるものの、独立な固有値・固有ベクトルを順次得ることができる。
【００２４】
第３の例として、第１の例の並列化ＦＤＴＤを中断・再開する方法を、図３、図１を用いて説明する。まず、図１の回収統合処理部１０５が中断要求を常時監視３１０するようにする。シグナルやスレッドを用いれば、ＦＤＴＤの計算を妨げずに中断要求を監視することができる。中断要求を検出した場合、部分計算同期処理３２０にうつる。ここではまず回収統合処理部１０５が各部分計算処理部１０６に中断要求を転送する。ｉ番目の部分計算処理部１０６は、この中断要求を検出した場合、ＦＤＴＤの計算をきりのよいところまで進めてから一時停止し、処理済のシミュレーション時刻ｔｉを回収統合処理部１０５に報告する。回収統合処理部１０５はこれらのシミュレーション時刻の最大値ｔｍａｘを計算し、各部分計算処理部１０６に通知する。各部分計算処理部１０６は、シミュレーション時刻ｔｍａｘまでＦＤＴＤの処理を通常どおりに進める。ここで音場回収統合処理３３０にうつる。各部分計算処理部１０６は、時刻ｔｍａｘにおいて保持している部分計算領域の音場（スカラー場）ｐおよび媒質速度場（ベクトル場）νを回収統合処理部１０５に出力する。回収統合処理部１０５は各部分計算処理部１０６からの上記データを回収・統合し、時刻ｔｍａｘにおける全体の音場ｐおよび媒質速度場νを作成する。ここで音場出力処理３４０にうつる。ここでは、上記の音場ｐ、媒質速度場ν、時刻ｔｍａｘの各データを所定のファイルに保存する。これで中断処理は完了である。再開の場合は、保存した音場ｐ、媒質速度場ν、時刻ｔｍａｘの各データを初期化処理部１０１で読み込めばよい。
【００２５】
【発明の効果】
以上のように、本発明により各種音響モデルにもとづき音場を高速に計算できる。
【図面の簡単な説明】
【図１】並列ＦＤＴＤ法による音場計算装置の構成図。
【図２】並列ＦＤＴＤ法による部分計算処理部間の境界データ交換の説明図。
【図３】並列ＦＤＴＤ法による音場計算装置の中断処理の説明図。
【図４】並列ノーマルモード法による音場計算装置の構成図。
【図５】音響モデルを記述する座標系を示す図。
【符号の説明】
１００…並列ＦＤＴＤ法による音場計算装置の構成およびデータの流れ、１０１…初期化処理部、１０２…分割処理部、１０３…配布処理部、１０４…部分計算処理部の集合、１０５…回収統合処理部、１０６…部分計算処理を担当する個々のＦＤＴＤサーバ、１１１…音源を規定するデータ、１１２…海洋などの媒質の場を規定するデータ（音速プロファイルや深度プロファイルなど）、１１３…差分化用の格子（範囲と刻み幅など）、１１４…トポロジー（部分計算における領域分割方法、部分計算領域と計算機とのマッピングなど）、１１５…差分スキーム（差分化方程式の係数）、１１６…初期音場（および初期速度場）、１１７…部分計算領域ごとに分割された差分スキーム、１１８…部分計算領域ごとに分割された初期音場（および初期速度場）、１１９…部分計算領域ごとに分割された差分化用格子、１２０…部分計算処理および回収統合処理の起動データ、１２１…全体の音場、２００…境界データ交換の説明、２１０…部分計算領域、２２０…データ交換領域、２３０…データ交換領域の厚さ、３００…音場計算中断の処理、３１０…中断要求検出処理、３２０…部分計算領域同期処理、３３０…音場回収統合処理、３４０…音場出力処理、３５０…中断時の音場データ、４００…並列ノーマルモード法による音場計算装置の構成およびデータの流れ、４０１…初期化処理部、４０２…固有値・固有ベクトル算出処理部、４０３…回収統合処理部、４１１…音源を規定するデータ、４１２…海洋などの媒質の場を規定するデータ（音速プロファイルなど）、４１３…差分化用の格子（範囲と刻み幅など）、４１４…差分スキーム（差分化方程式の係数）、４１５…音場。
【図１】
2005024354.tif 000013
【図２】
2005024354.tif 000014
【図３】
2005024354.tif 000015
【図４】
2005024354.tif 000016
【図５】
2005024354.tif 000017