JPH11514471A

JPH11514471A - 超並列計算機による傾斜ムーブアウト解析方法

Info

Publication number: JPH11514471A
Application number: JP9514397A
Authority: JP
Inventors: トーマスエイディッケンス; ウィリアムエイジュニアシュナイダー; デニスイーウィーレム
Original assignee: エクソンプロダクションリサーチカンパニー
Priority date: 1995-10-06
Filing date: 1996-10-02
Publication date: 1999-12-07
Also published as: CA2223849A1; AU7253096A; EP0856179A4; WO1997013213A1; NO981476L; DE69625369D1; AU705347B2; NO981476D0; CA2223849C; DE856179T1; EP0856179B1; EP0856179A1

Abstract

(57)【要約】本発明は並列プロセッサ（１０６、１０８、１１０）、好ましくは超並列プロセッサ（ＭＰＰ）（１００）上で地震デバイス（１０２、１０４）を処理する方法に関する。プロセスへの入力は、複数の地震ライン（２）の１つまたはそれ以上からの地震データである。本発明の最終成果は、地下のＤＭＯ修正済ゼロオフセット地震イメージである。異なる入力オフセットに対して本方法を繰り返すことによって、入射角の関数としての反射体の反射強度の変化を調べることができる。本方法は、２段階並列処理を含んでいる。第１に、到来トレースを並列処理して各到来トレース上の各点毎の１組のＤＭＯ修正済部分イメージを生成する。第２に、出力位置を並列処理して部分イメージを累積して組合わせ、出力トレースを発生させる。

Description

【発明の詳細な説明】超並列計算機による傾斜ムーブアウト解析方法発明の分野本発明は物理探鉱の分野に関する。詳しく述べれば、本発明は並列プロセッサを使用してプリスタック地震データの部分的にマイグレートされたイメージを求める方法を含む。発明の背景地下炭化水素鉱床の探索は、典型的にはデータ取得、解析、及び解釈手順のシーケンスを含む。データ取得相は、エネルギ源の使用して信号を生成し、この信号を土中へ伝播させ、種々の地下地質構造から反射させることを含む。反射した信号は、地表上の、または地表付近の、またはその上に横たわる水のボディ中の多数の受信機によって記録される。地震トレースと呼ばれることが多い受信信号は、時間、受信機位置、及び源位置の関数として変化し、そして最も重要なことには、信号を反射させた構造の物理特性の関数として変化する音響エネルギの振幅からなる。データ解析者は、これらのトレースを（地球）物理モデルと共に使用して地下の地質構造のイメージ（画像）を開発する。解析相は、調査中の地質構造の性質、及びデータセット自体の特性に依存して変化する手順を含む。しかしながら、一般的に言えば、典型的な地震データ処理努力の目的は、記録されたデータから地質構造のイメージを発生させることである。このイメージは、信号が土中へ伝送され、地下層によって減衰され、そして地質構造から反射される技法の理論的及び経験的モデルを使用して開発される。データ処理シーケンスの最終成果の質は、これらの解析手順の正確さに大きく依存する。最終相は、解析結果の解釈である。詳しく述べれば、解釈者の仕事は、地下の炭化水素鉱床が存在する範囲にアクセスし、付加的に試掘することが適正であるかどうか、または炭化水素回収シナリオを最適化するのは何かを決定するのを援助することである。この査定におけるイメージの解釈は、種々の異なる努力を含む。例えば、解釈者は、イメージ化された結果を検討してその地域の地下の地質の知識を入手することが多い。これは断層、向斜、背斜のような主構造特色をマークすることを含むかも知れない。次に、層準の予備輪郭画定（コンタリング）を遂行することができる。種々の垂直断面を横切る層準を、解釈された断層の相関と共に連続的に追跡するその後のステップも遂行することができる。当分野においては公知のように、地震シーケンスのデータ解析ステップの結果の質及び正確さは、この翻訳相の結果の正確さ及び有用さに重大なインパクトを与える。原則として、地震イメージは、地震波伝播の三次元物理モデルを使用して開発することができ、それによってデータ内の全ての反射の正確な深さ及び方位角のスケーリングが容易になる。解釈は、含まれている地質構造の性質に絞ることができ、イメージの正確さに関心は向けられないので、反射が正確に特定されていればデータ解釈は大幅に簡易化される。不幸にも、三次元物理モデルは耐えられない程長い計算時間を必要とすることが多く、解析時間及び費用の両者の負担を軽減するためには、地震解析はデータ処理努力を可能な限り大幅に簡易化することを余儀なくされる。解析は、３Ｄ計算に挑戦する他に、ボリュームの処理に挑戦することにも当面している。例えば、典型的なデータ取得を実行するためには、数百乃至数十万の源位置を含む可能性があり、各源位置は数百の受信機位置を有している。各源・受信機対は所望の出力イメージに対して価値のある貢献をなし得るから、データ処理負荷（即ち、入力／出力データ転送要求）は、計算負担から独立した、それ自体における負担であり得る。図１は、物理イメージが望まれる土地のある領域２０の斜視図である。地表１８には、地震データを取得する複数のショットライン（爆破線：shot line）２が示されている。図１Ａに示すように、ショットライン２は一連の位置からなり、これらの位置に地震源３が配置され、これらの源から地震信号５が土中に送信される。各線に沿って配置された受信機４は、各源からの信号を種々の地下反射物６から反射された後に受信する。上述した地震データ負担を管理する第１の方法は、データを取得した領域を注意深く限定することを含む。詳しく説明すれば、何等かの使用可能な予備的な地質及び物理情報を使用することにより、地震データを取得する必要がある表面領域を最小にすることを容易にすることができる。このような最小化により、最終的に取得されるデータの量が直接的に減少する。更に、ショットライン間の間隔を同様に注意深く計画することによりデータのボリュームが減少し、それによって解析努力が最適化される。そして最後に、使用される源及び受信機の数、及び相隣る源及び受信機位置間の間隔を最適化することも、データ解析に利益をもたらす。これらの努力は、ペナルティ無しに達成することはできない。例えば、ショットライン間の、または源と受信機との間の間隔を比較的広くすると、計算される地震イメージの解像力（分解能）が減少して解釈を行うことがより困難になる。更に、比較的狭い間隔を用いることなく複雑な地質特色を分解することはできない。そして最後に、比較的未踏の領域におけるような若干のデータ取得においては、データを取得する表面領域を最適化することはできない。その結果、データ取得計画を通して、データ処理負担を完全に排除することはできない。近年使用可能になった超並列プロセッサ（Massively Parallel Procssor：ＭＰＰ）は、地震データ解析に重要な機会を提供している。超並列プロセッサは、計算を同時に遂行できる複数の中央処理ユニット（ＣＰＵ）を有している。これらのＣＰＵを効率的に使用することによって、従来は数週間または数カ月を要した複雑な解析を、数日、もしくは多分数時間まで短縮することができる。しかしながら、この重要な長所は、効率的な計算アルゴリズムがＭＰＰソフトウェア内にエンコードされている場合に限って実現可能である。従って、ＭＰＰに地震データ解析を遂行させることは、多重ＣＰＵを利用する適当な計算アルゴリズムに対する挑戦でもある。この挑戦は、計算アルゴリズムが、最も一般的に既存地震解析ルーチンのために書かれている手法を考えることによって容易に検討することができる。最近までは、コンピュータは逐次計算と呼ぶ演算モードに頼っていた。逐次計算は、任意の時間に単一の手順だけを遂行する、即ち、多分データまたはイメージの単一のサブセットに絞られた解析ルーチンの使用を含んでいる。これは１つだけのＣＰＵを有するコンピュータの直接的な結果である。この理由から、単一のＣＰＵ上で使用できる唯一の最適化手順は、その手順またはサブセットの処理効率を向上させることである。しかしながら、最終的には全ての計算をその単一のＣＰＵによって遂行させなければならないから、高い性能を得るためのオプションは本質的に制限される。一方、ＭＰＰの多重ＣＰＵ能力は、同時計算という明白な利点を提供する。この利点は、作業を細分して種々のＣＰＵに遂行させることによって計算問題を解くのに必要な合計時間を短縮できることである。但し、この細分によって各ＣＰＵが有用に作業を遂行している間、他のＣＰＵも作業を遂行することを条件とする。不幸にも、多重ＣＰＵハードウェアの欠点は、ソフトウェア開発に長い間使用されてきた逐次処理方法を、より適切な並列化された計算方法に置換しなければならないことである。簡単に述べれば、ＭＰＰは、多重ＣＰＵを効率的に使用する処理方法を開発することを必要とするのである。理想的には、これらの方法は、プロセッサの間に作業を比較的均等に分散させ、全てのプロセッサが、他のプロセッサからの中間結果を待つことなく常時必要な計算を遂行しているように編成すべきである。並列処理方法を限定する、及びこれらの一旦限定された並列処理方法を最適化することへの挑戦は、地震データ処理分野において特に過酷である。一般に地震データは大量の個々のトレースからなり、それぞれ他のトレースからある程度独立して記録されている。ある時点に解析を単一の計算に集中することが要求される論理的には充分な逐次計算方法は、これらの独立トレースの解析に良好に適合する。これは、たとえ計算上の隘路が存在していても然りである。例えば、解析シーケンスのある部分は他の部分よりも比較的多くの計算時間を要求したり、他の計算が進行する前に完了させなければならなかったり、または例えば走時のような他のトレースと類似の入力データに頼ったりする。逐次処理では同時計算は行われないから、これらの隘路が、必要とされる合計処理時間を除いて、単一ＣＰＵの計算効率を低下させることはない。この合計時間要求に関することを除いて、解析者にとって、このような計算隘路が存在することが問題になることはない。しかしながら、ＭＰＰ計算能力の全ての長所を利用するために、目的が全てのＣＰＵに同時処理を遂行させることであることから、これらの隘路を排除することによって地震解析相を最適化するための方法を開発しなければならない。ＭＰＰのこの長所は、計算時間が単一のＣＰＵコンピュータ内のイメージ領域サイズ上に賦課する制限を考えることによって明白になる。例えば、図１の立方体２０のサイズを拡張することによって、または例えばショットライン２に付加的な源３及び受信機４を追加することによってイメージのサイズを増加させると、合計計算が増加する。計算時間に対するこの直接的なインパクトは、特にたとえ小さいイメージ領域であっても、そしてたとえ最高速度の逐次処理コンピュータを使用しても、数週間の計算時間を必要とするような重い負担をイメージサイズを最適化するための地震解析に賦課する。これとは対照的に、256程度、またはそれより多くの個々のＣＰＵを有することができるＭＰＰ上での効率的な処理は、より大きい領域によって必要となる付加的な作業のほんの一部分、例えば1/ 256だけを各ＣＰＵが引受けるので、計算時間が長くなる分は最小で済む。イメージ領域及びイメージ生成に必要な作業負荷のスケーラビリティがＭＰＰの主便益であり、もしこのような作業負荷スケーラビリティが得られる並列化地震処理方法が開発されれば、そしてその場合に限って、便益を得ることができる。効率的な並列化地震処理方法を決定するための基本的な考え方は、上述した地震解析プロセスを再検討することによって明白になる。上述したように、地震解析の目的は、測定された地震データを物理モデルを使用して解析し、地下のイメージを開発することである。従って、３つの各主要処理成分、即ちデータ、モデル、及びイメージは、ＭＰＰ内の種々のプロセッサの間に分散させるべき計算作業の候補と考えることができる。プロセッサ間に作業を分散させる１つのオプションは、入力地震トレースデータの異なるグループを異なるプロセッサに割当てることである。例えば、トレースを源位置によってグループ化し、異なるプロセッサに異なるグループを割当てることができる。同様に、出力イメージを細分して異なるプロセッサに割当てることができる。最後に、出力イメージを生成するために使用される物理モデルを、種々のプロセッサに割当てることができるグルーピングに細分することも可能である。（一般的にこのモデルは、処理努力の対象である数学的モデルが要求する演算に組み入れられていると考えられる。例えば、地震解析における数学的モデルは、波動方程式に基づくことが多い。）例えば、データを周波数ドメインに変換し、個々の周波数を個々のプロセッサに割当てることができる。また、これらのアプローチの組合わせを開発することも可能である。例えば、あるプロセッサのグループにイメージ内の特定の層準位置に対する個々の責を負わせながら、モデル内の特定周波数及びイメージ内の全ての深さに対して集団責任を割当てることができる。地震データ解析者にとっての挑戦は、地震データ、モデル、及びイメージをＭＰＰ内の個々のプロセッサに割当てることがきる成分に細分し、それによって各プロセッサにおける計算を他のプロセッサとは無関係に遂行できるようにする方法を決定することである。地震データ解析を個々の成分に細分することを、一般に地震分解と称している。「単一命令・多重データ」ストリーム（ＳＩＭＤ）ＭＰＰにおいては、処理要素は、典型的に複数のデータストリームに対して同一の演算を遂行する。一つの例は、Thinking Machines Corporation製のＣＭ２である。これらの種類のマシンは、典型的に共用メモリを有していない。即ち、各プロセッサはそれ自体の個別メモリユニットを有しているが、そのメモリ内の情報に他のプロセッサが直接アクセスすることはできない。個々のプロセッサは、典型的に制限された計算能力及びメモリを有している。多数の処理要素が存在しているにも拘わらず、共用メモリが存在していないことから、処理要素間のデータ転送がマシンの能力を効率的に使用する上での主要隘路になっている。たとえ超立方体配列におけるような精緻な相互接続技術を使用したとしても、プロセッサ間のデータの転送が、プログラムの走行時間の主要要因である。他のコンピュータは、数十または数百のアレイ内に遙かにパワフルな要素を有している。Cray Research Corporation製のＴ３Ｄが、この種類のマシンの例である。ＣＭ２におけるよりも遙かにパワフルな個々の処理要素を有している他にも、Ｔ３Ｄはより少ない要素と、物理的に分散され論理的に共用されるメモリとを有している。この複数命令・複数データストリーム（ＭＩＭＤ）マシンは、データの異なる部分に対して異なる演算を同時に遂行する異なる要素を有している。処理要素の数が少ないということは、ＳＩＭＤマシン程多くデータを転送しなくともよいことを意味している。個々の要素の精緻さ及びコストが増加しているので、またそれらの数が少ないので、効率的に使用するためには処理要素の負荷を平衡させる必要がある。付加的な要因は、各処理要素を総合データボリュームのより大きいサブセットに順応させなければならないことであり、データの格納を含む計算はより複雑になる。ムーアヘッド(Moorhead)の米国特許第5,404,296号に開示されているＭＰＰ上でのマイグレーション方法では、予め選択された規則的なパターンに配列された多数の処理要素が存在している。データは、データセット内の受信機位置の数とブロック内のショットの数との積が、使用可能な処理要素の合計数に等しくなるように始めにショット分類され、ショットのブロックに分割される。３Ｄ地震データボリュームは典型的には数百のショットと受信機位置とを含むので、上記開示はＳＩＭＤマシンに制限される。当分野においては公知のように、プリスタックマイグレーションプロセスは極めて大量の計算を必要とする。プリスタックマイグレーションの代替として、プリスタック部分マイグレーションを行い、次いでポストスタック時間または深さマイグレーションを遂行することができる。多くの場合、これによって精度を大きく落とすことなく計算時間が大幅に短縮される。プリスタック部分マイグレーションの一つの一般的方法は、傾斜(dip)ムーブメント（ＤＭＯ）プロセスである。ＤＭＯ修正のプロセスに関しては、デレゴフスキー及びロッカ(Deregows ki ，S.M．and Rocca，F.)の“Geometrical Optics and Wave Theory of Constant Offset Sections in Layered Media”，29 Geoph．Prosp．374（1981）に詳述されている。ＤＭＯオペレータは、２段階の結果として視覚化することができる。第１段階は、所与の源及び受信機位置に関する等時線（アイソクロン）を見出すことである。アイソクロンとは、源から各点まで、及びそれから受信機へ戻る走時が、観測された走時と等しいような点の集合を意味する。当分野に精通していれば理解されるように、定速度媒体の場合には、アイソクロンはショット及び受信機位置に焦点を有する楕円面である。この楕円面上の何れかの点の法線は地表上の位置と、法線によって限定される射線(ray)に対応するゼロオフセット走時とを限定する。ＤＭＯ演算は、データを始めの観測点から考え得る１組の表面位置まで移動させ、対応するゼロオフセット走時を決定することからなる。定速度媒体の場合には、表面位置は走時楕円面の長軸に沿う。キルヒホッフＤＭＯ手順は、１）各入力トレース毎に、複数の出力表面位置における１組の置換トレースを生成すること、それに続いて２）各出力位置において、各入力位置からの置換トレースの合計を求めること、そして３）部分的にマイグレートされたイメージを得ることからなる。デレゴフスキー及びロッカは、共通オフセット区分のための速度独立ＤＭＯオペレータを、ｔ₀＝ｔ_n｛１−（ｘ²/ｈ²）｝^1/2 （１）として表すことができることを示した。但し、ｔ₀はゼロオフセット時間、ｔ_nは正規ムーブアウト（ＮＭＯ）修正済時間、ｈは源・受信機オフセットの半分、そしてｘは中心点からのマイグレーション距離、即ち源・受信機対の中心点と等価ゼロオフセット置換トレースの位置との間の距離である。ゼロオフセット時間は、置換トレースに沿って測定された時間である。ＮＭＯ修正済時間は、平均自乗根（ｒｍｓ）速度を使用して通常のＮＭＯ手順によって求められる。当分野においては公知のように、傾斜修正済ｒｍｓ速度を使用する共通中心点（ＣＤＰ）の集まりを普通にスタッキングしたのでは、異なる傾斜の反射体を扱うことはできない。一方、ＤＭＯによって求められたゼロオフセット区分は、異なる傾斜からの反射を正確に表す。式（１）は、計算時間を大幅に短縮しながら、源及び受信機を通る垂直面内にエネルギを移動させ、ポストスタック・ゼロオフセットマイグレーション法によって処理できるイメージを発生させる２ＤＤＭＯオペレータを定義している。媒体が定速度を有していない場合のＤＭＯオペレータはより複雑である。定速度勾配を有する媒体の問題のための部分解については、アートリー及びヘール(A rtley，C．and Hale，D．)の“Dip-Moveoust Processing for Depth Variable V elocity”59 Geophysics 610(April 1994)に記載されている。しかしながら、アートリー及びヘールはＤＭＯオペレータのインライン成分だけを考えている。実際、垂直速度勾配を有する媒体の場合には、一般に「ゼロオフセットへのマイグレーション(Migration to Zero Offset ：ＭＺＯ)」オペレータと呼ばれるＤＭＯオペレータは、ショット及び受信機位置を通る垂直面内に限定されることはない。定速度勾配の場合のＭＺＯの問題の詳細な解析がディートリヒ及びコーエン(Dietrich，M．and Cohen，J.K．)の“Migration to Zero Offset( DMO)for a Constant Vertical Velocity Gradient：An Analytical Form ulatio n”41 Geophysical Prospecting 621，(1993)に記載されている。この場合のアイソクロンはｘ、ｙ、及びｚに関する４次多項式であり、ゼロオフセットへのマイグレーションオペレータは多値である。オペレータの１つの枝は、小さめの傾斜に関連するサドル状のオペレータであり、速度勾配がない場合の普通のＤＭＯオペレータにほぼ対応する。オペレータの第２の枝は、大きめの傾斜または転倒反射体から発生される。ＤＭＯオペレータのこれらの枝は、共にクロスライン成分を有している。ディートリヒ及びコーエンの作業は幾つかの限界を有している。第１に、これらのオペレータは射線理論的な考察に基づいている。その結果、オペレータの２つの枝が一緒になる場所で、オペレータの一次枝がかなり強い振幅で突然に終端する。更に、尖点終端から離れた場所においてさえ、ディートリヒ及びコーエンによって与えられる振幅は、ＤＭＯオペレータの振幅の問題に関連する従来技術の参照によって与えられる振幅と一致しなくなる。ゴンザレス(Gonzalez)らの米国特許第5,285,422号には、地下に垂直速度変化が存在する場合の３ＤＤＭＯを遂行する方法が開示されている。この方法では、定速度のＤＭＯオペレータが適用され、次いで第１のステップからの楕円面を圧搾（スクィーズ）する効果を有する残留オペレータが適用される。しかしながら、この方法は多値オペレータの一次枝に対処するだけである。シャイマン(Sch eiman)の米国特許第4,974,212号は、出力位置をある地面格子上の離散した組のビンに押しやる演算によってもたらされるスミアリングの問題を扱っている。しかしながら、この方法は垂直異質性の問題には触れていない。これらの限界に加えて、ゴンザレス及びシャイマンの両方法は共に単一ＣＰＵ向きの方法であり、ＭＰＰ上で効率的に実現することはできない。共用メモリ並列コンピュータ上での３ＤＤＭＯプロセスの実現が、ルー、ベル、及びリム(Lu，L.，Bell L.，and Lim，K．)の“Parallel Implementation o f 3-D DMO on Shared Memory Systems”64 SEG Extended Abstracts 214，(1994 )に記載されている。ルーらは、入力ショットの集まりに関してだけタスクを並列処理しているが、ＤＭＯ計算の出力相においては並列処理していない。しかしながら、この開示は単一オフセットデータのＤＭＯに制限され、加えて総合プロセスの効率を低下させる幾つかの限界を有している。詳述すれば、多重プロセッサはそれらの出力をＭＰＰの同一メモリ位置に書き込もうとする（または追加しようとする）ので、「メモリロッキング」手順がプロセス内に組み込まれている。プロセッサは所望メモリ位置へのアクセスを遅延させられることが多いために、このメモリロッキング手順によりプロセスの低速化がもたらされる。更に、ショットの各グループの処理の後に、類似の「ファイルロッキング」を使用してディスクファイルが更新される。これも、この方法における隘路である。ＣＭ２（またはThinking Machine ＣＭ５コンピュータ）上で３ＤＤＭＯプロセスをＭＰＰで実現することが、ビオンディ及びムーアヘッド(Biondi，B．an d Moorhead，W.D．)の“Data Parallel Algorithm for Kirchhoff 3-D DMO”，6 2 SEG Extended Abstracts 330，( 1992）に記載されている。この例は入力ショットの集まり、及び出力スウォス(swaths)を並列計算する。連続ショットの集まりは、全てのプロセッサが一杯になるまで連続するプロセッサ上にロードされ、各プロセッサはその入力の一部として１つまたはそれ以上の集まりを有している。各プロセッサは、複数の出力位置に対しても責を負っている。計算相中に各プロセッサは先ず、その入力トレースのための、及び出力タスクをも有している位置だけのためのＤＭＯ置換トレースを計算する。全てのプロセッサが初期入力トレースのための計算を完了すると、各プロセッサはその入力トレースを次のプロセッサへ送り、これらのＤＭＯ演算が繰り返される。この「パイプライン化」は、全ての出力位置のための出力トレースが計算されてしまうまで、必要に応じて複数回繰り返される。ＣＭ２またはＣＭ５の個々のプロセッサ上の使用可能なメモリが制限されているので、パイプライン化が必要である。またプロセッサは、最も近いオフセットから開始して入力トレースを増分的にドロップしてデータ転送を減少させる。この方法は、ショットの集まりが隣接し、且つ物理的に順番である、即ちプロセッサ間の負荷平衡が貧弱であることを必要とする。更に各プロセッサは出力の隣接部分に対して責を負っているから、複雑な速度モデルではより多くの計算が必要となるので、モデルの若干の部分において隘路が発生し得る。最後に、一般的な場合におけるように、若干のトレースが入力データセットから失われる場合に、このような計画における負荷平衡が劣化するようになる。ゴンザレス、ルーら、及びビオンディ及びムーアヘッドは、ＤＭＯ演算のカイネティックスのみを扱っている。即ち、ＤＭＯ計算は、置換トレースを生成するための直接的な加算を含んでいる。これらは、地震反射における振幅及び位相情報を適切に処理しない。ＤＭＯプロセスの一部として振幅及び位相情報を保存する問題に関しては、ブラック、シュライヒャー、及びツァング(Black，J.L.，Sc hleicher，K.L．and Zhang，L．)の“True-Amplitude Imaging and Dip Moveout ”58 Geophysics 47(January 1993)に記載されている。ブラックらは、適切な時点に種々の入力トレースからの貢献を単に加算するだけでは誤解を与えるような地下のイメージを発生する恐れがあると指摘している。詳述すれば、鋭く傾斜した反射体の振幅は、水平の反射体に対して減衰する。当分野に精通していれば理解されるように、振幅及び位相情報は地下の岩質及び流体含有の解釈に重要な役割を果している。ブラックらは、入力トレース間の距離が大き過ぎることに起因して、鋭い傾斜が不適切にイメージされるようなエイリアシングの問題をも検討している。しかしながら、ブラックらは、この問題に対する解決法を示唆していない。従って、ＤＭＯ計算の負荷を効率的に平衡させ、入力及び出力を並列処理し、複数のオフセットを処理することができ、垂直速度変化を処理することができ、ＭＺＯオペレータの両枝を使用して正確なイメージを発生し、そしてデータが隣接していることを要求せずに、ＭＰＰ上でＭＺＯプロセスを実現するような発明に対する要望が存在している。更に、この発明は地震データ内の振幅及び位相情報を従来技術よりも正確に保存すべきであり、またエイリアシング誤差を回避すべきである。本発明はこの要望を満足する。発明の概要本発明は並列プロセッサ、好ましくは超並列プロセッサ（ＭＰＰ）上で地震データを処理する方法である。プロセスへの入力は、複数の地震ラインの１つまたはそれ以上からの地震データである。本発明の最終成果は、地下のＤＭＯ修正されたゼロオフセット地震イメージである。異なる入力オフセットについて本方法を繰り返すことによって、入射角の関数としての反射体の反射強度の変化を調べることができる。地震データの振幅対オフセット解析の長所は、当分野に精通していれば理解されよう。本方法は２段階平行処理を含む。第１に、到来トレースを平行処理することによって、各到来トレース上の各点毎に１組のＤＭＯ修正された部分的イメージが生成される。第２に、出力位置を平行処理することによって、これらの部分的イメージが累積されて組合わされ、出力トレースが発生する。若干の実施例において、コンピュータへのデータ入力をどのような順序ででも格納できるようにしたことが本発明の特色であるが、この順序はシーケンス内の共通オフセット分類であることが好ましい。若干の実施例において、もし解析を遂行するプロセッサ上のメモリ割当てが一杯になれば、ディスク上にデータを書き込む他のプロセッサが存在して「解析」プロセッサ上のメモリ位置を空にするようにしたことが本発明の特色である。若干の実施例において、モデルを通してインタリーブされた出力位置に対する責を全ての解析プロセッサに負わせることによって、メモリ使用及びプロセッサ間通信を平衡させるようにしたことが本発明の特色である。このインタリービングは、個々の解析処理要素（ＰＥ）のメモリ要求を平衡させるので、特定のＰＥがそのメモリの限界に到達したことによりデータがディスク上に書き込まれる時の遊休時間を少なくする。若干の実施例において、ある所与の入力トレースのための全ての処理を、単一のプロセッサによって遂行させるようにしたことが本発明の特色である。これにより、入力トレースを１つのプロセッサから別のプロセッサへ引き渡す必要が排除され、別の通信負担も排除される。若干の実施例において、イメージングプロセスにおいて従来の方法よりも正確に振幅及び位相情報が保存されるようにしたことが本発明の特色である。イメージングプロセスは、垂直速度勾配が地下に存在する場合に使用することができ、鋭い傾斜を認識する能力を失うことなくデータ内の最大可能な周波数を使用する反エイリアシング特色を有している。本発明は、ＭＺＯ内の尖点の存在に起因する出力イメージ内の不連続を排除するが、これが従来技術に対する改善である。また、垂直速度勾配が存在する場合ＭＺＯの両枝をイメージングプロセスに使用することも本発明の特色である。これは、ＭＺＯオペレータの二次枝を発生させる鋭い傾斜及び転倒(overturned)反射体をイメージするのに本発明を有用ならしめる。図面の簡単な説明本発明及びその長所は、以下の添付図面を参照しての詳細な説明からより一層容易に理解されよう。図１及び１Ａは、３Ｄ地震探鉱に典型的に使用されるデータ取得構成を示す概要図である。図２は、フェザーリングの現象及びデータ内へのクロスライン成分の導入を示す図である。図３は、ＤＭＯによって遂行される部分的マイグレーション原理を示す図である。図４は、本方法を使用して地震データ処理を遂行することができる並列処理システムの解析相の簡易機能ブロック線図である。図５は、本発明の第１の実施例においてＤＭＯ修正に先立って各トレースに遂行される予備位相シフティング及びバンドパスフィルタリングのブロック線図である。図６は、本発明の第１の実施例においてイメージングのために１組のＤＭＯ修正済置換トレースを発生する際に各プロセッサによって遂行されるステップを示す図である。図７は、出力相のための並列処理システムの簡易機能ブロック線図である。図８は、個々の処理要素への出力位置のインタリーブされた割当てを示す図である。図９は、本発明の第２の実施例においてＤＭＯ修正に先立って各トレースに遂行される予備位相シフティング及びバンドパスフィルタリングのブロック線図である。図１０は、第２の実施例においてイメージングのためにＤＭＯ修正済トレースの予備セットを発生する際に各プロセッサによって遂行されるステップを示す図である。発明の詳細な説明本発明は、地震データを並列プロセッサ上で処理する方法である。本方法は、解析、及び／または、入力／出力タスクを遂行するために使用可能な１つより多いプロセッサを有するどのような型の並列処理コンピュータ上で使用することもできる。例えば、８つまでのプロセッサを有するCray Ｙ−ＭＰを使用することができる。ネットワークされたパーソナルコンピュータまたはワークステーションも使用することができる。好ましい実施例では、本方法は16またはそれ以上の処理要素を有する超並列プロセッサ（ＭＰＰ）上で使用される。このようなＭＰＰの一例は、Cray Research Corporation製のＴ３Ｄである。限定する意図はないが、便宜上以下にＭＰＰという略語を使用するが、これは本方法に適する何等かの並列処理コンピュータを意味するものとする。しかしながら、このようなＭＰＰの使用方法が、単に超並列であると商業市場で考えられているコンピュータへの本発明の応用を如何ようにも限定するものではない。公知のように、地震データ処理方法は、一般的に地表のある部分からのデータから開始される。図１は、物理イメージが欲しい土地の、ある領域の斜視図である。地表１８上には、地震データを取得するための複数のショットライン２が示されている。図１Ａに示してあるように、ショットライン２は一連の位置からなり、これらの位置に地震源３が配置され、地震源３から地震信号５が地中へ送信される。各ラインに沿って配置された受信機４は、各源位置から伝播され、種々の地下反射体６から反射した後の信号を受信する。プロセスへの入力は、地下２０の３Ｄデータボリュームのイメージを求めるために一連の地震ライン２から取得した地震データである。本発明の最終成果は、ＤＭＯ修正された、ゼロオフセット地震イメージである。更に、３Ｄボリュームの一部分に関する地震データは、「ケーブルフェザーリング」現象の結果として、単一のラインによる海洋データ取得からもたらされ得る。これを図２に示す。図２は、ケーブル４２を曳航している地震収集船４０の平面図であり、受信機４及び船４０の直後の源３がケーブル内にあることを示している。横流Ｖ_cが存在する場合、ケーブルの向きが船の運動方向Ｖ_vと一直線に整列しないことは明らかである。ＤＭＯ演算はケーブル上の源及び各受信機４からのデータを用いて、置換トレースを生成するＤＭＯ演算を遂行する。図２には、ＤＭＯの後にトレースが割当てられる位置を識別するビンの矩形格子（ｘ，ｙ）も示されている。ビンは地震データ収集の公称方向（即ち、ｘ軸）に沿って等間隔である。ＤＭＯ演算は、Ｖ_cが０でないことを条件として、単一の源位置から、源の各公称インラインｘ位置毎に異なるクロスラインｙ位置におけるデータを発生することが理解されよう。他の源位置３’についてもこれを繰り返すことによって、各ｘ位置における多くの異なるｙ位置のためのデータが得られる。データ収集中のケーブルの向きを測定すれば、このようなケーブルがフェザーリングした時のデータの３Ｄイメージを本方法によって求めることができる。図２には示してないが、ビンの番号付けは格子の左下隅から開始され、先ずｘ方向に増加して行き、次いで次のｙの値へ進んでｘ方向に値が増加する、等々になっている。これは、以下に図４を参照して説明する中間Ｉ／Ｏ手順を容易にするための便宜上そのようにしたまでである。ＤＭＯの原理を図３に示す。図３には、一定のオフセット地震区分上に単一の傾斜反射イベント６’が示されている。また、反射体６’上の反射点５０、及び各反射点を通るＤＭＯオペレータ４５も示されている。反射体６’のＤＭＯイメージ６”は、ＤＭＯオペレータ４５への接線として、即ち反射体点５０の等価ゼロオフセットイメージを表す点５０’を通る線として求められる。キルヒホッフＤＭＯは、個々の反射体６’上の孤立した点の代わりに、地震区分が多くの空間点及び時間点（ｘ，ｔ）にデータを有しているような共通（地球）物理データセットを扱う。これらの状況においては、イメージは各入力点を通るＤＭＯオペレータを「描く」ことによって求められる。当分野に精通していれば理解されるように、オペレータは、オペレータが傾斜に正接する点から離れているフレネルゾーンより小さいような領域の各反射体のある傾斜に沿って構造的に付加される。このようにして、各出力位置における結果的な貢献が合計されてＤＭＯイメージが得られる。本発明の特色の１つは、ＤＭＯタスクの並列処理が、出力段と比べて入力段（即ち、解析段）に関して異なることである。本方法の好ましい実施例の入力段に使用されるＭＰＰ並列処理システムは、図４を参照すると理解し易い。ＤＭＯ修正が適用されるプリスタックデータはテープ１０２上にあっても、またはディスクまたは他の入力データ記憶デバイス１０４上にあっても差し支えない。データは何等かの特定の順序にある必要はないが、もしデータが共通オフセットパネル内に分類され、パネル内でショット順にされていれば性能が改善される。解析者は、源・受信機中心点位置の範囲、及びＤＭＯが遂行される源・受信機オフセットを指定しよう。データセットは、これらがテープ１０２またはデバイス１０４上のデータトレースだけであるように前処理されていることが好ましい。当分野に精通していれば、この分類が好ましいことが理解されよう。ＭＰＰ１００は複数の処理要素（ＰＥ）からなり、その１つ１０６が制御ＰＥに指定される。好ましい実施例では、制御ＰＥが始めに入力データ記憶デバイスから地震トレースを読み出し、これらのトレースをＰＥのグループ１０８内のＰＥへ転送するタスクを遂行する。データを、ＰＥのグループ１０８内のデータ解析に使用可能な次のＰＥへ伝送することが、制御ＰＥのタスクの１つである。便宜上、グループ１０８内のＰＥを解析ＰＥと称するが、これらのＰＥの役割は本発明の後段では変化する。制御ＰＥは、そのメモリが許す限り多くのトレースを入力デバイスから予め読み出して、これらのトレースを非同期的に送り出す。勿論、始めはグループ１０８内の全てのＰＥは空であり、制御ＰＥはデータトレースをそれらに順番に送る。グループ１０８内の各ＰＥのタスクの１つは、ＤＭＯ演算を遂行することである。この演算の効果は、各入力トレースが、複数の出力位置に置換トレースを作り出すことである。各置換トレースは、各出力位置に出力イメージを作成する責を負っているＰＥに伝送される。プロセッサを出力位置に割当てることに関しては図７及び８を参照して後述する。プロセッサグループ１１０は、入力段中に入力／出力（Ｉ／Ｏ）機能を遂行する（後述）。図２を参照して上述したように、解析者は入力としてデータの３Ｄボリュームの出力（ｘ，ｙ）ビン位置を指定しよう（ここに、ｘはインライン位置であり、ｙはクロスライン位置である）。データが「ケーブルがフェザーリング」している海洋探査から得られたものであっても、標準のデータ格子から得られたものであっても、この指定は必要である。この指定に基づいて、グループ１０８内の各ＰＥは（ｘ，ｙ，ｔ）空間における複数のｘスライスに対して責を負う（但し、ｔはＰＥが出力イメージを発生する時間である）。一般的には本発明はある時点に地震データの単一のラインを処理し、従って（ｘ，ｙ）ビンは通常図１のライン２上に中心を有するスウォス内にあることに注目されたい。しかしながら、本発明が単一ラインのデータ処理だけに限定されるものではないことは、当分野に精通していれば理解されよう。置換トレースの伝送は相互の動作であり、グループ１０８内の各ＰＥは、新たに作成された置換トレースに関連するある出力位置に対して責を負うグループ１０８内の他の各ＰＥへデータを伝送するようになっている。ＰＥが、それ自身責を負っている出力位置のための置換トレースを作成するような場合も存在する。図４には、このプロセッサ間通信を、ＰＥのグループ１０８内の個々の処理要素間の破線相互接続によって示してある。ＤＭＯ修正及び置換トレースの作成のプロセスを、図５及び６を参照して説明する。この手順は定速度ＤＭＯ解析に関係しており、更に、層にした地球モデルのための近似解にも関係している。ＰＥのグループ１０８内の全てのＰＥは、制御ＰＥ１０６からグループへ送られる入力トレースのために、これらのタスクを並列に遂行する。しかしながら、従来技術に対してなされた本発明の改良は、入力トレースがグループ１０８内のＰＥのメモリ内へ読み込まれてしまうと、その特定の入力トレースの全ての処理がそのＰＥによってなされることである。当分野に精通していれば理解されるように、図５のステップ１２０における入力トレースは、デコンボリューション、スプレッディング補償、及びＮＭＯのようなルーチンステップを適用するために前処理されていることが好ましい。ＮＭＯ修正は、中心点（ｘ，ｙ）位置においてｒｍｓ速度関数を用いて遂行されているべきである。スプレッディング修正は、ＤＭＯ修正されたトレースを真の振幅で発生するために、本発明の能力を利用して適用さるべきである。もしトレースが前処理されていなければ、これは当分野においては公知の方法を使用して、ＰＥによって遂行される。ステップ１２１において、各ＰＥはそのトレースに関して源・受信機オフセット（ｈを源または受信機の中心点オフセットとして、典型的には代数表現２ｈによって定義される）が最小値ｄ_minよりも大きいか否かを調べる。もし大きければ、ステップ１２２においてフーリエ変換を適用した後に、ステップ１２４において（−ｉω）^1/2を乗算する（但し、ωは時相周波数）。この乗算は時間ドメインにおける半導関数に等価である。当分野においては公知のように、その結果は、45°の位相シフトと、ラインディフラクタ(line diffractor)からのデータの散乱を適切に斟酌するのに必要な高周波数のブーストとである。次いでトレースは、ステップ１２６において逆ＦＦＴによって時間ドメインへ戻すように変換される。当分野においては公知のように、源・受信機オフセットが小さい場合には、ＤＭＯオペレータは入力トレースの処理を最小に留めるべきである。もしその入力トレースの源・受信機分離が０であれば、ＤＭＯオペレータは入力トレースを同一位置へそのままコピーすべきであり、他の位置には何等の置換トレースをも発生すべきではない。何故ならば、一致源・受信機ジオメトリのためのアイソクロンは球形であり、アイソクロンへの対応法線がアイソクロン上の各点を源（及び受信機）位置へ戻して写像するからである。ステップ１２８及び１３０は、これらの小さいオフセット状況に対処することを意味している。もし２ｈが、ある最小量ｄ_minよりも小さければ、（−ｉω）^1/2演算に含まれる45°位相シフトは不適切である。ｄ_minの大きさは経験に基づいて選択するが、1800ｆｔの値が典型的な３Ｄデータセットにとって良好な結果を与えることを見出した。ｄ_minを指定する際に含まれる考察は、当分野に精通していれば理解されよう。もし源・受信機分離が０とｄ_minとの間に存在していれば、入力トレースに周波数の重み付けと、時間変化する位相シフトとが適用される（ステップ１２８及び１３０）。以下の計算が必要である。ａ）フーリエ変換、ｂ）√ωの乗算、ｃ）負の周波数における振幅のゼロイング、及び０及びナイキスト周波数における値を除く正の周波数における振幅の２倍化、ｄ）実数部分及び虚数部分の両者の√ωスケーリングを用いて解析信号を求めるために、ｃ）の出力の逆フーリエ変換、これはそのヒルベルト（または直交(q uadrature)変換）である。当分野に精通していれば理解されるように、ステップａ）及びｂ）は単なる√ω フィルタであり、一方ステップｃ）及びｄ）は濾波されたトレースの解析信号を与える。次に、ステップ１３０において、時間変化位相シフトが点毎に解析信号に適用される。当分野においては公知のように、瞬時位相は、直交（ヒルベルト変換された）成分と実数成分との比の逆正接（アークタンジェント）によって与えられる。この瞬時位相は、時間０からカットオフ時間ｔ_cまでは−π／４だけシフトされる。ｔ_cより大きい時間の場合には、−（πｐ）／２の位相シフトが適用される。但しｐ＝0.5min（２ｈ／ｄ_min，1.0）である。２ｈ／ｄ_mimは０と１との間にあるからｐは０と0.5との間にあり、従って位相シフトは０と−π／４との間で連続的に変化する。実際には、時間ｔ_cに中心を持つ線形遷移ゾーンが存在する。経験的に、この遷移ゾーンの幅は400ｍｓにセットされる。カットオフ時間ｔ_cの値は、ウィンボウ及びトランサム(Winbow，G.A．and Tra ntham，E.C．)の1994年の“Nonaliased amplitude-preserving DMO”，64th Ame rican International Meeting，Society Exploration Geophysicists，Expanded Abstracts，719-721niよって定義されている静止位相解析妥当性条件に基づいて求められる。この条件は、 ωｔ_n／（ｈ／Δｘ）²≪２（２）と書くことができる。ここに、ｔ_nは正規ムーブアウト時間、ｈは源・受信機分離の半分、そしてΔｘは入力トレース間隔である。経験的に、カットオフ時間ｔ_c＝0.8（ｈ／Δｘ）²／２πｆまでは静止位相条件が満足される。本発明では、地震データについて、より大きい時間の場合に入力データが源と受信機との間の中心点から導出されたものとして、または運動ＤＭＯ条件を必要としないものと等価で処理されるように、典型的な25Ｈｚのピーク周波数ｆを使用して時間ｔ_cを決定する。図５に示す手順は、ｄ_minより大きいかまたは等しいオフセットにおける位相シフトを常に−π／４に留める効果を与える。ゼロオフセットの場合、位相シフトは常に０である。オフセットが０乃至ｄ_minである場合には、位相シフトは小さい時間に対しては−π／４であり、時間が大きくなるにつれて０に接近する。時間変化位相シフトをこの特定例に限定することを意図するものではない。当分野に精通していれば、他の方法も明白であろう。オペレータの勾配（傾斜）が充分に大きく、隣接する置換トレース上の地震ウェーブレットが半波長より大きい相対シフトを受ける程になるとＤＭＯオペレータエイリアシングが発生する。その結果、高周波数入力トレースは低周波数トレースにおけるよりも低い傾斜においてオペレータエイリアシングを生じさせる。ＤＭＯプロセスにおけるエイリアシングを回避する１つの方法は、グレイ(Gray ，S．)がFrequency-Selective Design of the Kirchhoff Migration Operator， 40 Geophysical Prospecting 565，(1992)において提唱した計画を実施することによっている。この計画では、トレースをローパスフィルタの時間ドメイン表現を用いて畳込むことによって、周波数濾波済で、位相シフト済の各トレースの複数のローパス濾波済バージョンが作成される。カットオフ周波数は、等間隔セットを形成するように典型的に選択される（例えば、一実施例では周波数は５Ｈｚ乃至100Ｈｚの範囲にわたることができる）。当分野に精通していれば、使用するローパスフィルタの範囲をどのように選択するかは理解されよう。この計画は図５のステップ１３４に示されており、ステップ１３６において１組の濾波済トレースが得られる。以上によりＰＥは、図６の手順に従って１組の置換トレースを発生する準備が整う。上述したように、入力トレース上の各時間毎に、置換トレースの位置に依存するゼロオフセット時間が存在する。各ＰＥは、その入力トレースの源・受信機オフセットｈ、及びＤＭＯ速度νに基づいて、ＤＭＯオペレータのためのアパーチャを決定する。オペレータのアパーチャとは、置換トレースがゼロオフセット時間の関数として生成されるトレースの共通中心点からの最大距離のことである。源・受信機オフセット及び速度はアイソクロンを決定し、また当分野に精通していれば理解されるように、アイソクロンの最大水平広がりを超えたところから発するゼロオフセット反射は存在し得ない。図６のステップ１４０において、各ＰＥは、置換トレースの出力位置を選択する。始めは、これは、解析者が指定する置換トレースを所望する第１の所望出力位置であろう。ＰＥは、この出力トレース上の出力時間サンプルを選択する（ステップ１４２）。始めは、これは解析者が指定する出力トレース上の第１の所望出力時間サンプルであろう。ＤＭＯ速度が一定であるような実施例では、置換トレース上の時間ｔ₀へ移動させられるＮＭＯ修正済入力トレース上の時間サンプルｔ_nが計算される（図６のステップ１４４）。ｔ₀とｔ_nとの間の関係は式（１）によって与えられる。この計算を遂行する上で、ｘは入力トレースに関する源位置と受信機位置との中心点から出力トレースの位置までの距離であり、一方ｈは入力トレース上の源位置と受信機位置との間の距離の半分である。ステップ１４４においては、出力位置θ_dにおけるＤＭＯの勾配も、式（１）のｔ₀をｘについて微分することによって決定される。次に、周波数ｆ_maxが、 Δｘ tanθ_d＝（ν／２ｆ_max）｛１−（ｘ²／ｈ²）｝^1/2 （３）として計算される。好ましい実施例では、式（３）のΔｘは、どちらが大きくとも置換トレース間の間隔であり、νは伝播速度である。速度が一定の場合には、ＤＭＯオペレータ tanθ_dの勾配は式（１）から容易に求められる。置換トレースは、源と受信機とを結ぶ直線上に位置していよう。ｆ_maxの計算に続き、図５のステップ１３６の濾波済トレースの組からｆ_maxを超えない最高周波数を有する濾波済トレースを使用して、ＤＭＯオペレータの出力位置のための出力時間サンプルにおける置換トレース値が求められる（図６のステップ１４６）。当分野においては公知のように、所与の出力時間に対応する入力時間ｔ_nは、サンプルが存在する時間に対応しないかも知れない（即ち、入力トレース上の２つのサンプルされた値の間に存在するかも知れない）。このような場合、出力時間における値は、入力時間ｔ_nを直接ブラケットしている入力トレース上の値の間を線形補間することによって求める。当分野では公知のように、他の補間方法も使用することができる。本発明の特色は、補間された値に正しい振幅重み付けを適用して、鋭い傾斜反射が、それ程鋭くない傾斜反射に対して減衰し始めないようにすることである。振幅倍率（スケーリングファクタ）の決定はブラックらによって示されており、ここでは繰り返さない。補間及び置換トレースサンプルの振幅修正のプロセスはステップ１４８において遂行される。ステップ１５０において、出力位置における出力サンプルのための振幅修正済値がＰＥのメモリ内に格納される。次に、ステップ１５２において、それ以上の出力サンプルが存在するか否かを調べる。もし存在すれば、次のサンプルに関して処理を続けるために、ステップ１５４を通してステップ１４２へ戻る。もし出力トレースのための全てのサンプルが生成されていれば、他の出力位置において何等かの置換トレースを要求しているか否かを調べる（ステップ１５８）。もし要求していれば、次の出力位置に関して処理を続けるために、ステップ１６０を通してステップ１４０へ戻る。もしこの入力トレースのための全ての出力位置に関する置換トレースが生成されていれば、ＰＥは、置換トレースを生成した出力位置に対して責を負うＰＥに信号を送って、置換トレースが使用可能であることを指示する（ステップ１６２）。以下に説明するように、全ての解析ＰＥが、これらの置換トレースを受け取ると、ＰＥは、制御ＰＥによってそれに送られるべき次の入力トレースを待機し（ステップ１６８）、上述したようにステップ１４０から始まるその入力トレースのための置換トレースを生成し始める。解析ＰＥグループ１０８の別のタスクは、グループ内の他のＰＥから置換トレースを読み出すことである。図６のステップ１６４において、ＰＥが、置換トレースが他のＰＥの何れかによって生成されたことの信号を受け取ると、そのＰＥは出力トレースを生成するタスクを開始する。置換トレースは、ＰＥが責を負っている出力位置のための置換トレースを生成した他のＰＥのメモリから読み取られる（ステップ１６６）。もしその出力位置のためのトレースが現在存在していなければ、それは単に置換トレースを受け取って部分出力トレースを作成するだけである。しかしながら、もしその出力位置のための部分トレースが既に存在していれば、到来置換トレースが既存部分出力トレースに付加される。このようにして、ＰＥはキルヒホッフＤＭＯを完了させるために必要な加算ステップを遂行する。本方法の第１段階の結果は、各ＰＥが２つのタスクを遂行することである。第１に、各ＰＥは複数の出力位置のためのＤＭＯ修正済置換トレースを計算し、第２に、各ＰＥは他のＰＥのメモリから置換トレースを読み取り、キルヒホッフ加算をも遂行する。この段階中、個々のＰＥのメモリが一杯になることがあり得る。このようなことが発生すると、ＰＥのグループ１０８内の、メモリが限界に到達した何れかのＰＥは、入力／出力（Ｉ／Ｏ）プロセッサ（図４の１１０）として指定されたプロセッサのグループの１つへ信号を送る。Ｉ／ＯＰＥの数は、典型的にはＤＭＯイメージが計算されるオフセットの数に等しい。解析ＰＥは互いに通信し合って各オフセット毎にメモリ内の最小スタッキングビン数を決定し、次いでこの値をＩ／ＯＰＥに通信する。スタッキングビンの数は、図２を参照して説明済みである。Ｉ／ＯＰＥは、先ず、ディスク１１２（または解析者が選択した他のデータ記憶デバイス）からスタッキングビン内のデータを読み出す。この読み出しは、効率の理由から、大きいデータのブロックに対する要求を使用して最小スタッキングビン番号から開始される。次いでＩ／ＯＰＥはこれらを使用してスタッキングビンの内容を更新し、それらを共通オフセット形状でディスク１１２（図４の１１２）へ書き直し、プロセッサグループ１０８内のＰＥのメモリ位置を完全に空にする。Ｉ／Ｏグループ１１０内の各プロセッサはディスクの異なる部分にアクセスしてグループ１１０内では相互作用は存在しないから、図４のグループ１１０内には破線による相互作用を図示してない。最小スタッキングビン番号を決定し、ディスクから読み出し、更新し、そしてディスクへ書き直すこのプロセスは、解析ＰＥのメモリの置換トレースがフラッシュされてしまうまで繰り返される。これは、全ての解析ＰＥがそれらの解析タスクを完了するまで続けられる。この段階が終わると、スタッキングビンの内容は（ｔ，ｘ，ｙ）順にディスク１１２上にあり、共通オフセットデータボリュームに分離されている。各ＰＥが、広範囲のマイグレーション距離について置換トレースを生成し、フィルタリング及び補間を独立的に遂行し、そして出力位置に割当てられた置換トレースを合計するようになっている本発明は、固有のＰＥ間負荷平衡特性を有している。図８を参照して後述するように、ＰＥに出力位置を割当てることによってメモリの使用及びプロセッサ間の通信の平衡も得られる。本発明の出力段に使用される並列処理システムを図７に示す。ＭＰＰ１００は出力制御プロセッサ１０６を有しているが、プロセッサグループ１０８及び１１０はＰＥの単一のグループ１０９に置換されている。グループ１０９内のどのプロセッサも互いに通信し合わない。そのため、個々のＰＥ間の相互通信は示してない。図７のグループ１０９内の破線は、図示はしてないがより多くのＰＥが存在することを示しているのであって、図示したＰＥ間の相互通信を示しているのではない。更に、グループ１０９内のＰＥは同じタスクを遂行するので、解析ＰＥとＩ／ＯＰＥとの間に区別はない。グループ１０９が遂行するタスクは、ディスク１１２のようなデータ記憶デバイスからデータを読み出し、データを出力制御ＰＥ１０６へ伝送することである。ディスクからの読み出しは、ＰＥが個々にディスク１１２の異なる部分へアクセスして並列に行われる。出力制御ＰＥ１０６はディスク１１２からヘッダレコードを読み出し、これらをＰＥ１０９から受け取った地震データに付加して出力トレースを形成し、出力トレースをテープ１７２、または他の大容量記憶デバイス１７０のような出力データ記憶装置へ書き込む。出力書き込みは、出力媒体上に適切なトレース順を維持するために順次に遂行される。本発明の特色の１つは、処理相を解析ＰＥの間に分散させることによって負荷を平衡させることである。この負荷平衡特性の１つの面は、スタッキングビンの（ｘ，ｙ）格子から、ｘスライスに対する責任を解析ＰＥ間に分散させることである。図８に示す第２の負荷平衡特色では、出力データボリューム内の出力スライスの位置が、出力イメージングのために使用される４つのＰＥだけで示されているが、これに限定されるものではなく単に例示の目的に過ぎない。番号１９０乃至１９３で示されている個々のＰＥが循環的に繰り返されており、各ＰＥはデータボリューム全体を通して規則的な間隔で離間されている出力の部分に対する責を負っている。これにより、各ＰＥにモデルの局部的な部分に対して責を負わせるよりは良好なメモリ使用及びプロセッサ間通信の平衡が得られる。このインタリービングを用いなければ、即ち各プロセッサが３Ｄデータボリュームの異なる局部的部分に対する責を負っていれば、処理の異なる段階においてメモリ使用及びプロセッサ間通信に顕著な変化を生じ得る。例えば、入力データは通常はショットシーケンス順であり、インタリービングを用いない場合には、全てのＰＥは少数のＰＥに割当てられた出力位置のための置換トレースを計算することになろう。これらの少数のＰＥのメモリは急速に一杯になり、全てのプロセッサ間通信はこれらのＰＥに向けられることになろう。従って、異なるＰＥ上の負荷を平衡させることが一層困難になる。当分野に精通していれば理解されるように、このプロセスは、共通オフセットパネル内の入力トレースがショットシーケンス順になっていればより効率的に動作する。データがショットシーケンス順である場合には、ＤＭＯオペレータのアパーチャが制限されているために、所与の出力位置は、その位置に責を負っているＰＥのメモリが一杯になるまで置換トレースが数回合計されるようになる。従って、メモリ位置が数回再使用され、ディスクＩ／Ｏへの必要性が減少し、プロセス効率が増大する。本発明の第２の実施例は、地下の速度が深さと共に線形に、または任意に変化するような状況を扱う。当分野では公知のように、これはより実際的な状況である。ディートリヒ及びコーエンが示したように、ＭＺＯオペレータはクロスライン成分であるから、考え得る出力位置は最早源位置と受信機位置とを結ぶラインに限定されなくなる。更に、垂直速度勾配が存在する場合には、ＭＺＯオペレータは多値になる。一次枝上の第１の値は小さめの傾斜に対応し、一方ＭＺＯオペレータの二次枝上の第２の値は大きめの傾斜及び転倒反射体に対応する。結果は、ケーブルフェザーリングがない場合であっても、トレースが単一の垂直面に限定されずに、単一の地震ラインからの（ｘ，ｙ）格子上の置換トレースになるであろう。しかしながら、ディートリヒ及びコーエンは、クロスライン貢献は強く減衰するから、殆どの実際的な状況では、源と受信機とを結ぶライン上に中心を持つ幅0.4ｈのスウォスにオペレータを制限することで充分であることを示した。以下に説明する実施例はオペレータのインライン部分だけを実施し、クロスライン変位に関する置換トレースは計算しない。後述するように、クロスライン変位を組み入れた他の実施も本発明の範囲内にあるので、これが本発明を限定するものではない。本発明は、ＭＺＯオペレータの多値をも斟酌する。第２の実施例に含まれる計算を図９及び１０に示してあるが、これらは図５及び６に類似している。上述した計算と同じように、プロセスは入力トレースから開始される（ステップ１２０）。図９のステップ２００、２０２、２０４、２０６、２０８、及び２１０は、それぞれ図５のステップ１２１、１２２、１２４、１２６、１２８、及び１３０に関して説明したものと同一の演算を遂行する。換言すれば、大きいオフセットに対しては一定の位相シフトが適用され、小さいオフセットに対しては時間変化位相シフトが適用される。結果は上述したものと同一である。ｄ_minより大きいかまたは等しいオフセットにおいては位相シフトは常に−π／４に留まり、ゼロオフセットの場合には位相シフトは常に０であり、そして０乃至ｄ_minのオフセットに対する位相シフトは小さい時間の場合には−π／４であり、時間が大きくなるにつれて０に接近する。本発明は、時間変化位相シフトのこの特定例に限定されるものではない。当分野に精通していれば、他の方法も明らかであろう。定速度ＤＭＯの場合のように、濾波済トレースの１組のバンドパス済バージョンを計算することによって、エイリアシングは回避される。しかしながら、ＭＺＯオペレータは、濾波済トレース及びそれらの直交成分を必要とする。トレースの「実数」部分は、１組のローパスフィルタを使用してＭＺＯオペレータの一次枝のために使用され（ステップ２１４）、出力トレースが与えられる（ステップ２２０）。ローパスフィルタのヒルベルト変換済みの組を使用して（ステップ２１６）１組の直交トレースが計算される（ステップ２３０）。これらの直交トレースは、ＭＺＯオペレータの二次枝のために使用される。図１０を参照してＭＺＯ修正済トレースの生成を説明する。第１の実施例におけるように、出力トレースの生成は、出力位置及び出力トレース上の時間サンプルを参照して遂行される。しかしながらＭＺＯオペレータは、ＤＭＯオペレータのように簡単な解析式を有していない。処理を簡易化するために、ＭＺＯオペレータは、１組の表値として各解析ＰＥメモリ内に格納されている。この表は、図１０に示すステップに先立って行われる２段階において生成される。第１に、所与の入力位置及び時間に関して、ディートリヒ及びコーエンの開示に基づいて所望出力位置における対応ＭＺＯオペレータ時間が計算される。上述したように、導出された振幅に修正が適用される。当分野においては公知の補間プロセスによって、出力時間から入力時間への「逆方向」表が生成される。これらのオペレータのために使用されるメモリを少なくするために、処理されるオフセットの範囲の中の３つのオフセットのためのマッピングが計算される。これらの計算は、典型的には所望の出力サンプルレートよりも粗い時間サンプリング間隔で遂行される。典型的には、４ミリ秒（ｍｓ）でサンプルされたデータの場合、16ｍｓサンプリングレートでオペレータ値を計算すれば充分であることを見出した。ＭＺＯオペレータ自体に加えて、オペレータの勾配（アンチエイリアス・フィルタリングのために必要）及び振幅倍率（出力トレースの振幅を適切に処理するために必要）が計算され、格納される。ＭＺＯオペレータを３つのオフセットについて計算することに限定されるものではない。明らかに、もし入力データが１つのオフセットしか有していなければ、オペレータはそのオフセットについてのみ計算るだけでよい。図１０は、２つの類似した、しかし順次の計算を含んでいる。第１の計算はＭＺＯオペレータの一次枝のためのものであり（ステップ２４０−２６０）、それに続く第２の計算は二次枝のためのものである（ステップ２７０−２９０）。第１の実施例と同様に、全てのＰＥは、割当てられた入力トレースを使用してこのシーケンスを並列に遂行する。第１の出力位置が選択され（ステップ２４２）、出力サンプル数の初期選択がなされる（ステップ２４４）。ＭＺＯオペレータが本質的には多値であるので、当分野に精通していれば理解されるように、各入力時間毎に２つの出力時間が存在し得る。反対に、各出力時間毎に、入力トレース上に２つの時間が存在し得るが、これらは特定の出力位置に貢献しているプロセッサによって処理される。これらの一方はＭＺＯオペレータの一次枝に対応し、他方は二次枝に対応しよう。所望の出力時間に対応するＭＺＯオペレータが決定される（ステップ２４６）。当分野においては公知のように、これは予め計算された表値をルックアップし、補間を遂行することによって遂行される。この表は一次枝と二次枝とでは異なっていよう。次に、第１の実施例において遂行されたものと同じ手法で、アンチエイリアシング考察に基づいて入力トレースのローパス濾波済バージョンが選択される（ステップ２４８）。ＭＺＯオペレータの一次枝の場合、図９のステップ２２０において、即ちヒルベルト変換を行わずに、濾波済トレースの組から適切なトレースが選択される。ＭＺＯオペレータの二次枝の場合、図９のステップ２３０において、即ちヒルベルト変換が適用されて、濾波済トレースの組から適切なトレースが選択される。再度表をルックアップして振幅倍率が適用され、補間が遂行される（ステップ２５０）。当分野に精通していれば理解されるように、二次枝における勾配及び振幅倍率は、一次枝におけるそれらとは全て異なる。最初に計算されたＭＺＯオペレータの枝からの値は、ＰＥメモリ内に格納される（ステップ２５４）。第２の実施例における他のステップは、図６を参照して説明済の実施例のそれらに類似している。それ以上の出力サンプルが存在するか否かが調べられる（ステップ２５２）。もし存在すれば、処理はステップ２５８を通してステップ２４４へ戻り、次の出力サンプルへ進む。もし存在しなければ、トレースが生成されるそれ以上の出力位置が残っているか否かが調べられる（ステップ２６０）。もし残っていればそれらが発生され（ステップ２６２）、ステップ２４２へ戻る。このシーケンス全体は、二次枝（ステップ２７０乃至２９０）について繰り返され、その結果がステップ２５４において一次枝結果に加えられる。両枝が完了すると、ＰＥは、置換トレースを作成した出力位置に対して責を負っているＰＥに信号を送る（ステップ３０２）。ＰＥは、置換トレースが他のＰＥによって生成されたことを指示する信号をも受け取り（ステップ３０４）、そのＰＥが責を負っている出力位置のための置換トレースを受け取る（ステップ３０６）。次いでＰＥは、それに送られる次の入力トレースを待つ（ステップ３０８）。当分野に精通していれば理解されるように、ＭＺＯオペレータの一次及び二次枝が一緒になる個所で、オペレータの急速な終端が存在する。これは、ディートリヒ及びコーエンの解析が、オペレータの非鏡状(nonspecular)部分（波動理論的には、鏡状(specular)終端点の縁を通って伸びる）を予測しない射線理論の原理に基づいているからである。この点を尖点と呼ぶ。発生するかも知れない問題は、尖点における振幅が大きくなり、それによって存在していない不連続を結果的なイメージ内に導入することである。波動理論的に考察すると、尖点におけるオペレータの急速な終端は存在し得ない。本発明は、オペレータをほぼ非鏡状成分まで人工的に延長することによってこの問題に対処する。当分野に精通していれば、このような延長は波動理論的原理で予測されることが理解されよう。尖点を超えて0.05ｈ乃至0.10ｈの距離だけ延長すると、満足できる結果が得られることを経験的に見出した。この延長は、オペレータの一次枝の二次テイラー級数展開を使用して遂行される（図１０のステップ２４８）。この延長上では、振幅はCos²重み付けによってテーパー付けされ、自然に０にされる。前述したように、ディートリヒ及びコーエンによって導出された振幅は、他の先行技術結果と一致しない。モデル生成された合成地震動記録上では、ディートリヒ及びコーエンによって与えられたＭＺＯ振幅は誤った結果を与える。ＭＺＯオペレータの構造的補強が存在するフレネルゾーンの幅に基づいて、本発明の第２の実施例は、ディートリヒ及びコーエンの振幅ファクタに（ｔ₀）^1/2／ｈを乗算する。これは、ＭＺＯオペレータが時間の増加と共に狭くなるのを、及びオペレータが源・受信機分離の増加と共に広くなるのを補償する。この乗数は、ブラックらのｘ＝０スケールファクタに対応し、このスケールファクタはＭＺＯ解析に関しては厳密には正しくないが、この乗数を使用するとディートリヒ及びコーエンの結果が改善されることを経験的に見出した。ＭＺＯオペレータのクロスライン成分も含まれるような場合への拡張は、本発明の範囲内にある。この拡張により、図１０のステップ２４２（出力位置の選択）において考えるより多くの位置が存在するようになろう。クロスライン成分の生成が、たとえ源と受信機との間のライン上に中心を持つ幅0.4ｈのスウォスに制限されているとしても、典型的な海洋地震応用（240チャネル以上のレコーディングを有する）においては、位置の数は２桁まで増加するようになる。図４のＩ／Ｏプロセッサ１１０は、第２の実施例でも第１の実施例と同一の役割を果たす。同様に、第２の実施例における出力トレースの生成は、図７及び８の第１の実施例に使用したプロセスと同一である。当分野に精通していれば理解されるように、ＰＥの数が増加すると、各ＰＥによって処理されるデータの量が減少する。我々の経験では、もし到来データボリュームを細分し、各々が16または32のＰＥを使用して多重ジョブを走らせれば（ 64乃至128のＰＥを用いてデータボリューム全体のための単一のジョブを走らせるのではなく）最良の総合効率が得られる。後者の場合には、並行処理を発生させる他の効率より通信要求の方が重要になる。ゼロオフセット計算を遂行する多くの方法が存在し、また本発明の範囲がこれらの種々の方法を含み、それらに適用可能であることは理解されよう。

───────────────────────────────────────────────────── フロントページの続き (72)発明者シュナイダーウィリアムエイジュニアアメリカ合衆国テキサス州 77040 ヒューストンテイホーアベニュー 16210 (72)発明者ウィーレムデニスイーアメリカ合衆国テキサス州 77090 ヒューストンルワーレーン 1003

Claims

【特許請求の範囲】１．多重処理要素を有するコンピュータシステム上で地震データを処理するための２段並列処理方法において、ａ）ｉ）少なくとも１つの上記処理要素を入力制御タスクに割当て、 ii）少なくとも複数の上記処理要素を解析タスクに割当て、そして iii）少なくとも１つの上記処理要素を入力／出力タスクに割当てて、入力並列処理段を使用して上記地震データを処理するステップと、ｂ）ｉ）少なくとも１つの上記処理要素を出力制御タスクに割当て、そして ii）少なくとも複数の上記処理要素を読み出し／伝送タスクに割当てて、上記処理済地震データを出力並列処理段を使用して出力フォーマットに変換するステップと、を備えていることを特徴とする方法。２．上記地震データを入力データ記憶デバイスからダウンロードする入力制御タスクを更に備えている請求項１に記載の方法。３．上記地震データを上記解析処理要素へ伝送する入力制御タスクを更に備えている請求項１に記載の方法。４．出力部分を上記各解析処理要素に割り当てる入力制御タスクを更に備えている請求項１に記載の方法。５．上記割当てられる部分は、三次元データボリューム内のインタリーブ済スライスに対応している請求項４に記載の方法。６．上記地震データを前処理する解析タスクを更に備えている請求項１に記載の方法。７．上記各解析処理要素は、上記処理済地震データを少なくとも１つの他の上記解析処理要素へ伝送するようになっている請求項１に記載の方法。８．上記地震データを、ゼロオフセット等価フォーマットに処理する解析タスクを更に備えている請求項１に記載の方法。９．上記各解析処理要素は、ａ）上記地震データを上記制御処理要素から受け取り、ｂ）ゼロオフセット等価フォーマット部分貢献を決定し、そしてｃ）上記部分貢献と、少なくとも１つの他の上記解析処理要素から受け取った上記部分貢献とを合計することによって上記ゼロオフセット等価フォーマットを発生することによって、上記ゼロオフセット等価フォーマットを生成するようになっている請求項８に記載の方法。 10．上記解析処理要素から処理済地震データを受け取り、上記処理済地震データをデータ記憶デバイスに書き込む入力／出力タスクを更に備えている請求項１に記載の方法。 11．上記各入力／出力処理要素は、ａ）上記データ記憶デバイスから上記処理済地震データを読み出し、ｂ）上記読み出したデータと、少なくとも１つの解析処理要素から受け取ったゼロオフセット等価フォーマット部分貢献と組合わせ、そしてｃ）上記組合わされたデータを上記データ記憶デバイスに書き込むことによって、ゼロオフセット等価フォーマットデータを生成するようになっている請求項１０に記載の方法。 12．上記処理済地震データをデータ記憶デバイスから読み出し、上記処理済地震データを上記出力制御処理要素に伝送する読み出し／伝送タスクを更に備えている請求項１に記載の方法。 13．ａ）出力データヘッダレコードを読み出し、ｂ）上記出力データヘッダレコードと、上記処理済地震データとを組合わせて出力フォーマットデータを生成し、そしてｃ）上記出力フォーマットデータを出力データ記憶装置に書き込む出力制御タスクを更に備えている請求項１に記載の方法。