[PYTHON] Lösen Sie die Verzögerung der Interferometerbeobachtung

Einführung

Die Verzögerung $ \ tau $ hat eine Beziehung von $ \ phi = \ phi_0 + 2 \ pi \ nu \ tau $ zwischen der Phase und der Winkelfrequenz $ 2 \ pi \ nu $. Die Verzögerung $ \ hat {\ tau} $, die sich aus der Sichtbarkeit für jede Basislinie ergibt, gilt für die Verzögerung $ \ tau $ für jede Antenne.

\hat{\tau}_k = \tau_j - \tau_i \tag{3.1}

Es wird vertreten durch. $ i, j, k $ folgen kanonischer Reihenfolge. Diese Beziehung gilt auch dann, wenn allen Antennen ein gemeinsamer Verzögerungsversatz zugewiesen wird, sodass das Setzen der Verzögerung der Referenzantenne (Refant) auf 0 nicht ihre Allgemeingültigkeit verliert.

Minimum-Quadrat-Methode

Die Gleichung, die die Verzögerung $ \ tau_i $ für jede Antenne aus der Verzögerung $ \ hat {\ tau} _k $ löst, die für jede Basislinie nach der Minimum-Square-Methode gemessen wurde, lautet

\left(
\begin{array}{c}
\hat{\tau}_0  \\
\hat{\tau}_1  \\
\hat{\tau}_2  \\
\hat{\tau}_3  \\
\hat{\tau}_4  \\
\hat{\tau}_5  \\
\vdots
\end{array}
\right) = \left( \begin{array}{cccc}
 1 & 0 & 0 & \cdots \\
 0 & 1 & 0 & \cdots \\
-1 & 1 & 0 & \cdots \\
 0 & 0 & 1 & \cdots \\
-1 & 0 & 1 & \cdots \\
 0 & -1 & 1 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{array} \right)
\left(
\begin{array}{c}
\tau_1  \\
\tau_2  \\
\tau_3  \\
\vdots
\end{array}
\right) \tag{3.2}

Es wird so sein. Da wir $ \ tau_0 = 0 $ setzen, beträgt die Anzahl der Unbekannten $ N_a -1 $. Die mittlere $ (N_a - 1) \ mal N_b $ Größenmatrix ist $ P $, der beobachtete Betrag $ \ hat {\ tau} _k $ Vektor ist $ \ boldsymbol {Y} $, unbekannte Zahl $ \ tau_i $ Wenn der Vektor $ \ boldsymbol {X} $ ist, sind die simultanen Gleichungen, die die Verzögerung für jede Antenne im Sinne des minimalen Quadrats erhalten

(P^TP) \boldsymbol{X} = P^T \boldsymbol{Y} \tag{3.3}

Es wird sein. Daher ist $ \ boldsymbol {X} = (P ^ TP) ^ {-1} P ^ T \ boldsymbol {Y} $, so dass $ \ boldsymbol {X} $ erhalten werden kann.

Arbeitsersparnis bei der Berechnung

Wenn die Anzahl der Antennen und die Anzahl der Basislinien zunimmt, wird der Rechenaufwand für $ P $ enorm. ALMA verfügt über maximal 64 Antennen und maximal 2016 Basislinien, daher ist es schwierig, dies zu lösen. Glücklicherweise wird die Berechnung von $ P ^ TP $ bei Verwendung der kanonischen Reihenfolge wie folgt vereinfacht.

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{3.4}

Kurz gesagt, die diagonale Komponente ist $ N_a-1 $ und die anderen Komponenten sind $ -1 $. Da dies eine echte symmetrische Matrix ist, kann sie unter Verwendung der Cholesky-Zerlegung gelöst werden, um die untere Dreiecksmatrix $ L $ zu erhalten, die als $ LL ^ T = P $ ausgedrückt werden kann. Darüber hinaus kann die inverse Matrix von $ (P ^ TP) $ leicht direkt berechnet werden.

(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) \tag{3.5}

Es wird sein. Kurz gesagt, eine $ (N_a - 1) \ times (N_a - 1) $ - Größenmatrix, deren diagonale Komponenten $ \ frac {2} {N_a} $ und die anderen Komponenten $ \ frac {1} {N_a} $ sind. ist. Wenn Sie es in Python-Code schreiben, können Sie antNum als Anzahl der Antennen in einer Zeile verwenden.

import numpy as np
PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum

Alles was bleibt ist die Berechnung von $ P ^ T \ boldsymbol {Y} $. Selbst wenn wir vermeiden, $ P ^ T $ zu berechnen, was eine sehr spärliche Matrix ist, möchten wir Schleifen mit der Anzahl der Basislinien vermeiden, also führen wir Vektoroperationen innerhalb der Schleife mit der Anzahl der Antennen aus, wie unten gezeigt.

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

Die Listen ANT0 und ANT1 werden im Voraus unter Bezugnahme auf diesen Artikel Interferometer-Basisliniennummerierung durch kanonische Bestellung erstellt. index0 hat "ant_index als Basisendpunkt" Wenn man das innere Produkt von $ (P ^ TP) ^ {-1} $ und $ P ^ TY $ nimmt, das auf diese Weise erzeugt wurde, kann ein unbekannter Vektor $ X $ erhalten werden.

Delay_ant = [0.0] + (PTP_inv.dot(PTY)).tolist()

0.0 wird am Anfang der Liste eingefügt, da es sich um eine Liste mit Refant handelt. Auf diese Weise wird die Verzögerung für jede Basislinie

Delay_bl = Delay_ant[ANT1] - Delay_ant[ANT0]

Sie können es so bekommen.

Zusammenfassungscode

Unten finden Sie eine Python-Funktion, die eine Antennenverzögerung von einer Basislinienverzögerung erhält. Das Argument bl_delay ist eine Reihe von auf Baseline basierenden Verzögerungen, die in Canonical Ordering aufgereiht werden müssen. Der Rückgabewert ist eine Antennenverzögerung, die 0 ist, da sie mit Refant beginnt.

def cldelay_solve(bl_delay):
    blNum = len(bl_delay); antNum = Bl2Ant(blNum)[0]
    ant0, ant1 = np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    PTP_inv = (np.diag(np.ones(antNum - 1)) + 1.0) / antNum
    PTY = np.zeros(antNum - 1)
    for ant_index in range(1, antNum):
        index0 = np.where(ant0 == ant_index)[0].tolist()
        index1 = np.where(ant1 == ant_index)[0].tolist()
        PTY[ant_index - 1] += np.sum(bl_delay[index0]) - np.sum(bl_delay[index1])
    #
    return np.array( [0.0] + (PTP_inv.dot(PTY)).tolist() )
#

das ist alles.

Zurück zu Mokuji