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

Löse die Amplitude

In diesem Kapitel wird beschrieben, wie nur die Amplitude des Antennengewinns unter Verwendung der Amplitude der Sichtbarkeit gelöst wird, die durch die Interferometerbeobachtung erhalten wird. Anstatt es tatsächlich zu verwenden, empfehlen wir jedoch Lösen der komplexen Verstärkung durch Anpassen von Amplitude und Phase. Dies liegt daran, dass die Amplitude eher der Reisverteilung als der Normalverteilung folgt. Wenn also das Signal-Rausch-Verhältnis klein ist, neigt die Amplitude dazu, stark vorgespannt zu sein (während die komplexe Verstärkung normalverteilt ist). Wenn $ SNR> 5 $ ist, liegt dies nahe an einer Normalverteilung, sodass Sie die Methode in diesem Kapitel beschreiben können. Die Bedeutung dieses Kapitels besteht darin, das Verständnis [Lösen des Antennengewinns als komplexe Zahl] zu erleichtern (http://qiita.com/kamenoseiji/items/e0621a6a26ceecbbaac2).

Amplitude der Antennenverstärkung|G|Und basenbasierte Sichtbarkeitsamplitude|V|Die Beziehung zu ist wie folgt.

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

Auf der rechten Seite|V|_kIst die wahre Sichtbarkeitsamplitude auf der linken Seite|\hat{V}|_kIst die beobachtete Sichtbarkeitsamplitude. Die wahre Sichtbarkeitsamplitude ist erst bekannt, wenn sie kalibriert und bestimmt wurde.|V|_kGeteilt durch|\hat{V}|_k / |\hat{V}|_kDie Behandlung als linke Seite verliert nicht ihre Allgemeinheit. in diesem Fall,|G|Ist ein entsprechend skalierter Wert und hat nur eine Bedeutung als relative Amplitude. In diesem Fall die Beobachtungsgleichung(5.1)Ist das Folgende(5.2)Es sieht aus wie.

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

Einfache Lösung durch Logarithmus

Nehmen Sie den Logarithmus von Gleichung (5.2)

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

Und kann linearisiert werden. Das ist,

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

ist. Die mittlere ProzessionPToki auf der linken Seite\log |V|Vektor\boldsymbol{V},Auf der rechten Seite\log |G|Vektor\boldsymbol{G}Wenn du sagst\boldsymbol{V} = P \boldsymbol{G}Damit\boldsymbol{G}Lösen\boldsymbol{G} = (P^T P)^{-1} P^T \boldsymbol{V}ist.

Arbeitsersparnis bei der Berechnung

Berechnen Sie $ P ^ T P $ von Hand.

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}

Und seine inverse Matrix

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

Es wird sein. Kurz gesagt, die diagonale Komponente ist $ \ frac {2N_a -3} {2 (N_a-1) (N_a - 2)} $, und die nicht diagonale Komponente ist $ - \ frac {1} {2 (N_a-1) (N_a-) 2)} $. Diese Matrix kann in Python auch in einer Zeile codiert werden.

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

Die Berechnung von $ P ^ T \ boldsymbol {V} $ kann auf die gleiche Weise wie der Code beim Lösen der Verzögerung erfolgen.

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

Nichtlineare Lösung

Wenn Gleichung (5.2) ohne Logarithmus gelöst wird, wird sie zu einer nichtlinearen Gleichung. Unter der Annahme, dass der Anfangswert der Antennenverstärkungsamplitude $ G ^ 0 $ ist, ist der basenbasierte Rest $ r_k = \ hat {V} _k - G ^ 0_j G_i ^ 0 $. Der Korrekturvektor $ \ boldsymbol {c} $, der $ G ^ 0 $ gegeben wird, um die Restnorm $ \ displaystyle \ sum_k r_kr ^ * _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}{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}

Erfüllt die Gleichung. Die mittlere Matrix ist nicht linear, da das $ P $ -Element $ G $ enthält, das gelöst werden soll. Das Lösen von Gleichung (5.7) für $ c $ ergibt $ \ boldsymbol {c} = (P ^ T P) ^ {-1} P ^ T \ boldsymbol {r} $. Das so erhaltene $ \ boldsymbol {c} $ wird verwendet, um $ \ boldsymbol {G} $ zu verbessern.

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

Wiederholen Sie diesen Vorgang, bis $ \ boldsymbol {G} $ konvergiert, um eine Lösung zu erhalten. Die Konvergenz kann bestimmt werden, indem die Norm von $ \ boldsymbol {c} $ klein genug gemacht wird.

Arbeitsersparnis bei der Berechnung

Wie üblich wird $ P ^ T P $ manuell berechnet.

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}

Implementieren Sie dies in 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
#

Die inverse Matrix $ (P ^ T P) ^ {-1} $ kann nicht in einfacher Form geschrieben werden. Da es sich um eine echte symmetrische Matrix handelt, verwenden wir die Cholesky-Zerlegung, um die Gleichung zu lösen. Berechnen Sie $ \ boldsymbol {y} = P ^ {T} \ boldsymbol {r} $ im Voraus.

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

Verwenden Sie dann die Holesky-Zerlegung, um $ P ^ T P \ boldsymbol {c} = \ boldsymbol {y} $ für $ \ boldsymbol {c} $ zu lösen.

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

Addiere das so erhaltene $ \ boldsymbol {c} $ (Korrektur) zu $ \ boldsymbol {G} $ (antGain), verbessere es und wiederhole es, bis es konvergiert.

Zusammenfassungscode

Zusammenfassend wird der Python-Code zum Erhalten der Antennen-basierten Verstärkungsamplitude aus der Basislinien-basierten Sichtbarkeitsamplitude unten gezeigt. Die Funktion logamp_solve () ist eine logarithmisch linearisierte Lösung, die einen Vektor der Sichtbarkeitsamplitude als Argument verwendet und einen Vektor der Antennenverstärkungsamplitude zurückgibt. Die Funktion clamp_solve () ist eine nichtlineare Lösung und hat dieselben Spezifikationen wie logamp_solve () für ihre Argumente und ihren Rückgabewert. Intern wird logamp_solve (bl_amp) aufgerufen, um den Anfangswert der Antennenverstärkungsamplitude zu erhalten. Es ruft auch die Funktion PTPmatrix () auf, um $ (P ^ T P) ^ {-1} $ zu berechnen.

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
#

das ist alles. Tatsächlich ist es besser, Amplitude und Phase gleichzeitig als komplexe Zahlen zu lösen, als nur die Amplitude zu lösen. Zurück zu Mokuji