LoginSignup
2

More than 5 years have passed since last update.

干渉計観測の位相を解く

Last updated at Posted at 2017-01-13

位相を解く

位相$\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
#

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

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
2