二元干涉仪仿真

# 二元加法干涉仪仿真

# 点源推导输出功率

假设我们生活在一维空间，即一张纸上；射电源为点源，且离观测者足够远，到达观测者的电磁波近似于平面波。两天线固定指向，当射电源偏离天线指向的角度为$\theta$ 时，两路电磁波入射到天线上，走过的几何长度之差：$\Delta L=bsin\theta$

两路电磁波时延之差：$\Delta L=\frac{\Delta L}{c}=\frac{b}{c} sin\theta$

入射电磁波可以设为平面正弦波，两路同轴线长度、低噪放等配置一致，则有

天线 A 感应到的电压：$V_A=V_0sin[2\pi(\frac{x}{\lambda}-ft)]=V_0sin(kx-\omega t)$

天线 A 感应到电压的复数形式：$V_A=V_0e^{i(kx-\omega t)}$

天线 B 感应到电压的复数形式：$V_B=V_0e^{i[kx-\omega (t-\tau)]}$

经过合路器相加，有总输出电压：

$A=V_A+V_B=V_0e^{i(kx-\omega t)}(1+ e^{i\omega \tau})$

输出总功率：

$P=AA^*=V_0^2(1+e^{i\omega \tau})(1+e^{-i\omega\tau})\\=2A_0^2[1+cos(\omega \tau)]=2A_0^2[1+cos(2\pi fbsin\theta/c)] (干涉仪公式1)$

由于天体在天空中是运动的，因此固定天线指向时，$\theta$ 是随时间均匀变化的，有$\theta=\dot{\theta}t$，一天转 360°，相当于每秒转 360°/(24*3600)=0.00417°/s，即$\dot{\theta}=0.00417°/s$ ，用弧度表示即有$\dot{\theta}=7.27 \times 10^{-5}rad/s$。以天体恰好运行到天线指向处的时间为开始时间（$t=0$），则有$cos(2\pi fbsin\theta/c)=cos(2\pi fbsin(7.27 \times 10^{-5}t)/c)$，可以看出，接收机输出的总功率会随时间而上下震荡，这种随震荡被称为条纹（类比于光学条纹干涉实验中明暗相间的条纹，这里是时域上的条纹，考虑到$\theta=\dot{\theta}t$，$\Delta L=bsin\theta$，实际上是$\Delta L$ 上的条纹，也相当于空间上的条纹分布了）。

不考虑底噪$2A_0^2$ 项（即扣除），分别用长度为 3m、9m、15m 的基线，在 1425MHz 处观测射电源时，在射电源运行到天线指向处的前后半小时内，输出总功率随时间的变化图如下：

# 考虑源的角直径对输出功率的影响

实际上的射电源不是电源，而是有一定亮度分布，夸张来讲是类似于太阳和月亮这种面源。考虑面源的角直径为$\alpha$ 且亮度是均匀分布的，把这个因素加入到干涉仪公式 1 则有

$P=AA^*=2A_0^2\alpha^2[1+sinc(\pi f \alpha b/c)cos(2\pi fbsin\theta/c)] (干涉仪公式2)$

其中$sinc(\pi f \alpha b/c)$ 描述了射电源大小，$cos(2\pi fbsin\theta/c)$ 描述了条纹变化率

注：本文的 sinc (x) 定义为 sinx/x，下文都将使用这种定义（代码中可能例外，可能使用 sinc (x)=sin (pi x)/(pi x) 这种定义）

太阳的角直径约为 0.5°，换算为弧度即有$\alpha = 0.5 \times \frac{\pi}{180} \approx 0.00873 \text{ rad}$。不考虑底噪$2A_0^2\alpha^2$ 项（即扣除），分别用长度为 3m、9m、15m 的基线，在 1425MHz 处观测太阳，在射电源运行到天线指向处的前后半小时内，输出总功率随时间的变化图如下：

# 考虑天线方向图对输出功率的影响

类似于人的耳朵，天线对不同方向的响应（或者说灵敏度）不一样，相当于对不同方向会赋予不同的权重。假设天线对不同方向的权重呈高斯分布，即高斯波束。

（此处可略过，如果感觉复杂的话）对于面天线而言，当用口径为$D$ 的面天线观测波长为$\lambda$ 的电磁波时，其半功率波束宽度近似为$HPBW=70\lambda/D$，增益近似为$G=10log_{10}(D^2/\lambda^2)$，即有$D/\lambda=\sqrt{10^{(G/10)}}$，$HPBW=70/\sqrt{10^{(G/10)}}$。对于在 1425MHz 有 15dBi 增益的八木天线（换算为线性增益为 31.62），可以得到其半功率波束宽度约为$70/ \sqrt{31.6}=12.45°$。即设置高斯函数时，其在$12.45°/2=6.225°$（半功率波束宽度的一半）处的权重相对于 0° 处的权重已经衰减了约一半，对应的弧度为 0.109 rad。

对于高斯函数而言，可以将标准差设置为 0.109，即可保证在 0.109 处衰减约一半，即有

$W(\theta)=e^{-\frac{(\theta-\mu)^2}{2\sigma}}=e^{-\frac{\theta^2}{2*0.109}}$，如下图所示

把$W(\theta)$ 加入到干涉条纹公式 2，可以得到

$P=AA^*=2A_0^2\alpha^2e^{-\frac{\theta^2}{2*0.109}}[1+sinc(\pi f \alpha b/c)cos(2\pi fbsin\theta/c)] (干涉仪公式3)$

不考虑底噪$2A_0^2\alpha^2$ 项（即扣除），分别用长度为 3m、9m、15m 的基线，在 1425MHz 处观测太阳，在射电源运行到天线指向处的前后半小时内，输出总功率随时间的变化图如下：

# 二元加法干涉仪仿真代码

import numpy as np

import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']

plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像时负号 '-' 显示为方块的问题

# 1、2、3 代表考虑（点源）、（太阳）、（太阳 + 天线方向图）的条纹仿真

simulator = 3

# 仿真的基线长度列表

b_list = [3,9,15]

# 仿真的时间范围

t = np.linspace(-1800, 1800, 1000)

# 仿真的其他参数

f = 1425e6  # 频率（赫兹）

b =  3 # 基线长度（米）

c = 299792458  # 光速 (米 / 秒)

alpha = 0.00873 # 太阳的角直径（弧度）

# 绘制图形

plt.figure(figsize=(10, 8))

for b in b_list:

    # 条纹率项

    rate_term = np.cos(2 * np.pi * f * b * np.sin(7.27e-5 * t) / c)

    # 射电源尺寸因素项

    def sinc(x):

        return np.where(x == 0, 1, np.sin(x) / x) # 定义 sinc 函数

    size_term = sinc(np.pi * f * alpha * b / c)

    # 定义高斯方向图

    def gaussian(x, mean, stddev):

        return np.exp(-0.5 * ((x - mean) / stddev) ** 2)

    gain_term = gaussian(7.27e-5 * t,0,0.109)

    if str(simulator) == '1':

        # 计算函数值

        y = rate_term

    elif str(simulator) == '2':

        # 计算函数值

        y = rate_term*size_term

    elif str(simulator) == '3':

        y = rate_term*size_term*gain_term

    plt.plot(t, y,label=f'b={b} m')

plt.title('不同基线长度下的条纹率',fontsize=20)

plt.xlabel('时间 (s)',fontsize=18)

plt.ylabel('扣除底噪的总功率',fontsize=18)

plt.xticks(fontsize=15)

plt.yticks(fontsize=15)

plt.legend(fontsize=15)

plt.xlim(-1800, 1800)

plt.ylim(-1.5, 1.5)

plt.show()

# 二元干涉仪综合孔径成像

# 两天线接收到的总电压值

两天线配置一致，射电源的强度分布为$I(\theta)$，天线对不同方向的响应函数为$W(\theta)$，再考虑不同方向的电磁波入射到两天线上会有几何时延差导致的相位差，导致的幅度变化可归为$e^{i2\pi f\frac{b}{c}sin\theta}$ 项，则天线指向射电源中心角度$\theta_c$ 时，输出功率为

$P(\theta_c)=[W(\theta_1)*I(\theta_1)*e^{i2\pi f\frac{b}{c} sin\theta_1}+W(\theta_2)*I(\theta_2)*e^{i2\pi f\frac{b}{c} sin\theta_2}+\\W(\theta_3)*I(\theta_3)*e^{i2\pi f\frac{b}{c} sin\theta_3}+\dots+W(\theta_n)*I(\theta_n)e^{i2\pi f\frac{b}{c} sin\theta_n}]*\Delta\theta=\\=\int_{\theta_c-\pi}^{\theta_c+\pi}{W(\theta)I(\theta)e^{i2\pi f\frac{b}{c}sin\theta}d\theta} (干涉仪公式4)$

# 宽天线波束观测均匀面源得到总电压值表达式

设射电源为理想的强度均匀分布的面源，天线指向射电源的中心位置，此时角度为$\theta_c$，其角直径为$\alpha$，当$|\theta-\theta_c|\le \alpha/2$ 时，则有$I(\theta)=I_0$，当$|\theta-\theta_c|>\alpha/2$ 时，有$I(\theta)=0$。

如果天线的增益不高，相当于$W(\theta)$ 在$\theta_c$ 附近足够平坦，即当$|\theta-\theta_c|\le \alpha/2$ 时，可以认为$W(\theta)$ 为常数，即可认为是$A$，带入干涉仪公式 4 可以推导出

注：设置新变量为$\theta^{'}=\theta-\theta_c$

$P(\theta_c)=\int_{\theta_c-\alpha/2}^{\theta_c+\alpha/2}{W(\theta)I(\theta)e^{i2\pi f\frac{b}{c}sin\theta}d\theta}\\=\int_{-\alpha/2}^{\alpha/2}{W(\theta^{'})I(\theta^{'})e^{i2\pi f\frac{b}{c}sin(\theta'+\theta_c)}d\theta^{'}}\\=AI_0\int_{-\alpha/2}^{\alpha/2}{e^{i2\pi f\frac{b}{c}sin(\theta'+\theta_c)}d\theta^{'}}(干涉仪公式5)$

在$-\alpha/2$ 到$\alpha/2$ 的范围内（太阳的角直径$\alpha$ 仅为 0.5°），$\theta^{'}$ 足够小，利用$sin(\theta^{'}+\theta_c)=sin\theta^{'}cos\theta_c+cos\theta^{'}sin\theta_c=\theta^{'}cos\theta_c+sin\theta_c$，则有

$P(\theta_c)=AI_0e^{i2\pi f\frac{b}{c}sin\theta_c} \int_{-\alpha/2}^{\alpha/2}{e^{i2\pi f\frac{b}{c}\theta^{'}cos\theta_c}d\theta^{'}}\\=AI_0\alpha e^{i2\pi f\frac{b}{c}sin\theta_c}\frac{sin(\pi f\frac{b}{c}\alpha cos \theta_c)}{\pi f\frac{b}{c}\alpha cos \theta_c}\\=AI_0\alpha e^{i2\pi f\frac{b}{c}sin\theta_c} sinc(\pi f\frac{b}{c}\alpha cos \theta_c)(干涉仪公式6)$

如果射电源运行到最高处，即在中天位置附近， 即$\theta_c$ 足够小，$cos(\theta_c)$ 此时接近 1，于是不考虑底噪项时，干涉仪公式 6 等价于公式 3（只取复数的实部）。也可以看出$bcos\theta_c$ 是基线长度$b$ 在垂直于射电源指向方向的投影，是实际上的等效基线长度。

# 射电源强度分布和天线输出电压是一对傅里叶变化

由干涉仪公式 5 也可见（仅认为$W(\theta)$ 为常数$A$，强度分布$I(\theta)$ 不认为是常数）：

$P(\theta_c)=A\int_{-\alpha/2}^{\alpha/2}{I(\theta^{'})e^{i2\pi f\frac{b}{c}sin(\theta^{'}+\theta_c)}d\theta^{'}}$

由于射电源角直径非常小，因此$sin(\theta^{'}+\theta_c)=\theta^{'}cos\theta_c+sin\theta_c$

$P(\theta_c)=A\int_{-\alpha/2}^{\alpha/2}{I(\theta^{'})e^{i2\pi f\frac{b}{c}(\theta^{'} cos\theta_c+sin\theta_c)} d\theta^{'}}\\=Ae^{i2\pi f\frac{b}{c}sin\theta_c} \int_{-\alpha/2}^{\alpha/2}{I(\theta^{'})e^{i2\pi f\frac{b}{c}\theta^{'} cos\theta_c} d\theta^{'}}$

最后是为了获得强度的分布，而不是相位分布，因此把$P(\theta_c)e^{-i2\pi f\frac{b}{c}sin\theta_c}$ 记为$P'(\theta_c)$，显然$P(\theta_c)$ 和$P'(\theta_c)$ 的绝对值是一样的（而且我们已知$b$ 和$\theta_c$，也可以直接把这个相位项消除），即有：

$P'(\theta_c)=P(\theta_c)e^{-i2\pi f\frac{b}{c}sin\theta_c}=\int_{-\alpha/2}^{\alpha/2}{I(\theta^{'})e^{i2\pi f\frac{b}{c}\theta^{'} cos\theta_c} d\theta^{'}}$

其中$b cos(\theta_c)$ 为基线相对于射电源方向的投影，因此可将$f\frac{b}{c}cos(\theta_c)=\frac{b}{\lambda}cos(\theta_c)$ 记为$b_\lambda$，则有：

$P'(b_\lambda)=A\int_{-\alpha/2}^{\alpha/2}{I(\theta^{'})e^{i2\pi b_{\lambda} \theta^{'}}d\theta^{'}}(干涉仪公式8)$

可见干涉仪输出的电压值$E$ 和射电源的强度分布$I$ 是一对傅里叶变化，干涉仪公式 8 的逆变换为：

$I(\theta') = \int_{-\infty}^{\infty} P'(b_\lambda) e^{-i2\pi \theta^{'} b_\lambda} d b_\lambda(干涉仪公式9)$

因此我们可以利用地球自转（$\theta_c$ 变动）或者主动改变基线长度（$b$ 变动），从而改变$b_\lambda$，测量在不同$b_{\lambda}$ 下的输出电压，然后进行傅里叶级数叠加，即可得到射电源的强度分布，例如利用长度为 3m、6m、9m... 24m 一系列基线在 1425MHz（21cm）每天测量一次中天的射电源，可以得到

$I(\theta_c)=(P_{3m}e^{i2\pi*3\theta^{'}/0.21}+ P_{6m}e^{i2\pi*6\theta^{'}/0.21} +P_{9m}e^{i2\pi*9\theta^{'}/0.21}+...+P_{24m}e^{i2\pi*24\theta^{'}/0.21})*\Delta b_\lambda \\=(P_{3m}e^{i2\pi*3\theta^{'}/0.21}+ P_{6m}e^{i2\pi*6\theta^{'}/0.21} +P_{9m}e^{i2\pi*9\theta^{'}/0.21}+...+P_{24m}e^{i2\pi*24\theta^{'}/0.21})*3/0.21$

# 综合孔径实际观测考量

干涉仪公式 8 和公式 9 认为天线的方向图在指向射电源时是平坦的，即对输出功率毫无影响，然而如果利用地球自转而不跟踪射电源的话，射电源运行到不同方位时，天线对电磁波的响应权重肯定也变了，因此需要了解天线的详细方向图才能进行没法直接测量方向图。

然而业余观测时往往不知道方向图，因此放弃靠地球自转来进行综合孔径成像。可以固定天线一直指向射电源中天的位置，只观测中天的强度，由于天线方向图对此方向的响应权重肯定是不变的，因而无需考虑天线方向图的影响，视为常数$A$，在此情况下可利用干涉仪公式 8 和公式 9。

干涉仪的最小分辨角约为$\theta_{min}=\frac{\lambda}{b}$，即和观测波长成正比，和基线长度成反比。太阳的角直径为 0.5°，也就是 0.00873 rad，为使得二元干涉仪能够观测太阳，一系列基线长度之间的间距不能超过$\Delta b_{\lambda}=\lambda/0.00873=115\lambda$，比如在 1420MHz（21cm）观测，则间距不能超过 25m。我们把这个间距设为 3m，利用长度为 3m、6m、9m... 24m（手动改变基线长度）一系列基线每天测量一次中天的射电源强度，然后得到

$I(\theta_c)=(P_{3m}e^{i2\pi*3\theta^{'}/0.21}+ P_{6m}e^{i2\pi*6\theta^{'}/0.21} +P_{9m}e^{i2\pi*9\theta^{'}/0.21}+...+P_{24m}e^{i2\pi*24\theta^{'}/0.21})*3/0.21$

取其实部，即可得到太阳的一维强度分布，也就是所谓的综合孔径合成，在此未使用 uv 覆盖的概念，实际上是等价的。

当天线指向射电源中天位置进行观测时，干涉仪公式 6 中的$sinc(\pi f\frac{b}{c}\alpha cos \theta_c)=sinc(\pi f\frac{b}{c}\alpha)$ 即在理论上描述了使用不同长度基线观测到的幅度变化，可以此作为权重，编写代码如下，对在 1425MHz，用长度为 3m、6m、9m... 24m 长度的基线对太阳（角直径为 0.00873 rad）进行综合孔径结果进行仿真。结果如下图所示，可见其半功率衰减角为 0.00353*2=0.00706 rad，和理论的太阳大小相近。

为了进一步验证理论和程序的可靠性，尝试仿真 0.001 rad 角直径的射电源，仍然在 1425MHz 进行观测，一系列基线长度之间的间距不能超过$\Delta b_{\lambda}=\lambda/0.0001=1000\lambda$，即 210m，我们把这个间距设为 12m，利用长度为 12m、24m、36m... 210m（手动改变基线长度）一系列基线每天测量一次中天的射电源强度，然后得到的结果如下图，可见其半功率衰减角为 0.000415*2=0.000830 rad，和理论的射电源大小相近

# 综合孔径仿真代码

import numpy as np

import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']

plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像时负号 '-' 显示为方块的问题

# 定义参数

alpha = 0.00873 #角分辨率

#alpha = 0.001 #角分辨率

f = 1425e6  # 频率

c=3e8

lam = c/f # 波长

# 定义 b 值，除以波长

b_lambda_list = np.arange(3, 24, 3)/lam  # b = 3, 6, 9, ..., 24

#b_lambda_list = np.arange(12, 210, 12)/lam  # b = 12, 24, 36, ..., 210

# 计算不同 b 下的权重

weights = []

for b_lambda in b_lambda_list:

    weights.append(np.sinc(b_lambda * alpha))

weights = np.array(weights) / np.max(weights)  # 归一化权重

# 定义 theta' 的取值范围

theta_prime = np.linspace(-0.005, 0.005, 1000)

# 计算傅里叶级数和

result = np.zeros_like(theta_prime, dtype=complex)

for i, b_lambda in enumerate(b_lambda_list):

    result += weights[i] * np.exp(1j * 2 * np.pi * b_lambda * theta_prime)

result *= b_lambda_list[1]-b_lambda_list[0]  # 乘以常数

result = result/max(result)

# 绘图

plt.figure(figsize=(12, 8))

# 创建第一个子图

plt.subplot(2, 1, 1) 

plt.plot(b_lambda_list*lam,weights, marker='.',markersize=10, label='强度的角分布')

plt.xticks(fontsize=15)

plt.yticks(fontsize=15)

plt.xlabel('基线长度 (m)',fontsize=18)

plt.ylabel('归一化强度值',fontsize=18)

plt.legend(fontsize=15)

plt.title(f'各长度基线观测射电源输出功率的变化',fontsize=20)

plt.grid()

# 创建第二个子图

plt.subplot(2, 1, 2) 

plt.plot(theta_prime, np.real(result), label='强度的角分布')

# plt.plot(theta_prime, np.imag(result), label='Imag Part')

# plt.plot(theta_prime, np.imag(result), label='Imaginary Part', linestyle='dashed')

plt.title(f'{alpha} rad 角直径的射电源',fontsize=20)

plt.xlabel('角度 (rad)',fontsize=18)

plt.ylabel('归一化强度值',fontsize=18)

plt.xticks(fontsize=15)

plt.yticks(fontsize=15)

plt.legend(fontsize=15)

plt.xlim(-alpha*3, alpha*3)

# plt.ylim(-1.5, 1.5)

# 找到最接近 0.5 的值及其索引

target_value = 0.5

index = (np.abs(np.real(result) - target_value)).argmin()  # 找到最接近 0.5 的值的索引

closest_value = np.real(result)[index]  # 接近 0.5 的最小值

corresponding_theta = theta_prime[index]  # 对应的 theta_prime 值

plt.text(corresponding_theta, 0.5, f'({corresponding_theta:.6f},0.5)',fontsize=20)

plt.scatter(corresponding_theta, 0.5,color='r')

# plt.axhline(0.5, color='black', linewidth=0.5, linestyle='--')

# plt.axvline(0, color='black', linewidth=0.5, linestyle='--')

plt.grid()

plt.legend(fontsize=15)

plt.tight_layout()

plt.show()