Résolvez la phase

La phase $ \ phi $ est similaire à Delay

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

Cependant, il y a deux points à noter.

La phase est indéfinie $ 2n \ pi $. Il n'est pas évident que la phase basée sur la ligne de base soit, par exemple, $ 2 \ pi / 3 $ ou $ -4 \ pi / 3 $. Donc, $ \ hat {\ phi} _k = \ phi_j- \ phi_i \ pm 2n \ pi $, ce qui est problématique à résoudre ainsi que le délai.
Cette relation ne peut être utilisée que lorsque la phase de clôture est 0. Si le corps céleste est une source ponctuelle, la phase de fermeture sera 0, mais cette expression relationnelle ne vaut pas pour un corps céleste large avec une structure compliquée.

Concernant 2., nous allons y remédier en sélectionnant une source ponctuelle comme corps céleste d'étalonnage de phase. Concernant 1., dans le cas où la phase de chaque ligne de base est à l'intérieur de $ \ pm \ pi / 2 $, elle peut être traitée de la même manière que la solution de retard. Sinon, il est généralement préférable de résoudre avec une solution complexe.

Résoudre comme un nombre complexe

Le terme de phase du gain d'antenne est

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

Il est exprimé comme. Il satisfait l'équation (4.1) et n'a pas d'indétermination de $ 2n \ pi $. Cependant, puisqu'il s'agit d'une équation non linéaire, il est nécessaire de la répéter par approximation séquentielle. Dans ce qui suit, $ E_j = e ^ {i \ phi_j} $ est utilisé pour simplifier la notation. Supposons que la valeur initiale $ E ^ 0_i $ soit suffisamment proche de la solution optimale. Le résidu basé sur la ligne de base est $ r_k = E_k --E ^ 0_j E_i ^ {0 *} $. Le vecteur de correction $ \ boldsymbol {c} $ donné à $ \ phi ^ 0_j $ pour minimiser la norme résiduelle $ \ sum_k r_k r ^ * _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}{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}

Cela est illustré par l'équation. Notez que la phase de refant est mise à 0 et que la généralité n'est pas perdue, donc $ c_0 $ est toujours 0, pas inconnu. $ \ frac {\ partial} {\ partial \ phi_0} $ vaut également 0. Quand j'écris les éléments de la matrice du milieu $ P $,

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}

On dirait. Le gain de refant est $ E_0 = E ^ * _0 = 1 $. L'équation (4.3) est

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

Il est donné par les équations simultanées, donc si vous résolvez ceci pour $ \ boldsymbol {c} $

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

Ce sera. Ajoutez le $ \ boldsymbol {c} $ obtenu de cette manière à la valeur de phase initiale $ \ boldsymbol {\ phi} ^ 0 $, et ajoutez $ \ boldsymbol {\ phi} ^ 1 \ leftarrow \ boldsymbol {\ phi} ^ Améliorez la phase avec 0 + \ boldsymbol {c} $. Si vous répétez ceci, il finira par converger. Le jugement de convergence peut être jugé par la norme du vecteur de correction $ \ boldsymbol {c} $.

Économie de travail dans le calcul

Calculez manuellement $ P ^ {* T} P $ à l'avance.

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)

Et la même matrice que pour le delay. Par conséquent, la matrice inverse

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

Vous pouvez donc utiliser la même matrice que lorsque vous obtenez la solution pour le retard. Le calcul de $ P ^ {* T} \ boldsymbol {r} $ est un peu délicat, mais vous pouvez l'obtenir avec le code Python suivant. Il est inévitable de se retrouver dans une boucle avec le nombre d'antennes, mais une boucle avec le nombre de lignes de base peut être évitée.

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

Code récapitulatif

Pour résumer ce qui précède, nous montrons la fonction clphase_solve () qui résout la phase basée sur l'antenne à partir de la visibilité basée sur la ligne de base. L'argument Vis est la visibilité, qui est un vecteur dans lequel des nombres complexes de lignes de base sont alignés dans l'ordre canonique. IterNum est le nombre d'itérations, qui est fixé à 2 par défaut (convergence suffisante si le rapport S / N est bon). La valeur de retour est la phase basée sur l'antenne.

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
#

c'est tout. Retour à Mokuji

[PYTHON] Résolution de la phase des observations de l'interféromètre

Résolvez la phase

Résoudre comme un nombre complexe

Économie de travail dans le calcul

Code récapitulatif