pSolve全局解算中天球参考架的无净旋转约束

2023-10-06 · 预计阅读 15 分钟

编写 pSolve 全局解算的配置文件时，需要说明如何确定解算所得天球参考架的坐标轴指向。射电源坐标、台站位置和地球定向参数同时参与估计时，若整体旋转的自由度缺少约束，法方程就会在相应方向上亏秩。这里结合 pSolve 的代码，说明如何通过无净旋转约束确定天球参考基准。

一、天球无净旋转约束

1.1 基本概念

天球参考架的无净旋转约束（no-net-rotation，NNR）是对一组选定射电源的坐标改正量施加三条约束，使这组源相对于参考星表的整体旋转接近零，从而确定解算所得天球参考架的坐标轴指向。参考源通常选自某代 ICRF 的定义源。每个源的坐标仍然可以改变，约束的是整组坐标改正量中的旋转分量。在 pSolve 中，这三条约束方程并非等式约束，而是具有有限的标准差。

设第 $i$ 个源的先验坐标为 $(\alpha_i,\delta_i)$，待估坐标改正量为 $\Delta\boldsymbol{x}_i=(\Delta\alpha_i,\Delta\delta_i)^\mathsf T$。这里用的是原始赤经改正量 $\Delta\alpha_i$，两个坐标分量都以弧度表示。设 $\boldsymbol s_i$ 是源的方向单位矢量，$\boldsymbol\epsilon=(\epsilon_x,\epsilon_y,\epsilon_z)^\mathsf T$ 是微小旋转向量。采用主动旋转约定，保留一阶项，源方向及其坐标的变化为

\[\begin{aligned} \Delta\boldsymbol s_i &= \boldsymbol\epsilon\times\boldsymbol s_i,\\ \Delta\boldsymbol x_i &= \underbrace{\begin{pmatrix} -\cos\alpha_i\tan\delta_i & -\sin\alpha_i\tan\delta_i & 1\\ \sin\alpha_i & -\cos\alpha_i & 0 \end{pmatrix}}_{A_i}\boldsymbol\epsilon. \end{aligned}\]

1.2 等权情况（均匀权）

这里的等权指在拟合整体旋转时，给各源的赤经、赤纬改正量相同的权重。将各源的改正量 $\Delta\boldsymbol x_i$ 和矩阵 $A_i$ 按相同的源顺序分别纵向叠放，得到 $\Delta\boldsymbol x$ 和 $A$。若 $A^\mathsf T A$ 可逆，整体旋转的等权最小二乘估计为

\[\widehat{\boldsymbol\epsilon}=(A^\mathsf T A)^{-1}A^\mathsf T\Delta\boldsymbol x.\]

如果以先验星表为天球参考基准，就要求按上式拟合得到的整体旋转接近零，因此对坐标改正量施加约束 $A^\mathsf T\Delta\boldsymbol x\simeq\boldsymbol0$。更一般地，把三条约束方程的右端项记为 $\boldsymbol r=(r_1,r_2,r_3)^\mathsf T$，则有

\[\begin{aligned} \sum_i\left[-\cos\alpha_i\tan\delta_i\,\Delta\alpha_i+\sin\alpha_i\,\Delta\delta_i\right] &\simeq r_1,\\ \sum_i\left[-\sin\alpha_i\tan\delta_i\,\Delta\alpha_i-\cos\alpha_i\,\Delta\delta_i\right] &\simeq r_2,\\ \sum_i\Delta\alpha_i &\simeq r_3. \end{aligned}\]

若从先验星表到目标星表的坐标变化是一个主动旋转 $\boldsymbol\epsilon_0$，将相应的坐标变化 $A\boldsymbol\epsilon_0$ 代入三式，便得到右端项 $\boldsymbol r=A^\mathsf T A\,\boldsymbol\epsilon_0$。例如，绕 $z$ 轴旋转 $\epsilon_z$ 时，各源的赤经都增加 $\epsilon_z$。若有 $n$ 颗源的赤经参与约束，第三式的右端项就是 $r_3=n\epsilon_z$。这里的 $\simeq$ 表示约束方程允许有残差，第 1.4 节讨论约束标准差如何控制这些残差的权重。

1.3 非等权情况

各射电源的观测次数和观测几何有所不同，即使是同一颗源，赤经和赤纬的估计精度也可能不同。因此，在确定参考架坐标轴指向时，可以对不同坐标分量分别赋权。

设 $c_{kj}^{(0)}$ 为等权情况下第 $k$ 条约束方程中第 $j$ 个坐标参数的系数。pSolve 取当前法矩阵中与该参数对应的对角元素 $N_{jj}$ 作为权重，将原系数乘以 $N_{jj}$，再对第 $k$ 行的系数作归一化。对于非零的约束行，最终系数为

\[c_{kj}=\frac{c^{(0)}_{kj}N_{jj}}{\sqrt{\displaystyle\sum_\ell\left(c^{(0)}_{k\ell}N_{\ell\ell}\right)^2}}.\]

分母是加权后第 $k$ 行系数向量的欧氏范数，求和遍及该行涉及的所有坐标分量。

采用 $N_{jj}$ 加权，可以从最小二乘目标函数来理解。若固定其他参数，只让第 $j$ 个参数改变 $t$，目标函数中相应的二次项为 $N_{jj}t^2$。$N_{jj}$ 越大，目标函数沿这个参数方向的曲率就越大。联立估计所有参数时，最终的参数方差还取决于参数之间的相关性。

1.4 约束强度

将前面三条约束方程的系数按行排列，得到约束矩阵 $C$。每行对应一条约束方程，每列对应一个待估坐标分量。等权时，$C=A^\mathsf T$。非等权时，$C$ 的元素是第 1.3 节给出的加权系数 $c_{kj}$。三条方程采用相同的标准差 $\sigma$，约束及其对总 $\chi^2$ 的贡献可以写为

\[\begin{aligned} C\Delta\boldsymbol x &\simeq \boldsymbol r,\\ \chi^2_{\mathrm{constraint}} &= \frac{1}{\sigma^2}\left\|C\Delta\boldsymbol x-\boldsymbol r\right\|_2^2. \end{aligned}\]

$\sigma$ 越小，同样的偏差对 $\chi^2$ 的贡献越大，约束就越紧。$\sigma$ 的取值取决于研究目的。例如，Titov & Lambert（2013）在研究河外射电源视自行场中的银河系光行差信号时，先在各次观测实验的解算中对射电源坐标施加较松的无净旋转约束，再由所得的坐标时间序列估计各源的视自行。此前，Titov et al.（2011）的对比试验表明，收紧这一约束后，拟合得到的银河系光行差向点偏离了银心方向。

1.5 约束的结果

讨论全局解算所得天球参考架的坐标轴指向时，需要区分星表的三种作用。以下各坐标向量均按相同的顺序排列所选参考源的坐标分量。

先验星表提供解算使用的先验坐标 $\boldsymbol x_{\mathrm{apr}}$。
本次全局解算得到的星表（后验星表）给出坐标估计值 $\boldsymbol x_{\mathrm{new}}$。
参考星表提供确定天球参考基准所用的坐标 $\boldsymbol x_{\mathrm{ref}}$，例如某代 ICRF 的定义源坐标。

在全局解中，我们估计的是相对于先验坐标的改正量 $\Delta\boldsymbol x=\boldsymbol x_{\mathrm{new}}-\boldsymbol x_{\mathrm{apr}}$。若以先验星表为天球参考基准，就要求这组改正量中的整体旋转分量接近零，此时取 $\boldsymbol r=\boldsymbol0$，即

\[C\left(\boldsymbol x_{\mathrm{new}}-\boldsymbol x_{\mathrm{apr}}\right)\simeq\boldsymbol0.\]

若采用另一份参考星表作为基准，则先根据参考坐标与先验坐标之差计算右端项，再在全局解中施加相应的约束。两步分别写为

\[\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. \end{aligned}\]

将第一式代入第二式并消去先验坐标，得到

\[\boxed{C\left(\boldsymbol x_{\mathrm{new}}-\boldsymbol x_{\mathrm{ref}}\right)\simeq\boldsymbol0}.\]

因此，非零右端项补偿了参考星表与先验星表之间的整体旋转，使新解相对于参考星表满足无净旋转条件。 上述推导采用同一个线性化后的约束矩阵 $C$。不过在实际工作中，生成右端项时用参考坐标计算系数，解算时用当前先验坐标计算系数。在一阶近似下，两者可以视为相同。

在实际解算中，先验星表与参考星表可以不同。在传统的 Calc/Solve 流程中，理论时延可以预先算好。pSolve 则可以调用 VTD，在解算过程中按需计算理论时延（on the fly）。我喜欢这种做法，因为换用新的先验源表后，可以在本次解算中重新计算理论时延及其偏导数，只是解算用时也会相应增加。例如，计算理论时延时可以采用较新的源坐标，确定天球参考基准时仍可以使用 ICRF3 的定义源。

在这种情况，后验星表的天球参考架约束问题有两种处理途径。一种是在解算中相对于先验星表施加无净旋转约束，解算后利用共同参考源拟合整体旋转，再校准输出星表的参考架坐标轴指向。另一种是先完成一次解算，再用参考坐标减先验坐标计算所选源的约束右端项，重新运行全局解算，使新解以这份参考星表为天球参考基准。

二、控制文件中的配置

在控制文件的 $CONSTRAINTS 段中，用 NO_NET_ROTATION_SOURCE 设置射电源位置的无净旋转约束。可以用下列语句关闭这项约束。

NO_NET_ROTATION_SOURCE NO

设置约束的作用范围、权重、强度和参考源时，使用以下语法。

NO_NET_ROTATION_SOURCE [GLOBAL] [LOCAL] [UNIFORM / WEIGHTED] [SIGMA <sigma>] [RIGHT_PART <r1> <r2> <r3>] [YES / NO] [EXCEPT <names>]

方括号表示可选项，斜线分隔备选项。

取值或子项	作用与限制
`YES [EXCEPT <names>]`	选用所有符合条件的源，排除例外名单中的源。
`NO EXCEPT <names>`	只选用名单中符合条件的源。
`UNIFORM`	对入选的坐标分量等权，采用第 1.2 节的系数。省略权重选项时采用此模式。
`WEIGHTED`	按当前法矩阵的对角元素对坐标分量加权，再逐行归一化，见第 1.3 节。
`SIGMA <sigma>`	三条约束方程共用的标准差，单位为弧度，默认值为 $10^{-11}\ \mathrm{rad}$（约 $2.06\ \mu\mathrm{as}$）。
`RIGHT_PART <r1> <r2> <r3>`	三条约束方程的右端项，单位为弧度，默认为零。
`GLOBAL`	作用于作为全局参数估计的源坐标。
`LOCAL`	作用于作为局部参数估计的源坐标，可与 `GLOBAL` 并用。两者均省略时，默认使用 `GLOBAL`。

例如，Leonid Petrov 在 RFC 2026c 的示例文件 rfc_2026c_zero_net_sou_cns.cnt 中使用了下面的设置：

  NO_NET_ROTATION_SOURCE         SIGMA    1.00D-10          UNIFORM  \ * nnr_src
    RIGHT_PART  5.667530D-09 -1.214608D-09  2.850700D-08  NO EXCEPT  \ * nnr_src
    ... (后接ICRF1定义源名称列表)

这份配置用 NO EXCEPT 选取 ICRF1 的 212 颗定义源。由于省略了 GLOBAL 和 LOCAL，约束作用于作为全局参数估计的源坐标。名单中的源实际贡献哪些坐标分量，取决于本次解算估计了哪些参数，具体判断见第 3.3 节。UNIFORM 对应第 1.2 节的三条等权方程。RIGHT_PART 给出它们的目标值，其中第三个数是名单内各源的参考赤经减先验赤经之和。三条方程共用的标准差为 $10^{-10}\ \mathrm{rad}$，约合 $20.63\ \mu\mathrm{as}$。

三、程序实现

3.1 源代码与程序

这里主要涉及三个源文件。

源文件	作用
`psolve/progs/solve/batch/gconst.f`	读取 `NO_NET_ROTATION_SOURCE` 的配置
`psolve/progs/solve/norml/nnr_src.f`	选择源并构造三条约束方程
`psolve/progs/solve/norml/apply_cnstr.f`	将约束贡献加入法方程矩阵和右端向量

控制文件读取阶段的调用顺序是 BATCH → CTRLS → CTRLFL → GCONST → GYNEXL。实际计算时，局部参数和全局参数分别通过以下路径处理。

配置范围	实际调用路径
`GLOBAL`	`PRCES → GLO_NORM → NORML_MAIN → NNR_SRC`
`LOCAL`	`PRCES → GLO_FORW → PROC_DO → PSTLP → CNSTR → NNR_SRC`

NNR_SRC 根据输入参数 FGLOBAL_L4 区分当前处理的是全局参数（.TRUE.）还是局部参数（.FALSE.），再结合配置开关 NNR_SRC_GLO、NNR_SRC_LOC 判断是否建立约束。在独立解算模式（INDEPENDENT）下，GLOBAL 设置也可通过局部处理路径生效。两条路径均在建立约束后调用 APPLY_CNSTR，把约束加入法方程。

3.2 配置读入

以第二部分的配置为例，GCONST 将各项设置保存到下列变量中。

配置内容	内部变量
省略 `GLOBAL/LOCAL`	`NNR_SRC_GLO = .TRUE.`
`UNIFORM`	`KSIG_SCALE_NNRS = 'UN'`
`SIGMA 1.0D-10`	`NNR_SRC_SIGMA = 1.0D-10`
`RIGHT_PART`	`NNR_SRC_RTP(1:3)` 按配置中的顺序存放三个右端项
`NO`	`DEFSRC` 第 7 位为 0，即默认不参加
`EXCEPT` 后的源名单	源名存入共享源名数组，在 `SOUSUP(:,7)` 中设置相应的选择位
本条约束的名单范围	`IND_NNR_SRC(1:2)`

其中，GYNEXL 先读取限定词，遇到 YES 或 NO 后再处理例外名单。因此，GLOBAL、UNIFORM、SIGMA 等限定词应写在 YES/NO 前面。

3.3 建立约束方程

NNR_SRC 分两步选源。第一步，根据默认状态和例外名单判断是否选中。

KINCL = KBIT(DEFSRC, INT2(7))
...
IF (SRCSUP(J) .EQ. ISTRN_CHR(ISTAR)) THEN
   IF (KBIT4(SOUSUP(1,7), J)) KINCL = .NOT. KINCL
END IF

程序先取默认选择状态。若源名出现在本条约束的例外名单中，且相应选择位已设置，就翻转该源的状态。这里直接比较 SRCSUP(J) 与 ISTRN_CHR(ISTAR)，因此配置名单要使用与程序内部一致的源名。

第二步，检查该源在当前处理阶段是否有待估坐标分量。

IF (.NOT. KBIT4(LSTAR(1,1),ISTAR) .AND. &
    .NOT. KBIT4(LSTAR(1,2),ISTAR)) KINCL = .FALSE.

赤经和赤纬都未参与估计时，程序直接排除该源。其余源还要按源名和参数名（RIGHT ASCEN 或 DECLINATION）匹配当前参数表，只对匹配到的坐标分量写入系数。

在局部处理路径中，CNSTR 复制完整的参数名称表，再把副本中属于全局参数的名称置空，保留原来的列号。例如，第 10 列是全局参数，第 11 列是局部源的赤经，程序就清空第 10 列的名称。NNR_SRC 匹配到第 11 列时，写入的系数仍对应法方程的第 11 列。

CVEC 保存各条方程的系数，IXREF_PAR 记录系数对应的参数列号。对于赤经，NNR_SRC 写入以下三个系数：

CVEC(NC(1),1) = -COS(RA)*TAN(DEC)
CVEC(NC(2),2) = -SIN(RA)*TAN(DEC)
CVEC(NC(3),3) =  1.0D0

赤纬只在前两条方程中有非零系数：

CVEC(NC(1),1) =  SIN(RA)
CVEC(NC(2),2) = -COS(RA)

在 UNIFORM 模式下，这两组赋值直接给出第 1.2 节三条约束方程的系数。注意这里等权的是原始赤经、赤纬改正量。如果改为对切平面位移等权，坐标改正量与旋转设计矩阵都要改写。令 $T_i=\operatorname{diag}(\cos\delta_i,1)$，纯旋转引起的切平面位移为

\[\begin{aligned} \boldsymbol u_i &= T_i\Delta\boldsymbol x_i =\begin{pmatrix}\Delta\alpha_i\cos\delta_i\\ \Delta\delta_i\end{pmatrix}=B_i\boldsymbol\epsilon,\\ B_i &= T_i A_i =\begin{pmatrix} -\cos\alpha_i\sin\delta_i & -\sin\alpha_i\sin\delta_i & \cos\delta_i\\ \sin\alpha_i & -\cos\alpha_i & 0 \end{pmatrix}. \end{aligned}\]

对 $\boldsymbol u_i$ 等权时，零右端约束为 $\sum_i B_i^\mathsf T\boldsymbol u_i\simeq\boldsymbol0$。写回原始坐标后，第三条方程变为 $\sum_i\cos^2\delta_i\,\Delta\alpha_i\simeq0$，而 pSolve 执行的第三式中，每个赤经改正量的系数都是 $1$。源码第 194–196 行的注释恰好给出了切平面方程写回原始坐标后的三个赤经系数。

! -cos(dec)*sin(dec)*cos(ra)
! -cos(dec)*sin(dec)*sin(ra)
!  cos(dec)**2

当前代码执行的是含 TAN(DEC) 的公式。这一点可能是 pSolve 与其他SOLVE软件在天球参考架的无净旋转约束的不同。

在 WEIGHTED 模式下，NNR_SRC 用 IXREF_PAR 中的参数列号找到相应的法矩阵对角元素。以 NOR_MAT 的存储方式为例，程序先对行内每个参数执行以下运算。

IND = INDX8(IXREF_PAR(I,ICON),IXREF_PAR(I,ICON))
CVEC(I,ICON) = CVEC(I,ICON)*NOR_MAT(IND)

一行系数全部加权后，再调用 NORMALIZE_VECTOR4，得到第 1.3 节给出的归一化系数。

CALL NORMALIZE_VECTOR4(CVEC(1,ICON), NC(ICON))

即使入选参数的对角元素都相同，WEIGHTED 也会把每个非零行缩放为单位范数，UNIFORM 则保留原始系数。两种模式使用同一个 SIGMA 时，约束强度仍可能不同。从第 1.4 节的 $\chi^2$ 表达式可见，若要将一条约束等价地写成单位范数形式，就要把该行系数、右端项和标准差同时除以行范数。程序在这一步只归一化系数。使用 WEIGHTED 且目标非零时，应按归一化后的 $C$ 计算 RIGHT_PART。

系数确定后，NNR_SRC 将三条方程写入约束对象 CNSTROBJ。

CALL ADDCNS_NAM(..., NNR_SRC_RTP(ICON), NNR_SRC_SIGMA, ...)
CALL ADDCNS_EQU(..., IXREF_PAR(IPARM,ICON), CVEC(IPARM,ICON), ...)

ADDCNS_NAM 记录右端项和标准差，ADDCNS_EQU 将非零系数和对应的参数列号写入约束对象。APPLY_CNSTR 随后把约束加入法方程。用矩阵表示时，$C$ 按完整法方程的参数列号排列，其他参数对应的列为零。令 $W_c=\sigma^{-2}I_3$，法矩阵 $N$ 和右端向量 $\boldsymbol b$ 的更新为 $\begin{aligned} N_{\mathrm{new}} &= N+C^\mathsf T W_c C,\\ \boldsymbol b_{\mathrm{new}} &= \boldsymbol b+C^\mathsf T W_c\boldsymbol r. \end{aligned}$ 取默认值 $\boldsymbol r=\boldsymbol0$ 时，约束只更新法矩阵。第二部分的 RFC 示例给出了非零右端项，对应的贡献也要加到法方程的右端向量中。

从配置读入到输出报告，有四处值得检查。

检查 RIGHT_PART 的三个数是否写齐。 少写分量时，程序可能沿用零值。建议在读取配置时指出缺少哪个分量。
检查 SIGMA 后面是否有有效正数。 单写 SIGMA 可能沿用默认值，零值还可能进入 1/SIGMA**2 的计算。建议在读取控制文件时就报错。
检查三条约束能否确定三个旋转自由度。 仅有源参加还不足以保证三条方程独立。建议依据实际入选的源和坐标分量在源坐标参数空间中检查 $CA$ 是否满秩，并在秩不足时报告未受约束的旋转方向。
按实际参加约束的源生成报告名单。 目前只在处理赤经时登记源名，因此只估计赤纬的源即使参与了前两条约束，也可能在报告中漏列。建议任一坐标分量进入约束时就登记源名，最后去重输出。

开启约束名单输出后，可以在解算输出报告中查找以下两行。

NNR_SRC: ... sources participated in NNR_SRC constraints:
NNR_SRC sigma: ... right hand sides: ... rad

第一行给出程序登记的源数，随后的各行列出源名。NNR_SRC sigma: 一行给出实际使用的标准差和三个右端项。核对约束时，还应结合控制文件中的 GLOBAL/LOCAL 和 UNIFORM/WEIGHTED 设置。