[PYTHON] Lösen der Phase von Interferometer-Beobachtungen

Löse die Phase

Phase $ \ phi $ ähnelt Verzögerung

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

Es sind jedoch zwei Punkte zu beachten.

1. Die Phase hat $ 2n \ pi $ Unbestimmtheit. Es ist nicht offensichtlich, ob die basenbasierte Phase beispielsweise $ 2 \ pi / 3 $ oder $ -4 \ pi / 3 $ ist. Also $ \ hat {\ phi} _k = \ phi_j- \ phi_i \ pm 2n \ pi $, was sowohl problematisch zu lösen als auch zu verzögern ist.
2. Diese Beziehung kann nur verwendet werden, wenn die Abschlussphase 0 ist. Wenn der Himmelskörper eine Punktquelle ist, ist die Abschlussphase 0, aber dieser relationale Ausdruck gilt nicht für einen breiten Himmelskörper mit einer komplizierten Struktur.

In Bezug auf 2. werden wir uns damit befassen, indem wir eine Punktquelle als Himmelskörper für die Phasenkalibrierung auswählen. In Bezug auf 1. kann in dem Fall, in dem die Phase jeder Basislinie innerhalb von $ \ pm \ pi / 2 $ liegt, sie auf die gleiche Weise wie die Verzögerungslösung behandelt werden. Ansonsten ist es im Allgemeinen besser, mit einer komplexen Lösung zu lösen.

Löse als komplexe Zahl

Der Phasenterm des Antennengewinns ist

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

Es wird ausgedrückt als. Es erfüllt Gleichung (4.1) und hat keine Unbestimmtheit von $ 2n \ pi $. Da es sich jedoch um eine nichtlineare Gleichung handelt, muss sie durch sequentielle Approximation wiederholt werden. Im Folgenden wird $ E_j = e ^ {i \ phi_j} $ verwendet, um die Notation zu vereinfachen. Angenommen, der Anfangswert $ E ^ 0_i $ liegt nahe genug an der optimalen Lösung. Das auf der Basislinie basierende Residuum ist $ r_k = E_k - E ^ 0_j E_i ^ {0 *} $. Der Korrekturvektor $ \ boldsymbol {c} $, der $ \ phi ^ 0_j $ gegeben wird, um die Restnorm $ \ sum_k r_k r ^ * _k $ zu minimieren, ist

\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}

Es wird durch die Gleichung gezeigt. Beachten Sie, dass die Refant-Phase auf 0 gesetzt ist und die Allgemeinheit nicht verloren geht, sodass $ c_0 $ immer 0 ist und nicht unbekannt. $ \ frac {\ partiell} {\ partiell \ phi_0} $ ist ebenfalls 0. Wenn ich die Elemente der mittleren Matrix $ P $ schreibe,

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}

Es sieht aus wie. Der Refant-Gewinn beträgt $ E_0 = E ^ * _0 = 1 $. Gleichung (4.3) ist

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

Es ist durch die simultanen Gleichungen gegeben. Wenn Sie dies also für $ \ boldsymbol {c} $ lösen

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

Es wird sein. Addieren Sie das auf diese Weise erhaltene $ \ boldsymbol {c} $ zum anfänglichen Phasenwert $ \ boldsymbol {\ phi} ^ 0 $ und fügen Sie $ \ boldsymbol {\ phi} ^ 1 \ leftarrow \ boldsymbol {\ phi} ^ hinzu Verbessere die Phase als 0 + \ boldsymbol {c} $. Wenn Sie dies wiederholen, wird es schließlich konvergieren. Die Konvergenzbeurteilung kann anhand der Norm des Korrekturvektors $ \ boldsymbol {c} $ beurteilt werden.

Arbeitsersparnis bei der Berechnung

Berechnen Sie $ P ^ {* T} P $ im Voraus von Hand.

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)

Und die gleiche Matrix wie für die Verzögerung. Daher die inverse Matrix

(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)

Sie können also dieselbe Matrix verwenden, in der Sie die Lösung für die Verzögerung erhalten haben. Die Berechnung von $ P ^ {* T} \ boldsymbol {r} $ ist etwas schwierig, aber Sie können sie mit dem folgenden Python-Code erhalten. Es ist unvermeidlich, in einer Schleife mit der Anzahl der Antennen zu enden, aber eine Schleife mit der Anzahl der Basislinien kann vermieden werden.

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)
#

Zusammenfassungscode

Zusammenfassend zeigen wir die Funktion clphase_solve (), die die Antennenphase aus der Basisliniensichtbarkeit löst. Das Argument Vis ist Sichtbarkeit, ein Vektor, in dem eine komplexe Anzahl von Basislinien in kanonialer Reihenfolge angeordnet ist. IterNum ist die Anzahl der Iterationen, die standardmäßig auf 2 gesetzt ist (sie konvergiert ausreichend, wenn das S / N-Verhältnis gut ist). Der Rückgabewert ist die Antennenphase.

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
#

das ist alles. Zurück zu Mokuji