技術研究報告第 45 号 2019.11

戸田建設株式会社

微動アレイ探査における波の位相速度推定手法の概説
OVERVIEW OF METHODS FOR ESTIMATING MICROTREMOR DISPERSION CURVE

成田 修英*1, 保井 美敏*2, 小阪 宏之*1
Nobuhide NARITA, Mitoshi YASUI and Hiroyuki KOSAKA
Recently, the needs for information about soil are increasing in Japan (MLIT notification, 2016). We are being asked for more detailed and more extensive
soil information. In order to meet this need, we have to streamline soil investigation, and in order to do so, we want to use microtremor array exploration.
The microtremor exploration is not a very new method. However, the method has not been applied to exploration of surface soil as a bearing stratum forbuilding foundations until recent years (e.g. IMAI, 2013, IGARASHI, 2016). Therefore, we consider that adequate knowledge of the applicability of thismethod to surface soil has not been gathered. In this report, we overview the methods for microtremor array exploration (especially the methods for
estimating microtremor dispersion curve) and compare and discuss the applicability of these methods to surface soil.
Keywords : Microtremor Array Exploration, Dispersion Curve, Frequency – Wavenumber Method, SPAC Method, Simultaneous Multi – Inputs Wave Field
常時微動アレイ探査，分散曲線，F-K 法，SPAC 法，同時多入力波動場

1. はじめに1.1 背景近年では，杭の支持層確認に対する社会的な要請が厳しくなってきており 1) ，非常に高密度の地盤情報が必要な物件もある．また長周期地震動の検討では，基礎の支持層や工学的基盤よりも，ずっと深い地層の地盤情報が必要となる場合が多い．このように拡大する
地盤情報のニーズに対応するためには，地盤調査の合理化は必須である．
地盤調査の合理化を実現する技術の一つとして，「微動アレイ探査」が考えられる．これは地盤を掘削することに無しに，地表面の常時微動の振動データをパターン解析することにより地下構造を

図１ 常時微動の主な振動源と卓越周波数帯および利用目的推定する技術であり，この技術の活用によって，地盤調査の大幅な合理化が期待される．しかしながら，この技術が一般的な建築工事

浅い地層の調査に微動アレイ探査を適用するためには，同時多入のための地盤調査へ適用され始めたのは比較的最近のことであり 2),

力波動場を想定した振動測定およびパターン解析の計画を立てる

3) ，建築工事のための地盤調査の対象となるような比較的浅い地層必要があるが，前述の通り，そのための知見は十分に蓄積されていへの適用に関して，十分な知見が蓄積されているとは言い難い．

ない．そこで本報では，現在微動アレイ探査で主に使われているパ

一般に，低周波数の微動は深い地層の情報を含み，高周波数の微ターン解析の手法であるF-K 法 4) と SPAC 法(Spatial AutoCorrelation

動は浅い地層の構造を反映する．建築基礎設計で主な関心の対象と

Method)5) について，その基礎となる理論を解説し，同時多入力波動

なる数 10m 以浅を探査の対象とする場合，パターン解析の対象と場への適用性について比較する．

するのは概ね数 Hz 以上の微動であり，その帯域では交通振動など

なお微動アレイ探査におけるパターン解析は，微動の鉛直成分をの人間の活動が主な振動源となるが，自然現象が主な振動源となる用いる場合と水平成分を用いる場合とで解析手順がやや異なるが，低周波の微動と比べると，振動源が安定しない傾向がある（図１）．

ここではより手順が簡便で適用事例の多い，鉛直成分を用いる場合

建築工事の多い都市部では特にその傾向が顕著である．ここで「振のみを取り扱う．

動源が安定しない」とは，主要な振動の発生源が一つに絞れず，時間とともに振動源が移動したり，複数の振動源から同時に微動が発

1.2 微動アレイ探査とは生したりする状況を意味する．以下，このような微動波動場を同時初学者のため，微動アレイ探査がどのようなものかについて，簡多入力波動場と呼ぶ．また逆に，単一の振動源の寄与が卓越してい単に述べておく．

る微動波動場を単一入力波動場と呼ぶ．

図 2 は微動アレイ探査における微動の測定状況の例である．微動アレイ探査では，同図のように複数の振動センサ（加速度計もしく*1 戸田建設㈱技術開発センター 修士（工学）Research and Development Center, TODA CORPORATION, M.Eng.

*2 戸田建設㈱技術開発センター 工学修士

Research and Development Center, TODA CORPORATION, M.Eng.

9-1

微動アレイ探査における波の位相速度推定手法の概説

は速度計）を地表面上に配置し，各センサ間の相関性の分析から，波数とは，波長の逆数に2𝜋をかけたもので，角周波数𝜔における地表面を伝わる微動の速度を推定する．この速度が，前節で「パター波の位相速度𝑐(𝜔)と𝑐(𝜔) = 𝜔/‖𝒌‖の関係にある（ここに，‖𝒙‖はン」と呼んでいたものの正体である．周波数ごとの微動の速度から

ベクトル𝒙のノルム（ベクトルの各要素の 2 乗和平方根））．各周波

地下の弾性波速度構造を推定するのが微動アレイ探査の最終的な数において，F-K スペクトルの振幅がピーク値をとる波数に対応し目標であるが，本報ではそこまで立ち入らずに，その前段階の速度

た速度が，微動の位相速度である．

推定に焦点を当てる．この「周波数ごとの微動の速度」を位相速度したがって，𝑢(𝑡, 𝒓)を測定して(1)式のフーリエ変換を実行できれと呼び，周波数と位相速度の対応関係が作る曲線を分散曲線と呼ぶ．

ば，即座に位相速度が求まることになるが，現在の測定技術で直接的に𝑢(𝑡, 𝒓)を得ることは不可能である．アレイから得らえる離散かつ有限の微動データ𝑢(𝑡, 𝒓𝑖 )（𝒓𝑖 センサ𝑖の位置ベクトルとして）から位相速度推定する方法を考える必要がある．
(2) 位相速度推定のための変形(1)式のフーリエ変換は時間と空間に関するフーリエ変換である．
前節では，このフーリエ変換を実施するのに十分な測定を実施する
ことは現在の技術では不可能だと述べたが，それは空間的な測定点の数と密度を確保できないためで，時間的には十分な測定点の数（測定時間）と密度（サンプリングレート）が確保できる．つまり時間に関しては十分に精度の高いフーリエ変換が可能なので，

図２ 常時微動アレイ探査における微動の測定状況

𝑢(𝑡, 𝒓)を時間に関してのみフーリエ変換した結果を𝑢(𝜔, 𝒓)とおく（以下，このスペクトルを周波数-位置スペクトルと呼ぶ）と，これと F-K スペクトルの関係は，(1)式より，次式の通りに表せる．

2. 位相速度推定の原理∞

図 3 に同時多入力波動場の模式図を示す．同図は𝑁 + 1個のセン𝑢(𝜔, 𝒌) = ∫ 𝑢(𝜔, 𝒓)exp(i𝒌𝒓)𝑑𝒓サから成るアレイに𝑀個の振動源から波が到来状況を示している．
位相速度推定の手法には，この図のように円周上＋円の中心のセン

(2)

−∞

また，この逆変換は次式である．
∞

サ配置を要求するもの（代表例として SPAC 法）と，センサ配置の𝑢(𝜔, 𝒓) = ∫ 𝑢(𝜔, 𝒌)exp(−i𝒌𝒓)𝑑𝒌

(3)

−∞

制限は特に無いもの（代表例として F-K 法）とがある．

微動アレイ探査では，アレイから得らえた微動のデータを時間に関してフーリエ変換し，上記の関係を利用した近似計算によりピー
ク波数を推定し位相速度を求める．F-K 法では，卓越波数の推定に
主に(2)式の関係を利用し，SPAC 法は(3)式の関係を利用する．
2.2 F-K 法(1) センサ配置の制限
F-K 法にはセンサ配置の制限は無い．ただし，センサ配置によって位相速度の推定精度に差が生じるため，でたらめなセンサ配置で良い訳でもない．確保できるセンサの数とサイトの状況が許す限り，
アレイレスポンス（インパルス入力がアレイに入った場合の F-K パ

図３ 同時多入力波動場の模式図

ワスペクトル推定値）6) が少しでもデルタ関数に近くなるようにセンサを配置すべきである．
「アレイレスポンスがデルタ関数に近く

2.1 微動の数学的表現なるよう」だと分かりづらいかも知れないが，直感的にはどの方向(1) 基本形から微動入力が来ても，微動の到達時間差が検知できるようなセン

地表面の位置𝒓 ≔ (𝑥, 𝑦)における時刻𝑡の微動振幅を𝑢(𝑡, 𝒓)と置

サ配置を考えればよい．実用的には，奇数の頂点を持つ正多角形，

くと，この微動の時間 1 次元空間 2 次元に関する 3 重フーリエ変換

特に正 3 角形もしくは正 5 角形の配置が多用される．

は次式の𝑢(𝜔, 𝒌)のようになる．このスペクトルは，周波数-波数ス(2) 位相速度推定方法

ペクトルもしくは F-K スペクトルと呼ばれる（以下，F-K スペクト

F-K 法では，(2)式の積分をそのまま数値計算することによって，ルで呼称を統一する）．
∞

F-K スペクトルの推定値 𝑢̂(𝜔, 𝒌)を計算する．F-K 法における F-K∞𝑢(𝜔, 𝒌) = ∫ ∫ 𝑢(𝑡, 𝒓) exp(−i(𝜔𝑡 − 𝒌𝒓)) 𝑑𝑡𝑑𝒓

スペクトルの推定値は次式で定義される．

(1)

𝑁

−∞ −∞

𝑢̂(𝜔, 𝒌) ≔ ∑ 𝑤𝑖 𝑢(𝜔, 𝒓𝑖 )exp(i𝒌𝒓𝑖 )ここに，𝜔：角周波数，𝒌：波数ベクトル（𝒌𝒓はベクトル𝒌と𝒓の内𝑖=1積で𝒌 = (𝑘𝑥 , 𝑘𝑦 )のとき，𝒌𝒓 = 𝑘𝑥 𝑥 + 𝑘𝑦 𝑦）

ここに，𝑤𝑖 ：数値積分の重み係数，𝑁：センサの総数

9-2

(4)

技術研究報告第 45 号 2019.11

戸田建設株式会社

重み係数𝑤𝑖 を適当に決めてやって，総当たり式の探索アルゴリズ

トルの値を，波数平面上でコンター図として示したものである．F-ム（グリッドサーチ）で|𝑢̂(𝜔, 𝒌)|2 のピークを与える波数探すのが，K パワスペクトルの計算結果は，同図のようなコンター図として示F-K 法の卓越波数推定である．
このとき重み係数の決め方によって，すのが通例である．

F-K 法の中でも更にBFM（Beam Forming Method）とMLM（Maximum

F-K 法では，周波数毎に図５のような結果が得られるので，各周Likelihood Method）とに手法が細分化される．

波数におけるピーク波数から，各周波数での位相速度を求める．

BFM は全ての測定点で𝑤𝑖 = 1とする方法である．前述のアレイ
レスポンスは BFM におけるアレイのインパルス応答であり，これを𝑊(𝒌)とおくと，(4)式より𝑊(𝒌)は次のように表される．
𝑁

𝑊(𝒌) = ∑ exp(i𝒌𝒓𝑖 )

(5)

𝑖=1

BFM は非常に分かり易いが，通常実現可能な数点～１０点前後のアレイでは，あまり良い積分精度は得られない．分かり易いため
F-K 法の説明では良く登場するが，実際に用いられるのは，ほぼ

図５ F-K パワスペクトルの例

MLM のみである．

(4) MLM に用いるクロススペクトルの行列について

MLM の導出の詳細については，本報の主旨からからやや外れる

本報では,(6)式の行列𝑿は微動データのクロススペクトルを並べため付録 3 に示すことにして，ここでは結論だけ述べる．MLM はたものとして説明しているが，行列𝑿を作成するのに，クロススペ

F-K スペクトルの推定値𝑢̂(𝜔, 𝒌)に含まれる波数𝒌以外の波数成分のクトルでなくクロススペクトルをその振幅で割った複素コヒーレ

寄与，つまり波数軸に関するスペクトル漏れの影響を最小化するよンスを用いるべきだとする文献もある．これは MLM の発明者であうな𝑤𝑖 を数値最適化の手法により決定して，F-K スペクトルの推定

る Capon が推奨している手順 4) を踏襲しているものと考えられるを行う手法である（図 4）．

が，理論的には振幅で割る理由はない（付録 3 参照，また Capon も同文献の中で「振幅で割る手順は省いてもよい」と記している）．
Capon はクロススペクトルの正規化を行う理由として，「センサ相互の相対的出力誤差の影響を取り除くため（to remove the effects ofimproper sensor equalization）」としているが，これは単一の振動源の寄与が卓越する波動場を暗に仮定しており，今我々が検討の対象としている同時多入力波動場は想定されていない．同時多入力波動場
では，各振動源から生じる波相互の干渉によって測定点ごとに微動の振幅に違いが出ることが正常であり，正規化によってこの情報を除去してしまうと，逆に位相速度推定の精度が低下する 7) ．また現図４ MLM の考え方

在では，Capon の MLM の論文が発表された時代に比べセンサの精度や扱いやすさが格段に進歩しており，単一の振動源の寄与が卓越

MLM による F-K パワスペクトルは次式で与えられる．
|𝑢̂(𝜔, 𝒌)|2 =1𝒆∗ 𝑿−1 𝒆するような場合にも，クロススペクトルの正規化を行うメリットは(6)

それほど無い．

ここに，上添え字*はベクトルもしくは行列の共役転置（スカラの場

振幅の正規化によって位相速度推定に有用な情報が失われる可合は単なる共役複素数）を表し，𝒆, 𝑿はそれぞれ以下の式で定義さ能性がある以上，クロススペクトルの正規化は行うべきでないと思れる．なおベクトル𝒆はビームステアリングベクトルと呼ばれる．

われる．センサの個体差は気になる場合は，事前にハドルテストを行って適宜キャリブレーション等実施すれば良い．

exp(i𝒌𝒓1 )exp(i𝒌𝒓2 )𝒆≔{}⋮exp(i𝒌𝒓𝑁 )2.3 SPAC 法
F-K 法では，微動の位相速度推定に F-K パワスペクトルのピーク

𝑿𝑢(𝜔, 𝒓1 )𝑢∗ (𝜔, 𝒓1 ) 𝑢(𝜔, 𝒓2 )𝑢∗ (𝜔, 𝒓1 ) ⋯ 𝑢(𝜔, 𝒓𝑁 )𝑢∗ (𝜔, 𝒓1 )𝑢(𝜔, 𝒓1 )𝑢∗(𝜔, 𝒓2 ) 𝑢(𝜔, 𝒓2 )𝑢∗ (𝜔, 𝒓2 ) ⋯ 𝑢(𝜔, 𝒓𝑁 )𝑢∗ (𝜔, 𝒓2 )
≔[
]
⋮
⋮
⋱
⋮
∗
∗
∗
𝑢(𝜔, 𝒓1 )𝑢 (𝜔, 𝒓𝑁 ) 𝑢(𝜔, 𝒓2 )𝑢 (𝜔, 𝒓𝑁 ) ⋯ 𝑢(𝜔, 𝒓𝑁 )𝑢 (𝜔, 𝒓𝑁 )波数を用いたが，SPAC 法では周波数-位置スペクトル𝑢(ω, 𝒓)の空

また，これらの記号を使うと BFM による F-K パワスペクトルはには興味が無い場合が多いので，自己相関関数に着目したアプロー間に関する自己相関関数の卓越波数を用いる．位相速度推定で問題になるのは主に寄与が最大の波数成分であり，それ以外の波数成分
チは自然な発想と言える．建築のシステム同定の分野で言うと，RD

次のように簡潔な形で表現できる．
2∗|𝑢̂(𝜔, 𝒌)| = 𝒆 𝑿𝒆

法や自己相関関数法 8) に対応する手法である．

(7)(3) F-K パワスペクトルの例
図 5 に MLM による F-K パワスペクトル推定値の例を示す．同図は F-K 法のグリッドサーチにおける各格子点での F-K パワスペク

9-3

微動アレイ探査における波の位相速度推定手法の概説(1) センサ配置の制限せで構成されていることが分かっており，その中でも特に基本モー

基本的に，円周上と円の中心にセンサを配置する必要がある．たドと呼ばれる特定の位相速度を持つ波の寄与が卓越する場合が多だしイレギュラーなセンサ配置の研究もされており 9), 10) ，本報でもいので，この仮定（以下，単一速度の仮定と呼ぶ）は多くの微動波

後で半円アレイ 11) に SPAC 法を適用する．

動場で上手く機能する．

(2) 理想的な SPAC 係数

SPAC 法の卓越波数推定は，実測データの空間自己相関より計算される SPAC 係数と整合する𝐽0 (𝑘̃‖𝜸‖)を求めることにより行う．しSPAC 係数とは正規化された空間自己相関係数の方位平均である．
(3)式より，ある位置𝒓とそこから𝜸だけ離れた点の周波数-位置スペ

かし現実には，アレイから得られるデータを使って(9)式の積分を実クトルのクロススペクトル𝑥(𝜔, 𝒓, 𝜸)は次式のように表せるので，行し，空間自己相関関数を求めるのは困難である．F-K 法のような

𝑥(𝜔, 𝒓, 𝜸) ≔ 𝑢(𝜔, 𝒓 + 𝜸)𝑢∗ (𝜔, 𝒓)

アプローチで数値積分を実施する方法も考えられるが，SPAC 法で∞∞は方位平均をとる操作（(10)式）も必要であるため，(9)式の積分も−∞ −∞数値積分で解決しようとすると，必要なセンサ数が多くなりすぎて=∫∫ 𝑢(𝜔, 𝒌1 )𝑢∗ (𝜔, 𝒌2 ) exp(−i(𝒌1 (𝒓 + 𝜸) − 𝒌2 𝒓)) 𝑑𝒌1 𝑑𝒌2(8)

現実的ではない．そこで SPAC 法では，クロススペクトルと空間自これを用いて，空間自己相関関数𝐶(𝜔, 𝜸)は次式のように書ける．

己相関係数の関係を利用する．

∞(3) アレイを用いた SPAC 係数推定𝐶(𝜔, 𝜸) ≔ ∫ 𝑥(𝜔, 𝒓, 𝜸)𝑑𝒓(8)式と(9)式を組み合わせると，次式の関係が得られる．

−∞
∞

=∫

∞

𝑥(𝜔, 𝒓, 𝜸) = 𝐶(𝜔, 𝜸) +∫ 𝑢(𝜔, 𝒌1 )𝑢∗ (𝜔, 𝒌2 ) exp(−i𝒌1 𝜸) 𝛿(𝒌2 − 𝒌1 )𝑑𝒌1 𝑑𝒌2

−∞ −∞

∬

∞

= ∫ |𝑢(𝜔, 𝒌)|2 exp(−i𝒌𝜸) 𝑑𝒌𝑢(𝜔, 𝒌1 )𝑢∗ (𝜔, 𝒌2 ) exp(−i(𝒌1 (𝒓 + 𝜸) − 𝒌2 𝒓)) 𝑑𝒌1 𝑑𝒌2

𝒌1 ≠𝒌2

(9)

(15)

−∞

ここでも単一速度の仮定をおけば，異なる波数を持つ F-K スペク

空間自己相関関数を用いて，SPAC 係数は次のように定義される．
2𝜋
1𝜌SPAC =∫ 𝐶(𝜔, 𝜸)𝑑𝜃
2𝜋𝐶(𝜔, 0) 0

トルは全て異なる方向から来ていることになる．つまり異なる波数

2𝜋
1 ∞
=
∫ 𝑝̃(𝜔, 𝒌) (∫ exp(i‖𝒌‖‖𝜸‖ cos(𝜃 − 𝜑 + 𝜋))) 𝑑𝒌
2𝜋 −∞0ら発生していると考えられるので，(15)式の𝑢(𝜔, 𝒌1 )と𝑢(𝜔, 𝒌2 )はを持つ F-K スペクトル同士は，全て相異なる全く無関係な振動源か
無相関だと考えて良い．したがって，次式が成り立つ．
𝐸[𝑢(𝜔, 𝒌1 )𝑢∗ (𝜔, 𝒌2 )] = 𝐸[𝑢(𝜔, 𝒌1 )]𝐸[𝑢∗ (𝜔, 𝒌2 )]

(10)

(16)

ここに，𝑝̃(𝜔, 𝒌) ≔ |𝑢(𝜔, 𝒌)|2 /𝐶(𝜔, 0)，𝜃：座標原点から見た観測

ここに，𝐸[𝑥]：確率変数𝑥の期待値（「確率変数」に馴染みの無い読点の方位角（𝜸 = (‖𝜸‖ cos 𝜃 , ‖𝜸‖ sin 𝜃 )），𝜑：波数𝒌を持つ F-K ス
ペクトルの到来方向（𝒌 = (‖𝒌‖ cos 𝜑 , ‖𝒌‖ sin 𝜑 )）者は，単にランダムに変動する変数，乱数のことだと考えて良い）(10)式に Hansen の積分表示（付録１参照）の０次（𝑚 = 0）の場

トルの期待値もゼロであり𝐸[𝑢(𝜔, 𝒌1 )] = 𝐸[𝑢∗ (𝜔, 𝒌2 )] = 0である．

一般に微動の振幅の期待値はゼロだと考えられるので，F-K スペク

合を適用すると，積分が一つ消えて次の形となる．

よって(16)式より次式が成り立つ．
𝐸[𝑢(𝜔, 𝒌1 )𝑢∗ (𝜔, 𝒌2 )] = 0∞

̃ (𝜔, 𝒌)𝐽0 (‖𝒌‖‖𝜸‖)𝑑𝒌
𝜌SPAC = ∫ 𝑝

(11)

(17)

ここで(15)式の右左辺の期待値を考えると，右辺の下段の項（以下，−∞

ここに，𝐽𝑚 ：第一種𝑚次の Bessel 関数

この項を無相関項と呼び，𝑁(𝜔, 𝒓, 𝜸)と表記する）がゼロになるのここで(11)式を極座標変換し，微動の性質に関して仮定を一つ付けで，次式が成り立つ．

加えて式を簡略化する．まず極座標変換した(11)式は次のように表𝐸[𝑥(𝜔, 𝒓, 𝜸)] = 𝐸[𝐶(𝜔, 𝜸)](18)(18)式は，理想的な微動波動場においてクロススペクトルの期待せる．
2𝜋

𝜌SPAC = ∫0

値が空間自己相関関数の期待値に等しいことを示している．しかし∞̃(𝜔, ‖𝒌‖, 𝜑)𝐽0 (‖𝒌‖‖𝜸‖)𝑑‖𝒌‖ 𝑑𝜑∫ ‖𝒌‖𝑞(12)

このことは，現実の微動測定データの解析において，クロススペク−∞

ここに，𝑞̃(𝜔, ‖𝒌‖, 𝜑)：極座標変換した 𝑝̃(𝜔, 𝒌)

トルに平均化処理を施した値（異なる時間区間の算術平均（以下，さらに，「ある周波数における微動の伝播は単一の速度で表せる」ものと仮定して，その速度に対応する波数ベクトルのノルムを 𝑘̃ とす

区間平均と呼ぶ）でも周波数平滑化でも良い）と空間自己相関関数

れば，数平滑化などの算術平均で期待値を得るためには無限長のデータ𝛿(‖𝒌‖ − 𝑘̃)𝑞̃(𝜔, ‖𝒌‖, 𝜑) =2𝜋‖𝒌‖の期待値が一致することを意味するものではない．区間平均や周波が必要なのに対し，現実に得られるデータは有限だからである．し(13)

たがって，次式で定義される SPAC 係数の推定値𝜌̂SPAC は，無相関と書けるので，これを(12)式に代入すれば，SPAC 係数は次の簡潔な項の寄与分の誤差を含んだものとなる．
𝑁形で表せることになる．
̃‖𝜸‖)
𝜌SPAC = 𝐽0 (𝑘𝜌̂SPAC ≔(14)

複数の速度を考慮する試みも研究レベルでは行われているが 12), 13) ，実際，地表面を伝わる微動は幾つかの限られた種類の波の重ね合わ1
𝑆[𝑥(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )]∑𝑁𝑆[𝑥(𝜔, 𝒓0 , 0)]𝑖=1𝑁1𝑆[𝐶(𝜔, 𝒓𝑖 − 𝒓0 )] + 𝑆[𝑁(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )]= ∑𝑁
𝑆[𝐶(𝜔, 0)] + 𝑆[𝑁(𝜔, 𝒓0 , 0)]
𝑖=1

9-4

(19)

技術研究報告第 45 号 2019.11

戸田建設株式会社

ここに，𝑆[𝑥]：変数𝑥の有限サンプルによる算術平均，𝒓0 ：円形ア4
ここに，𝑃は振動源の総数を𝑀として，𝑃 ≔ ∑𝑀
𝑙=1|𝑢(𝜔, 𝒌𝑙 )| （振動レイ中心点（図 3 参照）の位置ベクトル

源𝑙に対応する波数ベクトルを𝒌𝑙 とする）
．また，𝜀1 , 𝜀2 , 𝜀3 はそれぞ

実データを用いた SPAC 係数の推定では，𝑆[𝐶(𝜔, 𝑟𝑖 − 𝑟0 )] ≫れ平均化処理によってゼロに近づく項である．区間平均か周波数平𝑆[𝑁(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )]となるよう，十分な数のサンプリングを行わな

滑化あるいはその両方の手段によって𝐴個のサンプルの算術平均かければならない．

ら SPAC 係数の推定値を計算しているものとすると，𝜀1 , 𝜀2 , 𝜀3 のオーダーは Landau のビッグオーを用いて以下の通りに表される．

算術平均に用いるサンプル数が十分であれば，𝑆[𝐶(𝜔, 𝑟𝑖 − 𝑟0 )] ≃𝐸[𝐶(𝜔, 𝑟𝑖 − 𝑟0 )]かつ𝑆[𝑁(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )] ≃ 0であるので，次式が成り立つ．
𝑁𝜌̂SPAC ≃1
𝐸[𝐶(𝜔, 𝒓𝑖 − 𝒓0 )]∑𝑁
𝐸[𝐶(𝜔, 0)]𝑉[|𝐶𝛾 |] + 𝐸[|𝐶0 |2 ]𝜀1 = 𝒪 (√
)
𝐴

(27)

𝑉[|𝐶0 |] + 𝐸[|𝐶0 |2 ]𝜀2 = 𝒪 (√
)
𝐴

(28)

(20)

𝑖=1

さらにアレイ円周上のセンサの数が十分に多ければ，次式も成り立
つので，
𝑁

1
𝐸[𝐶(𝜔, 𝒓𝑖 − 𝒓0 )]1 2𝜋 𝐸[𝐶(𝜔, 𝜸)]
∑
≃
∫
𝑑𝜃
𝑁
𝐸[𝐶(𝜔, 0)]
2𝜋 0 𝐸[𝐶(𝛾, 0)]

(21)

2 2

2

𝑖=1

𝜀3 = 𝒪

このとき𝜌̂SPAC は理論値𝜌SPAC に近似する．
(19)式が理論と整合したオリジナルの SPAC 係数推定式になるが，一般には次の形の推定式のいずれかが用いられることが多い．

√

𝑉 [|𝐶0 |2 + |𝐶𝛾 | ] + 𝐸 [(|𝐶0 |2 + |𝐶𝛾 | ) ]
𝐴

(29)

(
)
単一入力波動場では，空間自己相関関数の定義(9)式より，𝐶0 =2|𝐶𝛾 |かつ𝑃 = 𝐶02 = |𝐶𝛾 | なので，3 つの SPAC 係数推定値（𝜌̂SPAC ,𝑁

1
𝑆[𝑥(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )]𝜌̃SPAC ≃ ∑𝑁
𝑆[|𝑥(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )|]

(22)

𝑁
1
𝑆[𝑥(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )]−𝜌̃SPAC ≃ ∑|𝑆[𝑥(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )]|
𝑁
𝑖=1

(23)

−
𝜌̃SPAC , 𝜌̃SPAC）は全て近似値をとる．一方，同時多入力波動場では，
2

𝐶0 ≠ |𝐶𝛾 |かつ𝑃 ≠ 𝐶02 ≠ |𝐶𝛾 | なので，3 つの推定値は全て異なる値と

𝑖=1なる．このとき，SPAC 係数の定義（(10)式）に整合するのは𝜌̂SPAC−のみであるため，
前述の通り，𝜌̃SPAC と𝜌̃SPACの使用は不適切である．

また単一入力波動場では𝜀1 , 𝜀2 , 𝜀3 は常にゼロとなる（付録 4 も参(22)式， (23)式は，暗に単一入力波動場を仮定している．単一入照のこと）が，同時多入力波動場では，求めたい値である𝐶0 や𝐶𝛾 と力波動場では|𝑥(𝜔, 𝒓0 , 𝒓0 )| = |𝑥(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )|となるため，SPAC

同程度の絶対値を持ち得る．この影響を軽減するためには，十分な−係数の推定に𝜌̃SPAC や𝜌̃SPACを用いても良く，むしろこれらの式で計サンプル数（(27)式から(29)式における𝐴）を確保すること）を用い算することにより，センサ相互の相対的出力誤差が打ち消されてた平均化が必須となる．

SPAC 係数の推定精度が向上する 14) ．しかし同時多入力波動場に(5) SPAC 係数の例

おいては，各入力成分同士の干渉効果により，一般に|𝑥(𝜔, 𝒓0 , 𝒓0 )|≠ |𝑥(𝜔, 𝒓0 , 𝒓𝑖 − 𝒓0 )|となるため，同式の使用は不適切である（「3.

式より，このような結果が得られるので，各周波数について SPAC

同時多入力波動における F-K 法と SPAC 法の適用性」も参照のこ

係数を第一種 0 次の Bessel 関数に当てはめた時の引数を逆算し，位と）．

相速度を求めるのが SPAC 法の位相速度推定である．

図 6 に SPAC 係数（ 𝜌̃SPAC ）の例を示す．(19)式もしくは(22), (23)(4) SPAC 係数推定値と空間自己相関関数の関係
単一速度の仮定に従えば，異なる波数を持つ F-K スペクトル同士
は，全て相異なる無相関な振動源から発生しているものと考えられる．このことと，確率統計の基本的な知識を組み合わせると，SPAC
係数の推定値と空間自己相関関数の関係は以下の通りに整理できる（導出の詳細は付録４に記載する）．以下では，表記の簡略化のため空間自己相関関数𝐶(𝜔, 𝜸)を𝐶𝛾 と書く．
𝑁

𝜌̂SPAC ≃𝐸[𝐶𝛾 ] + 𝜀1
1∑
𝑁
𝐸[𝐶0 ] + 𝜀2

図６ SPAC 係数の例

(24)

𝑖=1
𝑁

1
𝜌̃SPAC ≃ ∑𝑁

𝐸[𝐶𝛾 ] + 𝜀12𝑖=1 √𝐸 [|𝐶𝛾 | ] + 𝐸[|𝐶0 |2 ] − 𝐸[𝑃] + 𝜀3

3. 同時多入力波動場における F-K 法と SPAC 法の適用性
実測データと数値シミュレーションによる模擬微動データを対(25)

象に，F-K 法（MLM）により推定される位相速度と，𝜌̂SPAC , 𝜌̃SPAC ,−𝜌̃SPACより推定される位相速度とを比較する．

𝑁

𝐸[𝐶𝛾 ] + 𝜀1
1−
𝜌̃SPAC≃ ∑𝑁|𝐸[𝐶𝛾 ] + 𝜀2 |

(26)

𝑖=1

9-5

微動アレイ探査における波の位相速度推定手法の概説

3.1 実測例
まず実測において，単一入力波動場に近いと考えられる例と同時
多入力波動場に近いと考えられる例を示す．現実の微動波動場においては，振動源の情報は得られないが，F-K 法より得らえる F-K パ
ワスペクトルのコンター形状より，一方は単一入力波動場に近く，他方は同時多入力波動場に近いと推測した．
測定はどちらの例も 500Hz サンプリングで 20 分間行い，20 分間(a) 14Hz分のデータ全長を 16.348 秒間の小区間に分割，各区間で計算したク(b) 19Hz

図８ 単一入力波動場に近いと考えられるサイトのF-K パワスペクトルロススペクトル（もしくはクロススペクトルの絶対値）の区間平均をとった結果を，F-K 法（MLM），SPAC 法（3 種）に適用した．この際，各小区間のクロススペクトルは，区間平均をとる前に２Hz 幅の Parzen 窓で周波数平滑化を行った．
また，アレイの形状はどちらの例も 2 重正三角形で，
F-K 法
（MLM）
は全点のデータを使用し，SPAC 法は半円形アレイの考え方 11) を適用して小さい方の正三角形プラス１点のみのデータを使用した（図7）．

図９ 単一入力波動場に近いと考えられるサイトの分散曲線(a) F-K 法(b) SPAC 法

図７ 微動の実測に用いたアレイ(a) 12Hz(b) 16Hz

図 10 同時多入力波動場に近いと考えられるサイトの F-K パワスペクトル(1) 単一入力波動場に近いと考えられる例
これは長崎県某所で測定した例であり，小さい方の正三角形の１辺は𝑙 =1.25m である．
図 8 にこのサイトでの F-K スペクトルの例，図 9 に各手法より推定された分散曲線を示す．図 8 より，このサイトでは単一の波数成分の寄与が卓越しており，このことから単一入力波動場に近い条件になっていることが推測される．また図 9 より，このサイトでは FK 法と SPAC 法（３種）の分散曲線は概ね整合する．
(2) 同時多入力波動場に近いと考えられる例

図１１ 同時多入力波動場に近いと考えられるサイトの分散曲線

これは神奈川県某所で測定した例であり，小さい方の正三角形の１辺は𝑙 =1.15m である．
図 10 にこのサイトでの F-K スペクトルの例，図 11 に各手法より

3.2 数値シミュレーション

推定された分散曲線を示す．図 10 より，このサイトでは多方向に(1) 模擬微動データの作成方法

スペクトルのピークが分布しており，このことから同時多入力波動

微動の性質について，以下の仮定をおく．① 微動波動場は𝑀個の場に近い条件になっていることが推測される．また図 11 より，こ有限な波数成分の重ね合わせとして表される．② 各波数成分におのサイトでは F-K 法，SPAC 法（3 種）の位相速度はあまり整合しけるセンサ間の位相差は，各成分の持つ波数とセンサの位置により−ていない． 𝜌̃SPAC
は特に目立って他の手法と合わない傾向にあるが，確定的に定まる．③ 各波数成分は時間軸についてホワイトノイズ

F-K 法と SPAC 法の残りの２種もそれほど合っていない．PS 検層やと見做すことが出来る．

表面波探査など他の物理探査は用いていないため，どの手法が最もこの仮定に基づけば，基準となる位置の微動として𝑀個の独立な実地盤の特徴を最も良く反映しているのか決めることは出来ない

ホワイトノイズを作成し，それに各成分の波数から定まる基準点とが，10Hz 以上の帯域で F-K 法と SPAC 法（3 種）を比較すると，F-

センサ位置との位相差を付加することにより，所要の波数成分を持K 法は比較的値が安定しているのに対し，SPAC 法の位相速度は値

つ模擬微動データを作成できる．図 12 は模擬微動データ作成のフの変動が大きく，地表面を伝わる波の位相速度としては F-K 法の結ローである．

果の方が自然に思われる．

9-6

技術研究報告第 45 号 2019.11

戸田建設株式会社

ここに，𝜀：𝑿の各要素の絶対値に対して十分に小さい正の実数（ダンピングファクターと呼ばれる）, 𝑰：単位行列
𝑿の各行において，（非対角項）／（対角項）の比は各センサ間のコヒーレンスを表しているので，(30)式左辺の計算は，𝜀の分だけ人為的に各センサ間の相関を下げる操作になっている．したがって，𝜀の値を大きくしすぎると計算の精度が低下するデメリットはある．
ただし，𝑿でなく𝑿を用いることで各センサ間の相関が過小評価されることの２次的な効果として，det(𝑿) ≃ 0の場合でも確実にdet(𝑿) ≫ 0となるため，MLM の計算は安定する．デメリットを最小限に抑えるため，𝜀の値は計算結果が安定する最小限に設定する場合が多い．
今回のシミュレーションでは，クロススペクトル行列𝑿の各要素の絶対値の平均を𝜇𝑋 として，𝜀 = 𝜇𝑋 × 10−5 を用いる．
(4) 模擬微動データを用いた位相速度推定
図14 に，
作成した模擬微動データより推定した位相速度を示す．
同図(a)より単一入力の場合は全ての方法で正しい位相速度推定が出来ているが，同図(b)より同時多入力の場合は F-K 法と 𝜌̂SPAC を用いた SPAC 法のみが設定した位相速度 100m/s を正しく推定できる結果となった．
SPAC 法の推定値の中で 𝜌̂SPAC を用いた推定値のみ正解を与える
結果となったのは，「2.3 (3) アレイを用いた SPAC 係数推定」で述
べた通り，同時多入力波動場では 𝜌̂SPAC のみが SPAC 法の理論に適

図 12 模擬微動データ作成のフロー−合するためである．また 𝜌̃SPAC を用いた結果と𝜌̃SPACを用いた結果ここに，𝒓𝑐：基準点の位置，𝑎𝑙： 𝑙番目の波数成分の周波数-位置ス

で比べると，𝜌̃SPAC を用いた結果の方が比較的正解に近い．これは，ペクトルの絶対値，𝜙：一様乱数（0 < 𝜙 < 1），𝜔𝑛𝑦𝑞：ナイキスト−(25)式と(26)式の比較より， 𝜌̃SPAC の方が𝜌̃SPAC
よりも 𝜌̂SPAC に近い

角周波数，𝛥𝜔: FFT の周波数刻み，ℱ −1 [𝑥]：変数𝑥の逆フーリエ変ためと考えられる．

換

F-K 法は基本的に単一入力波動場でも同時多入力波動場でも問題(2) シミュレーションの条件

ないが，同時多入力波動場では低周波数側で位相速度が過大評価さ単一入力波動場と同時多入力波動場，それぞれ１ケースずつシれる傾向にある．これは複数の入力に対応した複数の F-K パワスペ

ミュレーションを実施する．図 13 にシミュレーション条件の模式

クトルのピークの裾が重なり合うことにより，重なり合った裾の振図を示す．位相速度は一律 100m/s とし，同時多入力波動場の波数幅がピークの振幅を上回って偽のピークが発生し，位相速度の推定成分の数は２つ，振幅が等しく進行方向が逆向きの成分とする．

に誤差が生じることに起因する．この現象は F-K スペクトルの縮重

また，アレイは半径 1ｍの正５角形＋中心の配置とし，模擬デーと呼ばれる 16) ．低周波数側でスペクトルの縮重が生じる（生じ易タのサンプリングレートは 1kHz, データ長は 16.384 秒とする．

い）のは，低周波数側では波数が小さく各測定点間での位相差が小さいため相対的に SN 比が低下することに起因する．このシミュ
レーションに関しては模擬データそのものにはノイズは入っていないが，ダンピングファクター（(付 30)式の𝜀）由来のノイズ成分によりスペクトルの縮重が発生している．

(a) ケース１(b) ケース２

図 13 シミュレーション条件(3) F-K 法（MLM)の適用方法
このシミュレーションでは模擬微動にノイズは混入させないので，図 13 の条件下では，クロススペクトルの行列𝑿の行列式がゼロに非常に近くなってしまい，
定義通りの MLM の計算は難しい 15) ．

(a) ケース１

このような場合には，𝑿の代わりに次のような行列𝑿を用いて MLM

図 14 模擬微動データより推定された分散曲線の計算を行う 15) ．
𝑿 ≔ 𝑿 + 𝜀𝑰

(30)

9-7

(b) ケース２

微動アレイ探査における波の位相速度推定手法の概説考えられる．しかし実測データを用いた検討では， 𝜌̂SPAC による結2) 今井博, 青木智幸, 石井裕康 「微動チェーンアレーによる浅部探査 十数メートル程度までの深度に関する微動探査の適用性」，大成建設技術セ
ンター報, pp. 26_1－26_7, 2013
3) 五十嵐治人, 伊藤仁 「微動チェーンアレー探査による支持地盤の推定」,
錢高組技報, pp. 17－22, 2016

果は 𝜌̃SPAC による結果よりも悪くなっているように見える．これは4)

3.3 適用性検討のまとめ
数値シミュレーションの結果より（理論的にも），SPAC 法の中で同時多入力波動場に適用可能なのは 𝜌̂SPAC を用いた方法だけだと

J. CAPON “High-resolution frequency-wavenumber spectrum analysis”,
Proceedings of the IEEE, Vol. 52, No. 8, pp. 1408－1418, 1968

𝜌̂SPAC が実データに含まれる種々のノイズ，特にセンサ相互の相対的出力誤差の影響を受けやすい傾向にあることに起因していると5)−考えられる．この弱点の改善のために 𝜌̃SPAC や 𝜌̃SPACを用いた方法

Keiiti AKI “Space and time spectra of stationary stochastic waves, with specialreference to microtremors”, Bulletin of the Earthquake Research Institute, Vol.
35, pp. 415－456, 1957

が提案された訳だが 14) ，上に示した通り，同時多入力波動場で−は，𝜌̃SPAC も𝜌̃SPAC
も使うことはできない．したがって現行の SPAC

6)

Masanori HORIKE “Inversion of phase velocity of long-period microtremors

法の計算は，同時多入力波動場に対する適用には向いていないと考

to the s-wavevelocity structure down to the basement in urbanized areas”,

えられる．

Journal of Physics of the Earth, Vol. 33, pp. 59－96, 1985

一方 F-K 法の方は，同時多入力波動場でも比較的安定した結果が7)

成⽥修英, 保井美敏, 小阪宏之 「常時微動アレイ探査による傾斜基盤
サイトの支持層推定」, 第 53 回地盤工学研究発表会梗概集, pp. 1967－

得られる．ただしスペクトルの縮重には注意する必要がある．スペ

1968, 2018

クトルの縮重が疑われる場合には，アレイのサイズを拡大する（(付11)式の𝑟𝑚𝑎𝑥 を大きくする），測定時間を長くしてデータの SN 比を

8)

日本建築学会 「建築物の減衰」, 日本建築学会, pp. 65－96, 2000

向上させる（(付 11)式の𝜙𝑚𝑖𝑛 を小さくする）などして F-K 法の波数9)Michael W. ASTEN “On bias and noise in passive seismic data from finite
circular array data processed using SPAC methods”, GEOPHYSICS, Vol. 71,

分解能を上げる必要がある．

No. 6, pp. V153－V162, 2006

4. おわりに10)

白石英孝, 松岡達郎 「Lamb の問題に基づくレーリー波複素コヒーレ
ンス関数の離散定式とその応用-空間自己相関法の新しい解釈-」, 物理

4.1 得られた知見

探査, Vol. 58, pp. 137－146, 2005

本報では，F-K 法と SPAC 法の原理について説明し，同時多入力
波動場への適用性について検討した．同時多入力波動場への適用性

11)岡⽥広, 松岡達郎, 白石英孝, 八戸昭一 「微動アレー探査のための空間自己相関法：半円形アレーの適用について」, 物理探査学会学術講演

検討では，以下の知見を得た．

梗概集, pp. 183－186, 2003

まず SPAC 法は，本報で検討対象とした３種類の計算方法の何れについても，同時多入力波動場への適用には向いていない．一方，

12)

Kohji TOKIMATSU, Shuji TAMURA, and Hisaya KOJIMA “Effects ofF-K 法は，低周波数側でスペクトルの縮重が生じる場合があるもの

multiple modes on Rayleigh wave dispersion characteristics”, Journal ofの，高周波数側では同時多入力波動場でも安定した結果が得られる

Geotechnical Engineering, Vol. 118, No. 10, pp. 1529－1543, 1992

ため，同時多入力波動場への適用性は比較的良い．したがって，浅13)

池⽥達紀, 松岡俊文, 辻健, 林宏一 「SPAC 法における異なる相関距
離を考慮したマルチモード解析」, 物理探査, Vol. 64, pp. 127－138, 2011

い地層の探査など同時多入力波動場が想定される探査では F-K 法を採用すべきだと考えられる．

14)

岡⽥広, 松島健, 日高英二 「長周期微動に含まれる表面波の位相速度推定について: 空間自己相関法と周波数波数法の比較」, 北海道大学地
球物理学研究報告, Vol. 49, pp. 53－62, 1987

4.2 今後の課題
F-K 法のピーク波数探索アルゴリズムはグリッドサーチであるた

15)Michael W. ASTEN and J. D. HENSTRIDGE “Array estimators and the use

め，波数の探索範囲や探索刻みの設定に地盤に関する先験的な情報

of microseisms for reconnaissance of sedimentary basins”, GEOPHYSICS, Vol.

もしくは解析者の試行錯誤が要求される．また探索対象の格子点の49, No. 11, pp. 1829－1837, 1984

数だけデータが出力されるため，解析後の後処理も SPAC 法に比べ16)

岡⽥広, 凌甦群, 佐藤洋, 笹谷務, 宮腰研, 斎藤誠治, 石川顕, 南雲秀て手間がかかる．この問題の解決のため，以下の検討課題が考えら

樹 「微動利用の地下構造推定における周波数・波数スペクトル法と空れる．

間自己相関法の比較検討」, 物理探査学会学術講演論文集, pp. 105－109,
1995① より効率的に F-K 法を実施するアルゴリズムを開発する．②
単一入力波動場では SPAC 法を利用するため，単一入力か同時多入17)

成⽥修英, 保井美敏, 山本健史, 小阪宏之 「FK パワースペクトル・
ピーク波数の直接推定法」, 第 52 回地盤工学研究発表会梗概集, pp.

力か簡便に判定する手法を開発する．このうち①については，すで1775－1776, 2017に検討に着手しており一定の成果を得ている 17), 18) が, 実用的なアルゴリズムの完成には未だ至っていない．

18)

成⽥修英, 保井美敏, 小阪宏之 「常時微動アレイ観測記録に対する主成分分析の適用性検討」, 第 15 回日本地震工学シンポジウム, pp. 206－215, 2018

参考文献
1) 国土交通省「基礎ぐい工事の適正な施工を確保するために講ずべき措置の実施に向けて」, aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
http://www.mlit.go.jp/totikensangyo/const/totikensangyo_const_tk1_000116.htm
l（2019 – 7 - 24 参照）

19)

例えば 森口繁一, 宇⽥川銈久, 一松信 「岩波数学公式３ 特殊関数」,
岩波書店, 1987

9-8

技術研究報告第 45 号 2019.11

20)戸田建設株式会社

例えば A. N. KOLMOGOROV, I.G. ZHUBENKO and A. V. PROKHOROV

𝑆[𝑥(𝜇𝑥 , 𝜎𝑥 )] = 𝜇𝑥 +（丸山哲郎, 馬場良和 訳）「コルモゴロフの確率論入門」, 森北出版株

21)𝑥 ′ (0, 1)(付 9)式会社, 2003

ここに，𝑥 ′ (0, 1)：(付 8)式の𝑥(0, 1)とは別の平均ゼロ，分散１の確

宮腰研 「微動探査法による表面波位相速度推定の基礎的研究」, 北海

率変数(付 9)式より，シグナルが𝜇𝑥 ，ノイズの分散が𝜎𝑥2 のデータ𝐴個の算道大学大学院理学研究科博士論文, 1998
22)

𝜎𝑥
√𝐴

術平均は，Landau のビッグオーを用いて次のように書ける．
𝜎𝑥
𝑆[𝑥(𝜇𝑥 , 𝜎𝑥 )] = 𝜇𝑥 + 𝒪 ( )(付 10)√𝐴

Hiroshi ARAI and Kohji TOKIMATSU “S-wave velocity profiling by joint
inversion of microtremor dispersion curve and horizontal-to-vertical (h/v)
spectrum”, Bulletin of the Seismological Society of America, Vol. 95, No. 5, pp.
1766－1778, 200523)

付録 2. 適用性の補足：ピーク波数検出限界

R. T. LACOSS “Data adaptive spectral analysis method”, GEOPHYSICS, Vol.

位相差の検出分解能には限界があるので小さすぎる波数は検出

36, No. 4, pp. 661－675, 1971
24)

25)

Ikuo CHO, Taku TADA, and Yuzo SHINOZAKI “A new method to determine

できないし，波数の検出に周期関数を利用しているので最小センサ

phase velocities of Rayleigh waves from microseisms”, GEOPHYSICS, Vol. 69,

間隔間の位相差が１サイクルを超えるような波数も基本的に検出

No. 6, pp. 1535－1551, 2004

不能である（この問題は空間エイリアシング 15） と呼ばれている）．

先名重樹 「地震防災のための地盤情報と今後の展開」, MSS 技報, Vol.

つまり単一のセンサ配置で検出できるピーク波数の範囲，すなわち24, 2014

位相速度を推定できる周波数帯域には一定の限界がある．したがって微動アレイ探査の実施においては，探査対象周波数（深度）と探付録 1. 本報で用いた数学の基礎知識査可能な周波数が一致するようにセンサ配置を計画する必要があ

付 1.1 Hansen の積分表示

る．単一のセンサ配置では探査対象周波数を全てカバーできないこ

第 1 種 Bessel 関数は指数関数を用いて次のように表すことが出とも多いが，そのような場合には，センサの配置を変えながら複数来る．これは Hansen の積分表示と呼ばれる 19) ．
1 2𝜋
𝜋
𝐽𝑚 (𝑧) =∫ exp (i (𝑧 cos 𝜃 − 𝑚 (𝜃 + ))) 𝑑𝜃2𝜋 02回の測定を実施する．
(付 1)
付 2.1 小波数（長波長）側の限界
アレイの最大センサ間隔（SPAC 法の場合は解析に使用するアレ付 1.2 確率変数の基礎知識

イ半径）を𝑟𝑚𝑎𝑥 とし，アレイと測定対象サイトの波動場の性質およ

確率変数の性質として本報で用いるのは，以下の性質のみである．
以下の性質に関する説明，証明については，確率論の教科書を参照

び測定条件（測定時間，サンプリングレート等）によって定まる２測定点間の位相差の検出下限値を𝜙𝑚𝑖𝑛 とする．このときアレイにされたい 20) ．

よって検出できる最小波数𝑘𝑚𝑖𝑛 は次式で表される．

(1) 独立な確率変数の性質
𝑘𝑚𝑖𝑛 ≥

２つの確率変数𝑥, 𝑦が独立であり，確率変数𝑥の期待値を𝐸[𝑥]，分散を𝑉[𝑥]で表すとき，以下の関係が成り立つ（(付 2)式については，

𝜙𝑚𝑖𝑛
𝑟𝑚𝑎𝑥

(付 11)

𝜙𝑚𝑖𝑛 の一つの目安として，F-K 法で2π⁄3 ～2𝜋/5程度，SPAC 法

独立でなくても成り立つ）．
𝐸[𝑥 + 𝑦] = 𝐸[𝑥] + 𝐸[𝑦](付 2)

でπ⁄3 ～𝜋/5程度という値が知られている 21), 22) ．ただし，𝜙𝑚𝑖𝑛 は
基本的にデータの SN 比によって定まる値であり，この値は測定に𝐸[𝑥𝑦] = 𝐸[𝑥]𝐸[𝑦](付 3)用いた機器の特性や測定対象の波動場の性質によって大きく変動

𝑉[𝑥] = 𝐸[𝑥 2 ] − (𝐸[𝑥])2(付 4)

し得るので，文献値は測定計画を立てる際などに用いる参考値と考𝑉[𝑥 + 𝑦] = 𝑉[𝑥] + 𝑉[𝑦](付 5)

えるべきである．

また，確率変数𝑥の定数倍𝑎𝑥について，次式が成り立つ．
𝑉[𝑎𝑥] = 𝑎2 𝑉[𝑥](付 6)付 2.2 大波数（短波長）側の限界

さらに分散が全て𝑉𝑥 である互いに独立な確率変数𝑥1 , 𝑥2 . ⋯ . 𝑥𝑁 があ

空間エイリアシングを避けるためには，
F-K 法では複素正弦波が，るとき，その算術平均の分散は，(付 5)式と(付 6)式を組み合わせてSPAC 法では第１種０次の Bessel 関数が，それぞれ１価関数として

次のように書ける．

扱える波数の範囲内にのみ限定して解析を行う必要がある．

𝑥1 + 𝑥2 + ⋯ + 𝑥𝑁
𝑁𝑉𝑥 𝑉𝑥𝑉[]= 2 =𝑁𝑁
𝑁

したがってアレイの最小センサ間隔（SPAC 法の場合は解析に使(付 7)用するアレイ半径）が𝑟𝑚𝑖𝑛 のとき，アレイによって検出できる最大(2) 標準化と算術平均

波数𝑘𝑚𝑎𝑥 は次式で表される．

期待値𝜇, 分散𝜎 2 の確率変数を𝑥(𝜇, 𝜎)と書くとき，任意の𝑥(𝜇, 𝜎)は，平均ゼロ，分散 1 の確率変数𝑥(0, 1)を使って次のように表せる．
𝑥(𝜇𝑥 , 𝜎𝑥 ) = 𝜇𝑥 + 𝜎𝑥 𝑥(0, 1)(付 8)

𝐴個のサンプルを使った𝑥(𝜇𝑥 , 𝜎𝑥 )の算術平均を考えると，(付 7)式，(付 8)式より，次式が成り立つ．

9-9𝑘𝑚𝑎𝑥 ≤

𝜙𝑚𝑎𝑥
𝑟𝑚𝑖𝑛

2𝜋 (F − K 法)
𝜙𝑚𝑎𝑥 = {𝑗1,1 ≃ 1.22𝜋 (SPAC 法)(付 12)(付 13)

微動アレイ探査における波の位相速度推定手法の概説

ここに，𝜙𝑚𝑎𝑥：各手法の評価関数が１価関数として扱える最大の引ターゲット項のパワを保持しつつ，漏れ項の影響を低減するような数，𝑗𝑛,𝑘 ：第１種 n 次の Bessel 関数の第𝑘零点

𝑤𝑖 を決める必要がある．

ただし(付 13)式は最小センサ間隔が𝑟𝑚𝑖𝑛 である以外は完全な測定

そこでターゲット項のパワを保持するための条件を考える．これ

データが得られている場合を想定した理想的な上限であり，実現可は，（ターゲット項）= 𝑢′ (𝜔, 𝜿)∑𝑤𝑖 であることより，∑𝑤𝑖 が定数で能な有限のアレイより得られるデータでは，限界値はこれよりも小あれば良い．この条件，つまり∑𝑤𝑖 =定数（ここでは∑𝑤𝑖 = 1とすさくなる．実アレイにおける限界は用いるアレイと評価対象の波動る）の条件下で，F-K パワスペクトルを最小化すると，ターゲット

場の性質によって変化するが，例えば二重正三角形アレイでは，次項は保持されて漏れ項の影響のみが最小化される．したがって，F-に示す通りの値となることが知られている 21) ．

K パワスペクトルに関する次式のような条件付き最小化問題を解2𝜋けば，求める最尤の𝑤𝑖 が得られる．

(F − K 法)(付 14)

𝜙𝑚𝑎𝑥 ≃ {√3
𝜋 (SPAC 法)

𝑁

最小化

𝐺 ≔ |𝑢̂ (𝜔, 𝜿)|2 + 𝜆 (∑ 𝑤𝑖 − 1)

(付 19)

𝑖=1

付録 3. MLM の導出

ここに，𝜆：Lagrange の未定乗数

MLM すなわち Maximum Likelihood Method とは，確率過程のパ

これを整理して解きやすくするために，次のベクトルを定義する．

ラメータ推定における最尤法の事を指し，元々は地表面の微動を確

𝑤1 exp(i𝜿𝒓1 )𝑤2 exp(i𝜿𝒓2 )𝒗≔{}⋮
𝑤𝑁 exp(i𝜿𝒓𝑁 )率過程として捉える立場から出てきた手法であるが，Lacoss の解釈
23) を応用すれば，本文(4)式の重み係数𝑤 を決定する手順の一つと𝑖しての解釈も成り立つ．この解釈に基づく MLM の導出が，恐らく(付 20)

最も予備知識を必要としない導出手順になるので，以下ではその説

これとクロススペクトル行列𝑿およびビームステアリングベクトル

明を行う．

𝒆を組み合わせると，F-K パワスペクトル推定値と∑𝑤𝑖 は以下のよ

ある特定の波数𝜿における F-K スペクトルを評価することを考えうにも書ける．

る（以下， 𝜿をターゲット波数と呼ぶ）．ターゲット波数𝜿における|𝑢̂(𝜔, 𝜿)|2 = 𝒗∗ 𝑿𝒗

最適な重み係数を決定するために，まず周波数-位置スペクトルを𝑁

ターゲット波数𝜿が寄与する部分とそれ以外とに分離する．波数平∑ 𝑤𝑖 = 𝒆∗ 𝒗

面上でターゲット波数𝜿近傍の微小領域を除く領域を𝑅とすると，

𝑖=1

本文(3)式より，周波数位置スペクトルは次のように分離できる．

𝐺 = 𝒗∗ 𝑿𝒗 + 𝜆(𝒆∗ 𝒗 − 1)

(付 15)(付 23)(付 23)式の𝐺を最小化するベクトル𝒗が決定されると，(付 20)式より，𝑤𝑖 も自動的に決定されるので，𝐺を最小化する𝒗をまず求めるこ𝑅

ここで表記の簡略化のため，とにする．𝐺を最小化する𝒗は，𝐺の𝒗に対する勾配がゼロになる点

𝑢′ (𝜔, 𝜿) ≔ 𝑢(𝜔, 𝜿)𝑑𝒌

(付 16)

𝑅(𝜔, 𝜿) ≔ ∫ 𝑢(𝜔, 𝒌) exp(−i𝒌𝒓) 𝑑𝒌(付 17)の𝒗であるので，次式を満たす．
𝜕𝐺= 𝒗∗ 𝑿 + 𝜆𝒆∗ = 0𝜕𝒗

𝑅とおき，(付 15)式を F-K スペクトルの推定式（本文(4)式）に代入す(付 24)

これを整理すると，未定乗数𝜆について以下の２つの解が得られる．

ると，ターゲット波数の寄与とそれ以外の部分に分離された F-K ス

𝜆 = −𝒗∗ 𝑿𝒆, −ペクトルの推定式が次のように得られる．

𝒗∗ 𝒆𝒆∗ 𝑿−1 𝒆(付 25)

このうちひとつめの解は𝒗 = 0となる解であり，𝑤𝑖 = 0になってし𝑁

𝑢̂(𝜔, 𝜿) = 𝑢′ (𝜔, 𝜿) ∑ 𝑤𝑖
𝑖=1

(付 22)

これを用いると，最小化の目的関数𝐺 は次式のように整理できる．

𝑢(𝜔, 𝒓) = 𝑢(𝜔, 𝜿)𝑑𝒌 exp(−i𝜿𝒓)+ ∫ 𝑢(𝜔, 𝒌) exp(−i𝒌𝒓) 𝑑𝒌(付 21)

まうので意味が無い．２つめの解について考えると，この解の分子𝒗∗ 𝒆はベクトル𝒗, 𝒆の定義より𝒗∗ 𝒆 = ∑𝑤𝑖 なので，制約条件∑𝑤𝑖 = 1(付 18)

より１でなければならない．つまり，次式が成り立つ．

𝑁+𝑅(𝜔, 𝜿) ∑ 𝑤𝑖 exp(i𝜿𝒓𝑖 )

𝜆=−

𝑖=1

(付 18)式下段の項（以下，この項を漏れ項と呼び，上段右辺の項

1(付 26)𝒆∗ 𝑿−1 𝒆

これを(付 24)式に代入すると，𝐺を最小化する𝒗が求まり，それをさらに(付 21)式に代入すると|𝑢̂(𝜔, 𝜿)|2 が定まる．

をターゲット項と呼ぶ）が𝑢̂(𝜔, 𝜿)に与える影響を最小化するような𝑤𝑖 が上手く設定できれば，F-K スペクトル推定値のスペクトル漏𝒗𝑜𝑝𝑡 =れ（本文図 4 参照）も最小化され，最尤の F-K スペクトル推定が実
現する．

𝑿−1 𝒆

|𝑢̂(𝜔, 𝜿)|2 =

漏れ項を最小化する𝑤𝑖 を単純に考えるとすぐに𝑤𝑖 = 0が思いつ
くが，これだとターゲット項もゼロになってしまうので意味が無い．

9-10(付 27)𝒆∗ 𝑿−1 𝒆
1𝒆∗ 𝑿−1 𝒆

ここに，𝒗opt ：𝐺を最小化する𝒗(付 28)

技術研究報告第 45 号 2019.11

戸田建設株式会社

以上より，MLM による F-K パワスペクトル推定値が導かれた．
また(付 27)式および𝒗の定義（(付 20)式）より，重み係数𝑤𝑖 は次式のようになる．
∑𝑁
𝑗=1 𝑞𝑖𝑗 exp(i𝜿(𝒓𝑗 − 𝒓𝑖 ))𝑤𝑖 = 𝑁∑𝑖=1 ∑𝑁
𝑗=1 𝑞𝑖𝑗 exp(i𝜿(𝒓𝑗 − 𝒓𝑖 ))∗𝑁(𝜔, 𝒓, 𝜸) = ∑ 𝑢𝑙 𝑢𝑚exp(−i(𝒌𝑙 (𝒓 + 𝜸) − 𝒌𝑚 𝒓))(付 34)𝑙≠𝑚(2) F-K スペクトルの累乗の期待値に関する仮定(付 29)

クロススペクトルの２乗の期待値およびクロススペクトルの４(付 29)式に測定データに依存する項𝑞𝑖𝑗 が含まれていることから分乗の期待値の評価をするため，次の仮定をおく．この仮定が成り立かるように，MLM における重み係数𝑤𝑖 は測定データに応じて特性つのであれば，評価の手順は大幅に簡略化される．
𝐸[|𝑢𝑙 |𝑚 (𝑢𝑙 )𝑛 ] = 0が変化する適応フィルタとしての性質を持つ．

(𝑚 ≥ 0, 𝑛 ≥ 1)

(付 35)(付 35)式の仮定は，単一速度の仮定が成り立ち，異なる波数をもつ

付録 4. SPAC 法に現れるクロススペクトルの算術平均評価

F-K スペクトル同士が全て互いに独立した別の振動源から発生してここではクロススペクトルおよびその２乗の算術平均と空間自いるのであれば，自然に成り立つと考えられるものである．以下，

己相関関数の関係について考える．以下，表記の簡略化のため，ク
ロススペクトル𝑥(𝜔, 𝒓, 𝜸)を𝑥𝛾 , F-K スペクトル𝑢(𝜔, 𝒌𝑙 )を𝑢𝑙 と書く．
(付 10)式および本文(18)式より，クロススペクトルと空間自己相

その説明を示す．
まず F-K スペクトルは次のように書くことが出来る．
(付 36)

𝑢𝑙 = |𝑢𝑙 | exp(i𝜙𝑙 )関関数の関係に対する無相関項の影響は次式のように書ける．

ここに，𝜙𝑙 ：F-K スペクトル𝑢𝑙 の位相
ここで，我々は振動源に関する情報は何も持っていないので，位相𝑉[𝑥𝛾 ]
𝑆[𝑥𝛾 ] = 𝐸[𝐶𝛾 ] + 𝒪 (√)
𝐴

(付 30)

𝜙𝑙 は0 ≤ 𝜙𝑙 ≤ 2𝜋の範囲で一様にランダムな値を取り得ると仮定
するしかない．このとき次式が成り立つ．
(付 37)

𝐸[exp(i𝜙𝑙 )] = 0

また，本文(22)式のクロススペクトルの正規化において，割り算の分母の項に含まれる推定誤差について検討するため，クロススペ

この複素指数関数の𝑛乗を考えると，(exp(i𝜙𝑙 ))𝑛 = exp(i𝑛𝜙𝑙 )なの

クトルの絶対値の２乗について考えると，この平均は次のように書で，この関数の偏角は0 ≤ 𝑛𝜙𝑙 ≤ 2𝜋の範囲で一様にランダムな値ける．なお本文(22)式は２乗した値ではなく，絶対値そのものが分をとると考えられ，その期待値はやはりゼロである（次式）．

母になっているが，絶対値そのままでは取り扱いが困難なため，こ𝐸[(exp(i𝜙𝑙 ))𝑛 ] = 0 (𝑛 ≥ 1)
(付 38)
また，|𝑢𝑙 |は𝑢𝑙 に対応する振動源における加振力の大きさから，𝜙𝑙2

1/2

こではS[|𝑥𝛾 |] ≃ (𝑆 [|𝑥𝛾 | ])となることを仮定して，絶対値の２乗

は加振力の位相から，それぞれ独立に定まると考えられる．した

平均について考えるものとする．

がって，次式が成り立つと仮定できる．
𝐸[|𝑢𝑙 |𝑚 (𝑢𝑙 )𝑛 ] = 𝐸[|𝑢𝑙 |𝑚+𝑛 (exp(i𝜙𝑙 ))𝑛 ]2

𝑉 [|𝑥𝛾 | ]𝑆 [|𝑥𝛾 | ] = 𝐸 [|𝑥𝛾 | ] + 𝒪 (√
)
𝐴
2

2

(付 31)= 𝐸[|𝑢𝑙 |𝑚+𝑛 ]𝐸[(exp(i𝜙𝑙 ))𝑛 ](付 39)(付 38)式を(付 39)式に代入すると，(付 35)式が導ける．

(付 30)，(付 31)式より，クロススペクトルおよびその絶対値の２(3) クロススペクトルの絶対値の２乗の期待値

乗の算術平均と空間自己相関関数の関係を得るためには，クロスス

クロススペクトルの絶対値の２乗は次の通りに書ける．

ペクトルの絶対値の２乗の期待値，クロススペクトルの絶対値の４乗の期待値と空間自己相関関数の関係を知る必要がある（(付 4)式

𝑀

𝑀

𝑀

𝑀

2

∗ 𝑢∗ 𝑢|𝑥𝛾 | = ∑ ∑ ∑ ∑ (𝑢𝑙 𝑢𝑚
𝑛 𝑜より誤差項に含まれれる分散はこれらの期待値より計算できる）．

𝑙=1 𝑚=1 𝑛=1 𝑜=1

付 4.1 準備× exp (−i((𝒓 + 𝜸)(𝒌𝑙 − 𝒌𝑛 )(1) F-K スペクトルと振動源の関係
振動源の総数を𝑀とすると，本文(15)式より，位置𝒓の測定点とそ− 𝒓(𝒌𝑚 − 𝒌𝑜 ))))

こから𝜸ずれた点とのクロススペクトルは次式で表せる．

(付 40)

𝑀

以下，表記の簡略化のため，(付 40)式右辺の被和分関数を𝑓𝑙𝑚𝑛𝑜 と書

𝑥𝛾 = ∑|𝑢𝑙 |2 exp(−i𝒌𝑙 𝜸)
𝑙=1

く．(付 35)式より，𝑓𝑙𝑚𝑛𝑜 の期待値がゼロとならないのは，以下の２(付 32)

つのパターンのみである．
𝑙 = 𝑚かつ𝑛 = 𝑜：

∗ exp(−i(𝒌 (𝒓 + 𝜸) − 𝒌 𝒓))+ ∑ 𝑢𝑙 𝑢𝑚
𝑙𝑚

このパターンに該当する(𝑙, 𝑚, 𝑛, 𝑜)の集合を𝐴とすると，𝑙≠𝑚𝑀

またこのとき，空間自己相関関数と無相関項は以下の通りに表され∑る．

(𝑙,𝑚,𝑛,𝑜)∈𝐴𝑀

𝑓𝑙𝑚𝑛𝑜 = ∑ ∑|𝑢𝑙 |2 |𝑢𝑛 |2 exp(−i𝜸(𝒌𝑙 − 𝒌𝑚 ))𝑙=1 𝑛=1𝑀𝐶𝛾 = ∑|𝑢𝑙 |2 exp(−i𝒌𝑙 𝜸)(付 33)= |𝐶𝛾 |

𝑙=1

9-11

2

(付 41)

微動アレイ探査における波の位相速度推定手法の概説

𝑙 = 𝑛かつ𝑚 = 𝑜：𝑆[𝑥𝛾 ] = 𝐸[𝐶𝛾 ] + 𝒪 (√このパターンに該当する(𝑙, 𝑚, 𝑛, 𝑜)の集合を𝐵とすると，
𝑀(付 42)

2

𝑆 [|𝑥𝛾 | ]𝑙=1 𝑚=1

また集合𝐴と集合𝐵に重複して含まれる(𝑙, 𝑚, 𝑛, 𝑜)の組み合わせ

2= 𝐸 [|𝐶𝛾 | ] + 𝐸[|𝐶0 |2 ] − 𝐸[𝑃]として，𝑙 = 𝑚 = 𝑛 = 𝑜があり，このとき𝑓𝑙𝑚𝑛𝑜 = |𝑢𝑙 |4 である．よっ(付 50)て，クロススペクトルの絶対値の２乗の期待値は次のようになる．
2(付 49)

𝑀

𝑓𝑙𝑚𝑛𝑜 = ∑ ∑ |𝑢𝑙 |2 |𝑢𝑚 |2 = |𝐶0 |2∑(𝑙,𝑚,𝑛,𝑜)∈𝐵

𝑉[𝐶𝛾 ] + 𝐸[|𝐶0 |2 ])𝐴

2𝐸 [|𝑥𝛾 | ] = 𝐸 [|𝐶𝛾 | ] + 𝐸[|𝐶0 |2 ] − 𝐸[𝑃](付 43)

𝑉[𝑥𝛾 ] = 𝑉[𝐶𝛾 ] + 𝐸[|𝐶0 |2 ] − 𝐸[𝑃]√+𝒪

(付 43)式および本文(18)式（𝐸[𝑥𝛾 ] = 𝐸[𝐶𝛾 ]）より次式が成り立つ．

2 2

2

𝑉 [|𝐶0 |2 + |𝐶𝛾 | ] + 𝐸 [(|𝐶0 |2 + |𝐶𝛾 | ) ]𝐴()(付 44)

付録 5. 比較的近年の位相速度推定手法: 長らの方法 24)

ここで，|𝐶0 |2 = ∑|𝑢𝑙 |2 |𝑢𝑚 |2 = 𝑃 + ∑𝑙≠𝑚|𝑢𝑙 |2 |𝑢𝑚 |2 より，𝑃 ≪|𝐶0 |2 とすれば，次式も成り立つ．

長らの提案したこれらの方法は，SPAC 法の拡張と捉えることが𝒪(𝑉[𝑥𝛾 ]) = 𝒪(𝑉[𝐶𝛾 ] + 𝐸[|𝐶0 |2 ])(付 45)

できる．原理は SPAC 法に類似しているが，一般的な SPAC 法に比(4) クロススペクトルの絶対値の４乗の期待値

べると，センサアレイの中心に測定点を必要とせず低周波数側に高クロススペクトルの絶対値の４乗は次の通りに書ける．

い感度を持つ手法などSPAC 法の欠点を補うことのできる手法を含
むため，近年適用が進みつつある 25) ．

4∗ 𝑢∗ 𝑢 𝑢∗ 𝑢 𝑢 𝑢∗|𝑥𝛾 | = ∑ (𝑢𝑙 𝑢𝑚
𝑛 𝑜 𝑝 𝑞 𝑠 𝑣

原理をかいつまんで説明すると以下のようになる．長らの方法では，まず複素正弦関数で重み付けした微動データの方位積分を以下

× exp (−i ((𝜸 + 𝒓)(𝒌𝑙 − 𝒌𝑛 − 𝒌𝑝 + 𝒌𝑠 )のように定義する．
2𝜋− 𝒓(𝒌𝑚 − 𝒌𝑜 − 𝒌𝑝 + 𝒌𝑣 ))))𝑧𝑚 (𝑡, 𝒓) ≔ ∫ 𝑢(𝑡, 𝒓) exp(−i𝑚𝜃) 𝑑𝜃

(付 51)

0(付 46)

単一の振動源の寄与が卓越し，かつ単一速度の仮定が成立するもの

2 乗の場合と同様に，期待値がゼロになる(𝑙, 𝑚, 𝑛, 𝑜, 𝑝, 𝑞, 𝑠, 𝑣)の組として，本文(3)式のフーリエ変換と Hansen の積分表示（(付 1)式）

み合わせと，そうならない組み合わせを考える．期待値がゼロになを組み合わせると，𝑧𝑚 (𝑡, 𝒓)を時間軸に関してフーリエ変換した結らない組み合わせのパターンは全部で 24 通りあり，被和分関数の果は，以下のように書ける．

4

和が|𝐶𝛾 | となるのが４パターン，|𝐶0 |4 となるのが４パターン，|𝐶𝛾 | |𝐶0 |2 となるのが１６パターンで，全てのパターンで𝑙 = 𝑚 =𝜋𝑧𝑚 (𝜔, 𝒓) = 2𝜋 exp (i𝑚 ) 𝐽𝑚 (𝑘̃‖𝜸‖)2𝑛 = 𝑜 = 𝑝 = 𝑞 = 𝑠 = 𝑣となるパターンが重複している．以上まとこれを用いて，長らは次のような関係を導いた．

2

めると，クロススペクトルの絶対値の４乗の期待値は次のようになる．

2
𝐽0 (𝑘̃‖𝜸‖)
𝑧0 (𝜔, 𝒓)𝑧0∗ (𝜔, 𝒓)=()
𝑧1 (𝜔, 𝒓)𝑧1∗ (𝜔, 𝒓)
𝐽1 (𝑘̃‖𝜸‖)(付 53)

𝜌𝐻0 =

2𝑧0 (𝜔, 𝒓)𝑧0∗ (𝜔, 𝒓)= (𝐽0 (𝑘̃‖𝜸‖))∗𝑧1 (𝜔, 0)𝑧1 (𝜔, 0)(付 54)

𝜌𝐻1 =

2𝑧1 (𝜔, 𝒓)𝑧1∗ (𝜔, 𝒓)= (𝐽1 (𝑘̃‖𝜸‖))∗
𝑧0 (𝜔, 0)𝑧0 (𝜔, 0)(付 55)𝐽0 (𝑘̃‖𝜸‖)
𝑧0 (𝜔, 0)𝑧0∗ (𝜔, 𝒓)=∗2
𝑧1 (𝜔, 𝒓)𝑧1 (𝜔, 𝒓)(𝐽1 (𝑘̃‖𝜸‖))(付 56)𝜌𝐶𝐶𝐴 =
4

4

𝐸 [|𝑥𝛾 | ] = 4𝐸 [|𝐶𝛾 | ] + 4[|𝐶0 |4 ](付 47)
2

2+ 16𝐸 [|𝐶𝛾 | |𝐶0 | ] − 23𝑄
8
ここに，𝑄 ≔ ∑𝑀
𝑙=1 𝐸[|𝑢𝑙 | ]ここで，𝑃 2 = 𝑄 + ∑𝑙≠𝑚|𝑢𝑙 |4 |𝑢𝑚 |4 より，𝑄 ≪ 𝑃 2 ≪ |𝐶0 |4と仮定

𝜌𝑉 =すると，次式も成り立つ．
2

𝒪 (𝑉 [|𝑥𝛾 | ])2

2 2

(付 52)

以上の式において，中辺は微動の実測データより（近似的に）計算できるので，その値と整合するような𝑘̃ を右辺の式より求めれば，(付 48)= 𝒪 (𝑉 [|𝐶0 |2 + |𝐶𝛾 | ] + 𝐸 [(|𝐶0 |2 + |𝐶𝛾 | ) ])微動波動場の卓越波数が得られる．特に(付 53)式を用いる CCA 法
は，センサアレイ中心に測定点を必要としないため，SPAC 法の適用が難しいようなサイトでも CCA 法なら適用が可能である場合が付 4.2 算術平均の評価
以上計算結果を(付 30)，(付 31)式に代入して整理すると，クロスあり，長らの方法の中でも代表的な方法とされている．

ススペクトルおよびその絶対値の２乗の算術平均は以下の通りに表される．

9-12