拟合一个真实 EHT 闭合三角形

先固定观测，再开始拟合

文件 SR1_M87_2017_100_lo_hops_netcal_StokesI.csv 保存 2017 年 4 月 10 日 M87 的真实 EHT-HOPS SR1 Stokes I 观测，采用低频段，文件头中心频率为 227.0707 GHz。“网络标定”是数据发布时经历的一步处理，并不表示各站残余大气相位已经变成已知天空相位。

原始文件为 194347 字节。程序先检查固定源文件的 SHA-256，再检查摘出的三行是否逐字节等于原文件第 1272、1274、1283 行，同时保留两行文件头。三条记录都是 UTC 04.90972209 小时。数据仓库以 PDDL 1.0 发布这些测量，教程读取器和拟合程序独立编写。

python3 fit-eht-triangle.py --data-dir . --output reproduced-eht.json
# 比较 reproduced-eht.json 与 eht-triangle-results.json
# 原始 CSV 的 SHA-256：
# 1dc51e3e31a840ec34de5b5ff4d8a6af0cbb73cd8c6db85e917366ed3ec7924a
# 三行摘录的 SHA-256：
# 62ffea8dadaae26aa7fb3788c0456c4625d0e13eaad61e57529ba20311d86502

对应阅读：EHT 2019-D01-01 M87 data release, pinned commit 27ae76a

读清同一时刻的三个站点

AA 是 ALMA，AZ 是亚利桑那的亚毫米望远镜，LM 是墨西哥的大毫米望远镜。下面各列依次为 UTC 小时、第一站、第二站、东向空间频率、北向空间频率、振幅、相位和热噪声振幅不确定度。空间频率已经用波长表示，振幅和不确定度为央斯基，相位用度。

发布者依次保存 AA 到 AZ、AA 到 LM、AZ 到 LM。要沿 AA 到 AZ 到 LM 再回 AA，就得把 AA 到 LM 反向，取复共轭：振幅不变，相位变号。导出的空间频率相加只差几千个波长，远小于十亿波长的基线。保留这点舍入差，不必修改原始数值。

04.90972209,AA,AZ,3028965632.0000,-4496652800.0000,0.10473071,179.9880,0.00501226
04.90972209,AA,LM,2296739584.0000,-3319742720.0000,0.07996017,59.4402,0.00309308
04.90972209,AZ,LM,-732228224.0000,1176914944.0000,0.51672293,-62.7827,0.04022401

对应阅读：EHT 2019-D01-01 M87 data release, pinned commit 27ae76a · EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

计算相位和热噪声近似

把 AA–AZ 与 AZ–LM 相位相加，再减去保存的 AA–LM 相位，折回一个周期，得到 $57.7651^{\circ}$ 。程序还独立构造三个复数，取有向乘积的辐角，两个算法给出相同结果。

信噪比较高时，一点垂直于信号的复噪声，会带来大约“噪声除以振幅”的相位扰动。假设三条热误差独立，把相位方差相加。这里最低振幅信噪比为 $12.85$ ，闭合相位不确定度约 $5.68545^{\circ}$ 。这项估计没有纳入振幅标定误差、相关误差和不能闭合的系统效应。

\Phi=\operatorname{wrap}(179.9880^{\circ}-62.7827^{\circ}-59.4402^{\circ})=57.7651^{\circ}

最后一项取减号，因为三角形经过这条边的方向相反，与它本身正负无关。

\sigma_\Phi\simeq\left[\sum_{e=1}^3\left(\frac{\sigma_e}{A_e}\right)^2\right]^{1/2}=0.09923\,\mathrm{rad}

传播误差时先用弧度，最后再换成度。这是当前一个三角形的高信噪比热噪声近似。

对应阅读：EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

选一个能表示不对称的经验模型族

用两个圆高斯分量，半高全宽都为 $w$ 。主分量放原点，次分量与主分量通量比为 $r$ ，位移为 $(d_l,d_m)$ 。两者间距记为 $d$ ，角度 $\psi$ 从东向北量。高斯标准差 $s=w/\sqrt{8\ln2}$ ，换成弧度后放进指数。

归一化圆高斯的傅里叶变换仍是高斯。给双分量相量乘上这个包络，便得到下式。总通量 $F$ 乘所有振幅，却从闭合相位里抵消。整体位置也不能由闭合相位测出，故把主分量固定在原点，省去这里无法辨认的平移。

V(u,v)=F e^{-2\pi^2s^2(u^2+v^2)}\frac{1+r e^{-2\pi i(ud_l+vd_m)}}{1+r},\qquad(d_l,d_m)=d(\cos\psi,\sin\psi)

共有五个拟合量：间距、方向角、通量比、共同宽度和总通量。它是用来考察信息量的经验模型，不据此断言 M87 真的由两个这样的发光体组成。

对应阅读：NRAO Essential Radio Astronomy, Chapter 3, aperture synthesis

通量先解析求解，再搜索形状

给定一个候选形状，用 $q_e$ 表示第 $e$ 条边在单位总通量下的可见度振幅。振幅残差是 $A_e-Fq_e$ 。除以 $\sigma_e$ 、平方、求和，对 $F$ 求导，便得到一条线性方程和下面的加权解。数值搜索只需要处理四个形状参数。

程序先搜索十到一百微角秒的间距、每五度一个方向、六种通量比和五种初始宽度。选出八个彼此分开的起点，再逐次缩小坐标步长。输出保存了搜索设置，可以原样重跑。有限搜索足以展示多个好解，不保证找到全局最优。

\widehat F=\frac{\sum_e A_eq_e/\sigma_e^2}{\sum_e q_e^2/\sigma_e^2}

候选信号非零时，分母为正。这个通量随当前源模型一起拟合；独立零基线测量会提供另一项信息。

Q=\sum_e\frac{(A_e-\widehat Fq_e)^2}{\sigma_e^2}+\frac{\operatorname{wrap}(\Phi_{\mathrm{model}}-\Phi_{\mathrm{obs}})^2}{\sigma_\Phi^2}

相位差要先折回一周，避免把周期两端本来很近的角度当成很远。目标值采用上面明确说明的近似热误差。

比较真正保存下来的拟合结果

一个细化结果的间距为 $60.16\,\mu\mathrm{as}$ ，方向 $240.00^{\circ}$ ，通量比 $0.596$ ，宽度 $23.56\,\mu\mathrm{as}$ ，总通量 $0.67307\,\mathrm{Jy}$ 。按原始行顺序，预测振幅为 $(0.104679,0.080008,0.516147)\,\mathrm{Jy}$ ，闭合相位 $57.7359^{\circ}$ ，近似目标值 $0.000575$ 。

另一个解的间距却是 $97.44\,\mu\mathrm{as}$ ，方向 $230.00^{\circ}$ ，通量比 $0.592$ ，宽度 $24.24\,\mu\mathrm{as}$ 。目标值仍只有 $0.000681$ 。两者都远比热误差要求得更贴近观测。三条振幅加一个闭合角，却拟合五个自由量，小残差当然会留下很多选择。

若强制闭合相位为零，仅相位一项便约为 $103.2$ 。这个角确实携带很强的不对称信息，同时容许许多不同的不对称形状。这里不报告约化卡方，因为模型参数数目已经超过四项数据条件。

把实测三角形与合成成像分清楚

实测计算没有把发布的原始站点相位直接做逆变换，再宣称那就是天空图像。它拟合的是能避开未知站点相位的振幅和有向闭合相位。真正的二维成像还需要更多三角形、更多标定信息，并处理测量之间共用的误差。

交互阵列则采用已知的有限像素合成天空，另外保存独立理想标定源和总功率值。这个受控实验允许我们比较脏图、两种正则化反演以及单独保存的已知源。两部分相互补充：一部分读懂发布者实际测到了什么，另一部分在真值已知时研究稀疏采样怎样影响重建。

对应阅读：EHT 2019-D01-01 M87 data release, pinned commit 27ae76a · EHT Collaboration, First M87 Results IV, arXiv:1906.11241v1

换一个三角形时，留下哪些东西

保存发布版本、原始行、站点顺序、时间定义、频率和误差列。先画有向路径，再加相位；确认三条边来自同次积分，按明确的误差模型传播不确定度。比较不止一种候选形状，并把预测信号与拟合参数一起保存。下一条基线就能检验这些候选到底在哪里不同。