[PYTHON] Lösung des komplexen Gewinns von Interferometer-Beobachtungen

Passen Sie die Amplitude und Phase an und lösen Sie den Antennengewinn als komplexe Zahl. Es gibt 2 Punkte:

Komplexe Zahlen, dh "reelle" anstelle von "Amplitude" und "Phase" Löse als "Teil und Imaginärteil". Da erwartet wird, dass der Real- und Imaginärteil der Sichtbarkeit normal verteilt ist, wird verhindert, dass die Methode der kleinsten Quadrate verzerrt wird (Amplitude). Daher ist die Methode der kleinsten Quadrate voreingenommen.
Die Refant-Phase wird auf 0 gesetzt und verliert ihre Allgemeinheit nicht. Mit anderen Worten, der Imaginärteil des Refant-Gewinns $ \ boldsymbol {G} _0 $ ist auf 0 festgelegt.

Wie im Fall von "Lösen der Amplitude" verwendet die Amplitude eine geeignete Skalierung.

\hat{V}_k = \boldsymbol{G}_j \boldsymbol{G}^*_i \tag{6.1}

Lösen wir nach $ \ boldsymbol {G} $. Der unbekannte Vektor ist $ 2N_a, wie $ \ boldsymbol {G} = (reG_0, reG_1, reG_2, \ dots, reG_ {Na-1}, imG_1, imG_2, \ dots, imG_ {Na-1}) $. - Dargestellt durch 1 $ reelle Zahl. Schreiben wir den Real- und Imaginärteil von $ als $ R bzw. I $, wie z. B. $ \ boldsymbol {G} = R + iI. Da die beobachtete Gleichung eine nichtlineare Gleichung ist, die Unbekannte in der Matrix $ P $ enthält, wird sie durch die iterative Methode gelöst. Wenn der Anfangswert $ \ boldsymbol {G ^ 0} $ ist, ist der Restvektor $ r_k = \ hat {V} _k - G ^ 0_j G_i ^ {* 0} $ und der Korrekturbetragvektor für dieses $ \ boldsymbol {c } $ Erfüllt $ \ boldsymbol {r} = P \ boldsymbol {c} $ und wird durch $ \ boldsymbol {c} = (P ^ TP) ^ {-1} P ^ T \ boldsymbol {r} $ erhalten .. Wenn Sie die Komponenten der Matrix $ P $ ausschreiben,

P = \left( \begin{array}{ccccccccc}
 R_1 &  R_0 & 0 & 0 & \cdots &  I_0 & 0 & 0 & \cdots \\
 R_2 &  0 & R_0 & 0 & \cdots &  0 & I_0 & 0 & \cdots \\
 0 & R_2 & R_1 & 0 & \cdots &  I_2 & I_0 & 0 & \cdots \\
 R_3 & 0 & 0 & R_0 & \cdots &  0 & 0 & I_0 & \cdots \\
 0 & R_3 & 0 & R_1 & \cdots &  I_3 & 0 & I_1 & \cdots \\
 0 & 0 & R_3 & R_2 & \cdots &  0 & I_3 & I_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots \\
 I_1 &  -I_0 & 0 & 0 & \cdots &  R_0 & 0 & 0 & \cdots \\
 I_2 &  0 & -I_0 & 0 & \cdots &  0 & R_0 & 0 & \cdots \\
 0 & I_2 & -I_1 & 0 & \cdots &  -R_2 & R_0 & 0 & \cdots \\
 I_3 & 0 & 0 & -I_0 & \cdots &  0 & 0 & R_0 & \cdots \\
 0 & I_3 & 0 & -I_1 & \cdots &  -R_3 & 0 & R_1 & \cdots \\
 0 & 0 & I_3 & -I_2 & \cdots &  0 & -R_3 & R_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots \\
\end{array} \right) = \left( \begin{array}{cc}
A & B \\
C & D \\
\end{array} \right) \tag{6.2}

Und wird durch vier Submatrix $ A, B, C, D $ dargestellt.

A =  \left( \begin{array}{ccccc}
 R_1 &  R_0 & 0 & 0 & \cdots \\
 R_2 &  0 & R_0 & 0 & \cdots \\
 0 & R_2 & R_1 & 0 & \cdots  \\
 R_3 & 0 & 0 & R_0 & \cdots \\
 0 & R_3 & 0 & R_1 & \cdots  \\
 0 & 0 & R_3 & R_2 & \cdots  \\
\vdots & \vdots & \vdots & \vdots & \ddots  \\
\end{array} \right) \\
B = \left( \begin{array}{cccc}
  I_0 & 0 & 0 & \cdots \\
  0 & I_0 & 0 & \cdots \\
 I_2 & I_0 & 0 & \cdots \\
  0 & 0 & I_0 & \cdots \\
 I_3 & 0 & I_1 & \cdots \\
 0 & I_3 & I_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
C = \left( \begin{array}{ccccc}
 I_1 &  -I_0 & 0 & 0 & \cdots  \\
 I_2 &  0 & -I_0 & 0 & \cdots  \\
 0 & I_2 & -I_1 & 0 & \cdots  \\
 I_3 & 0 & 0 & -I_0 & \cdots \\
 0 & I_3 & 0 & -I_1 & \cdots \\
 0 & 0 & I_3 & -I_2 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\
D = \left( \begin{array}{cccc}
R_0 & 0 & 0 & \cdots \\
0 & R_0 & 0 & \cdots \\
-R_2 & R_0 & 0 & \cdots \\
0 & 0 & R_0 & \cdots \\
-R_3 & 0 & R_1 & \cdots \\
0 & -R_3 & R_2 & \cdots \\
\vdots & \vdots & \vdots & \ddots \\
\end{array} \right) \\

Arbeitsersparnis bei der Matrixberechnung

Es wird wie folgt als Submatrix ausgedrückt.

P^T P = \left( \begin{array}{cc}
A^T A + C^TC & A^T B + C^T D \\
B^T A + D^T C & B^TB + D^T D
\end{array} \right) \tag{6.3}

Die Berechnung von $ A ^ TA $ erfolgt nach der Formel 5.7 von "Lösen der Amplitude". Verwenden Sie sie also so, wie sie ist. $ C ^ TC $ ist $ A ^ TA $ sehr ähnlich

C^TC = \left( \begin{array}{ccccc}
\sum_k I^2_k - I^2_0 & -I_0 I_1 & -I_0 I_2 & -I_0 I_3 & \cdots \\
-I_0 I_1 & \sum_k I^2_k - I^2_1 &  -I_1 I_2 & -I_1 I_3 & \cdots \\
-I_0 I_2 & -I_1 I_2 & \sum_k I^2_k - I^2_2 & -I_2 I_3 & \cdots \\
-I_0 I_3 & -I_1 I_3 & -I_2 I_3 & \sum_k I^2_k - I^2_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.4}

Es wird sein. Die Berechnung von $ B ^ TA $ und $ D ^ TC $ ist wie folgt.

B^TA = \left( \begin{array}{ccccc}
R_1 I_0 & \boldsymbol{R}\cdot\boldsymbol{I} - R_1 I_1 & R_1 I_2 & R_1 I_3 & \cdots \\
R_2 I_0 & R_2 I_1  & \boldsymbol{R}\cdot\boldsymbol{I} - R_2 I_2 &  R_2 I_3 & \cdots \\
R_3 I_0 & R_3 I_1 & R_3 I_2 & \boldsymbol{R}\cdot\boldsymbol{I} - R_3 I_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.5}

D^TC = \left( \begin{array}{ccccc}
R_0 I_1 & -\boldsymbol{R}\cdot\boldsymbol{I} + R_1 I_1 & R_2 I_1 & R_3 I_1 & \cdots \\
R_0 I_2 & R_1 I_2  & -\boldsymbol{R}\cdot\boldsymbol{I} + R_2 I_2 &  R_3 I_2 & \cdots \\
R_0 I_3 & R_1 I_3 & R_2 I_3 & -\boldsymbol{R}\cdot\boldsymbol{I} + R_3 I_3 & \cdots \\
\vdots & \vdots & \vdots & \vdots &\ddots
\end{array} \right) \tag{6.6}

Beachten Sie, dass $ I_0 = 0 $. Es wird auch berechnet durch $ A ^ TB = (B ^ TA) ^ T $, $ C ^ TD = (D ^ TC) ^ T $. Zusammenfassend gesagt, wenn Sie eine Funktion schreiben, um $ P ^ TP $ in Python zu erhalten, ist dies wie folgt.

def ATAmatrix(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 CTCmatrix(Gain):  # Gain is a vector of antenna-based gain amplitude (imag)
    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 ATBmatrix(Gain):  # Gain is a vector of antenna-based complex gain
    antNum = len(Gain); normG = Gain.real.dot(Gain.imag)
    PTP = np.zeros([antNum, antNum]) + Gain.imag
    for ant_index in range(antNum):
        PTP[ant_index,:] = Gain.real[ant_index]* Gain.imag + Gain.imag[ant_index]* Gain.real
        PTP[ant_index, ant_index] = 0.0
    #
    return PTP[1:antNum]
#
def PMatrix(CompSol):
    antNum = len(CompSol); matSize = 2*antNum-1
    PM = np.zeros([matSize, matSize])
    PM[0:antNum][:,0:antNum]   = ATAmatrix(CompSol.real) + CTCmatrix(CompSol.imag) # ATA + CTC
    PM[antNum:matSize][:,antNum:matSize] = ATAmatrix(CompSol.imag)[1:antNum][:,1:antNum] + CTCmatrix(CompSol.real)[1:antNum][:,1:antNum] # BTB + DTD
    PM[antNum:matSize][:,0:antNum]   = ATBmatrix(CompSol) # ATB + CTD
    PM[0:antNum][:,antNum:matSize]   = PM[antNum:matSize][:,0:antNum].T # BTA + DTC
    return PM
#

Das Argument CompSol ist der Vektor des anfänglichen Antennengewinns (komplexe Zahl). Die Rückgabewertmatrix PM hat reale Komponenten und eine Größe von $ (2N_a - 1, 2N_a - 1) $.

Berechnen Sie dann $ P ^ T r $. Die folgende Python-Funktion PTdotR () verwendet den Anfangswertvektor CompSol des Antennengewinns (komplexe Zahl) und den Restvektor (komplexe Zahl) Cresid als Argumente und berechnet und gibt $ P ^ T r $ zurück.

def PTdotR(CompSol, Cresid):
    antNum = len(CompSol)
    blNum = antNum* (antNum-1) / 2
    ant0, ant1= np.array(ANT0[0:blNum]), np.array(ANT1[0:blNum])
    PTR = np.zeros(2*antNum)
    for ant_index in range(antNum):
        index0 = np.where(ant0 == ant_index)[0].tolist()
        index1 = np.where(ant1 == ant_index)[0].tolist()
        PTR[range(ant_index)]          += (CompSol[ant_index].real* Cresid[index0].real + CompSol[ant_index].imag* Cresid[index0].imag)
        PTR[range(ant_index+1,antNum)] += (CompSol[ant_index].real* Cresid[index1].real - CompSol[ant_index].imag* Cresid[index1].imag)
        PTR[range(antNum, antNum+ant_index)]    += (CompSol[ant_index].imag* Cresid[index0].real - CompSol[ant_index].real* Cresid[index0].imag)
        PTR[range(antNum+ant_index+1,2*antNum)] += (CompSol[ant_index].imag* Cresid[index1].real + CompSol[ant_index].real* Cresid[index1].imag)
    #
    return PTR[range(antNum) + range(antNum+1, 2*antNum)]
#

Zusammenfassungscode

Die Funktion gainComplex (), die die Verstärkung des Antennenkomplexes durch das iterative Verfahren unter Verwendung der obigen PMatrix () und PTdotR () löst, ist unten gezeigt. Das Argument bl_vis ist eine auf der Basislinie basierende Sichtbarkeit und folgt der kanonischen Reihenfolge. Der Niter ist die Anzahl der Iterationen, aber ungefähr zwei Mal reichen aus, um zu konvergieren. Da $ P ^ {T} P $ eine reelle symmetrische Matrix ist, lösen wir die Gleichung, indem wir die untere Dreiecksmatrix L durch Cholesky-Zerlegung erhalten.

def gainComplex( bl_vis, niter=2 ):
    blNum  =  len(bl_vis)
    antNum =  Bl2Ant(blNum)[0]
    ant0, ant1, kernelBL = ANT0[0:blNum], ANT1[0:blNum], KERNEL_BL[range(antNum-1)].tolist()
    CompSol = np.zeros(antNum, dtype=complex)
    #---- Initial solution
    CompSol[0] = sqrt(abs(bl_vis[0])) + 0j
    CompSol[1:antNum] = bl_vis[kernelBL] / CompSol[0]
    #----  Iteration
    for iter_index in range(niter):
        PTP        = PMatrix(CompSol)
        L          = np.linalg.cholesky(PTP)         # Cholesky decomposition
        Cresid     = bl_vis - CompSol[ant0]* CompSol[ant1].conjugate()
        t          = np.linalg.solve(L, PTdotR(CompSol, Cresid))
        correction = np.linalg.solve(L.T, t)
        CompSol    = CompSol + correction[range(antNum)] + 1.0j* np.append(0, correction[range(antNum, 2*antNum-1)])
    #
    return CompSol
#

das ist alles.

Zurück zu Mokuji