\section{多阵列随机信号 DoA 估计}
\label{sec:multiArrayStoDoaEst}

本节是多阵列随机 DoA 支线，用于说明多颗卫星或多个阵列如何通过几何映射共享同一组全局未知参数。它与后续多星 DoA--Doppler 主线在“多阵列几何约束”上相通，但不引入参考星 Doppler 参数 $\nu_{\mathrm{ref}}$，因此阅读主线时可将本节视为多阵列 DoA 的背景和对照。

\subsection{参考阵列角度参数化的多阵列随机 DoA 模型}
\label{subsec:projection_operator}

在多阵列（多卫星）协同 \gls{doa} 估计中，不同阵列由于位置与指向各不相同，即便面对同一组地面目标，其局部观测角度也不一致。因此，需要建立从某一参考阵列的角度参数到任意阵列角度参数的映射，确保所有阵列在 \gls{mle} 中共享一套全局未知参数。本节给出这一映射的几何构造，并定义角度映射算子。

\subsubsection{参考阵列与角度映射}

为便于参数表示，将第 $r$ 个阵列（例如 $r=0$）作为参考阵列，并将其局部来向集合记为
%
\begin{align}
  \bm{\Theta}_\mathrm{ref}
  \triangleq
  \left[\bm{\theta}_{\mathrm{ref},0},\ldots,\bm{\theta}_{\mathrm{ref},K-1}\right],
\end{align}
%
对于其余阵列 $\ell \neq r$，通过几何关系（卫星轨道、姿态以及坐标变换）构造从参考阵列角度到第 $\ell$ 个阵列局部角度的已知角度映射算子 $\mathcal T_\ell(\,\cdot\,)$，使得
%
\begin{align}
  \bm{\Theta}_\ell
  = \mathcal T_\ell(\bm{\Theta}_\mathrm{ref}),\qquad \ell=0,\ldots,L-1,
\end{align}
%
则第 $\ell$ 个阵列的导向矩阵可以写成
%
\begin{align}
  \mathbf{A}_\ell(\bm{\Theta})
  =
  \mathbf{A}_\ell\left(\bm{\Theta}_\ell\right)
  =
  \mathbf{A}_\ell\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right),
\end{align}
%
从而所有阵列的协方差矩阵 $\{\mathbf{R}_\ell\}$ 最终都只依赖于同一组参考角度参数 $\bm{\Theta}_\mathrm{ref}$。

\subsubsection{卫星--地面几何关系}

设地面目标、第 $\ell$ 个卫星阵列、参考阵列（例如卫星 $0$）的 \gls{ecef} 坐标分别为 $\mathbf{p}_g, \mathbf{p}_\ell, \mathbf{p}_{\mathrm{ref}} \in \mathbb{R}^3$，从卫星 $\ell$ 指向地面目标的全局方向向量定义为
\begin{align}
    \mathbf{u}_{\ell}^{\mathrm{global}}
    =
    \frac{\mathbf{p}_g - \mathbf{p}_\ell}{
          \|\mathbf{p}_g - \mathbf{p}_\ell\|}
    \in \mathbb{R}^3 .
\end{align}

几何上，$\mathbf{u}_{\ell}^{\mathrm{global}}$ 表示地面目标在 \gls{ecef} 坐标系下的方向，而 \gls{doa} 估计所依赖的角度参数必须定义在卫星阵列的局部坐标系内。因此需要构造卫星局部坐标系并将全局方向向量旋转至其坐标系下。

\subsubsection{局部坐标系的建立}

对于卫星 $\ell$，首先定义其局部坐标系在 \gls{ecef} 坐标系下的坐标轴全局方向向量 $\{\mathbf{x}_\ell, \mathbf{y}_\ell, \mathbf{z}_\ell\}$。不失一般性，假设阵列法向指向地心，也即取局部 $\mathbf{z}$ 轴为
%
\begin{align*}
  \mathbf{z}_\ell = -\frac{\mathbf{p}_\ell}{\|\mathbf{p}_\ell\|},
\end{align*}
%
然后选取一个与 $\mathbf{z}_\ell$ 不平行的向量
%
\begin{align*}
  \mathbf{v} = \begin{cases}
    [1,0,0]^\Tran, & |z_{\ell,x}| < 0.9,\\
    [0,1,0]^\Tran, & \text{otherwise},
  \end{cases}
\end{align*}
%
其中 $z_{\ell,x}$ 表示 $\mathbf{z}_\ell$ 在 \gls{ecef} 坐标线 $x$ 方向的分量，以此构造局部 $x$ 轴
%
\begin{align*}
  \mathbf{x}_\ell = \frac{\mathbf{v} \times \mathbf{z}_\ell}{ \|\mathbf{v} \times \mathbf{z}_\ell\| },
\end{align*}
%
最后，局部 $y$ 轴由右手法则给出：
%
 \begin{align*}
  \mathbf{y}_\ell = \mathbf{z}_\ell \times \mathbf{x}_\ell .
\end{align*}
%
由此，从\emph{局部坐标系}到\emph{全局坐标系}的旋转矩阵可写为
%
\begin{align}
  \boldsymbol{\Omega}_\ell
  = \begin{bmatrix}
      \mathbf{x}_\ell, \mathbf{y}_\ell, \mathbf{z}_\ell
    \end{bmatrix}
  \in \mathbb{R}^{3\times 3},
\end{align}
%
其中列向量为局部坐标轴在全局坐标系下的方向。

利用该旋转矩阵可将 \gls{ecef} 坐标系下的全局方向向量旋转到阵列的局部坐标系中
%
\begin{align}
  \mathbf{u}_{\ell}^{\mathrm{local}}
  = \boldsymbol{\Omega}_\ell^\Tran \mathbf{u}_{\ell}^{\mathrm{global}}
  = \begin{bmatrix}
    u_{\ell,x}, u_{\ell,y}, u_{\ell,z}
  \end{bmatrix}^\Tran.
\end{align}
%
由此得到以阵列为参考系的局部方位角与俯仰角分别为
\begin{align}
    \alpha_\ell^{\mathrm{loc}} = \operatorname{atan2}(u_{\ell,y}, u_{\ell,x}),\qquad
    \beta_\ell^{\mathrm{loc}}  = \arcsin(u_{\ell,z}).
\end{align}

\subsubsection{角度映射算子}

我们希望在似然函数中使用统一的一组未知方向参数来描述所有阵列的视线方向。为此，选定参考阵列（例如卫星 $0$），将其局部方向参数
%
\begin{align*}
  \bm{\theta}_{\mathrm{ref}} = (\alpha_{\mathrm{ref}}^{\mathrm{loc}}, \beta_{\mathrm{ref}}^{\mathrm{loc}})
\end{align*}
%
作为全局未知参数。对应的局部单位方向向量为
%
\begin{align}
  \mathbf{u}_{\mathrm{ref}}^{\mathrm{local}}(\bm{\theta}_{\mathrm{ref}})
  =
  \begin{bmatrix}
    \cos\beta_{\mathrm{ref}}^{\mathrm{loc}} \cos\alpha_{\mathrm{ref}}^{\mathrm{loc}} \\[1mm]
    \cos\beta_{\mathrm{ref}}^{\mathrm{loc}} \sin\alpha_{\mathrm{ref}}^{\mathrm{loc}} \\[1mm]
    \sin\beta_{\mathrm{ref}}^{\mathrm{loc}}
  \end{bmatrix} .
\end{align}
%
利用局部到全局的旋转矩阵 \(\boldsymbol{\Omega}_{\mathrm{ref}}\)，可得到参考阵列视线在 \gls{ecef} 坐标系下的方向
%
\begin{align}
  \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}})
  = \boldsymbol{\Omega}_{\mathrm{ref}}
    \mathbf{u}_{\mathrm{ref}}^{\mathrm{local}}(\bm{\theta}_{\mathrm{ref}}).
\end{align}
%
接下来即可根据参考阵列的 \gls{ecef} 坐标 $\mathbf{p}_{\mathrm{ref}}$ 与视线方向 $\mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}})$ 恢复地面目标的空间位置。

\begin{enumerate}
  \item 假定目标位于地球表面，并将地球近似为半径为 $\mathcal{R}_\mathrm{E}$ 的球体，则地面点 $\mathbf{p}_g$ 的坐标可表示为
  %
  \begin{align}
    \mathbf{p}_g
    = \mathbf{p}_{\mathrm{ref}}
      + r_{\mathrm{ref}}\,
        \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}}),
  \end{align}
  %
  其中距离参数 $r_{\mathrm{ref}} > 0 $由射线与球面 $\|\mathbf{p}_g\| = \mathcal{R}_\mathrm{E}$ 的交点条件确定。由于 $\mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}$ 为单位向量，有
  %
  \begin{align*}
    \left\|\mathbf{p}_{\mathrm{ref}} + r_{\mathrm{ref}}\mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}\right\|^2
    = \mathcal{R}_\mathrm{E}^2
    \quad\Longrightarrow\quad
    r_{\mathrm{ref}}^2
    + 2\,\mathbf{p}_{\mathrm{ref}}^\Tran \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}} r_{\mathrm{ref}}
    + \left\|\mathbf{p}_{\mathrm{ref}}\right\|^2 - \mathcal{R}_\mathrm{E}^2
    = 0,
  \end{align*}
  %
  从而
  %
  \begin{align}
    r_{\mathrm{ref}}
    =
    - \mathbf{p}_{\mathrm{ref}}^\Tran
      \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}})
    - \sqrt{
      \left(\mathbf{p}_{\mathrm{ref}}^\Tran
            \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}})\right)^2
      - \left(\|\mathbf{p}_{\mathrm{ref}}\|^2 - \mathcal{R}_\mathrm{E}^2\right)
    },
  \end{align}
  %
  其中取使地面点位于参考阵列视线前方的正根。

  \item 假定目标位于地球椭球表面，其 \gls{ecef} 坐标满足
  %
  \begin{align}
    \mathbf{p}_g^\Tran \mathbf{B}\,\mathbf{p}_g = 1,
  \end{align}
  %
  其中 $\mathbf{B}$ 为描述地球椭球的对称正定矩阵。  例如，对于赤道半轴为 $a$、极半轴为 $b$ 的轴对称椭球，有 $\mathbf{B} = \operatorname{diag}\bigl(1/a^2,\,1/a^2,\,1/b^2\bigr)$，则地面点可写为
  %
  \begin{align}
    \mathbf{p}_g
    = \mathbf{p}_{\mathrm{ref}}
      + r_{\mathrm{ref}}\,
        \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}}),
  \end{align}
  %
  其中距离参数 $r_{\mathrm{ref}} > 0$ 由射线与椭球的交点条件 $\mathbf{p}_g^\Tran \mathbf{B}\,\mathbf{p}_g = 1$ 确定。将上述射线方程代入可得关于 $r_{\mathrm{ref}}$ 的二次方程
  %
  \begin{align}
    A r_{\mathrm{ref}}^2 + B r_{\mathrm{ref}} + C = 0,
  \end{align}
  %
  其中 $A = (\mathbf{u}_{\mathrm{ref}}^{\mathrm{global}})^\Tran \mathbf{B} \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}$，$B = 2\mathbf{p}_{\mathrm{ref}}^\Tran \mathbf{B} \mathbf{u}_{\mathrm{ref}}^{\mathrm{global}}$，$C = \mathbf{p}_{\mathrm{ref}}^\Tran \mathbf{B} \mathbf{p}_{\mathrm{ref}} - 1$，从而
  %
  \begin{align}
    r_{\mathrm{ref}}
    = \frac{-B \pm \sqrt{B^2 - 4AC}}{2A},
  \end{align}
  其中取沿视线方向、对应地面点位于参考阵列视线前方的正根。

\end{enumerate}

在得到 $\mathbf{p}_g$ 后，第 $\ell$ 个阵列（卫星）指向该地面点的全局单位方向向量为
%
\begin{align}
    \mathbf{u}_\ell^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}})
    =
    \frac{
      \mathbf{p}_g - \mathbf{p}_\ell
    }{
      \left\|\mathbf{p}_g - \mathbf{p}_\ell\right\|
    }.
\end{align}
%
利用第 $\ell$ 个阵列的旋转矩阵 $\boldsymbol{\Omega}_\ell$ 可得到其局部坐标系中的方向
%
\begin{align}
  \mathbf{u}_\ell^{\mathrm{local}}(\bm{\theta}_{\mathrm{ref}})
  = \boldsymbol{\Omega}_\ell^\Tran
    \mathbf{u}_\ell^{\mathrm{global}}(\bm{\theta}_{\mathrm{ref}})
\end{align}
%
由此，第 $\ell$ 个阵列在自身局部坐标系下的方位角与俯仰角可写为
%
\begin{align}
    \alpha_\ell^{\mathrm{loc}}(\bm{\theta}_{\mathrm{ref}})
    = \operatorname{atan2} \left(
       \left[\mathbf{u}_\ell^{\mathrm{local}}(\bm{\theta}_{\mathrm{ref}})\right]_y,
       \left[\mathbf{u}_\ell^{\mathrm{local}}(\bm{\theta}_{\mathrm{ref}})\right]_x
       \right),\quad
    \beta_\ell^{\mathrm{loc}}(\bm{\theta}_{\mathrm{ref}})
    = \arcsin\left(
       \left[\mathbf{u}_\ell^{\mathrm{local}}(\bm{\theta}_{\mathrm{ref}})\right]_z
       \right).
\end{align}
%
于是可以定义从参考阵列角度到第 $\ell$ 个阵列角度的\emph{非线性}角度映射算子
%
\begin{align}
  \bm{\theta}_\ell
  = \mathcal T_\ell(\bm{\theta}_{\mathrm{ref}})
  \triangleq
  \Bigl(
    \alpha_\ell^{\mathrm{loc}}(\bm{\theta}_{\mathrm{ref}}),\,
    \beta_\ell^{\mathrm{loc}}(\bm{\theta}_{\mathrm{ref}})
  \Bigr).
\end{align}
%
可以看出，$\mathcal T_\ell(\,\cdot\,)$ 由参考阵列与第 $\ell$ 个阵列的 \gls{ecef} 坐标 $\mathbf{p}_{\mathrm{ref}}, \mathbf{p}_\ell$、对应的姿态矩阵 $\boldsymbol{\Omega}_{\mathrm{ref}}, \boldsymbol{\Omega}_\ell$，以及地球模型参数 $\mathcal{R}_\mathrm{E}$ 唯一确定。因此，在似然函数中可以仅以参考阵列的方向参数$\bm{\theta}_{\mathrm{ref}}$ 作为未知量，并通过几何映射 $\mathcal T_\ell$ 将其映射到各个阵列的局部角度，实现多阵列的统一参数化建模。

\subsection{参考阵列角度参数化下的多阵列似然函数求取}

假设信源服从圆对称复高斯分布 $\mathbf{x}_t \sim \mathcal{CN} \left( 0,\mathbf R_x \right)$，为简化讨论，假设各信源相互不相关，则 $\mathbf R_x$ 为对角矩阵，且其协方差矩阵非奇异 $\det \left( \mathbf R_x \right) \ne 0$，时域不相关
%
\begin{equation}
  \begin{cases}
    \mathbb{E} \left[\mathbf{x}_t \mathbf{x}_s^\Herm \right] = \mathbf R_x \delta(t-s)
    \\
    \mathbb{E} \left[ \mathbf{x}_t \mathbf{x}_s^\Tran \right] = 0
  \end{cases}, \quad\forall t,s.
\end{equation}
%
由于噪声 $\mathbf{n}_{\ell,t} \sim \mathcal{CN}(0,\sigma_{n,\ell}^2 \mathbf{I}_{M_\ell})$，且与 $\mathbf{x}_t$ 不相关。则在给定 $(\bm{\Theta},\mathbf R_x,\{\sigma_{n,\ell}^2\})$ 的条件下，$\{\mathbf{y}_{\ell,t}\}$ 为相互独立的圆对称复高斯向量，也即 $\mathbf{y}_{\ell,t} \sim \mathcal{CN}(\mathbf{0},\mathbf{R}_\ell)$，因而联合似然函数为
%
\begin{align}
  L\left(\{\mathbf{R}_\ell\},\bm{\Theta}\middle|\{\mathbf{y}_{\ell,t}\}\right)
  &= \prod_{\ell=0}^{L-1} \prod_{t=0}^{T-1}
  \frac{1}{\pi^{M_\ell}\det(\mathbf{R}_\ell)}
  \exp\!\left(-\mathbf{y}_{\ell,t}^\Herm \mathbf{R}_\ell^{-1}\mathbf{y}_{\ell,t}\right)
  \nonumber\\
  &= \prod_{\ell=0}^{L-1}
  \pi^{-M_\ell T} \left[\det(\mathbf{R}_\ell)\right]^{-T}
  \exp\!\left(
  - T\,\tr\left(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell\right)
  \right).
\end{align}
%
其中第 $\ell$ 个阵列的样本协方差矩阵为 $\hat{\mathbf{R}}_\ell \triangleq \displaystyle\frac{1}{T} \sum_{t=0}^{T-1} \mathbf{y}_{\ell,t}\mathbf{y}_{\ell,t}^\Herm$，对应的对数似然函数为
%
\begin{equation}
  \ln \mathcal{L}\left(\{\mathbf{R}_\ell\},\bm{\Theta}\right)
  = \mathrm{Const}
  - T \sum_{\ell=0}^{L-1} \log\det(\mathbf{R}_\ell)
  - T \sum_{\ell=0}^{L-1} \tr\left(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell\right).
\end{equation}

\subsubsection{原始问题建模}

引入参考局部来向集合 $\bm{\Theta}_\mathrm{ref}$ 将多阵列对数似然函数写为
%
\begin{equation}
  \ln \mathcal{L}\left(\{\mathbf{R}_\ell\},\bm{\Theta}_\mathrm{ref}\right)
  = \mathrm{Const}
  - T \sum_{\ell=0}^{L-1} \log\det(\mathbf{R}_\ell)
  - T \sum_{\ell=0}^{L-1} \tr\left(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell\right),
\end{equation}
%
其中 $\mathbf{R}_\ell
  =
  \mathbf{A}_\ell\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  \,\mathbf R_x\,
  \mathbf{A}_\ell^\Herm\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  + \sigma_{n,\ell}^2 \mathbf{I}_{M_\ell}.
$ 则多阵列的 \gls{sml}/\gls{uml} 优化问题可以表述为
%
\begin{alignat}{2}
  \label{eq:multi_arr_sto_ml_doa_opt}
  &\underset{\bm{\Theta}_\mathrm{ref},\,\mathbf R_x,\,\{\sigma_{n,\ell}^2\}}{\mathrm{minimize}} \quad
  && \sum_{\ell=0}^{L-1}
  \left[
  \log\det(\mathbf{R}_\ell)
  + \tr\left(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell\right)
  \right] \nonumber\\
  &\text{subject to} \quad
  && \mathbf{R}_\ell
  =
  \mathbf{A}_\ell\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  \,\mathbf R_x\,
  \mathbf{A}_\ell^\Herm\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  + \sigma_{n,\ell}^2 \mathbf{I}_{M_\ell},\ \ell=0,\ldots,L-1, \nonumber\\
  &&& \bm{\Theta}_\mathrm{ref} \in \varTheta_\mathrm{ref},\quad
  \mathbf R_x \succeq \mathbf{0},\quad
  \sigma_{n,\ell}^2 > 0.
\end{alignat}
%
可以看到，相比 $\eqref{eq:sto_ml_doa_opt}$ 的单阵列情形，多阵列的随机最大似然只是将代价函数在阵列索引 $\ell$ 上累加，同时通过映射算子 $\mathcal T_\ell(\,\cdot\,)$ 将所有阵列的导向矩阵统一为同一组参考角度参数 $\bm{\Theta}_\mathrm{ref}$，从而保证最终的参数估计结果只有一组角度。

\subsubsection{集中似然表达}

严格而言，多阵列随机模型中各阵列共享同一信号协方差矩阵 $\mathbf R_x$。为了得到与单阵列情形形式一致、便于实现的集中似然形式，本文采用一种常见的近似处理：对每个阵列分别应用单阵列的集中化步骤，得到局部的 $\widehat{\mathbf R}_{x,\ell}(\bm{\Theta}_\mathrm{ref})$ 与 $\widehat{\sigma}_{n,\ell}^2(\bm{\Theta}_\mathrm{ref})$，再在 $\ell$ 维度上进行联合搜索。


与单阵列情形类似，我们只关心来向参数 $\bm{\Theta}_\mathrm{ref}$，因此希望将对数似然函数中的其他参量（包括信号协方差 $\mathbf R_x$ 以及噪声方差 $\{\sigma_{n,\ell}^2\}$）用 $\bm{\Theta}_\mathrm{ref}$ 和观测 $\{\mathbf{y}_{\ell,t}\}$ 来表征。对 \eqref{eq:multi_arr_sto_ml_doa_opt} 中 $\mathbf R_x$ 和 $\{\sigma_{n,\ell}^2\}$ 分别求偏导并令其为零，可以得到与单阵列情形类似的极大似然估计分别为
%
\begin{equation}
  \widehat{\mathbf R}_{x,\ell}(\bm{\Theta}_\mathrm{ref}) = (\mathbf{A}_\ell^\Herm  \mathbf{A}_\ell)^{-1} \mathbf{A}_\ell^\Herm  \hat{\mathbf{R}}_\ell \mathbf{A}_\ell (\mathbf{A}_\ell^\Herm  \mathbf{A}_\ell)^{-1} - \widehat{\sigma}_{n,\ell}^2(\bm{\Theta}_\mathrm{ref}) (\mathbf{A}_\ell^\Herm \mathbf{A}_\ell)^{-1},
\end{equation}
%
\begin{equation}
  \widehat{\sigma}_{n,\ell}^2(\bm{\Theta}_\mathrm{ref}) = \frac{1}{M_\ell-K} \operatorname{tr}\left[ \left(\mathbf{I}_{M_\ell} - \mathbf{A}_\ell(\mathbf{A}_\ell^\Herm  \mathbf{A}_\ell)^{-1} \mathbf{A}_\ell^\Herm  \right) \hat{\mathbf{R}}_\ell \right],
\end{equation}
%
将这些估计回代到多阵列对数似然函数中，$\displaystyle\sum_{\ell}\tr(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell)$ 在最优点处变为常数项，可被消去，从而得到关于 $\bm{\Theta}_\mathrm{ref}$ 的集中随机最大似然估计
%
\begin{equation}
  \bm{\Theta}_{\mathrm{CSML}}
  = \argmin_{\bm{\Theta}_\mathrm{ref}}
  \sum_{\ell=0}^{L-1}
  \log\det\!\left[
  \mathbf{A}_\ell\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  \widehat{\mathbf R}_{x,\ell}(\bm{\Theta}_\mathrm{ref})
  \mathbf{A}_\ell^\Herm\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  + \widehat{\sigma}_{n,\ell}^2(\bm{\Theta}_\mathrm{ref})\mathbf{I}_{M_\ell}
  \right].
\end{equation}
%
换言之，多阵列的 \gls{csml} 估计可以通过最小化联合代价函数
%
\begin{equation}
  \mathcal{J}_\mathrm{joint}(\bm{\Theta}_\mathrm{ref})
  \triangleq
  \sum_{\ell=0}^{L-1}
  \log\det\!\left[
  \mathbf{A}_\ell\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  \widehat{\mathbf R}_{x,\ell}(\bm{\Theta}_\mathrm{ref})
  \mathbf{A}_\ell^\Herm\left(\mathcal T_\ell(\bm{\Theta}_\mathrm{ref})\right)
  + \widehat{\sigma}_{n,\ell}^2(\bm{\Theta}_\mathrm{ref})\mathbf{I}_{M_\ell}
  \right]
\end{equation}
%
而得到，其中 $\bm{\Theta}_\mathrm{ref}$ 是唯一的待估角度参数向量，其通过已知的角度映射算子 $\{\mathcal T_\ell(\,\cdot\,)\}$ 统一约束了所有阵列的局部来向，从而自然实现了多星观测在同一角度域下的联合最大似然融合。

\subsection{仿真结果}

下面基于多阵列角度映射进行仿真。算法的流程是：
%
\begin{enumerate}
  \item 选择其中一个阵列为参考阵列；
  \item 从参考阵列的视角选择初始角度 $\bm{\theta}_\mathrm{ref}^{(0)}$，计算得到代价函数；
  \item 从参考阵列的 \gls{ecef} 坐标 $\mathbf{p}_\mathrm{ref}$ 出发，沿着 $\bm{\theta}_\mathrm{ref}^{(0)}$ 方向做射线，求得射线和地表的交点（也即假设用户位于地表）；
  \item 由交点和第 $\ell$ 个阵列的 \gls{ecef} 坐标 $\mathbf{p}_g$ 和 $\mathbf{p}_\ell$，求得第 $\ell$ 个阵列视角下的信号来向 $\bm{\theta}_\ell^{(0)}$，计算得到代价函数；
  \item 将所有阵列的代价函数相加，得到总的代价函数；
  \item 遍历可行角度集，找到使代价函数最小的角度，即为目标角度。
\end{enumerate}

图 \ref{fig:doaEcefDualSatSingleSourceUpa20000Snap} 给出了基于若干算法得到的估计结果，分别是：
%
\begin{enumerate}
  \item 单阵列 \gls{music}，仅进行角度估计，作为估计性能的基线；
  \item 多阵列平均 \gls{music}，两个阵列独立进行 \gls{doa} 估计，第二个阵列估计得到的角度结合阵列坐标得到和地表的交点，再变换到第一个阵列的视角下，两个角度取平均作为估计结果；
  \item 单阵列 \gls{mle}：按照上述算法流程，但是仅包含一颗卫星，同样会进行角度 $\to$ 地表坐标 $\to$ 角度的变换；
  \item 多阵列平均 \gls{mle}：两颗卫星分别按照算法流程进行处理，第二个阵列的估计结果类似多阵列 \gls{music} 统一到第一个阵列的视角下；
  \item 多阵列联合 \gls{mle}：按照上述算法流程进行处理，包含两颗卫星。
\end{enumerate}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.8\textwidth]{fig/doaEcefDualSatSingleSourceUpa20000Snap.pdf}
  \caption{多阵列随机 DoA 支线：ECEF 坐标系下单/双星 UPA、单信源场景中 MUSIC 与 MLE 的 DoA MSE（20000 快拍）}
  \label{fig:doaEcefDualSatSingleSourceUpa20000Snap}
\end{figure}

从估计结果来看，\gls{mle} 的估计结果不理想，初步怀疑求射线和地表交点这一步会引入较大的非线性误差，下一步的仿真考虑遍历地面坐标，只进行坐标 $\to$ 角度的映射，看能否改善估计误差。

\subsection{地面坐标参数化的多阵列随机 DoA 模型}

假设信源服从圆对称复高斯分布 $\mathbf{x}_t \sim \mathcal{CN} \left( 0,\mathbf R_x \right)$，为简化讨论，假设各信源相互不相关，则 $\mathbf R_x$ 为对角矩阵，且协方差矩阵非奇异 $\det \left( \mathbf R_x \right) \ne 0$，时域不相关
%
\begin{equation}
  \begin{cases}
    \mathbb{E} \left[\mathbf{x}_t \mathbf{x}_s^\Herm \right] = \mathbf R_x \,\delta(t-s),\\[1mm]
    \mathbb{E} \left[ \mathbf{x}_t \mathbf{x}_s^\Tran \right] = 0
  \end{cases}, \quad\forall t,s.
\end{equation}
%
噪声满足 $\mathbf{n}_{\ell,t} \sim \mathcal{CN}(0,\sigma_{n,\ell}^2 \mathbf{I}_{M_\ell})$，且与 $\mathbf{x}_t$ 不相关。
在给定 $(\bm{\Xi},\mathbf R_x,\{\sigma_{n,\ell}^2\})$ 的条件下，$\{\mathbf{y}_{\ell,t}\}$ 为相互独立的圆对称复高斯向量，也即 $\mathbf{y}_{\ell,t} \sim \mathcal{CN}(\mathbf{0},\mathbf{R}_\ell)$，联合似然函数为
%
\begin{align}
  L\left(\{\mathbf{R}_\ell\},\bm{\Xi}\middle|\{\mathbf{y}_{\ell,t}\}\right)
  &= \prod_{\ell=0}^{L-1} \prod_{t=0}^{T-1}
  \frac{1}{\pi^{M_\ell}\det(\mathbf{R}_\ell)}
  \exp\!\left(-\mathbf{y}_{\ell,t}^\Herm \mathbf{R}_\ell^{-1}\mathbf{y}_{\ell,t}\right)
  \nonumber\\
  &= \prod_{\ell=0}^{L-1}
  \pi^{-M_\ell T} \left[\det(\mathbf{R}_\ell)\right]^{-T}
  \exp\!\left(
  - T\,\tr\left(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell\right)
  \right),
\end{align}
%
其中第 $\ell$ 个阵列的样本协方差矩阵为
%
\begin{equation}
  \hat{\mathbf{R}}_\ell
  \triangleq \frac{1}{T} \sum_{t=0}^{T-1} \mathbf{y}_{\ell,t}\mathbf{y}_{\ell,t}^\Herm,
\end{equation}
%
对应的对数似然函数为
%
\begin{equation}
  \ln \mathcal{L}\left(\{\mathbf{R}_\ell\},\bm{\Xi}\right)
  = \mathrm{Const}
  - T \sum_{\ell=0}^{L-1} \log\det(\mathbf{R}_\ell)
  - T \sum_{\ell=0}^{L-1} \tr\left(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell\right).
\end{equation}

\subsubsection{原始问题建模}

为建立多阵列协同估计中的统一参数化，本文以地面目标的空间位置作为全局未知量。
设第 $k$ 个地面目标在 \gls{ecef} 坐标系下的位置为
$\mathbf{p}_{g,k} \in \mathbb{R}^3$，将所有目标位置堆叠为
%
\begin{align}
  \bm{P}_g
  \triangleq
  \bigl[\mathbf{p}_{g,0},\ldots,\mathbf{p}_{g,K-1}\bigr]
  \in \mathbb{R}^{3\times K},
\end{align}
%
并将其矢量化后记为
%
\begin{align}
  \bm{\Xi} \triangleq \vect(\bm{P}_g) \in \mathbb{R}^{3K},
\end{align}
%
作为多阵列随机信号模型中的统一未知参数。

给定卫星阵列的 \gls{ecef} 坐标 $\{\mathbf{p}_\ell\}$ 与姿态矩阵 $\{\boldsymbol{\Omega}_\ell\}$，地面目标 $\mathbf{p}_{g,k}$ 相对于第 $\ell$ 个阵列的视线方向由
%
\begin{align}
  \mathbf{u}_{\ell,k}^{\mathrm{global}}(\bm{P}_g)
  &= \frac{\mathbf{p}_{g,k} - \mathbf{p}_\ell}
          {\|\mathbf{p}_{g,k} - \mathbf{p}_\ell\|},\\
  \mathbf{u}_{\ell,k}^{\mathrm{local}}(\bm{P}_g)
  &= \boldsymbol{\Omega}_\ell^\Tran \mathbf{u}_{\ell,k}^{\mathrm{global}}(\bm{P}_g),
\end{align}
%
从而得到第 $\ell$ 个阵列视角下第 $k$ 个目标的局部方位角与俯仰角
%
\begin{subequations}
\begin{align}
  \alpha_{\ell,k}^{\mathrm{loc}}(\bm{P}_g)
  &= \operatorname{atan2}\!\left(
     [\mathbf{u}_{\ell,k}^{\mathrm{local}}(\bm{P}_g)]_y,\,
     [\mathbf{u}_{\ell,k}^{\mathrm{local}}(\bm{P}_g)]_x
     \right),\\
  \beta_{\ell,k}^{\mathrm{loc}}(\bm{P}_g)
  &= \arcsin\!\left(
     [\mathbf{u}_{\ell,k}^{\mathrm{local}}(\bm{P}_g)]_z
     \right).
\end{align}
\end{subequations}
%
记第 $\ell$ 个阵列的局部来向集合为
%
\begin{align}
  \bm{\Theta}_\ell(\bm{P}_g)
  \triangleq
  \bigl[
   \bm{\theta}_{\ell,0}(\bm{P}_g),\ldots,
   \bm{\theta}_{\ell,K-1}(\bm{P}_g)
  \bigr],\quad
  \bm{\theta}_{\ell,k}(\bm{P}_g)
  \triangleq
  \bigl(\alpha_{\ell,k}^{\mathrm{loc}}(\bm{P}_g),\beta_{\ell,k}^{\mathrm{loc}}(\bm{P}_g)\bigr),
\end{align}
%
则第 $\ell$ 个阵列的导向矩阵可写为
%
\begin{align}
  \mathbf{A}_\ell(\bm{P}_g)
  \triangleq
  \mathbf{A}_\ell\bigl(\bm{\Theta}_\ell(\bm{P}_g)\bigr).
\end{align}
%
在此记号下，第 $\ell$ 个阵列的信号协方差矩阵为
%
\begin{align}
  \mathbf{R}_\ell(\bm{P}_g,\mathbf R_x,\sigma_{n,\ell}^2)
  =
  \mathbf{A}_\ell(\bm{P}_g)\,
  \mathbf R_x\,
  \mathbf{A}_\ell^\Herm(\bm{P}_g)
  + \sigma_{n,\ell}^2 \mathbf{I}_{M_\ell}.
\end{align}
%
因此，多阵列的随机最大似然/无约束最大似然优化问题可表述为
%
\begin{alignat}{2}
  \label{eq:multi_arr_sto_ml_geo_opt}
  &\underset{\bm{P}_g,\,\mathbf R_x,\,\{\sigma_{n,\ell}^2\}}{\mathrm{minimize}} \quad
  && \sum_{\ell=0}^{L-1}
  \left[
  \log\det\bigl(\mathbf{R}_\ell(\bm{P}_g,\mathbf R_x,\sigma_{n,\ell}^2)\bigr)
  + \tr\Bigl(\mathbf{R}_\ell^{-1}(\bm{P}_g,\mathbf R_x,\sigma_{n,\ell}^2)\,\hat{\mathbf{R}}_\ell\Bigr)
  \right] \nonumber\\
  &\text{subject to} \quad
  && \bm{P}_g \in \mathcal{G},\quad
     \mathbf R_x \succeq \mathbf{0},\quad
     \sigma_{n,\ell}^2 > 0,\ \ell=0,\ldots,L-1,
\end{alignat}
%
其中 $\mathcal{G}$ 表示地面目标位置的可行集合（例如 WGS-84 椭球表面或给定覆盖区域）。
可以看到，相比 \eqref{eq:sto_ml_doa_opt} 的单阵列情形，多阵列的随机最大似然只是将代价函数在阵列索引 $\ell$ 上累加，同时通过几何关系
$\bm{P}_g \mapsto \bm{\Theta}_\ell(\bm{P}_g)$
将所有阵列的导向矩阵统一为同一组地面目标坐标参数，从而保证最终的参数估计结果在几何上对应唯一的一组地面位置。

\subsubsection{集中似然表达}

严格而言，多阵列随机模型中各阵列共享同一信号协方差矩阵 $\mathbf R_x$。
为了得到与单阵列情形形式一致、便于实现的集中似然形式，本文采用一种常见的近似处理：对每个阵列分别应用单阵列的集中化步骤，得到局部的
$\widehat{\mathbf R}_{x,\ell}(\bm{P}_g)$ 与
$\widehat{\sigma}_{n,\ell}^2(\bm{P}_g)$，再在 $\ell$ 维度上进行联合搜索。

与单阵列情形类似，我们只关心地面目标坐标参数 $\bm{P}_g$，因此希望将对数似然函数中的其他参量（包括信号协方差 $\mathbf R_x$ 以及噪声方差 $\{\sigma_{n,\ell}^2\}$）用 $\bm{P}_g$ 和观测 $\{\mathbf{y}_{\ell,t}\}$ 来表征。
在上述近似框架下，对每个阵列的局部似然函数分别对 $\mathbf R_x$ 和 $\sigma_{n,\ell}^2$ 求偏导并令其为零，可得到单阵列形式的闭式极大似然估计
%
\begin{equation}
  \widehat{\mathbf R}_{x,\ell}(\bm{P}_g)
  = \bigl(\mathbf{A}_\ell^\Herm(\bm{P}_g)\,\mathbf{A}_\ell(\bm{P}_g)\bigr)^{-1}
    \mathbf{A}_\ell^\Herm(\bm{P}_g)\,\hat{\mathbf{R}}_\ell\,
    \mathbf{A}_\ell(\bm{P}_g)\,
    \bigl(\mathbf{A}_\ell^\Herm(\bm{P}_g)\,\mathbf{A}_\ell(\bm{P}_g)\bigr)^{-1}
    - \widehat{\sigma}_{n,\ell}^2(\bm{P}_g)\,
      \bigl(\mathbf{A}_\ell^\Herm(\bm{P}_g)\,\mathbf{A}_\ell(\bm{P}_g)\bigr)^{-1},
\end{equation}
\begin{equation}
  \widehat{\sigma}_{n,\ell}^2(\bm{P}_g)
  = \frac{1}{M_\ell-K} \operatorname{tr}\!\left[
      \Bigl(\mathbf{I}_{M_\ell}
      - \mathbf{A}_\ell(\bm{P}_g)
        \bigl(\mathbf{A}_\ell^\Herm(\bm{P}_g)\,\mathbf{A}_\ell(\bm{P}_g)\bigr)^{-1}
        \mathbf{A}_\ell^\Herm(\bm{P}_g)
      \Bigr)\,
      \hat{\mathbf{R}}_\ell
    \right].
\end{equation}
%
将这些局部极大似然估计回代到多阵列对数似然函数中，
$\sum_{\ell}\tr(\mathbf{R}_\ell^{-1}\hat{\mathbf{R}}_\ell)$
在最优点处变为与 $\bm{P}_g$ 无关的常数项，可被消去，从而得到关于地面目标坐标参数的集中随机最大似然估计
%
\begin{equation}
  \bm{P}_{g,\mathrm{CSML}}
  = \argmin_{\bm{P}_g}
  \sum_{\ell=0}^{L-1}
  \log\det\!\left[
  \mathbf{A}_\ell(\bm{P}_g)\,
  \widehat{\mathbf R}_{x,\ell}(\bm{P}_g)\,
  \mathbf{A}_\ell^\Herm(\bm{P}_g)
  + \widehat{\sigma}_{n,\ell}^2(\bm{P}_g)\,\mathbf{I}_{M_\ell}
  \right].
\end{equation}
%
换言之，多阵列的 \gls{csml} 估计可以通过最小化联合代价函数
%
\begin{equation}
  \mathcal{J}_\mathrm{joint}(\bm{P}_g)
  \triangleq
  \sum_{\ell=0}^{L-1}
  \log\det\!\left[
  \mathbf{A}_\ell(\bm{P}_g)\,
  \widehat{\mathbf R}_{x,\ell}(\bm{P}_g)\,
  \mathbf{A}_\ell^\Herm(\bm{P}_g)
  + \widehat{\sigma}_{n,\ell}^2(\bm{P}_g)\,\mathbf{I}_{M_\ell}
  \right]
\end{equation}
%
而得到，其中 $\bm{P}_g$ 是唯一的待估地面坐标参数，其通过几何映射
$\bm{P}_g \mapsto \{\bm{\Theta}_\ell(\bm{P}_g)\}$ 统一约束了所有阵列的局部来向，从而自然实现了基于地面位置参数的多星观测联合最大似然融合。

本节的关键作用是展示“多个局部阵列观测由同一个全局几何参数约束”的思想。后续多星静态 DoA--Doppler 主线会保留这一共享全局参数的结构，但把待估目标改写为 $\boldsymbol\eta_{\mathrm{stat}}=[\boldsymbol\xi^\Tran,\nu_{\mathrm{ref}}]^\Tran$，并显式利用各星链路 Doppler 的参考星参数化。