干渉計観測の位相を解く

#位相を解く
位相$\phi$は遅延と同様に

\hat{\phi}_k = \phi_j - \phi_i \tag{4.1}

という関係で表されますが、2点ほど注意が必要です。

1. 位相には$2n\pi$の不定性があります。基線ベースの位相が例えば$2\pi/3$なのか$-4\pi/3$なのか、自明ではありません。だから、$\hat{\phi}_k = \phi_j - \phi_i \pm 2n\pi$となってしまい、遅延と同様に解くのは問題があります。
2. この関係式が使えるのは、closure phaseが0になる場合だけです。天体が点源であればclosure phaseは0になりますが、複雑な構造をもつ広がった天体ではこの関係式は成り立ちません。

2.に関しては位相較正天体として点源を選ぶことで対処します。1.については、基線毎のの位相が$\pm \pi/2$に収まっているようなケースではそのまま遅延の解と同様に扱ってもよいです。そうでない一般には複素解で解くのがよいでしょう。

##複素数として解く
アンテナゲインの位相項は

e^{i\hat{\phi}_k} = e^{i\phi_j} e^{-i\phi_i} \tag{4.2}

のように表されます。これは式(4.1)を満たし、かつ$2n\pi$の不定性がありません。ただし非線形な方程式ですので逐次近似による反復が必要です。以下では表記を簡略にするために$E_j = e^{i\phi_j}$とおきます。
初期値$E^0_i$が最適解に十分近いとします。基線ベースの残差は$r_k = E_k - E^0_j E_i^{0*}$です。残差のノルム$\sum_k r_k r^*_k$ を最小にするために$\phi^0_j$に与える修正ベクトル$\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}{cccc}
 \frac{\partial {E_1 E^*_0}}{\partial \phi_1} &  \frac{\partial {E_1 E^*_0}}{\partial \phi_2} &  \frac{\partial {E_1 E^*_0}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_2 E^*_0}}{\partial \phi_1} &  \frac{\partial {E_2 E^*_0}}{\partial \phi_2} &  \frac{\partial {E_2 E^*_0}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_2 E^*_1}}{\partial \phi_1} &  \frac{\partial {E_2 E^*_1}}{\partial \phi_2} &  \frac{\partial {E_2 E^*_1}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_3 E^*_0}}{\partial \phi_1} &  \frac{\partial {E_3 E^*_0}}{\partial \phi_2} &  \frac{\partial {E_3 E^*_0}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_3 E^*_1}}{\partial \phi_1} &  \frac{\partial {E_3 E^*_1}}{\partial \phi_2} &  \frac{\partial {E_3 E^*_1}}{\partial \phi_3} & \cdots \\
 \frac{\partial {E_3 E^*_2}}{\partial \phi_1} &  \frac{\partial {E_3 E^*_2}}{\partial \phi_2} &  \frac{\partial {E_3 E^*_2}}{\partial \phi_3} & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
c_1  \\
c_2  \\
c_3  \\
\vdots
\end{array}
\right) \tag{4.3}

という方程式で示されます。なお、refantの位相は0と定めて一般性を失いませんので、$c_0$は未知数ではなく常に0です。$\frac{\partial}{\partial \phi_0}$も0です。
真ん中の行列$P$の要素を書き出してみると、

P = \left( \begin{array}{cccc}
i{E_1 E^*_0} &  0 & 0 & \cdots \\
0 &  i{E_2 E^*_0} &  0 & \cdots \\
-i{E_2 E^*_0} &  i{E_2 E^*_0} &  0\\
0 &  0 & i{E_3 E^*_0}  & \cdots \\
-i{E_3 E^*_1} &  0 & i{E_3 E^*_1}  & \cdots \\
0 &  -i{E_3 E^*_2} & i{E_3 E^*_2}  & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right) \tag{4.4}

のようになります。なおrefantのゲインについては$E_0 = E^*_0 = 1$です。式(4.3)は

\boldsymbol{r} = P \boldsymbol{c}

という連立方程式で与えられますので、これを$\boldsymbol{c}$について解くと

P^{*T}\boldsymbol{r} = P^{*T} P \boldsymbol{c} \\
\boldsymbol{c} = (P^{*T} P)^{-1} P^{*T}\boldsymbol{r}

となります。このようにして得られた$\boldsymbol{c}$を位相の初期値$\boldsymbol{\phi}^0$に加え、$\boldsymbol{\phi}^1 \leftarrow \boldsymbol{\phi}^0 + \boldsymbol{c}$として位相を改良します。これを反復すればいつかは収束するでしょう。収束判定は修正ベクトル$\boldsymbol{c}$のノルムで判定できます。

計算の省力化

$P^{*T} P$をあらかじめ手計算しておきます。

P^{*T} P = \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)

となり、遅延のときと同じ行列になります。従って逆行列も

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

ですので、遅延の解を得るときと同じ行列が使えます。
$P^{*T}\boldsymbol{r}$の計算は少し手間がかかりますが、以下のPythonコードで得られます。アンテナ数のループになってしまうのは仕方ないですが、基線数のループは避けられます。

PTY = np.zeros(antNum - 1, dtype=complex)
for ant_index in range(1,antNum):
  index0, index1 = np.where(ant0 == ant_index)[0].tolist(), np.where(ant1 == ant_index)[0].tolist()
  Y = np.zeros(blNum, dtype=complex)
  Y[index0] = -1.0j* antGain[ant_index].conjugate()* antGain[ant1[index0]]
  Y[index1] =  1.0j* antGain[ant_index]* antGain[ant0[index1]].conjugate()
  PTY[ant_index - 1] += Y.dot(resid)
#

まとめのコード

以上をまとめて、基線ベースのビジビリティからアンテナベースの位相を解く関数 clphase_solve() を示します。引数のVisはビジビリティで、Canonial orderで基線数の複素数が並んだベクトルです。iterNumは反復回数で、デフォルトで2にしています（S/N比がよければ十分に収束する）。返り値はアンテナベースの位相です。

def clphase_solve(Vis, iterNum = 2):
    # Initial setup
    Vis = Vis / abs(Vis)    # Normalization (discard amplitude)
    blNum = len(Vis); antNum = Bl2Ant(blNum)[0]
    ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    # Partial matrix
    PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum
    antGain = np.append(1.0 + 0.0j, Vis[KERNEL_BL[0:antNum-1]])
    #---- Loop for iteration
    for iter_index in range(iterNum):
        resid = Vis - (antGain[ant0] * antGain[ant1].conjugate())
        PTY = np.zeros(antNum - 1, dtype=complex)
        #---- Calculate P*T r
        for ant_index in range(1,antNum):
            index0 = np.where(ant0 == ant_index)[0].tolist()
            index1 = np.where(ant1 == ant_index)[0].tolist()
            Y = np.zeros(blNum, dtype=complex)
            Y[index0] = -1.0j* antGain[ant_index].conjugate()* antGain[ant1[index0]]
            Y[index1] =  1.0j* antGain[ant_index]* antGain[ant0[index1]].conjugate()
            PTY[ant_index - 1] += Y.dot(resid)
        #
        #---- Correction for the antenna-based phase
        antPhs  = np.append(0, np.angle(antGain[1:antNum]) + PTP_inv.dot(PTY.real))
        antGain = np.cos(antPhs) + 1.0j* np.sin(antPhs)
    #
    return antPhs
#

以上です。
もくじにもどる