LoginSignup
0
0

More than 5 years have passed since last update.

@kamenoseiji(Seiji Kameno)

干渉計観測の振幅を解く

Last updated at Posted at 2017-01-16

振幅を解く

この章では、干渉計観測で得られるビジビリティの振幅を使ってアンテナゲインの振幅だけについて解く方法を一応記します。ですが実際には使わずに、振幅と位相を合わせて複素ゲインを解く方をお奨めします。なぜなら、振幅は正規分布ではなくRice分布に従うため、signal-to-noise ratioが小さいときは振幅が大きめにバイアスされる傾向があるからです（一方、複素ゲインは正規分布する）。$SNR > 5$くらいなら正規分布に近くなるので、この章の方法を摘要してもいいでしょう。
この章の意義は、アンテナゲインを複素数として解く方法を理解しやすくするための前振りです。

アンテナゲインの振幅$|G|$と基線ベースのビジビリティ振幅$|V|$との関係は以下のようになります。

|\hat{V}|_k = |G|_j |G|_j |V|_k\tag{5.1}

右辺の$|V|_k$は真のビジビリティ振幅、左辺の$|\hat{V}|_k$は観測されたビジビリティ振幅です。真のビジビリティ振幅は較正して定めるまで分からないですから、いっそ$|V|_k$で割った$|\hat{V}|_k / |\hat{V}|_k$を左辺として扱っても一般性は失いません。この場合、$|G|$は適当にスケーリングされた値となり、相対的な振幅としての意味しか持ちません。この場合観測方程式(5.1)は、以下の(5.2)のようになります。

|\hat{V}|_k = |G|_j |G|_j \tag{5.2}

対数による簡易解

式(5.2)の対数をとると

\log |\hat{V}|_k = \log |G|_j + \log |G|_j \tag{5.3}

となり、線形化できます。すなわち

\left(
\begin{array}{c}
\log |\hat{V}_0|  \\
\log |\hat{V}_1|  \\
\log |\hat{V}_2|  \\
\log |\hat{V}_3|  \\
\log |\hat{V}_4|   \\
\log |\hat{V}_5|   \\
\vdots
\end{array}
\right) = \left( \begin{array}{ccccc}
1 & 1 & 0 & 0 & \cdots \\
1 & 0 & 1 & 0 & \cdots \\
0 & 1 & 1 & 0 & \cdots \\
1 & 0 & 0 & 1 & \cdots \\
0 & 1 & 0 & 1 & \cdots \\
0 & 0 & 1 & 1 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
\log |{G}_0|  \\
\log |{G}_1|  \\
\log |{G}_2|  \\
\log |{G}_3|  \\
\vdots
\end{array}
\right) \tag{5.4}

です。真ん中の行列を$P$とおき、左辺の$\log |V|$のベクトルを$\boldsymbol{V}$, 右辺の$\log |G|$のベクトルを$\boldsymbol{G}$とおくと、$\boldsymbol{V} = P \boldsymbol{G}$ですから、$\boldsymbol{G}$について解くと$\boldsymbol{G} = (P^T P)^{-1} P^T \boldsymbol{V}$ です。

計算の省力化

$P^T P$を手計算します。

P^TP = \left( \begin{array}{ccccc}
 N_a - 1 & 1 & 1 & 1 & \cdots \\
1 & N_a - 1 &  1 & 1 & \cdots \\
1 & 1 & N_a-1 & 1 & \cdots \\
1 & 1 & 1 & N_a-1 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{5.5}

となり、その逆行列は

(P^TP)^{-1} = \frac{1}{2(N_a-1)(N_a - 2)} \left( \begin{array}{ccccc}
2N_a -3 & -1 & -1 & -1 & \cdots \\
-1 & 2N_a -3 &  -1 & -1 & \cdots \\
-1 & -1 & 2N_a -3 & -1 & \cdots \\
-1 & -1 & -1 & 2N_a -3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{5.6}

となります。要するに対角成分が$\frac{2N_a -3}{2(N_a-1)(N_a - 2)}$, 非対角成分が$-\frac{1}{2(N_a-1)(N_a - 2)}$です。この行列もPythonなら1行でコードできます。

import numpy as np
PTP_inv = ((2.0* antNum - 2.0)* np.diag(np.ones(antNum)) - 1.0) / (2.0* (antNum - 1.0)* (antNum - 2.0))

$P^T \boldsymbol{V}$の計算も、遅延について解いたときのコードとほぼ同様にしてできます。

ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
PTV = np.zeros(antNum)
for ant_index in range(antNum):
  index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
  PTV[ant_index] += np.sum(log_bl[index0]) + np.sum(log_bl[index1])
#

非線形解

対数をとらずに式(5.2)を解く場合は非線形方程式になります。アンテナゲイン振幅初期値を$G^0$とすると、基線ベースの残差は$r_k = \hat{V}_k - G^0_j G_i^0$となります。残差のノルム$\displaystyle \sum_k r_kr^*_k$ を最小にするために$G^0$に与える修正ベクトル$\boldsymbol{c}$は、

\left(
\begin{array}{c}
r_0  \\
r_1  \\
r_2  \\
r_3  \\
r_4  \\
r_5  \\
\vdots
\end{array}
\right) = \left( \begin{array}{ccccc}
 G_1 &  G_0 & 0 & 0 & \cdots \\
 G_2 &  0 & G_0 & 0 & \cdots \\
 0 & G_2 & G_1 & 0 & \cdots \\
 G_3 & 0 & 0 & G_0 & \cdots \\
 0 & G_3 & 0 & G_1 & \cdots \\
 0 & 0 & G_3 & G_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
c_0  \\
c_1  \\
c_2  \\
c_3  \\
\vdots
\end{array}
\right) \tag{5.7}

という方程式を満たします。真ん中の行列を$P$の要素が解くべき$G$を含んでいるので、非線形になるわけです。式(5.7)を$c$について解くと$\boldsymbol{c} = (P^T P)^{-1} P^T \boldsymbol{r}$です。こうして得られた$\boldsymbol{c}$を用いて$\boldsymbol{G}$を改良します。

\boldsymbol{G} \leftarrow \boldsymbol{G} + \boldsymbol{c}

これを$\boldsymbol{G}$が収束するまで繰返せば、解が得られます。収束は$\boldsymbol{c}$のノルムが十分小さくなることで判定できます。

計算の省力化

例によって$P^T P$を手計算します。

P^TP = \left( \begin{array}{ccccc}
\sum_k G^2_k - G^2_0 & G_0 G_1 & G_0 G_2 & G_0 G_3 & \cdots \\
G_0 G_1 & \sum_k G^2_k - G^2_1 &  G_1 G_2 & G_1 G_3 & \cdots \\
G_0 G_2 & G_1 G_2 & \sum_k G^2_k - G^2_2 & G_2 G_3 & \cdots \\
G_0 G_3 & G_1 G_3 & G_2 G_3 & \sum_k G^2_k - G^2_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{5.8}

これをPythonで実装します。

def PTPmatrix(Gain):  # Gain is a vector of antenna-based gain amplitude (real)
  antNum = len(Gain); normG = Gain.dot(Gain)
  PTP = np.zeros([antNum, antNum]) + Gain
  for ant_index in range(antNum): 
    PTP[ant_index,:] *= Gain[ant_index]
    PTP[ant_index, ant_index] = normG - Gain[ant_index]**2
  #
  return PTP
#

逆行列$(P^T P)^{-1}$は簡単な形には書けません。実対称行列ですので、Cholesky分解を使って方程式を解くことにします。あらかじめ$\boldsymbol{y} = P^{T}\boldsymbol{r}$を計算しておきます。

y = np.zeros(antNum)
for ant_index in range(antNum):
  index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
  y[ant_index] += antGain[ant1[index0]].dot(resid[index0])
  y[ant_index] += antGain[ant0[index1]].dot(resid[index1])
#

次にCholesky分解を使って$P^T P \boldsymbol{c} = \boldsymbol{y}$を$\boldsymbol{c}$について解きます。

L = np.linalg.cholesky(PTP)
t = np.linalg.solve(L, y)
correction = np.linalg.solve(L.T, t)

このようにして得られた$\boldsymbol{c}$ (correction) を$\boldsymbol{G}$ (antGain) に加えて改良し、収束するまで反復します。

まとめのコード

以上をまとめた、基線ベースのビジビリティ振幅からアンテナベースのゲイン振幅を得るPythonコードを下記に示します。
関数logamp_solve() は対数をとって線形化した解法で、引数にビジビリティ振幅のベクトルをとり、アンテナゲイン振幅のベクトルを返します。
関数clamp_solve() は非線形解法で、引数と返り値はlogamp_solve()と同じ仕様です。内部でアンテナゲイン振幅の初期値を得るためにlogamp_solve(bl_amp)を呼び出しています。また、$(P^T P)^{-1}$の計算をするための関数 PTPmatrix() を呼び出します。

def logamp_solve(bl_amp):
    blnum  =  len(bl_amp); log_bl =  np.log(bl_amp)
    antNum =  Bl2Ant(blnum)[0]
    ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    PTP_inv = ((2.0* antNum - 2.0)* np.diag(np.ones(antNum)) - 1.0) / (2.0* (antNum - 1.0)* (antNum - 2.0))
    PTV = np.zeros(antNum)
    for ant_index in range(antNum):
        index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
        PTV[ant_index] += np.sum(log_bl[index0]) + np.sum(log_bl[index1])
    #
    return np.exp(PTP_inv.dot(PTV))
#
def PTPmatrix(Gain):  # Gain is a vector of antenna-based gain amplitude (real)
    antNum = len(Gain); normG = Gain.dot(Gain)
    PTP = np.zeros([antNum, antNum]) + Gain
    for ant_index in range(antNum):
        PTP[ant_index,:] *= Gain[ant_index]
        PTP[ant_index, ant_index] = normG - Gain[ant_index]**2
    #
    return PTP
#
def clamp_solve(bl_amp, niter=2):
    blnum  =  len(bl_amp)
    antGain = logamp_solve(bl_amp); antNum = len(antGain)
    ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    for iter_index in range(niter):
        resid = bl_amp - antGain[ant0]* antGain[ant1]
        y = np.zeros(antNum)
        for ant_index in range(antNum):
            index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
            y[ant_index] += antGain[ant1[index0]].dot(resid[index0])
            y[ant_index] += antGain[ant0[index1]].dot(resid[index1])
        #
        L = np.linalg.cholesky(PTPmatrix(antGain))
        t = np.linalg.solve(L, y)
        correction = np.linalg.solve(L.T, t)
        antGain += correction; antGain = abs(antGain)
    #
    return antGain
#

以上です。実際には振幅だけを解くより、振幅と位相を同時に複素数として解く方がいいでしょう。
もくじにもどる

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0