josh/doa_doppler_est
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

struct(sec): 正文文件重命名为 sec1/2/3/6，主线外的归档

Browse Source 
按 agent_instruction.md 阶段 4 的命名约定：
- detDoaEst → sec1_singleArrayDoa
- detDoaDopplerEst → sec2_singleArrayDoaDoppler
- multiArrayStatDetDoaDopplerEst → sec3_multiSatSingleUser
- multiArrayDynDetDoaDopplerEst → sec6_multiFrameCp

不进入主线的内容归档到 archive/：
- stoDoaEst（空文件）、stoDoaDopplerEst（仅注释）
- multiArrayStoDoaEst（随机信号模型不进入主线）
- multiArrayDetDoaDopplerDelayEst（与新论文方向不符）
- 98-future work、99-literature（前言已并入"后续方向"段）

主文件 \input 顺序对齐新命名，并为 sec4/sec5/sec7 留出注释占位。
abstract 同步更新以反映学习笔记定位与新章节结构。
本提交仅做结构搬迁，内容改写见后续 content commit。
This commit is contained in:
josh
2026-06-14 21:46:57 +08:00
parent 1be1379e53
commit 5d9a652b77
11 changed files with 11 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
archive/98-future work.tex
+82
View File
@@ -0,0 +1,82 @@
\section{后续研究方向}
在当前多星协作 DoA-Doppler 联合估计工作的基础上，后续可以从以下几个方向进一步扩展。
\subsection{扩展到定位场景：面向 DoA-Doppler-时延联合估计}
在当前多星协作 DoA-Doppler 联合估计工作的基础上，后续可进一步引入传播时延参数，研究多星场景下的 DoA-Doppler-时延联合估计问题，建立多星联合观测模型，对三类参数进行统一估计，并利用它们与用户位置、速度之间的几何关系，进一步实现定位或测速。
主要目标包括：
\begin{itemize}
\item 基于 MLE 的高精度联合估计；
\item 基于 空时二维/三维谱搜索的子空间方法（如空时 MUSIC）；
\item 面向高维参数空间的低复杂度迭代算法或交替优化算法；
\item 相应的 CRB 推导，用于刻画 DoA-Doppler-时延联合估计及定位的性能下界。
\end{itemize}
但是引入时延后，待估参数维度明显上升，计算复杂度会显著增加，会给算法设计带来更大挑战。
\subsection{考虑星间传输受限：分布式处理与低比特量化}
当前模型通常默认多星接收数据能够无损、实时地汇聚到中心节点进行全相干处理，这一假设在实际卫星系统中偏理想化。由于星间链路带宽、时延和计算资源受限，后续可研究在受限传输条件下的联合估计问题。
可以考虑分布式处理，也即各卫星本地完成部分信号处理，仅上传局部统计量、局部似然信息或低维特征，在中心节点进行融合；或者低比特量化，也即引入 1-bit ADC、混合精度 ADC 或压缩采样机制，以减少需要传输的原始观测数据量。
主要目标包括：
\begin{itemize}
\item 传输开销与估计性能之间的权衡关系；
\item 分布式/量化条件下 MLE、MUSIC 或迭代融合算法的设计；
\item 相应的量化 CRB、分布式 CRB 或性能损失分析。
\end{itemize}
\subsection{考虑模型失配与鲁棒联合估计}
当前分析通常建立在理想信号模型之上，但实际系统中往往存在模型失配。后续可以研究在模型不准确条件下的鲁棒 DoA-Doppler 联合估计。
可考虑的失配因素包括：
\begin{itemize}
\item 阵列流形失配，如阵元位置偏差、幅相误差、互耦等；
\item 卫星平台相关误差，如姿态误差、轨道参数误差、基线误差；
\item 时频同步误差，如星间时钟偏差、频偏；
\item 信道扰动，如相位抖动或短时不稳定效应。
\end{itemize}
这一方向通常可表述为带不确定集约束的优化问题，往往是非凸的，因此需要设计近似求解或迭代优化算法，并通过仿真验证鲁棒性提升。鲁棒算法适合形成篇幅较短但问题明确的论文，适合发 letter。
\subsection{面向超大规模阵列的随机矩阵理论分析}
当卫星端阵列规模、观测快拍数或协作卫星数进一步增大时，传统有限维分析可能难以准确刻画系统行为。此时可以引入随机矩阵理论，研究大系统条件下的渐近性能。
该方向可以分析的性能主要包括：
\begin{itemize}
\item \textbf{参数估计性能界}：例如大规模阵列下 DoA-Doppler 联合估计的渐近 CRB 或 Fisher 信息矩阵的确定性等价；
\item \textbf{子空间方法性能}：例如 MUSIC 谱峰分辨能力、噪声子空间泄漏、信号子空间与噪声子空间分离条件；
\item \textbf{协方差估计误差影响}：研究有限快拍条件下样本协方差偏差对参数估计的影响；
\item \textbf{可分辨性与相变现象}：分析两个相近信号在大系统条件下何时可分辨、何时出现分辨失效；
\item \textbf{量化或压缩条件下的性能极限}：例如 1-bit 量化、大规模阵列与有限快拍共同作用下的性能退化规律。
\end{itemize}
随机矩阵理论的优点是不仅可以解释仿真现象，还可以给出简洁的渐近表达式，用于揭示阵列规模、快拍数、信噪比和协作卫星数之间的基本规律；但对问题建模要求较高，需要首先明确所考虑的大系统极限形式。
\subsection{其他可进一步考虑的扩展方向}
除上述方向外，还可考虑以下扩展：
\begin{itemize}
\item \textbf{多用户场景}：由单用户扩展到多用户 DoA-Doppler 联合估计，并研究参数配对与用户分离问题；
\item \textbf{动态跟踪场景}：由静态估计扩展到时序跟踪，结合状态空间模型研究 DoA-Doppler 连续跟踪；
\item \textbf{宽带或近场模型}：针对超大孔径阵列或宽带信号，研究窄带远场假设失效时的联合估计问题；
\item \textbf{异步协作场景}：考虑多星间并非严格相干同步，而存在未知时延、相位差和频率偏差时的联合估计。
\end{itemize}
\subsection{与物理层安全的结合}
在上述方向基础上，需要进一步引入物理层安全（Physical Layer Security, PLS）的概念，使参数估计与安全传输形成结合。可能的切入方式包括：
\begin{itemize}
\item \textbf{基于参数估计的安全波束赋形}：利用多星估计得到的用户方向或位置参数，设计面向合法用户的协作波束，并抑制潜在窃听方向；
\item \textbf{鲁棒安全设计}：在 DoA、Doppler 或位置存在估计误差时，联合优化波束赋形与人工噪声，提高系统的保密性能；
\item \textbf{基于物理特征的认证}：利用多星观测到的 DoA-Doppler-Delay 组合作为物理层特征，实现合法用户识别与伪装检测；
\item \textbf{受限传输下的隐私保护融合}：在分布式估计或压缩传输场景中，同时考虑“估计性能--传输开销--信息泄露”之间的权衡。
\end{itemize}
archive/99-literature.tex
+37
View File
@@ -0,0 +1,37 @@
\section{初步文献调研}
\subsection{多星全相干合并建模分析}
在理想条件下，本文采用的多星协同 \gls{doa} 估计并不是在观测数据层面直接将多颗卫星阵列拼接成一个“超级阵列”，而是基于单星随机信号模型，通过共享一组全局角度参数，在似然函数中实现相干合并。多星最大似然估计的对数似然可写为各单星代价函数的加和，其形式与单星随机信号 \gls{mle} 基本一致，只是将本地角度替换为 $\mathcal{P}_\ell(\boldsymbol{\Theta}_\mathrm{ref})$。在实现上，可以先对每颗卫星独立构造协方差矩阵与本地阵列流形，再在全局角度域中将这些单星代价相干叠加以完成联合估计。
从统计结构来看，该方法相当于在单星 \gls{mle} 之外增加“全局到本地角度映射”与“多星求和”两个步骤，其核心估计机制仍与传统随机信号 \gls{mle} 保持一致。因此，在完全同步与完美标定的假设下，多星全相干融合的主要工作更多体现在几何关系的建模与联合性能分析，而非估计算法框架的改变。
这一现象与部分校准阵列和分布式子阵列的既有研究结果相吻合~\cite{pesaventoDirectionFindingPartly2002, seeDirectionofArrivalEstimationPartly2004, zhangDirectionFindingPartly2023, suDirectionFindingPartly2025}。当子阵间的相对相位、阵元位置及姿态信息均已校准时，多子阵协同 \gls{doa} 估计通常可以等价为单一大阵列的估计问题。本文的全相干建模与此类方法在本质上是一致的，因此整体表达形式也较为直接。
从性能角度来看，对于两颗结构相同且噪声功率一致的卫星阵列，其联合 Fisher 信息近似为单星的两倍，从而使单参数的 \gls{crb} 约降低一半，等效于获得约 $3$\,dB 的理论性能提升。此外，当多颗卫星处于不同轨道位置时，等效阵列孔径进一步增大，可在某些几何结构下获得额外的精度提升。因此，全相干模型在后续工作中可作为多星协同性能的理论上界。
\subsection{分布式阵列相关文献的调研结果}
虽然全相干模型提供了一个便于分析的理论基线，但在实际星座系统中完全满足相关假设并不容易。星间时钟偏移、载频漂移、姿态扰动和链路延迟等因素都会破坏严格相干性，使得直接采用全相干建模具有局限性。因此，有必要参考分布式阵列的相关研究，探讨在不完美条件下如何有效利用多星观测信息。
现阶段，与卫星通信直接结合的 \gls{doa} 估计研究仍然有限，更多的成果来自阵列信号处理和分布式子阵网络的研究体系。以下内容给出若干具有代表性的方向。
\subsubsection{基于部分校准阵列和弱同步模型的多星 DoA}
多星系统可视为时变的分布式阵列网络。现有关于部分标定阵列和分布式子阵列的研究~\cite{seeDirectionofArrivalEstimationPartly2004, zhangDOAEstimationExploiting2024, chowdhuryDirectionofArrivalEstimationExploiting2021, wangEnhancedDOAEstimation2024} 表明，当子阵间存在未知相位漂移或阵元位置误差时，通过对子阵内部与子阵间的约束结构进行建模，仍可在一定程度上恢复虚拟全孔径的角度估计能力。
基于这一思路，在卫星场景中可以构建“星内完全标定、星间弱同步”的模型，并进一步研究在星间信息交换受限的条件下如何实现稳定且性能接近虚拟全孔径的联合估计。
\subsubsection{基于角域稀疏性的 DoA 估计}
在密集 \gls{leo} 星座下，卫星可覆盖范围极大，但实际地面来向往往只有有限数量，因此角域呈现明显稀疏性。压缩感知与稀疏表示~\cite{chowdhuryDirectionofArrivalEstimationExploiting2021, liHighResolutionDOA2025} 提供了一类适合这类结构的 \gls{doa} 估计方法，尤其适用于快拍数有限、信噪比较低或仅有少量卫星参与协作的情况。
此外，多星在轨道运行中具有不断变化的相对几何结构，这些时变特征可进一步作为额外约束，引入到“多星多时刻”场景的时–空联合稀疏估计框架中。
\subsubsection{DoA、TDOA 与 FDOA 的联合参数估计}
现有卫星定位研究~\cite{haoInterferenceGeolocationSatellite2021, elgamoudiSurveyRecentTechniques2021} 主要依赖 \gls{tdoa} 和 \gls{fdoa} 测量，而基于 \gls{doa} 的方法通常仅在单星场景有所涉及。多星场景使三类观测量可以在同一几何框架下结合，从而获得更稳健的联合估计能力。
构建联合似然函数并推导联合 \gls{crb} 能够揭示不同观测量之间的互补关系。例如，\gls{doa} 估计可以降低对高精度时钟同步的依赖，而 \gls{tdoa}/\gls{fdoa} 提供了对目标位置与运动状态更敏感的测量；几何结构随时间变化也能为多参量融合带来额外增益。
在星间链路受限时，还可考虑分布式求解策略，即各卫星仅交换有限的中间信息，通过迭代的方式逼近集中式最优解。典型方法包括基于“共识（consensus）”思想的分布式最大似然方法~\cite{chenDistributedAlgorithmsArray2021}：各卫星先利用本地观测进行更新，再通过少量信息交换使整个系统在迭代中逐步趋于一致。这类方法通信开销可控，适用于未来大规模密集星座的协作定位任务。
archive/multiArrayDetDoaDopplerDelayEst.tex
+199
View File
@@ -0,0 +1,199 @@
\section{基于参考星链路参数化的多阵列确定信号 DoA--Doppler--时延估计扩展}
\label{sec:multiArrayDetDoaDopplerDelayEst}
本节在静态多星 DoA--Doppler 主线的基础上进一步加入参考星链路时延 $\tau_{\mathrm{ref}}$。因此，本节不是重新定义一套几何模型，而是在第~\ref{subsec:multiArrayCommonNotation} 小节的全局方向 $\boldsymbol\xi$、局部导向矢量 $\mathbf a_\ell(\boldsymbol\xi)$ 和参考星 Doppler 参数化基础上，额外引入各星链路时延映射。
\subsection{信号模型}
与前述静态 DoA--Doppler 模型一致，固定选择参考卫星 $\ell_{\mathrm{ref}}$。在加入时延后，DoA--Doppler--delay 扩展模型的主参数统一写为
\begin{equation}
\boldsymbol\eta_{\tau}
=
\begin{bmatrix}
\boldsymbol\xi\\
\tau_{\mathrm{ref}}\\
\nu_{\mathrm{ref}}
\end{bmatrix},
\end{equation}
其中 $\tau_{\mathrm{ref}}$ 与 $\nu_{\mathrm{ref}}$ 分别为参考卫星链路上的时延和归一化 Doppler。
第 $\ell$ 颗卫星的链路时延与归一化 Doppler 分别由参考星参数化写为
\begin{equation}
\tau_\ell
=
f_{\tau,\ell}(\boldsymbol\xi,\tau_{\mathrm{ref}}),
\qquad
\nu_\ell
=
f_{\nu,\ell}(\boldsymbol\xi,\nu_{\mathrm{ref}}),
\qquad
\ell=0,1,\ldots,L-1,
\end{equation}
并满足参考星约束
\begin{equation}
f_{\tau,\ell_{\mathrm{ref}}}(\boldsymbol\xi,\tau_{\mathrm{ref}})
=
\tau_{\mathrm{ref}},
\qquad
f_{\nu,\ell_{\mathrm{ref}}}(\boldsymbol\xi,\nu_{\mathrm{ref}})
=
\nu_{\mathrm{ref}}.
\end{equation}
对 $\ell\neq\ell_{\mathrm{ref}}$，上述映射由已知卫星位置、速度、姿态以及参考星链路参数诱导得到。
考虑已知训练序列 $s(t)$。加入链路时延后，第 $\ell$ 颗卫星在第 $t$ 个采样点的接收向量写为
\begin{equation}
\mathbf y_\ell(t)
=
\zeta_\ell
\mathbf a_\ell(\boldsymbol\xi)
s\!\bigl(t-\tau_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}})\bigr)
e^{j\nu_\ell(\boldsymbol\xi,\nu_{\mathrm{ref}})t}
+
\mathbf n_\ell(t),
\qquad t=0,1,\ldots,T-1,
\end{equation}
当 $\tau_\ell$ 非整数时，$s(t-\tau_\ell)$ 表示连续基带训练波形经过分数时延后在第 $t$ 个采样时刻的取值，等价地可由插值或分数时延滤波实现；若仅考虑整数采样时延，则应额外约束 $\tau_\ell\in\mathbb Z$。
其中为书写紧凑记
\begin{equation}
\tau_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}})
\triangleq
f_{\tau,\ell}(\boldsymbol\xi,\tau_{\mathrm{ref}}),
\qquad
\nu_\ell(\boldsymbol\xi,\nu_{\mathrm{ref}})
\triangleq
f_{\nu,\ell}(\boldsymbol\xi,\nu_{\mathrm{ref}}).
\end{equation}
这里 $\zeta_\ell\in\mathbb C$ 为未知复干扰系数，噪声满足
\begin{equation}
\mathbf n_\ell(t)
\sim
\mathcal{CN}(\mathbf 0,\sigma_{n,\ell}^2\mathbf I_{M_\ell}).
\end{equation}
定义包含时延和 Doppler 的时间模板
\begin{equation}
\mathbf q_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
\triangleq
\begin{bmatrix}
s\!\bigl(0-\tau_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}})\bigr)e^{j\nu_\ell(\boldsymbol\xi,\nu_{\mathrm{ref}})\cdot0}\\
s\!\bigl(1-\tau_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}})\bigr)e^{j\nu_\ell(\boldsymbol\xi,\nu_{\mathrm{ref}})\cdot1}\\
\vdots\\
s\!\bigl(T-1-\tau_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}})\bigr)e^{j\nu_\ell(\boldsymbol\xi,\nu_{\mathrm{ref}})(T-1)}
\end{bmatrix}.
\end{equation}
相应地，矩阵观测模型为
\begin{equation}
\mathbf Y_\ell
=
\zeta_\ell
\mathbf a_\ell(\boldsymbol\xi)
\mathbf q_\ell^\Tran(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
+
\mathbf N_\ell.
\end{equation}
向量化后，定义
\begin{equation}
\tilde{\mathbf y}_\ell
\triangleq
\vect(\mathbf Y_\ell),
\qquad
\tilde{\mathbf n}_\ell
\triangleq
\vect(\mathbf N_\ell),
\qquad
\mathbf g_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
\triangleq
\mathbf q_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
\otimes
\mathbf a_\ell(\boldsymbol\xi),
\end{equation}
可得
\begin{equation}
\tilde{\mathbf y}_\ell
=
\zeta_\ell\mathbf g_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
+
\tilde{\mathbf n}_\ell.
\end{equation}
本节中的干扰参数可写为
\begin{equation}
\boldsymbol\psi_{\tau}
=
\Big[
\{\zeta_\ell\}_{\ell=0}^{L-1},
\{\sigma_{n,\ell}^2\}_{\ell=0}^{L-1}
\Big]^\Tran.
\end{equation}
\subsection{联合似然与集中式估计}
在星间噪声独立、且第 $\ell$ 颗卫星噪声方差为 $\sigma_{n,\ell}^2$ 的条件下，省略与参数无关的常数项后，联合负对数似然为
\begin{equation}
\mathcal L_{\tau}(\boldsymbol\eta_{\tau},\boldsymbol\psi_{\tau})
=
\sum_{\ell=0}^{L-1}
\left[
M_\ell T\log\sigma_{n,\ell}^2
+
\frac{1}{\sigma_{n,\ell}^2}
\left\|
\tilde{\mathbf y}_\ell
-
\zeta_\ell\mathbf g_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
\right\|_2^2
\right].
\end{equation}
对于固定的 $(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})$ 和 $\{\sigma_{n,\ell}^2\}$，$\zeta_\ell$ 的最小二乘解为
\begin{equation}
\hat\zeta_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
=
\frac{
\mathbf g_\ell^\Herm(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})\tilde{\mathbf y}_\ell
}{
\|\mathbf g_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})\|_2^2
}.
\end{equation}
代回后，第 $\ell$ 颗卫星的集中残差为
\begin{equation}
\mathcal E_{\tau,\ell}(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})
\triangleq
\|\tilde{\mathbf y}_\ell\|_2^2
-
\frac{
\left|
\mathbf g_\ell^\Herm(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})\tilde{\mathbf y}_\ell
\right|^2
}{
\|\mathbf g_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})\|_2^2
}.
\end{equation}
若进一步对噪声方差进行轮廓化，则在星间异方差假设下，DoA--Doppler--delay 扩展模型的集中估计问题为
\begin{equation}
\widehat{\boldsymbol\eta}_{\tau}
=
\arg\min_{\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}}}
\sum_{\ell=0}^{L-1}
M_\ell T
\log \mathcal E_{\tau,\ell}(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}}).
\end{equation}
若各星噪声方差已知，则也可采用加权集中似然形式
\begin{equation}
\widehat{\boldsymbol\eta}_{\tau}
=
\arg\max_{\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}}}
\sum_{\ell=0}^{L-1}
\frac{1}{\sigma_{n,\ell}^2}
\frac{
\left|
\mathbf g_\ell^\Herm(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})\tilde{\mathbf y}_\ell
\right|^2
}{
\|\mathbf g_\ell(\boldsymbol\xi,\tau_{\mathrm{ref}},\nu_{\mathrm{ref}})\|_2^2
}.
\end{equation}
该扩展模型与静态 DoA--Doppler 模型的区别在于：$\boldsymbol\eta_{\tau}$ 在 $\boldsymbol\eta_{\mathrm{stat}}$ 的基础上显式加入 $\tau_{\mathrm{ref}}$，并通过 $f_{\tau,\ell}(\boldsymbol\xi,\tau_{\mathrm{ref}})$ 将各星时延耦合在同一个参考星参数化框架中。因此，delay 扩展不会改变多星 DoA--Doppler 主线的符号体系，只是在同一层级上增加一类可观测链路参数。
因此，阅读顺序上应先理解静态 DoA--Doppler 和动态连续相位主线，再把本节视为在同一参考星参数化框架下增加时延观测的一种横向扩展。若后续继续推导该模型的 CRB，也应沿用 $\boldsymbol\eta_{\tau}$ 与干扰参数分组，而不要重新引入一套独立的几何符号。
archive/multiArrayStoDoaEst.tex
+619
View File
@@ -0,0 +1,619 @@
\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 的参考星参数化。
archive/stoDoaDopplerEst.tex
+5
View File
@@ -0,0 +1,5 @@
\section{随机信号 DoA-Doppler 估计}
\subsection{CRB 推导}
\subsection{随机/无条件最大似然}
archive/stoDoaEst.tex
View File