[PYTHON] Résolution de l'amplitude des observations d'interféromètre

Résoudre l'amplitude

Dans ce chapitre, nous expliquerons comment résoudre uniquement l'amplitude du gain d'antenne en utilisant l'amplitude de la visibilité obtenue par l'observation de l'interféromètre. Cependant, au lieu de l'utiliser réellement, nous recommandons Résoudre le gain complexe en faisant correspondre l'amplitude et la phase. En effet, l'amplitude suit la distribution de Rice plutôt que la distribution normale, donc lorsque le rapport signal sur bruit est petit, l'amplitude a tendance à être biaisée dans une large mesure (alors que le gain complexe est normalement distribué). Si $ SNR> 5 $, il sera proche d'une distribution normale, vous pouvez donc décrire la méthode dans ce chapitre. L'importance de ce chapitre est de faciliter la compréhension Comment résoudre le gain d'antenne en tant que nombre complexe.

Amplitude du gain d'antenne|G|Et amplitude de visibilité basée sur la ligne de base|V|La relation avec est la suivante.

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

Sur le côté droit|V|_kEst la vraie amplitude de visibilité, sur le côté gauche|\hat{V}|_kEst l'amplitude de visibilité observée. L'amplitude de visibilité réelle n'est pas connue tant qu'elle n'est pas étalonnée et déterminée.|V|_kDivisé par|\hat{V}|_k / |\hat{V}|_kTraiter comme le côté gauche ne perd pas sa généralité. dans ce cas,|G|Est une valeur correctement mise à l'échelle et n'a qu'une signification en tant qu'amplitude relative. Dans ce cas, l'équation d'observation(5.1)Est le suivant(5.2)On dirait.

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

Solution simple par logarithme

Prenant la logarithmique de l'équation (5.2)

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

Et peut être linéarisé. C'est,

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

est. La procession du milieuPToki, sur le côté gauche\log |V|Vecteur\boldsymbol{V},Sur le côté droit\log |G|Vecteur\boldsymbol{G}Si tu le dis\boldsymbol{V} = P \boldsymbol{G}Alors\boldsymbol{G}Résoudre\boldsymbol{G} = (P^T P)^{-1} P^T \boldsymbol{V}est.

Économie de travail dans le calcul

Calculez $ P ^ T P $ à la main.

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}

Et sa matrice inverse

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

Ce sera. En bref, la composante diagonale est $ \ frac {2N_a -3} {2 (N_a-1) (N_a --2)} $, et la composante hors diagonale est $ - \ frac {1} {2 (N_a-1) (N_a-) 2)} $. Cette matrice peut également être codée sur une ligne en Python.

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

Le calcul de $ P ^ T \ boldsymbol {V} $ peut être fait de la même manière que le code lors de la résolution du délai.

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

Solution non linéaire

Lors de la résolution de l'équation (5.2) sans prendre le logarithme, elle devient une équation non linéaire. En supposant que la valeur initiale de l'amplitude du gain d'antenne est $ G ^ 0 $, le résidu basé sur la ligne de base est $ r_k = \ hat {V} _k --G ^ 0_j G_i ^ 0 $. Le vecteur de correction $ \ boldsymbol {c} $ donné à $ G ^ 0 $ pour minimiser la norme résiduelle $ \ displaystyle \ sum_k r_kr ^ * _k $ est

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

Satisfait l'équation. La matrice du milieu est non linéaire car l'élément $ P $ contient $ G $ à résoudre. La résolution de l'équation (5.7) pour $ c $ donne $ \ boldsymbol {c} = (P ^ T P) ^ {-1} P ^ T \ boldsymbol {r} $. Le $ \ boldsymbol {c} $ ainsi obtenu est utilisé pour améliorer $ \ boldsymbol {G} $.

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

Répétez ceci jusqu'à ce que $ \ boldsymbol {G} $ converge pour obtenir une solution. La convergence peut être déterminée en rendant la norme de $ \ boldsymbol {c} $ suffisamment petite.

Économie de travail dans le calcul

Comme d'habitude, $ P ^ T P $ est calculé manuellement.

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}

Implémentez ceci en 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
#

La matrice inverse $ (P ^ T P) ^ {-1} $ ne peut pas être écrite sous une forme simple. Puisqu'il s'agit d'une matrice symétrique réelle, nous utiliserons la décomposition de Cholesky pour résoudre l'équation. Calculez à l'avance $ \ boldsymbol {y} = P ^ {T} \ boldsymbol {r} $.

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

Utilisez ensuite la décomposition de Holesky pour résoudre $ P ^ T P \ boldsymbol {c} = \ boldsymbol {y} $ pour $ \ boldsymbol {c} $.

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

Ajoutez le $ \ boldsymbol {c} $ (correction) ainsi obtenu à $ \ boldsymbol {G} $ (antGain), améliorez-le et répétez jusqu'à ce qu'il converge.

Code récapitulatif

Pour résumer ce qui précède, le code Python permettant d'obtenir l'amplitude de gain basée sur l'antenne à partir de l'amplitude de visibilité basée sur la ligne de base est illustré ci-dessous. La fonction logamp_solve () est une solution logarithmique linéarisée qui prend un vecteur d'amplitude de visibilité comme argument et renvoie un vecteur d'amplitude de gain d'antenne. La fonction clamp_solve () est une solution non linéaire et a les mêmes spécifications que logamp_solve () pour ses arguments et sa valeur de retour. En interne, logamp_solve (bl_amp) est appelé pour obtenir la valeur initiale de l'amplitude du gain de l'antenne. Il appelle également la fonction PTPmatrix () pour calculer $ (P ^ T P) ^ {-1} $.

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
#

c'est tout. En fait, il vaut mieux résoudre l'amplitude et la phase en tant que nombres complexes en même temps plutôt que de résoudre uniquement l'amplitude. Retour à Mokuji