こんにちは．18のhiraです．
僕は水中ロボットを製作するアクア研で活動しています．そんな水中探査では観測により地形図の作製をすること(マッピング)が目標だったりするのですが，水中では電波が急速に減衰するためレーザー測量を使うことはできません．そのため超音波ソナーを使います．しかしこの音波というものは指向性に乏しく，あらゆる方向の反射波が重なってしまいます．そこで，今回はその反射波を数理的に解析し高度なマッピングを行うソナーの信号処理について書いていきます．
“フェイズドアレイ”って?
水中観測に欠かせないソナーですが，ここで多く使われるのがフェイズドアレイという手法です．戦艦や戦闘機関連の軍用のイメージが強いと思いますが，軍用迎撃レーダーや気象観測レーダー，電波望遠鏡，超音波エコー装置やソナーなど波の性質を使うあらゆる観測装置に使われている手法です．
専用の素子を使って対象に当たり跳ね返ってきた電波や超音波が観測される時間から距離を求めるという1つの素子の役割は変わらりません．しかしフェイズド(位相)アレイ(配列)という言葉通り，このシステムは数個~数百個の素子を並べ，素子同士の観測データの差異を解析することでより大規模な観測を行うことが目的です．(反射波と言っていますが電波望遠鏡など受動観測では入力波と読み替えてください)
素子が1つしかない観測の場合，波(電波や超音波)を発振してから t 時間後に反射波が到来したということが観測できます．そのため波の速度を v とすれば「 vt/2 の距離に障害物がある」ということが分かります．しかしその障害物が正面にあるのか，もっと左にあるのか，右にあるのかという「どの角度から来たのか」という情報を得ることはできません．現実的には市販の超音波センサだと正面から離れるにつれ波のエネルギーが弱まるので，観測可能な強い反射波は正面付近の障害物ぐらいからしか得られません．しかし正確な角度は分かりません．
これに対しフェイズドアレイは複数の素子を使い，ある角度から到来した反射波の素子ごとの到来時間のズレを計測(時間計測なら高精度です)して，波の到来角度に変換します．
ここでは一旦到来波とは反対にある角度に指向的に強い波を発振することを考えましょう．「特定の方向にだけ聞こえるスピーカー」を知っていますか?これにもフェイズドアレイが用いられており，特定の角度方向では波同士が強め合うため強い合成波が聞こえ，他の角度では合成波が上手くまとまらず相殺されてしまうため(場合によっては不本意に強まってしまうこともあるが)聞こえません．
波上の図は6つの素子を横1列に並べたものです．左の紫の素子から出た波は大きく広がっています．これは他の波よりも先に発振したためです．対して右の黄緑の素子の波はほとんど広がっていません．これは左の素子たちに比べ大きな遅延時間の後に発振されたからです．このように各素子に特定の遅延時間を与え発振タイミングをずらしています．
各波の位相が0の時の波(先頭)の合成波の波面が赤線で示されています．なんと，きれいに直線状になっています．これはある入射角度で波が強め合っているということを意味し，入射角度に指向的な波が作られたと言えます．このように1つ1つの素子はフツーの素子でも，遅延時間を制御することで方向まで制御した波を作るのがフェイズドアレイです．
具体的な遅延時間を求めてみましょう． M 個の素子でアレイが構成され，素子同士の間の距離を d とし，入射角度を \theta にしたいとすると，求める m 番目( m=1,2,...,M )の素子の1番目の素子に対する遅延時間を \tau_m(\theta) とすると，\begin{aligned}\tau_m(\theta) = (m-1)\frac{d\sin\theta}{v}\end{aligned}となります．
反対に到来波の角度を推定したい場合には，各素子ごとに到来波を観測し，観測したい角度に対応する遅延時間のズレを持つ波が観測されるかを求めます．あたかも好きな方向に向けたアンテナ(受信機)で観測したかのように波を再構成できるのです．
フェイズドアレイによる波の信号解析
フェイズドアレイにより，各素子の観測値を持ち寄ることで任意の方向に向いた受信機を計算により再現できることが分かりました．そしてこの計算手法が最も大事な要素です．計算方法を工夫することで，
・ 一度の計算で複数の方向をまとめて計算できる．
・ より厳密な角度の情報を得られる
・ アレイに使う素子の数を減らすことができる
というメリットがあります．ここでは到来方向推定に使われる「MUSIC法」と，それをソナーによるマッピングに活かす方法を紹介します．
MUSIC法
MUSIC法は MUltiple SIgnal Classification の略で波のより正確な到来角度推定ができる(角度分解能が高いという)手法です．
M 個の素子センサアレイで時間 \Delta t ごとにセンサの値を読み取る観測を行うと考えます．この時センサで観測される値を受信信号ベクトル {\bf x} とします．この観測値 {\bf x} を深堀りしてみましょう．この M 個の素子による受信信号ベクトルが方位 \theta からの到来波を含むとします．超音波など観測する波の波長を f とすると， m 番目の素子が1番目の素子に対して観測する波の位相差 a_m(\theta) は，
\begin{aligned} a_m(\theta) = \exp(-j2\pi f\tau_m(\theta)) \end{aligned}だけずれていると言えます．ここで角度 \theta_1 , \theta_2 の2方向から s_1 , s_2 の波形を持った波が観測されるとします．この s_1 , s_2 は適当な振幅や位相の正弦波と考えてください．時刻 k の時の観測値行列は M\times1 行列 {\bf x} (k)で表され，受信値は
\begin{aligned}{\bf x}(k) = {\bf a}(\theta_1)s_1(k)+{\bf a}(\theta_2)s_2(k)+{\bf n}(k)\end{aligned}と表せます． {\bf n}(k) は白色ノイズベクトルです．このようにある到来波がその角度による遅延分だけ素子ごとに位相のずれた振幅が検出されます．その波ごとの和が振幅の和が {\bf x}(k) となります．
そして近隣地がどれだけ似ているかを表す「空間自己相関」を {\bf x} に対して求めます．転置して虚数部の符号を反転させるエルミート転置を ...^H ，期待値を E[...] とすると，時刻 k 周辺のの {\bf x} の空間自己相関行列 {\bf R_x} は，
\begin{aligned} {\bf R_x} = E[{\bf x}^H(k){\bf x}(k)] \end{aligned}となります．定常な環境であれば期待値をとるには現在の時刻を K として，
\begin{aligned} E[{\bf x}^H(k){\bf x}(k)] = \frac{1}{K}\sum^K_{k=1}{\bf x}^H(k){\bf x}(k) \end{aligned}ともできます．波の反射源の数を L とし， {\bf R_x} の固有値のうち小さい M-L 個の固有値に対応する固有ベクトルを求めます．これは雑音固有ベクトルと呼ばれ，これを並べた M\times (M-L) 行列を {\bf E_N} とします．到来波があるか求めたい N 個の方位に対して次の式から角度スペクトラム P_{MU}(\theta_n), n=1,2,...,N を求めます．
\begin{aligned}P_{MU}(\theta)=\frac{{\bf a}^H(\theta_n){\bf a}(\theta_n)}{{\bf a}^H(\theta_n){\bf E_N}(\theta_n){\bf E_N}^H(\theta_n){\bf a}(\theta_n)}\end{aligned}スペクトラムとあるように， \theta と P_{MU}(\theta) をプロットすると \theta_n において強いスペクトラムとなります． \theta_n の到来波があればその角度のときの分母が0になります． P_{MU}(\theta) の値が無限大となるからです．
ここで， -60^\circ , -45^\circ , -15^\circ , 0^\circ , 10^\circ , 35^\circ , 50^\circ ,の7方位から来る到来波を検出できるか試してみましょう．素子数 M=16 ，サンプリング回数 K=10 としました．わかりやすさのため，角度スペクトラムを
\begin{aligned} \log_{10}\frac{P_{MU}(\theta)}{\max\bigl(P_{MU}(\theta)\bigr)} \end{aligned}とすると，推定結果は次のグラフのようになりました．
このように入射波のある角度に対しては強いスペクトラムが検出されます．閾値を設け，それより強い値が観測されれば波があると判断できます．
MUSIC法の問題点
さて，波の到来角度を鋭く推定できるMUSIC法ですがこれをソナーに使おうとすると問題があります．なぜならMUSIC法は時間をかけるほど正確・多数の波の推定ができるのに対し，ソナーではなるべく少ない時間間隔で区切る必要があるからです．言葉ではよく分かりませんね．式や図から見ていきましょう．
角度スペクトラムの計算において， {\bf R_x} は観測データ列 {\bf x} を用いました．この {\bf x} の1列は M 個のセンサがある時刻 k に計測した波の振幅の観測データです．先ほど {\bf R_x} は「今までの平均」としましたが時刻 k における {\bf R_x} の詳しい式としては，
\begin{aligned} {\bf R_x(k)} = \frac{1}{P}\sum^{P-1}_{p=0}{\bf x}^H(k-p){\bf x}(k-p) \end{aligned}となります．ここで新たに出てきた p は時刻 k から何個遡った観測値まで使うかを意味します．つまり {\bf R_x}(k) は時刻 k の直前の何個かの観測値 {\bf x}(k-P+1) ~ {\bf x}(k) の P 個の観測値 {\bf x} を使い作製します．
MUSIC法はこの P 個のデータの中に波が来たか，すなわち時間 P のどこかで波が来ているかを推定します．この P が大きいほど得られるスペクトルはより鋭く正確になります．
さらに波が一度に多数到来する場合，MUSIC法では一度のスペクトラム計算で見分けられる個数に制限があります．その個数はアレイを構成する素子の個数までで，さらにその個数すなわち M 個以上の P 個の {\bf x} を与えなけらばなりません．例えば素子の数が16個で，一度に12個の波が別々の方角から来るとします．この場合12個波のスペクトルが現れるには最低でも P=12 のデータ列 {\bf x} が必要です．これは {\bf R_x} のランクに依存し，アレイ素子の数だけの波を推定しようとすると {\bf R_x} をフルランクにする P \geq M の観測値が必要です．仮に P=4 程度でスペクトルを無理やり計算しようとすると，
とてもじゃないけれど波がどこにあったか分かりませんよね…ちなみに P=16 , P=50 でやってみると，
このように P を増やすほどきれいなスペクトラムが得られます．
ならば P を大きくすればいいのでは？となりますが，これはソナーには致命的です．なぜなら「 P 時間のどこかで波が来たのは分かるが，いつ来たかは分からない」からです．ソナーは反響音が帰ってくるまでの時間を計測することで音速x時間から距離を求めます．そのため「いつ帰って来たか」の時刻をピンポイントで求めることが不可欠です．ゆえにいつ帰ってきたかを正確に知りたいならば P を小さくする必要があります．
お察しの通り，これはMUSIC法の P を大きくしたいという要望と全く正反対です． P を大きくすればより正確な到来角度が分かり多くの波にも対応できるが，距離推定の誤差が大きくなる．逆に P を小さくすれば正確な到来角度は分かりにくく多くの波は判別できないけど，距離を正確に求められるという反対の関係です．
どうにかしてこのトレードオフを一挙解決する「 P を小さくして到来時刻を正確に求め，かつ正確な到来角度と多くの波も識別できる」魔法のような方法はないでしょうか?これを可能にするのが三菱重工が開発した「適応アレイ信号処理」によるMUSIC法の改良版[1]です．
適応アレイ信号処理
先のMUSIC法の問題はどうやって時間窓( P の値)を広げずに相関点数を稼いで {\bf R_x}(k) をフルランクにするかといえます．それを解決するのが空間平均法です．
空間平均法では並んだアレイ素子から何個か選び小さなグループ(サブアレイ)と見なします．このサブアレイでそれぞれ {\bf R_x}(k) を求め，最後各サブアレイで求めた {\bf R_x}(k) を合算します．これにより時間は同じでも相関点数の増えた {\bf R_x}(k) を得ることができます．
具体的な方法は，まず M 個の素子からなるアレイを下図のようにグループに分割し， M_S 個の素子からなるサブアレイを N_S 個作ります．
が求まります．元の標本相関行列の相関点数は P でしたが，新たに求めたこの {\bf R_x}(k) の標本相関点数は N_SP です．これにより相関点数が N_S 倍になりました．
さらに倍増させましょう．次はForward-Backward空間平均法を用います．これはアレイのデータをそのまま並べた行列と，アレイを逆順に並べたデータの複素共役行列の2つを使うことでより {\bf R_x}(k) の標本相関点数を増加させます．[2]
現在データ {\bf x} は
\begin{aligned} \begin{pmatrix} x_1(k-P+1) & x_1(k-P+2) & \cdots & x_1(k) \\ x_2(k-P+1) & x_2(k-P+2) & \cdots & x_2(k) \\ \vdots & \vdots & \ddots & \vdots\\ x_M(k-P+1) & x_M(k-P+2) & \cdots & x_M(k) \end{pmatrix} \end{aligned}と並んでいます．これが正順のForward行列です．ある要素ないし行列 a の複素共役を \bar a と表すと {\bf x} のBackward行列 {\bf x^b} は，
\begin{aligned} \begin{pmatrix} \overline{x_M(k-P+1)} & \overline{x_M(k-P+2)} & \cdots & \overline{x_M(k)} \\ \vdots & \vdots & \ddots & \vdots \\ \overline{x_2(k-P+1)} & \overline{x_2(k-P+2)} & \cdots & \overline{x_2(k)} \\ \overline{x_1(k-P+1)} & \overline{x_1(k-P+2)} & \cdots & \overline{x_1(k)} \end{pmatrix} \end{aligned}となります．具体的には
\begin{aligned}{\bf J}=\begin{pmatrix} 0 & \cdots & 0 & 1\\ 0 & \cdots & 1 & 0\\ \vdots & 1 & \ddots & \vdots \\ 1 & 0 & \cdots & 0 \end{pmatrix} \end{aligned}の {\bf J} を用いて
\begin{aligned} {\bf x^b} = J\overline {\bf x}\end{aligned}となります．これよりForwardと同じくBackwardの自己相関行列 {\bf R_x^b} を以下より求めます．
\begin{aligned} {\bf R_x^b} = E[{\bf x^b}^H(k){\bf x^b}(k)] \end{aligned}最後にForwardとBackwardの自己相関行列の和をとることで，Forward-Backward空間平均法による自己相関行列 {\bf R_x^{bf}} は，
\begin{aligned} {\bf R_x^{bf}} = \frac{{\bf R_x}+{\bf R_x^{b}}}{2}\end{aligned}となります．これをMUSIC法の角度スペクトラム式に利用することで相関点数がForwardだけの時( {\bf R_x} のとき)の2倍になります．先ほどのサブアレイによる空間平均法と合わせると P 時間の計測により 2N_SP 個の相関点を得ることができます．元の P 個の 2N_S 倍です．
またこれを用いて設定すべきサブアレイの個数を求めることができます． {\bf R_x^{bf}} をフルランクになるようにするには，
\begin{aligned} 2N_SP \geq M_S \end{aligned}を満足する必要があります．この時，各変数の関係は
\begin{aligned} N_S = M - M_S + 1 \end{aligned}であるといえるので，上記2式を満たす M_S の値は
\begin{aligned}M_S = \frac{M+1}{1+0.5P^{-1}}\end{aligned}となります．
これが本当に通用するのかさっきの例で試してみましょう．素子の数が16個で，一度に12個の波が別々の方角から来るとします．ここで P=4 とすると M_S = 15 ， N_S = 2 となります．つまり15個の素子をもつアレイが2つあります．この時相関点数は 2N_SP = 16 となります．波の個数12個よりは大きいようです．
計算結果が以下の図です．
おおっ！なんとたった P=4 で12個の波全てが検出されています！従来のMUSIC法では P=16 とか数倍も時間窓を広げていたので大きな時間分解能の向上です!これならばっちりソナーにも使えそうです．
一つのドでかいアレイもいいけれど，小さなアレイをたくさん使うことでより多くの情報を抜き出せるのは意外でした．また無秩序の代名詞のようなノイズでも「平均0，相関がない」という自然の秩序を持っていて，うまくあしらえば消すことができるというのも面白い発見でした．
そういえば今日(昨日)はクリスマスで，もういくつ寝るとお正月です．課題は24時間年末年始無休で猛威を振るってきそうですが，伊達巻ぐらいはつくってお正月を迎えたいです．
参考資料
・[1]藤島泰郎, 田崎智「適応アレイ信号処理の導入によるアクティブソーナーの高精度化」三菱重工技報 Vol.55 No.2 (2018), <http://www.mhi.co.jp/technology/review/pdf/552/552020.pdf>
・[2]伊東直樹「最大の電界強度を有する到来波の到来方向推定に関する研究」横浜国立大学, 平成18年2月28日, <http://www.arailab.dnj.ynu.ac.jp/thesis/h17/main-ito-naoki.pdf>