pSolve 配套工具 src_alignment 的使用与实现

· 预计阅读 9 分钟

src_alignment 是 pSolve 的一个配套工具程序。它读取先验星表、参考星表，以及可选的解算所得星表，比较它们的参考架坐标轴指向，并生成可用于后续全局解算的 NO_NET_ROTATION_SOURCE 配置。程序同时输出整体旋转参数的拟合结果和约束方程的右端项，两组数各有不同的用途。

一、程序结构

目录中有三个文件。src_alignment.f 包含主程序 SRC_ALIGNMENT_MAIN 和两个计算子程序 SRC_ALIGNMENT、ROTANG_SOU。src_alignment_version.i 设置版本号，makefile 定义编译和链接规则。程序直接链接 utils/getpar/read_sou.o，因此星表格式、源名和坐标误差的读取要结合 getpar 中的 READ_SOU 分析。

主程序先读取 pSolve 全局解算控制文件，从 VTD_CONF 找到 VTD 控制文件，再从其中的 SOURCE_COORDINATES: 找到先验星表。读取各份星表并确定选源方式后，程序进入 SRC_ALIGNMENT，匹配源、计算坐标偏差，最后生成约束文件。具体流程如下图所示。

flowchart TD
    M["SRC_ALIGNMENT_MAIN 解析命令行和控制文件"] --> R["READ_SOU 读取星表"]
    M --> L["根据观测数或已有源列表确定选源方式"]
    R --> S["SRC_ALIGNMENT 匹配源并累加约束右端项"]
    L --> S
    S --> D["ROTANG_SOU 拟合三组整体旋转参数"]
    D --> O["SRC_ALIGNMENT 填写比较结果和 RIGHT_PART"]
    O --> W["主程序调用 WR_TEXT 写出约束文件"]

SRC_ALIGNMENT 先累加约束右端项，再调用 ROTANG_SOU 做三次整体旋转参数拟合。虽然 RIGHT_PART 的数值独立于拟合结果，只要其中一次拟合失败，程序仍会提前返回，约束文件也不会写出。

二、命令行参数与输入输出

src_alignment generic_solution_name ref_sou_file numobs_min [ZERO] [ivrb]

第一个参数是解算文件的路径前缀，即解算标识。例如，给出 nju_2026b_zero 后，程序会寻找 nju_2026b_zero.cnt、nju_2026b_zero.sou 和 nju_2026b_zero.cns。第二个参数是参考星表的路径。三份星表分别承担以下作用：

  • 先验星表 CAT_APR。从 <前缀>.cnt 中的 VTD_CONF 找到 VTD 配置，再读取其 SOURCE_COORDINATES: 指定的星表。这份星表提供全局解估计坐标改正量的参考值。
  • 参考星表 CAT_REF。由命令行第二个参数指定，提供目标坐标。程序还用这份星表的赤经、赤纬计算旋转模型和约束方程的系数。
  • 解算所得星表 CAT_EST。来自 <前缀>.sou，通常由 getpar 从解算输出报告中提取，提供坐标估计值、各分量的形式误差以及实际采用的观测数 NOBS_USED。

第三个参数 numobs_min 控制选源方式。它取正整数时，程序要求存在 .sou，并从参考星表与解算所得星表的共同源中，选取 NOBS_USED >= numobs_min 的源。它取零时，程序先尝试从 .cnt 读取已有的约束源名单。如果 .cns 存在，就清空先前的名单，改从其中的 NNR_SRC: 段读取。最终名单为空时，程序回到 NOBS_USED >= 0 的筛选方式，逐个匹配参考源。

指定 ZERO 后，程序跳过 .sou，在内部令 CAT_EST = CAT_REF。未指定 ZERO 但 .sou 不存在时，程序在零阈值下也采用这一处理。正的观测数阈值需要解算所得星表，因此与 ZERO 同时使用会报错 7311。约束右端项仍由参考坐标与先验坐标之差计算，所以使用 ZERO 也可能得到非零的 RIGHT_PART。

ivrb 控制程序正常结束时的终端输出，默认值为 1。ivrb > 0 时，程序打印生成文件的前 20 行，再显示保存路径。ivrb <= 0 时，程序省略这段终端输出，仍计算并写出文件。当前实现只判断是否大于零，因此 1、2 等正整数的效果相同。错误信息由错误处理程序另外输出。未使用 ZERO 时，ivrb 放在第四个参数；使用 ZERO 时，放在第五个参数。例如，以下两条命令都显示正常结束时的输出。

src_alignment solution reference.src 0 1
src_alignment solution reference.src 0 ZERO 1

完整输出文件的路径形式为 /tmp/<SOL_ID>_net_sou_cns.cnt，其中 SOL_ID 是程序确定的解算标识。文件头记录命令、版本、输入星表、参考源数量，以及三组 Net rotation 和 Constraint 结果。其后是 NO_NET_ROTATION_SOURCE、RIGHT_PART 和参考源列表。

三、星表读取与定义源匹配

READ_SOU 根据星表首行识别格式，逐行调用 READ_SOU_LINE，把每条记录的 %IVS_NAME 放入名称数组。SRC_ALIGNMENT 从参考星表出发，在解算所得星表中找同名记录，再到先验星表中查找。三份星表之间采用射电源的 IVS 名称作等值比较。参考源在解算所得星表中找不到时，程序跳过该源。通过匹配和选源条件后，若在先验星表中找不到，程序报错 7351。若入选源在先验星表中出现多次，程序报错 7352。参考星表和解算所得星表的重名记录也需要检查，因为当前匹配循环可能将它们重复计入。

对成功匹配的源，程序形成三组坐标差。将各源的赤经、赤纬按相同顺序排列为坐标向量，可写为 \(\begin{aligned} \boldsymbol d_{\mathrm{RA}} &= \boldsymbol x_{\mathrm{ref}}-\boldsymbol x_{\mathrm{apr}},\\ \boldsymbol d_{\mathrm{EA}} &= \boldsymbol x_{\mathrm{est}}-\boldsymbol x_{\mathrm{apr}},\\ \boldsymbol d_{\mathrm{ER}} &= \boldsymbol x_{\mathrm{est}}-\boldsymbol x_{\mathrm{ref}}. \end{aligned}\) 变量后缀 RA、EA、ER 分别对应参考减先验、解算所得减先验、解算所得减参考。例如，D_RA_DEL 保存的是参考赤纬减先验赤纬。输出文件用 reference versus a priori、estimated versus a priori 和 estimated versus a reference 标记这三组比较。三组坐标差满足 $\boldsymbol d_{\mathrm{EA}}=\boldsymbol d_{\mathrm{RA}}+\boldsymbol d_{\mathrm{ER}}$。程序对它们使用相同的源、系数和拟合权重，因此输出的三组 Net rotation 和三组 Constraint 也分别满足这一相加关系，核对时需计入输出舍入的影响。

四、具体实现

设 $\boldsymbol d$ 表示上述任一组坐标差。对第 $i$ 个源，程序用参考星表中的 $(\alpha_i,\delta_i)$ 构造矩阵

\[A_i= \begin{pmatrix} -\cos\alpha_i\tan\delta_i & -\sin\alpha_i\tan\delta_i & 1\\ \sin\alpha_i & -\cos\alpha_i & 0 \end{pmatrix}.\]

这里的两个坐标分量是原始赤经差 $\Delta\alpha_i$ 和赤纬差 $\Delta\delta_i$，单位均为弧度。上述系数对应主动旋转约定 $\Delta\boldsymbol s=\boldsymbol\epsilon\times\boldsymbol s$。把各源的 $A_i$ 纵向排列为 $A$，程序分别计算整体旋转的加权最小二乘估计 $\widehat{\boldsymbol\epsilon}$ 和等权约束方程的三个线性组合 $\boldsymbol\psi$。

\[\begin{aligned} \widehat{\boldsymbol\epsilon} &= (A^\mathsf T W A)^{-1}A^\mathsf T W\boldsymbol d,\\ \boldsymbol\psi &= A^\mathsf T\boldsymbol d. \end{aligned}\]

第一式要求 $A^\mathsf T W A$ 可逆，其结果对应输出中的 Net rotation。第二式直接对坐标差按系数求和，其结果对应 Constraint。

4.1 数据加权

ROTANG_SOU 把每个源的赤经和赤纬作为两条方程，分别使用各自的形式误差。对第 $j$ 个坐标分量，法方程中的权重为

\[w_j=\frac{1}{s_j^2+f^2},\qquad f=0.10\ \mathrm{mas}.\]

$f$ 是旋转拟合的误差下限，在 SRC_ALIGNMENT 中直接赋值。

ERR_FLOOR = 0.10D0*MAS__TO__RAD

MAS__TO__RAD 将毫角秒转换为弧度。这一数值固定写在程序中，修改时需要改动源码并重新编译 src_alignment。

有 .sou 且未指定 ZERO 时，三次拟合全部使用 CAT_EST 的 ALP_ERR、DEL_ERR。即使拟合参考星表与先验星表之间的整体旋转，权重仍取自解算所得星表。各坐标分量的权重分别计算，$W$ 为对角矩阵，相当于将赤经、赤纬的误差视为不相关。以参考星表代替 CAT_EST 时，程序改用其中的分量误差。这里需要结合输入格式检查。READ_SOU_LINE 对旧式 SOU-MODFILE 等格式只读取综合位置误差，ALP_ERR 和 DEL_ERR 未被赋值。此时，程序会用尚未赋值的 ALP_ERR 和 DEL_ERR 计算旋转参数拟合的权重，因此得到的权重可能不可靠。

源码中的 WEI 为 $\sqrt{w_j}$。DIAD_CVT_S 接收 WEI**2 来累加法矩阵，NORVEC_UPD 在内部使用 WEI**2 来累加右端向量。程序还计算了旋转角的形式误差，即逆法矩阵三个对角元素的平方根，输出文件只记录旋转角的估计值。

4.2 约束右端项生成

生成配置时，程序写出 UNIFORM，把 PSI_RA 写入 RIGHT_PART，并用 NO EXCEPT 列出参考源。令 $A_{\mathrm{ref}}$ 表示用参考坐标构造的矩阵，右端向量及其第三个分量为

\[\begin{aligned} \boldsymbol r &= A_{\mathrm{ref}}^\mathsf T\left(\boldsymbol x_{\mathrm{ref}}-\boldsymbol x_{\mathrm{apr}}\right),\\ r_3 &= \sum_i\left(\alpha_{i,\mathrm{ref}}-\alpha_{i,\mathrm{apr}}\right). \end{aligned}\]

$r_3$ 是所选源的赤经差之和。生成的三条约束方程共用 SIGMA 1.00D-10，即标准差为 $10^{-10}\ \mathrm{rad}$。它控制后续全局解中的约束强度，而第 4.1 节的 $0.10\ \mathrm{mas}$ 用于本工具的整体旋转参数拟合，二者的含义不同。

Net rotation 和 RIGHT_PART 对应不同的线性运算。两份星表的实际坐标偏差通常不能完全用整体旋转来描述。将拟合出的整体旋转换算为各源的坐标变化，只能得到整体旋转对应的部分，无法还原包含拟合残差在内的全部坐标实际差。因此，用这些旋转引起的坐标变化生成右端项，也未必能还原程序直接用整组坐标偏差求得的结果。在后续全局解中，NNR_SRC 用先验坐标建立约束系数，并将读入的右端项施加到所选源的待估坐标分量上。当两份星表相近、参与的源和坐标分量一致时，一阶近似下可将生成右端项和施加约束时使用的矩阵视为同一个 $C$。记新解坐标为 $\boldsymbol x_{\mathrm{new}}$，则有

\[\begin{aligned} &\boldsymbol r = C\left(\boldsymbol x_{\mathrm{ref}}-\boldsymbol x_{\mathrm{apr}}\right),\\ &C\left(\boldsymbol x_{\mathrm{new}}-\boldsymbol x_{\mathrm{apr}}\right) \simeq \boldsymbol r,\\ &C\left(\boldsymbol x_{\mathrm{new}}-\boldsymbol x_{\mathrm{ref}}\right) \simeq \boldsymbol0. \end{aligned}\]

这里的 $\simeq$ 表示以有限标准差施加约束。这样生成的右端项，使新解算出的射电源表相对于参考星表满足无净旋转条件。