数学基础·

用概率密度函数（Probability Density Function, PDF）描述一组随机数据，即 p(x0,x1,⋯ ,xn−1;θ)p(x_0, x_1, \cdots, x_{n-1};\theta)，PDF 以未知量 θ\theta 为参数，以 n=1n=1，θ\theta 表示均值为例，数据的 PDF 为

p(x0;θ)=12πσ2exp⁡(−(x0−θ)22σ2)p(x_0;\theta) = \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp{\left( -\dfrac{(x_0-\theta)^2}{2\sigma^2} \right)}

然后就可以根据 x0x_0 的观测值推断 θ\theta 的值。

估计量性能评估·

将数据建模为

xn=A+wnx_n = A + w_n

其中 wnw_n 为加性高斯白噪声（Additive Gaussian White Noise, AWGN），即 wi∼N(0,σ2)w_i\sim \mathcal{N}(0, \sigma^2) 表示均值为 00，方差为 σ2\sigma^2 的高斯分布，并且所有样本是互不相关的。利用下式即数据的样本均值估计 AA

A^=1N∑n=0N−1xn\hat{A} = \dfrac{1}{N} \sum_{n = 0}^{N-1} x_n

满足

E(A^)=E(1N∑n=0N−1xn)=1N∑n=0N−1(xn)=Avar(A^)=var(1N∑n=0N−1xn)=1N2∑n=0N−1var(xn)=1N2Nσ2=σ2N\begin{aligned} E(\hat{A}) &= E\left(\dfrac{1}{N} \sum_{n = 0}^{N-1} x_n\right) = \dfrac{1}{N}\sum_{n = 0}^{N-1}(x_n) = A\\ \mathrm{var} (\hat{A}) &= \mathrm{var} \left(\dfrac{1}{N} \sum_{n = 0}^{N-1} x_n\right) = \dfrac{1}{N^2}\sum_{n = 0}^{N-1}\mathrm{var}(x_n) = \dfrac{1}{N^2} N\sigma^2 = \dfrac{\sigma^2}{N}\\ \end{aligned}

最小方差无偏估计·

无偏估计量·

无偏估计意味着估计量的平均值为未知参数 θ\theta 的真值，如果

E(θ^)=θa<θ<bE(\hat{\theta}) = \theta\qquad a < \theta < b

说明估计量是无偏的，其中 (a,b)(a, b) 表示 θ\theta 可能的取值范围。

对于同一个参数有多个估计可用的情况，即 {θ^0,θ^1,⋯ ,θ^n−1}\{\hat{\theta}_0, \hat{\theta}_1, \cdots, \hat{\theta}_{n-1}\}，对这些组合求平均，即

θ^=1n∑i=0n−1θ^i\hat{\theta} = \dfrac{1}{n}\sum_{i = 0}^{n-1} \hat{\theta}_i

假设每个估计量都是无偏的，方差相同且互不相关，则

E(θ^)=θ,var(θ^)=var(θ^0)nE(\hat{\theta}) = \theta, \mathrm{var}(\hat{\theta}) = \dfrac{\mathrm{var}(\hat{\theta}_0)}{n}

因此，求平均的估计值越多，方差越小，当 n→∞n \to \infty 时，θ^→θ\hat{\theta}\to \theta。

最小方差准则·

均方误差定义为

mse(θ^)=E[(θ^−θ)2]\mathrm{mse}(\hat{\theta}) = E [(\hat{\theta} - \theta)^2]

衡量估计值偏离真值的平方偏差的统计平均值。

mse(θ^)=E[(θ^−θ)2]=E{[(θ^−E(θ^))+(E(θ^)−θ)]2}\begin{aligned} \mathrm{mse}(\hat{\theta}) &= E [(\hat{\theta} - \theta)^2]\\ &= E\left\{ \left [ \left(\hat{\theta} - E(\hat{\theta}) \right) + \left( E(\hat{\theta}) - \theta \right) \right]^2 \right\}\\ \end{aligned}

其中第一部分是估计量围绕其数学期望的随机波动，第二部分是估计量的期望围绕真值的波动，展开后得到

mse(θ^)=var(θ^)+b2(θ)\mathrm{mse}{(\hat{\theta})} = \mathrm{var}(\hat{\theta}) + b^2(\theta)

其中 b(θ)=E(θ^)−θb(\theta) = E(\hat{\theta}) - \theta，表示估计量的偏差。上式表明 MSEMSE 是由 估计量的方差 和 偏差 引起的误差组成的。

对于有偏的估计量，MSEMSE 与参数 θ\theta 的真值有关，最小化 MSEMSE 可能导致不可实现的估计量，如果约束估计量为无偏的，然后最小化方差，得到的估计量就是最小方差无偏（Minimum Variance Unbiased, MVU）估计量。

最小方差无偏估计的存在性·

不总是存在一个估计量 θ^\hat{\theta}，对于所有的 θ\theta，其方差都小于其它无偏估计量！

Cramer-Rao 下限·

观测到单个样本

x0=A+w0,w0∼N(0,σ2)x_0 = A + w_0, \quad w_0\sim \mathcal{N}(0, \sigma^2)

无偏估计满足 A^=x0,var(A^)=σ2\hat{A}=x_0, \mathrm{var}(\hat{A})= \sigma^2.

考虑 PDF 的自然对数

ln⁡p(x0;A)=−ln⁡2πσ2−(x0−A)22σ2\ln p(x_0; A) = -\ln{\sqrt{2\pi\sigma^2}} - \dfrac{(x_0-A)^2}{2\sigma^2}

一阶导数为

∂ln⁡p(x0;A)∂θ=x0−Aσ2\dfrac{\partial{\ln p(x_0; A)}}{\partial{\theta}} = \dfrac{x_0-A}{\sigma^2}

负的二阶导数（对数似然函数的曲率）为

−∂2ln⁡p(x0;A)∂A2=1σ2-\dfrac{\partial^2{\ln p(x_0; A)}}{\partial{A}^2} = \dfrac{1}{\sigma^2}

随着 σ2\sigma^2 的减少而增加，并且已知 var(A^)=σ2\mathrm{var}(\hat{A}) = \sigma^2，则

var(A^)=1−∂2ln⁡p(x0;A)∂A2\mathrm{var}(\hat{A}) = \dfrac{1}{-\dfrac{\partial^2{\ln p(x_0; A)}}{\partial{A}^2}}

更一般的度量是

−E[∂2ln⁡p(x0;A)∂A2]-E\left [ \dfrac{\partial^2{\ln p(x_0; A)}}{\partial{A}^2} \right]

表示对数自然函数的平均曲率，值越大，表示估计量的方差越小。

直观理解：

PDF 的 xx 取固定值时，PDF 是参数 A 的（似然）函数，图像越尖锐，估计参数 A 的精度越高。用 负的二阶导数 定量描述 尖锐程度。

似然函数的对数曲线对参数越陡，表示数据对该参数越敏感，信息越多，可达的方差下界越小。

标量参数的 CRLB·

假设对于所有的参数 θ\theta，概率密度函数 p(x;θ)p(\mathbf{x};\theta) 满足正则条件

E[∂ln⁡p(x0;θ)∂θ]=0E\left [ \dfrac{\partial{\ln p(x_0;\theta)}}{\partial{\theta}} \right] = 0

那么，任何无偏估计量 θ^\hat{\theta} 的方差满足

var(θ^)≥1−E[∂2ln⁡p(x;θ)∂θ2]\mathrm{var}(\hat{\theta}) \geq \dfrac{1}{-E\left [\dfrac{\partial^2{\ln p(\mathbf{x};\theta)}}{\partial{\theta}^2}\right]}

当且仅当

∂ln⁡p(x0;θ)∂θ=I(θ)(g(x)−θ)\dfrac{\partial{\ln p(x_0;\theta)}}{\partial{\theta}} = I(\theta)(g(\mathbf{x}) - \theta)

时，对有所 θ\theta 达到下限的无偏估计量可以求得，估计量 θ^=g(x)\hat{\theta}= g(\mathbf{x}) 时 MVU 估计量，最小方差是 1I(θ)\dfrac{1}{I(\theta)}.

推导过程·

左边：

E[∂ln⁡p(x;θ)∂θ]=∫∂ln⁡p(x;θ)∂θp(x;θ)dx=∫∂p(x;θ)∂θdx\begin{aligned} E\left [ \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} \right] &= \int \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} p(x;\theta) \mathrm{d}x\\ &= \int \dfrac{\partial{p(x;\theta)}}{\partial{\theta}} \mathrm{d}x\\ \end{aligned}

右边：

0=∂1∂θ=∂∂θ∫p(x;θ)dx0 = \dfrac{\partial 1}{\partial \theta} = \dfrac{\partial}{\partial \theta} \int p(x;\theta) \mathrm{d}x

根据牛顿-莱布尼茨公式

F′(t)=ddt∫a(t)b(t)f(x,t)dx=f(b(t),t)b′(t)−f(a(t),t)a′(t)+∫a(t)b(t)∂∂tf′(x,t)dx\begin{aligned} F'(t) &= \dfrac{\mathrm{d}}{\mathrm{d}t} \int_{a(t)}^{b(t)} f(x, t) \mathrm{d}x\\ &= f(b(t), t)b'(t) - f(a(t), t) a'(t) + \int_{a(t)}^{b(t)} \dfrac{\partial}{\partial t} f'(x, t) \mathrm{d}x \end{aligned}

当 边界随参数变化带来的边界项 f(b(t),t)b′(t)−f(a(t),t)a′(t)=0f(b(t), t)b'(t) - f(a(t), t) a'(t)=0 时，求积分与求偏导运算可以交换。

假设正则条件满足，两式中的求积分和求偏导运算可以交换，说明 PDF 的非零边界和参数 θ\theta 无关。

下面推导标量参数 α=g(θ)\alpha = g(\theta) 的 CRLB，对于所有无偏估计量

E(α^)=∫α^p(x;θ)dx=α=g(θ)E(\hat{\alpha}) = \int \hat{\alpha} p(x;\theta)\mathrm{d}x = \alpha = g(\theta)

在满足正则条件的前提下，对等式两边求导，得到

∂g(θ)∂θ=∂∂θ∫α^p(x;θ)dx=∫α^∂p(x;θ)∂θdx=∫α^∂ln⁡p(x;θ)∂θp(x;θ)dx\dfrac{\partial g(\theta)}{\partial \theta} =\dfrac{\partial}{\partial \theta}\int \hat{\alpha} p(x;\theta)\mathrm{d}x =\int \hat{\alpha} \dfrac{\partial p(x;\theta)}{\partial \theta}\mathrm{d}x =\int \hat{\alpha} \dfrac{\partial \ln p(x;\theta)}{\partial \theta}p(x;\theta)\mathrm{d}x

根据正则条件 E[∂ln⁡p(x;θ)∂θ]=0E\left [ \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} \right] = 0，可以得到

∫α∂ln⁡p(x;θ)∂θp(x;θ)dx=αE[∂ln⁡p(x;θ)∂θ]=0\int \alpha \dfrac{\partial \ln p(x;\theta)}{\partial \theta}p(x;\theta)\mathrm{d}x = \alpha E\left [ \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} \right] = 0

两式作差得到

∫(α^−α)∂ln⁡p(x;θ)∂θp(x;θ)dx=∂g(θ)∂θ\int (\hat{\alpha} - \alpha) \dfrac{\partial \ln p(x;\theta)}{\partial \theta}p(x;\theta)\mathrm{d}x = \dfrac{\partial g(\theta)}{\partial \theta}

根据柯西-施瓦茨不等式，

(∂g(θ)∂θ)2=(∫(α^−α)∂ln⁡p(x;θ)∂θp(x;θ)dx)2≤∫(α^−α)2p(x;θ)dx⋅∫(∂ln⁡p(x;θ)∂θ)2p(x;θ)dx=var(α^)⋅E[(∂ln⁡p(x;θ)∂θ)2]var(α^)≥(∂g(θ)∂θ)2E[(∂ln⁡p(x;θ)∂θ)2]\begin{aligned} \left( \dfrac{\partial g(\theta)}{\partial \theta} \right)^2 &= \left( \int (\hat{\alpha} - \alpha) \dfrac{\partial \ln p(x;\theta)}{\partial \theta}p(x;\theta)\mathrm{d}x \right)^2\\ &\leq \int (\hat{\alpha} - \alpha)^2 p(x;\theta)\mathrm{d}x \cdot \int \left(\dfrac{\partial \ln p(x;\theta)}{\partial \theta}\right)^2p(x;\theta)\mathrm{d}x\\ &= \mathrm{var}(\hat{\alpha}) \cdot E\left [ \left(\dfrac{\partial \ln p(x;\theta)}{\partial \theta}\right)^2 \right]\\ \mathrm{var}(\hat{\alpha}) &\geq \dfrac{\left( \dfrac{\partial g(\theta)}{\partial \theta} \right)^2}{E\left [ \left(\dfrac{\partial \ln p(x;\theta)}{\partial \theta}\right)^2 \right]} \end{aligned}

由正则化条件，

E[∂ln⁡p(x;θ)∂θ]=0E\left [ \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} \right] = 0

得到

∫∂ln⁡p(x;θ)∂θp(x;θ)dx=0\int \dfrac{\partial \ln p(x;\theta)}{\partial \theta}p(x;\theta)\mathrm{d}x = 0

积分边界与参数 θ\theta 无关，于是

∂∂θ∫∂ln⁡p(x;θ)∂θp(x;θ)dx=0∫[∂2ln⁡p(x;θ)∂θ2p(x;θ)+∂ln⁡p(x;θ)∂θ∂p(x;θ)∂θ]dx=0∫[∂2ln⁡p(x;θ)∂θ2p(x;θ)+∂ln⁡p(x;θ)∂θ∂ln⁡p(x;θ)∂θp(x;θ)]dx=0−E[∂2ln⁡p(x;θ)∂θ2]=E[(∂ln⁡p(x;θ)∂θ)2]\begin{aligned} \dfrac{\partial}{\partial \theta}\int \dfrac{\partial \ln p(x;\theta)}{\partial \theta}p(x;\theta)\mathrm{d}x &= 0\\ \int \left[ \dfrac{\partial^2 \ln p(x;\theta)}{\partial \theta^2}p(x;\theta) + \dfrac{\partial \ln p(x;\theta)}{\partial \theta}\dfrac{\partial p(x;\theta)}{\partial \theta} \right]\mathrm{d}x &= 0\\ \int \left[ \dfrac{\partial^2 \ln p(x;\theta)}{\partial \theta^2}p(x;\theta) + \dfrac{\partial \ln p(x;\theta)}{\partial \theta}\dfrac{\partial \ln p(x;\theta)}{\partial \theta} p(x;\theta) \right]\mathrm{d}x &= 0\\ -E\left [ \dfrac{\partial^2 \ln p(x;\theta)}{\partial \theta^2} \right] &= E\left [ \left(\dfrac{\partial \ln p(x;\theta)}{\partial \theta}\right)^2 \right]\\ \end{aligned}

最终得到，

var(α^)≥(∂g(θ)∂θ)2−E[∂2ln⁡p(x;θ)∂θ2]\mathrm{var}(\hat{\alpha}) \geq \dfrac{\left( \dfrac{\partial g(\theta)}{\partial \theta} \right)^2}{-E\left [ \dfrac{\partial^2 \ln p(x;\theta)}{\partial \theta^2} \right]}

当且仅当无偏估计量 α^\hat{\alpha} 与对数似然函数的一阶偏导数呈线性关系，即

∂ln⁡p(x;θ)∂θ=1c(α^−α)\dfrac{\partial \ln p(x;\theta)}{\partial \theta} = \dfrac{1}{c}(\hat{\alpha} - \alpha)

时成立，其中 cc 与 xx 无关。

直观理解：

对数似然函数的一阶偏导数反映“秤”对真实值的敏感程度，取等条件表示估计值和 敏感程度 呈固定的比例，“秤”越敏感，估计的误差就按照这个比例调整，几步浪费精度，也不高估“秤”的能力，最终将误差压到最低值。

以 α=g(θ)=θ\alpha = g(\theta) = \theta 为例，达到 CRLB 时，

∂ln⁡p(x;θ)∂θ=1c(θ)(θ^−θ)∂2ln⁡p(x;θ)∂θ2=−1c(θ)+(θ^−θ)∂1c(θ)∂θ−E[∂2ln⁡p(x;θ)∂θ2]=1c(θ)\begin{aligned} \dfrac{\partial \ln p(x;\theta)}{\partial \theta} &= \dfrac{1}{c(\theta)}(\hat{\theta} - \theta)\\ \dfrac{\partial^2 \ln p(x;\theta)}{\partial \theta^2} &= -\dfrac{1}{c(\theta)} + (\hat{\theta}-\theta)\dfrac{\partial \frac{1}{c(\theta)}}{\partial \theta}\\ -E\left [ \dfrac{\partial^2 \ln p(x;\theta)}{\partial \theta^2} \right] &= \dfrac{1}{c(\theta)}\\ \end{aligned}

定义 Fisher 信息 I(θ)=−E[∂2ln⁡p(x;θ)∂θ2]I(\theta) = -E\left [ \dfrac{\partial^2 \ln p(x;\theta)}{\partial \theta^2} \right]，则

c(θ)=1I(θ)c(\theta) = \dfrac{1}{I(\theta)}

矢量参数的 CRLB·

现将前一部分的结果扩展到估计矢量参数 θ=[θ1θ2⋯θp]⊤\boldsymbol{\theta}=[\theta_1 \theta_2\cdots \theta_p]^{\top}，假定 θ^\hat{\boldsymbol{\theta}} 是无偏估计，矢量参数的 CRLB 允许对每隔元素的方差放置一个下限，即

var(θi^)≥[I−1(θ)]ii\mathrm{var}(\hat{\theta_i}) \geq [\mathbf{I}^{-1}(\boldsymbol{\theta})]_{ii}

其中 I(θ)\mathbf{I}(\boldsymbol{\theta}) 是 p×pp\times p 的 Fisher 信息矩阵，

[I(θ)]ij=−E[∂2ln⁡p(x;θ)∂θi∂θj]i=1,2,⋯ ,p;j=1,2,⋯ ,p[\mathbf{I}(\boldsymbol{\theta})]_{ij} = -E\left [ \dfrac{\partial^2 \ln p(\mathbf{x};\boldsymbol{\theta})}{\partial \theta_i\partial \theta_j} \right]\quad i = 1,2,\cdots, p; j = 1,2,\cdots, p

假设 PDF p(x;θ)p(\mathbf{x};\boldsymbol{\theta}) 满足正则条件

E[∂ln⁡p(x;θ)∂θ]=0E\left [ \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} \right] = 0

则任何无偏估计量 θ^\hat{\boldsymbol{\theta}} 的协方差矩阵满足

Cθ^≥I−1(θ^)\mathbf{C}_{\hat{\boldsymbol{\theta}}} \geq \mathbf{I}^{-1}(\hat{\boldsymbol{\theta}})

当且仅当

[I(θ)]ij=I(θ)(g(x)−θ)[\mathbf{I}(\boldsymbol{\theta})]_{ij} = \mathbf{I}(\boldsymbol{\theta})(\mathbf{g}(x)-\boldsymbol{\theta})

时可达下限。

推导过程·

下面推导矢量参数 α=g(θ)\boldsymbol{\alpha} = \mathbf{g}(\boldsymbol{\theta}) 的 CRLB，考虑无偏估计量

E(αi^)=αi=[g(θ)]ii=1,2,⋯ ,rE(\hat{\alpha_i}) = \alpha_i = [\mathbf{g}(\boldsymbol{\theta})]_i \quad i = 1,2,\cdots, r

根据正则条件 E[∂ln⁡p(x;θ)∂θ]=0E\left [ \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} \right] = 0，可以得到

∫α∂ln⁡p(x;θ)∂θp(x;θ)dx=αE[∂ln⁡p(x;θ)∂θ]=0\int \alpha \dfrac{\partial \ln p(x;\theta)}{\partial \theta}p(x;\theta)\mathrm{d}x = \alpha E\left [ \dfrac{\partial{\ln p(x;\theta)}}{\partial{\theta}} \right] = 0

两式作差得到

∫(α^i−αi)∂ln⁡p(x;θ)∂θip(x;θ)dx=∂[g(θ)]i∂θi\int (\hat{\alpha}_i - \alpha_i) \dfrac{\partial \ln p(x;\theta)}{\partial \theta_i}p(x;\theta)\mathrm{d}x = \dfrac{\partial [g(\theta)]_i}{\partial \theta_i}

当 i≠ji\not= j 时，

∫(α^i−αi)∂ln⁡p(x;θ)∂θjp(x;θ)dx=∫(α^i−αi)∂p(x;θ)∂θjdx=∂∂θj∫α^ip(x;θ)dx−αiE[∂ln⁡p(x;θ)∂θj]=∂αi∂θj=∂[g(θ)]i∂θj\begin{aligned} \int (\hat{\alpha}_i - \alpha_i) \dfrac{\partial \ln p(x;\theta)}{\partial \theta_j}p(x;\theta)\mathrm{d}x &= \int (\hat{\alpha}_i - \alpha_i) \dfrac{\partial p(x;\theta)}{\partial \theta_j}\mathrm{d}x\\ &= \dfrac{\partial}{\partial \theta_j}\int \hat{\alpha}_i p(x;\theta)\mathrm{d}x - \alpha_iE\left [ \dfrac{\partial\ln p(x;\theta)}{\partial \theta_j} \right]\\ &= \dfrac{\partial \alpha_i}{\partial \theta_j}\\ &= \dfrac{\partial [g(\theta)]_i}{\partial \theta_j}\\ \end{aligned}

组合成矩阵形式

∫(α^−α)∂ln⁡p(x;θ)⊤∂θjp(x;θ)dx=∂[g(θ)]∂θ\int (\hat{\boldsymbol{\alpha}} - \boldsymbol{\alpha}) \dfrac{\partial \ln p(\mathbf{x};\boldsymbol{\theta})^{\top}}{\partial \theta_j}p(\mathbf{x};\boldsymbol{\theta})\mathrm{d}x = \dfrac{\partial [\mathbf{g}(\boldsymbol{\theta})]}{\partial \boldsymbol{\theta}}

对于任意的 t×1t\times 1 矢量 a\mathbf{a} 和 p×1p\times 1 矢量 b\mathbf{b}，

∫a⊤(α^−α)∂ln⁡p(x;θ)⊤∂θjbp(x;θ)dx=a⊤∂[g(θ)]∂θb\int \mathbf{a}^{\top}(\hat{\boldsymbol{\alpha}} - \boldsymbol{\alpha}) \dfrac{\partial \ln p(\mathbf{x};\boldsymbol{\theta})^{\top}}{\partial \theta_j} \mathbf{b} p(\mathbf{x};\boldsymbol{\theta})\mathrm{d}x = \mathbf{a}^{\top}\dfrac{\partial [\mathbf{g}(\boldsymbol{\theta})]}{\partial \boldsymbol{\theta}}\mathbf{b}

由柯西-施瓦茨不等式

a⊤∂[g(θ)]∂θb≤∫a⊤(α^−α)(α^−α)⊤ap(x;θ)dx⋅∫b⊤∂ln⁡p(x;θ)∂θj∂ln⁡p(x;θ)⊤∂θjbp(x;θ)dx=a⊤Cα^abTI(θ)b\begin{aligned} \mathbf{a}^{\top}\dfrac{\partial [\mathbf{g}(\boldsymbol{\theta})]}{\partial \boldsymbol{\theta}}\mathbf{b} &\leq \int \mathbf{a}^{\top}(\hat{\boldsymbol{\alpha}} - \boldsymbol{\alpha})(\hat{\boldsymbol{\alpha}} - \boldsymbol{\alpha})^{\top}\mathbf{a}p(\mathbf{x};\boldsymbol{\theta})\mathrm{d}x \cdot \int \mathbf{b}^{\top} \dfrac{\partial \ln p(\mathbf{x};\boldsymbol{\theta})}{\partial \theta_j}\dfrac{\partial \ln p(\mathbf{x};\boldsymbol{\theta})^{\top}}{\partial \theta_j} \mathbf{b} p(\mathbf{x};\boldsymbol{\theta})\mathrm{d}x\\ &= \mathbf{a}^{\top} C_{\hat{\alpha}} \mathbf{a} \mathbf{b}^T I(\boldsymbol{\theta}) \mathbf{b}\\ \end{aligned}

并且

E[∂ln⁡p(x;θ)∂θi∂ln⁡p(x;θ)∂θj]=−E[∂2ln⁡p(x;θ)∂θi∂θj]=[I(θ)]ijE\left [ \frac{\partial \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \theta_i} \frac{\partial \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \theta_j} \right] = -E\left [ \frac{\partial^2 \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \theta_i \partial \theta_j} \right] = [\mathbf{I}(\boldsymbol{\theta})]_{ij}

令

b=I−1(θ)∂g(θ)T∂θa\mathbf{b} = \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{g}(\boldsymbol{\theta})^T}{\partial \boldsymbol{\theta}} \mathbf{a}

进一步得到

(aT∂g(θ)∂θI−1(θ)∂g(θ)T∂θa)2⩽aTCα^a(aT∂g(θ)∂θI−1(θ)∂g(θ)T∂θa)\left( \mathbf{a}^T \frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{g}(\boldsymbol{\theta})^T}{\partial \boldsymbol{\theta}} \mathbf{a} \right)^2 \leqslant \mathbf{a}^T \mathbf{C}_{\hat{\alpha}} \mathbf{a} \left( \mathbf{a}^T \frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{g}(\boldsymbol{\theta})^T}{\partial \boldsymbol{\theta}} \mathbf{a} \right)

由于 I(θ)\mathbf{I}(\boldsymbol{\theta}) 正定，得到

aT(Cα^−∂g(θ)∂θI−1(θ)∂g(θ)T∂θ)a⩾0\mathbf{a}^T \left( \mathbf{C}_{\hat{\alpha}} - \frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{g}(\boldsymbol{\theta})^T}{\partial \boldsymbol{\theta}} \right) \mathbf{a} \geqslant 0

aT(α^−α)=c∂ln⁡p(x;θ)T∂θb=c∂ln⁡p(x;θ)T∂θI−1(θ)∂g(θ)T∂θa\mathbf{a}^T (\hat{\alpha} - \alpha) = c \frac{\partial \ln p(\mathbf{x}; \boldsymbol{\theta})^T}{\partial \boldsymbol{\theta}} \mathbf{b} = c \frac{\partial \ln p(\mathbf{x}; \boldsymbol{\theta})^T}{\partial \boldsymbol{\theta}} \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{g}(\boldsymbol{\theta})^T}{\partial \boldsymbol{\theta}} \mathbf{a}

取等条件为

∂g(θ)∂θI−1(θ)∂ln⁡p(x;θ)∂θ=1c(α^−α)\frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \frac{1}{c} (\hat{\alpha} - \alpha)

与上一部分推导标量的类似，考虑 α=g(θ)=θ\boldsymbol{\alpha} = \mathbf{g}(\boldsymbol{\theta}) = \boldsymbol{\theta} 时，∂g(θ)∂θ=I\dfrac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \mathbf{I}

∂ln⁡p(x;θ)∂θ=1cI(θ)(θ^−θ)\frac{\partial \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \frac{1}{c} \mathbf{I}(\boldsymbol{\theta}) (\hat{\boldsymbol{\theta}} - \boldsymbol{\theta})

∂ln⁡p(x;θ)∂θi=∑k=1p[I(θ)]ikc(θ)(θ^k−θk)∂2ln⁡p(x;θ)∂θi∂θj=∑k=1p([I(θ)]ikc(θ)(−δkj)+∂([I(θ)]ikc(θ))∂θj(θ^k−θk))\begin{aligned} \frac{\partial \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \theta_i} &= \sum_{k = 1}^p \frac{[\mathbf{I}(\boldsymbol{\theta})]_{ik}}{c(\boldsymbol{\theta})} (\hat{\theta}_k - \theta_k)\\ \frac{\partial^2 \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \theta_i \partial \theta_j} &= \sum_{k = 1}^p \left( \frac{[\mathbf{I}(\boldsymbol{\theta})] _{ik}}{c(\boldsymbol{\theta})} (-\delta_{kj}) + \frac{\partial \left( \dfrac{[\mathbf{I}(\boldsymbol{\theta})]_{ik}}{c(\boldsymbol{\theta})} \right)}{\partial \theta_j} (\hat{\theta}_k - \theta_k) \right) \end{aligned}

[I(θ)]ij=−E[∂2ln⁡p(x;θ)∂θi∂θj]=[I(θ)]ijc(θ)[\mathbf{I}(\boldsymbol{\theta})] _{ij} = -E\left [ \frac{\partial^2 \ln p(\mathbf{x}; \boldsymbol{\theta})}{\partial \theta_i \partial \theta_j} \right] = \frac{[\mathbf{I}(\boldsymbol{\theta})]_{ij}}{c(\boldsymbol{\theta})}

附录·

下文推导基于以下假设：

信号到达方向对于各个阵元是一致的；信号扫过阵列的过程中只考虑相移，忽略包络的变化。

设阵元数为 MM，快拍数为 NN，观测数据为

x(t)=A(θ)s(t)+ω(t)t=1,2,⋯ ,Nx(t) = \mathbf{A}(\theta)\mathbf{s}(t) + \boldsymbol{\omega}(t) \quad t = 1,2,\cdots, N

其中

  • s(t)∈CM×1\mathbf{s}(t) \in \mathbb{C}^{M\times 1}，表示信号矢量；
  • ω(t)∈CM×1\boldsymbol{\omega}(t)\in \mathbb{C}^{M\times 1}，复加性高斯白噪声，各个阵元的噪声独立，且时间上互不相关，满足分布 ω(t)∼CN(0,σ2IM)\boldsymbol{\omega}(t)\sim \mathcal{CN}(0, \sigma^2\mathbf{I_M})；
  • x(t)∈CM×1\mathbf{x}(t) \in \mathbb{C}^{M\times 1}，表示含噪声的接收数据矢量，满足分布 x(t)∼CN(A(θ)s(t),σ2IM)\mathbf{x}(t)\sim \mathcal{CN}(\mathbf{A}(\theta)\mathbf{s}(t), \sigma^2\mathbf{I_M}).

复高斯随机变量概率密度函数的推导·

设复随机变量

Z=X+jYZ = X + jY

其中 XX 与 YY 均为实值随机变量，若 ZZ 是一个复原对称高斯随机变量，则

  • XX 与 YY 相互独立；

  • XX 与 YY 满足 X∼N(μ,σ22),Y∼N(μ,σ22)X\sim \mathcal{N}(\mu, \dfrac{\sigma^2}{2}), Y\sim \mathcal{N}(\mu, \dfrac{\sigma^2}{2})

  • ZZ 的均值表示为 μ=μX+jμY\mu = \mu_{X} + j\mu_{Y}；

下面推导其概率密度函数（PDF），本质上是其 实部与虚部的二维联合 PDF。

因为 XX 与 YY 独立，则联合密度为

fX,Y(x,y)=fX(x)fY(y)f_{X, Y}(x, y) = f_{X}(x) f_{Y}(y)

其中

fX(x)=1πσ2exp⁡(−(x−μX)2σ2)fY(y)=1πσ2exp⁡(−(y−μY)2σ2)\begin{aligned} f_{X}(x) = \dfrac{1}{\sqrt{\pi\sigma^2}}\exp{\left( -\dfrac{(x-\mu_{X})^2}{\sigma^2} \right)}\\ f_{Y}(y) = \dfrac{1}{\sqrt{\pi\sigma^2}}\exp{\left( -\dfrac{(y-\mu_{Y})^2}{\sigma^2} \right)}\\ \end{aligned}

因此

fX,Y(x,y)=fX(x)fY(y)=1πσ2exp⁡(−(x−μX)2+(y−μY)2σ2)\begin{aligned} f_{X, Y}(x, y) &= f_{X}(x) f_{Y}(y)\\ &= \dfrac{1}{\pi\sigma^2}\exp{\left( -\dfrac{(x-\mu_{X})^2 + (y-\mu_{Y})^2}{\sigma^2} \right)}\\ \end{aligned}

利用

Z=X+jY,μ=μX+jμYZ = X + jY, \mu = \mu_X + j\mu_Y

得到

∣z−μ∣2=∣(x−μX)+j(y−μY)∣2=(x−μX)2+(y−μY)2|z - \mu|^2 = |(x-\mu_X) + j(y-\mu_Y)|^2 = (x-\mu_X)^2 + (y-\mu_Y)^2

则

fZ(z)=fX,Y(ℜz,ℑz)=fX,Y(x,y)=1πσ2exp⁡(−(x−μX)2+(y−μY)2σ2)=1πσ2exp⁡(−∣z−μ∣2σ2)\begin{aligned} f_Z(z) &= f_{X, Y}(\Re{z}, \Im{z})\\ &= f_{X, Y}(x, y) = \dfrac{1}{\pi\sigma^2}\exp{\left( -\dfrac{(x-\mu_{X})^2 + (y-\mu_{Y})^2}{\sigma^2} \right)}\\ &= \dfrac{1}{\pi\sigma^2}\exp{\left(-\dfrac{|z-\mu|^2}{\sigma^2} \right)}\\ \end{aligned}

得到了一维 圆对称复高斯分布 的概率密度函数。

对于 mm 维独立的复高斯随机变量，其联合概率密度函数为

p(z)=∏i=1m1πσ2exp⁡(−∣z−μ∣2σ2)=∏i=1m1(πσ2)mexp⁡(−∑i=1m∣z−μ∣2σ2)=∏i=1m1(πσ2)mexp⁡(−∣∣z−μ∣∣2σ2)\begin{aligned} p(\mathbf{z}) &= \prod_{i = 1}^m \dfrac{1}{\pi\sigma^2}\exp{\left(-\dfrac{|z-\mu|^2}{\sigma^2} \right)}\\ &= \prod_{i = 1}^m \dfrac{1}{(\pi\sigma^{2})^m}\exp{\left(-\dfrac{\sum_{i = 1}^m |z-\mu|^2}{\sigma^2} \right)}\\ &= \prod_{i = 1}^m \dfrac{1}{(\pi\sigma^{2})^m}\exp{\left(-\dfrac{||\mathbf{z} - \boldsymbol{\mu}||^2}{\sigma^2} \right)}\\ \end{aligned}

对数似然函数的推导·

观测数据 x(t)\mathbf{x}(t) 满足分布 x(t)∼CN(A(θ)s(t),σ2IM)\mathbf{x}(t)\sim \mathcal{CN}(\mathbf{A}(\theta)\mathbf{s}(t), \sigma^2\mathbf{I_M})，均值向量为 μ(t)=A(θ)s(t)\boldsymbol{\mu}(t) = \mathbf{A}(\theta)\mathbf{s}(t)，带入上式，得到单个快拍的 PDF 为

p(x(t)∣θ)=1(πσ2)Mexp⁡(−1σ2∣∣x(t)−A(θ)s(t)∣∣2)p(\mathbf{x}(t)|\theta) = \dfrac{1}{(\pi \sigma^2)^{M}} \exp{\left(-\dfrac{1}{\sigma^2}||\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)||^2\right)}

由于 NN 个快拍在时间上相互独立，则似然函数为各快拍 PDF 的乘积，即

L(θ)=∏t=1Np(x(t)∣θ)=1(πσ2)MNexp⁡(−1σ2∑t=1N∣∣x(t)−A(θ)s(t)∣∣2)\begin{aligned} \mathcal{L}(\theta) &= \prod_{t = 1}^N p(\mathbf{x}(t)|\theta)\\ &= \dfrac{1}{(\pi \sigma^2)^{MN}} \exp{\left(-\dfrac{1}{\sigma^2}\sum_{t = 1}^{N}||\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)||^2\right)} \end{aligned}

最终的对数似然函数为

ln⁡L(θ)=−MNln⁡πσ2−1σ2∑t=1N∣∣x(t)−A(θ)s(t)∣∣2\ln \mathcal{L}(\theta) = -MN \ln{\pi\sigma^2} - \dfrac{1}{\sigma^2}\sum_{t = 1}^{N}||\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)||^2

Fisher 矩阵的推导·

对于单个参数 θ\theta，Fisher 矩阵为

I(θ)=−E[∂2∂θ2ln⁡L(θ)]\mathcal{I}(\theta) = -\mathrm{E}\left [\dfrac{\partial^2}{\partial\theta^2}\ln\mathcal{L}(\theta)\right]

对于矢量参数 θ∈Rp\boldsymbol{\theta}\in \mathbb{R}^p，Fisher 矩阵为

Iij(θ)=−E[∂2∂θi∂θjln⁡L(θ)]\mathcal{I}_{ij}(\boldsymbol{\theta}) = -\mathrm{E}\left [\dfrac{\partial^2}{\partial\theta_i\partial \theta_j}\ln\mathcal{L}(\boldsymbol{\theta})\right]

首先计算对数似然函数的一阶导数（忽略常数项）

ln⁡L(θ)=C−1σ2∑t=1N∣∣x(t)−A(θ)s(t)∣∣2=C−1σ2∑t=1N(x(t)−A(θ)s(t))†(x(t)−A(θ)s(t))\begin{aligned} \ln \mathcal{L}(\theta) &= C - \dfrac{1}{\sigma^2}\sum_{t = 1}^{N}||\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)||^2\\ &= C - \dfrac{1}{\sigma^2}\sum_{t = 1}^{N}\left(\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)\right)^{\dagger}\left(\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)\right)\\ \end{aligned}

令 z=x(t),z0(θ)=A(θ)s(t),e=z−z0z=\mathbf{x}(t),z_0(\theta) = \mathbf{A}(\theta)\mathbf{s}(t), e = z-z_0，关注和参数 θ\theta 有关的部分，

J(θ)=−1σ2∑t=1N(x(t)−A(θ)s(t))†(x(t)−A(θ)s(t))=−1σ2∑t=1Ne†e\begin{aligned} \mathcal{J}(\theta) &= - \dfrac{1}{\sigma^2}\sum_{t = 1}^{N}\left(\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)\right)^{\dagger}\left(\mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t)\right)\\ &= - \dfrac{1}{\sigma^2}\sum_{t = 1}^{N}e^{\dagger}e\\ \end{aligned}

则有

∂e†e∂θ=(−∂z0∂θ)†e+e†(−∂z0∂θ)=(−∂z0∂θ)†e+[(−∂z0∂θ)†e]∗=−2ℜ(∂z0∂θ)†e=−2ℜ(∂z0∂θ)†(z−z0)\begin{aligned} \dfrac{\partial e^{\dagger}e}{\partial \theta} &= \left(-\dfrac{\partial z_0}{\partial \theta}\right)^{\dagger}e + e^{\dagger}\left(-\dfrac{\partial z_0}{\partial \theta}\right)\\ &= \left(-\dfrac{\partial z_0}{\partial \theta}\right)^{\dagger}e + \left [\left(-\dfrac{\partial z_0}{\partial \theta}\right)^{\dagger}e\right]^*\\ &= -2\Re{\left(\dfrac{\partial z_0}{\partial \theta}\right)^{\dagger}e}\\ &= -2\Re{\left(\dfrac{\partial z_0}{\partial \theta}\right)^{\dagger}(z-z_0)}\\ \end{aligned}

令 dθ(t)=∂z0∂θ=∂A(θ)s(t)∂θ\mathbf{d}_{\theta}(t) = \dfrac{\partial z_0}{\partial \theta}= \dfrac{\partial \mathbf{A}(\theta)\mathbf{s}(t)}{\partial \theta} 代入到对数似然函数的导数中

ln⁡L(θ)∂θ=2σ2∑t=1Nℜdθ†(t)[x(t)−A(θ)s(t)]\dfrac{\ln\mathcal{L}(\theta)}{\partial \theta} = \dfrac{2}{\sigma^2} \sum_{t = 1}^{N} \Re{\mathbf{d}^{\dagger}_{\theta}(t)\left [ \mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t) \right]}

下面计算对数似然函数的二阶导数

已经得到一阶导数为

∂ln⁡L(θ)∂θ=2σ2∑t=1Nℜdθ†(t)[x(t)−A(θ)s(t)]=2σ2∑t=1Nℜdθ†(t)ω(t)\begin{aligned} \dfrac{\partial \ln\mathcal{L}(\theta)}{\partial \theta} &= \dfrac{2}{\sigma^2} \sum_{t = 1}^{N} \Re{\mathbf{d}^{\dagger}_{\theta}(t)\left [ \mathbf{x}(t) - \mathbf{A}(\theta)\mathbf{s}(t) \right]}\\ &= \dfrac{2}{\sigma^2} \sum_{t = 1}^{N} \Re{\mathbf{d}^{\dagger}_{\theta}(t)\boldsymbol{\omega}(t)}\\ \end{aligned}

对于复加性高斯白噪声，ω(t)∼CN(0,σ2IM)\boldsymbol{\omega}(t)\sim \mathcal{CN}(0, \sigma^2\mathbf{I_M})，因此

E(ω(t))=0, E(ω(t)ω(t)†)=σ2IM\mathrm{E}(\boldsymbol{\omega}(t)) = 0, \space \mathrm{E}(\boldsymbol{\omega}(t)\boldsymbol{\omega}(t)^{\dagger}) = \boldsymbol{\sigma}^2I_M

并且

E(dθ(t)†ω(t))=0, E(ω(t)†dθ(t))=0\mathrm{E}(\mathbf{d}_{\theta}(t)^{\dagger}\boldsymbol{\omega}(t)) = 0, \space \mathrm{E}(\boldsymbol{\omega}(t)^{\dagger}\mathbf{d}_{\theta}(t)) = 0

于是似然函数的二阶导为

∂2ln⁡L(θ)∂θ2=2σ2∑t=1Nℜ∂∂θ[dθ†(t)ω(t)]\dfrac{\partial^2 \ln\mathcal{L}(\theta)}{\partial \theta^2} = \dfrac{2}{\sigma^2} \sum_{t = 1}^{N} \Re{\dfrac{\partial}{\partial \theta}\left[\mathbf{d}^{\dagger}_{\theta}(t)\boldsymbol{\omega}(t)\right]}

其中

∂ω(t)∂θ=−dθ(t)∂∂θ[dθ†(t)ω(t)]=(∂dθ(t)∂θ)†ω(t)+dθ†(t)(−dθ(t))\begin{aligned} \dfrac{\partial \boldsymbol{\omega}(t)}{\partial \theta} &= -\mathbf{d}_{\theta}(t)\\ \dfrac{\partial}{\partial \theta}\left[\mathbf{d}^{\dagger}_{\theta}(t)\boldsymbol{\omega}(t)\right] &= \left(\dfrac{\partial \mathbf{d}_{\theta}(t)}{\partial \theta}\right)^{\dagger}\boldsymbol{\omega}(t) + \mathbf{d}^{\dagger}_{\theta}(t)(-\mathbf{d}_{\theta}(t)) \end{aligned}

代入对数似然函数二阶导数的表达式得到

∂2ln⁡L(θ)∂θ2=2σ2∑t=1Nℜ∂∂θ[dθ†(t)ω(t)]=2σ2∑t=1Nℜ(∂dθ(t)∂θ)†ω(t)−dθ†(t)dθ(t)\begin{aligned} \dfrac{\partial^2 \ln\mathcal{L}(\theta)}{\partial \theta^2} &= \dfrac{2}{\sigma^2} \sum_{t = 1}^{N} \Re{\dfrac{\partial}{\partial \theta}\left[\mathbf{d}^{\dagger}_{\theta}(t)\boldsymbol{\omega}(t)\right]}\\ &= \dfrac{2}{\sigma^2} \sum_{t = 1}^{N} \Re{\left(\dfrac{\partial \mathbf{d}_{\theta}(t)}{\partial \theta}\right)^{\dagger}\boldsymbol{\omega}(t) - \mathbf{d}^{\dagger}_{\theta}(t)\mathbf{d}_{\theta}(t)}\\ \end{aligned}

其中 dθ(t)=∂A(θ)s(t)∂θ\mathbf{d}_{\theta}(t) = \dfrac{\partial \mathbf{A}(\theta)\mathbf{s}(t)}{\partial \theta}.

下面取期望得到 Fisher 信息

I(θ)=−E[∂2∂θ2ln⁡L(θ)]\mathcal{I}(\theta) = -\mathrm{E}\left [\dfrac{\partial^2}{\partial\theta^2}\ln\mathcal{L}(\theta)\right]

对数似然函数的二阶表达式中

E[(∂dθ(t)∂θ)†ω(t)]=0E[dθ†(t)dθ(t)]=dθ†(t)dθ(t)\begin{aligned} \mathrm{E}\left[ \left(\dfrac{\partial \mathbf{d}_{\theta}(t)}{\partial \theta}\right)^{\dagger}\boldsymbol{\omega}(t) \right] &= 0\\ \mathrm{E}\left [ \mathbf{d}^{\dagger}_{\theta}(t)\mathbf{d}_{\theta}(t) \right] &= \mathbf{d}^{\dagger}_{\theta}(t)\mathbf{d}_{\theta}(t) \end{aligned}

那么

I(θ)=2σ2∑t=1N∣∣dθ(t)∣∣2=2σ2∑t=1N∣∣A′(θ)s(t)∣∣2\mathcal{I}(\theta) = \dfrac{2}{\sigma^2} \sum_{t = 1}^{N}||\mathbf{d}_{\theta}(t)||^2 = \dfrac{2}{\sigma^2} \sum_{t = 1}^{N}||\mathbf{A}'(\theta)\mathbf{s}(t)||^2

对于多参数的情况，即 θ=[θ1,θ2,⋯ ,θp]⊤\boldsymbol{\theta} = [\theta_1,\theta_2,\cdots, \theta_p]^{\top}，

Iij=2σ2∑t=1Nℜdθ†(t)dθ(t)=2σ2∑t=1Nℜs(t)†(Aθi′(t))†Aθj′(t)s(t)\mathcal{I}_{ij} = \dfrac{2}{\sigma^2} \sum_{t = 1}^{N}\Re{\mathbf{d}^{\dagger}_{\theta}(t)\mathbf{d}_{\theta}(t)} = \dfrac{2}{\sigma^2} \sum_{t = 1}^{N}\Re{\mathbf{s}(t)^{\dagger}\left( \mathbf{A}'_{\theta_i}(t) \right)^{\dagger}\mathbf{A}'_{\theta_j}(t)\mathbf{s}(t)}

CRLB 的推导·

对于实参数 θ\theta，其 Cramér–Rao 下界 为

Var(θ^)≥1I(θ)\mathrm{Var}(\hat{\theta}) \geq \dfrac{1}{\mathcal{I}(\theta)}

则

CRB(θ)=σ22∑t=1N∣∣A′(θ)s(t)∣∣2\mathrm{CRB}(\theta) = \dfrac{\sigma^2}{2\sum_{t = 1}^{N}||\mathbf{A}'(\theta)\mathbf{s}(t)||^2}

单信号源+均匀线阵·

下面考虑 单信号源，ULA 阵列，设阵元数为 MM，阵元间距 d=λ2d=\dfrac{\lambda}{2}，ULA 导向矢量为

a(θ)=[1e−j2πdsin⁡θλ⋮e−j2πd(M−1)sin⁡θλ]\mathbf{a}(\theta) = \begin{bmatrix} 1 \\ e^{-j\frac{2\pi d\sin{\theta}}{\lambda}} \\ \vdots \\ e^{-j\frac{2\pi d(M-1)\sin{\theta}}{\lambda}} \\ \end{bmatrix}

导向矢量对参数 θ\theta 的导数为

∂a(θ)∂θ=−j2πdλcos⁡(θ)[1e−j2πdsin⁡θλ⋮(M−1)e−j2πd(M−1)sin⁡θλ]\dfrac{\partial a(\theta)}{\partial \theta} = -j\dfrac{2\pi d}{\lambda}\cos(\theta) \begin{bmatrix} 1 \\ e^{-j\frac{2\pi d\sin{\theta}}{\lambda}} \\ \vdots \\ (M-1)e^{-j\frac{2\pi d(M-1)\sin{\theta}}{\lambda}} \\ \end{bmatrix}

2-范数平方为

∣∣a′(θ)∣∣2=(2πdλcos⁡(θ))2∑m=0M−1m2=(2πdλcos⁡(θ))2M(M−1)(2M−1)6\begin{aligned} ||a'(\theta)||^2 &= \left(\dfrac{2\pi d}{\lambda}\cos(\theta)\right)^2 \sum_{m = 0}^{M-1} m^2\\ &= \left(\dfrac{2\pi d}{\lambda}\cos(\theta)\right)^2 \dfrac{M(M-1)(2M-1)}{6} \end{aligned}

最终得到单信源， ULA 阵列的 Cramér–Rao 下界 为

CRB(θ)=σ22∑t=1N∣∣A′(θ)s(t)∣∣2=σ22N∣s(t)∣2∣∣a′(θ)∣∣2=3λ2N⋅SNR⋅(2πdcos⁡(θ))2M(M−1)(2M−1)\begin{aligned} \mathrm{CRB}(\theta) &= \dfrac{\sigma^2}{2\sum_{t = 1}^{N}||\mathbf{A}'(\theta)\mathbf{s}(t)||^2}\\ &= \dfrac{\sigma^2}{2N|\mathbf{s}(t)|^2||a'(\theta)||^2}\\ &= \dfrac{3\lambda^2}{N\cdot\mathrm{SNR}\cdot\left(2\pi d\cos(\theta)\right)^2 M(M-1)(2M-1)} \end{aligned}

单信号源+均匀面阵·

下面考虑 单信号源，UPA 阵列，设水平方向阵元数为 MM，垂直方向阵元数为 NN，则总阵元数为 MNMN，设入射角的俯仰角为 ϕ\phi，方位角为 θ\theta，则 UPA 的导向矢量表示为

a(θ,ϕ)=e−j2πdλ(msin⁡θcos⁡ϕ+nsin⁡θsin⁡ϕ)a(\theta, \phi) = e^{-j\frac{2\pi d}{\lambda}(m\sin\theta\cos\phi + n\sin\theta\sin\phi)}

其中，

{m=0,1,⋯ ,M−1n=0,1,⋯ ,N−1\left\{ \begin{aligned} m &= 0,1,\cdots, M - 1\\ n &= 0,1,\cdots, N - 1\\ \end{aligned} \right.

导向矢量对方位角和俯仰角求偏导

∂am,n∂θ=−j2πdλ(mcos⁡θcos⁡ϕ+ncos⁡θsin⁡ϕ)am,n∂am,n∂ϕ=−j2πdλ(−msin⁡θsin⁡ϕ+nsin⁡θcos⁡ϕ)am,n\begin{aligned} \dfrac{\partial a_{m, n}}{\partial \theta} &= -j\dfrac{2\pi d}{\lambda}(m\cos\theta\cos\phi + n\cos\theta\sin\phi) a_{m, n}\\ \dfrac{\partial a_{m, n}}{\partial \phi} &= -j\dfrac{2\pi d}{\lambda}(-m\sin\theta\sin\phi + n\sin\theta\cos\phi) a_{m, n}\\ \end{aligned}

下面求 Fisher 矩阵，对于二维参数 θ=[θ,ϕ]\boldsymbol{\theta} = [\theta, \phi]，

I=[IθθIθϕIϕθIϕϕ]\mathcal{I} = \begin{bmatrix} \mathcal{I}_{\theta\theta} & \mathcal{I}_{\theta\phi}\\ \mathcal{I}_{\phi\theta} & \mathcal{I}_{\phi\phi}\\ \end{bmatrix}

其中

Iij=2σ2∑t=1Nℜdθ†(t)dθ(t)=2N∣s(t)∣2σ2ℜ((ai′)†aj′\mathcal{I}_{ij} = \dfrac{2}{\sigma^2} \sum_{t = 1}^{N}\Re{\mathbf{d}^{\dagger}_{\theta}(t)\mathbf{d}_{\theta}(t)} = \dfrac{2N|s(t)|^2}{\sigma^2} \Re{((a_i')^{\dagger}a_j'}

对于方位角 θ\theta

Iθθ=2N∣s(t)∣2σ2ℜ(∂a∂θ)†(∂a∂θ)\mathcal{I}_{\theta\theta} = \dfrac{2N|s(t)|^2}{\sigma^2} \Re{\left(\dfrac{\partial \mathbf{a}}{\partial \theta}\right)^{\dagger}\left(\dfrac{\partial \mathbf{a}}{\partial \theta}\right)}

∣∂am,n∂θ∣2=(2πdλ)2cos⁡2θ(mcos⁡ϕ+nsin⁡ϕ)2\left| \dfrac{\partial a_{m, n}}{\partial \theta} \right|^2 = \left(\dfrac{2\pi d}{\lambda}\right)^2\cos^2\theta(m\cos\phi + n\sin\phi)^2

则

(∂a∂θ)†(∂a∂θ)=(2πdλ)2cos⁡2θ∑m=0M−1∑n=0N−1(mcos⁡ϕ+nsin⁡ϕ)2=(2πdλ)2cos⁡2θ∑m=0M−1∑n=0N−1(m2cos⁡2ϕ+n2sin⁡2ϕ+2mncos⁡ϕsin⁡ϕ)\begin{aligned} \left(\dfrac{\partial \mathbf{a}}{\partial \theta}\right)^{\dagger}\left(\dfrac{\partial \mathbf{a}}{\partial \theta}\right) &= \left(\dfrac{2\pi d}{\lambda}\right)^2\cos^2\theta \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1} (m\cos\phi + n\sin\phi)^2\\ &= \left(\dfrac{2\pi d}{\lambda}\right)^2\cos^2\theta \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1} (m^2\cos^2\phi + n^2\sin^2\phi + 2mn\cos\phi\sin\phi)\\ \end{aligned}

其中

∑m=0M−1∑n=0N−1m2cos⁡2ϕ=Ncos⁡2ϕM(M−1)(2M−1)6∑m=0M−1∑n=0N−1n2sin⁡2ϕ=Msin⁡2ϕN(N−1)(2N−1)6∑m=0M−1∑n=0N−12mncos⁡ϕsin⁡ϕ=2cos⁡ϕsin⁡ϕM(M−1)2N(N−1)2\begin{aligned} \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1}m^2\cos^2\phi &= N\cos^2\phi\dfrac{M(M-1)(2M-1)}{6}\\ \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1}n^2\sin^2\phi &= M\sin^2\phi\dfrac{N(N-1)(2N-1)}{6}\\ \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1}2mn\cos\phi\sin\phi &= 2\cos\phi\sin\phi\dfrac{M(M-1)}{2}\dfrac{N(N-1)}{2}\\ \end{aligned}

得到

Iθθ=2N∣s(t)∣2σ2ℜ(∂a∂θ)†(∂a∂θ)=2N∣s(t)∣2σ2{(2πdλ)2cos⁡2θ∑m=0M−1∑n=0N−1(m2cos⁡2ϕ+n2sin⁡2ϕ+2mncos⁡ϕsin⁡ϕ)}=2N∣s(t)∣2σ2(2πdλ)2cos⁡2θ(cos⁡2ϕMN(M−1)(2M−1)6+sin⁡2ϕMN(N−1)(2N−1)6+cos⁡ϕsin⁡ϕMN(M−1)(N−1)2)\begin{aligned} \mathcal{I}_{\theta\theta} &= \dfrac{2N|s(t)|^2}{\sigma^2} \Re{\left(\dfrac{\partial \mathbf{a}}{\partial \theta}\right)^{\dagger}\left(\dfrac{\partial \mathbf{a}}{\partial \theta}\right)}\\ &= \dfrac{2N|s(t)|^2}{\sigma^2} \left\{ \left(\dfrac{2\pi d}{\lambda}\right)^2\cos^2\theta \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1} (m^2\cos^2\phi + n^2\sin^2\phi + 2mn\cos\phi\sin\phi) \right\}\\ &= \dfrac{2N|s(t)|^2}{\sigma^2} \left(\dfrac{2\pi d}{\lambda}\right)^2\cos^2\theta \left( \cos^2\phi\dfrac{MN(M-1)(2M-1)}{6} + \sin^2\phi\dfrac{MN(N-1)(2N-1)}{6} + \cos\phi\sin\phi\dfrac{MN(M-1)(N-1)}{2} \right)\\ \end{aligned}

同理，可以求得

Iϕϕ=2N∣s(t)∣2σ2(2πdλ)2sin⁡2θ(sin⁡2ϕMN(M−1)(2M−1)6+cos⁡2ϕMN(N−1)(2N−1)6−cos⁡ϕsin⁡ϕMN(M−1)(N−1)2)Iθϕ=2N∣s(t)∣2σ2(2πdλ)2sin⁡2θ24⋅MN⋅(sin⁡2ϕ[(N−1)(2N−1)−(M−1)(2M−1)]+3cos⁡2ϕ(M−1)(N−1))\begin{aligned} \mathcal{I}_{\phi\phi} &= \dfrac{2N|s(t)|^2}{\sigma^2} \left(\dfrac{2\pi d}{\lambda}\right)^2\sin^2\theta \left( \sin^2\phi\dfrac{MN(M-1)(2M-1)}{6} + \cos^2\phi\dfrac{MN(N-1)(2N-1)}{6} - \cos\phi\sin\phi\dfrac{MN(M-1)(N-1)}{2} \right)\\ \mathcal{I}_{\theta\phi} &= \dfrac{2N|s(t)|^2}{\sigma^2} \left(\dfrac{2\pi d}{\lambda}\right)^2\dfrac{\sin2\theta}{24}\cdot MN\cdot \left(\sin2\phi \left [(N-1)(2N-1) - (M-1)(2M-1)\right] + 3\cos2\phi(M-1)(N-1) \right) \end{aligned}

最终的 Fisher 矩阵

I=[IθθIθϕIθϕIϕϕ]\mathcal{I} = \begin{bmatrix} \mathcal{I}_{\theta\theta} & \mathcal{I}_{\theta\phi}\\ \mathcal{I}_{\theta\phi} & \mathcal{I}_{\phi\phi}\\ \end{bmatrix}

最终的 CRLB 表示为

CRLBθ=[I−1]θθ=IϕϕIθθIϕϕ−Iθϕ2CRLBϕ=[I−1]ϕϕ=IθθIθθIϕϕ−Iθϕ2\begin{aligned} \mathrm{CRLB}_{\theta} &= [\mathcal{I}^{-1}]_{\theta\theta} = \dfrac{\mathcal{I}_{\phi\phi}}{\mathcal{I}_{\theta\theta}\mathcal{I}_{\phi\phi}-\mathcal{I}^2_{\theta\phi}}\\ \mathrm{CRLB}_{\phi} &= [\mathcal{I}^{-1}]_{\phi\phi} = \dfrac{\mathcal{I}_{\theta\theta}}{\mathcal{I}_{\theta\theta}\mathcal{I}_{\phi\phi}-\mathcal{I}^2_{\theta\phi}}\\ \end{aligned}