4. 位置估计的几何基础
物理源项反演要回答两个问题:源有多强、源在哪里。位置估计是后者的核心,而它本质上是一个几何问题——把若干传感器上的测量值(浓度、剂量率、到达时间、方位角等)翻译成空间中的约束,再求这些约束的交点或最优拟合点。本节围绕这条主线,讲清坐标系怎么选、测量模型怎么写、定位方程怎么解,以及误差如何被几何构型放大。
一、坐标系与空间变换
反演中常用的坐标系有三类。第一类是局部直角坐标系(ENU:东-北-天),以某个参考点为原点,适合描述厂区、实验厅内的源与探测器位置,单位通常为米。第二类是大地坐标系(经纬度高程,或ECEF地心地固坐标),适合大范围、跨区域的监测网络。第三类是传感器自身坐标系,例如探测器阵列的本地框架,用于描述阵列内部各单元的相对位置和朝向。
三者之间通过刚体变换联系。若已知传感器本地坐标 $\mathbf{p}_l$,要转到ENU坐标 $\mathbf{p}_e$,需要旋转矩阵 $\mathbf{R}$ 和平移向量 $\mathbf{t}$:
p_e = R · p_l + t
R = Rz(ψ) · Ry(θ) · Rx(φ) // 依次绕Z、Y、X轴旋转,ψ为航向角
实际工程中最容易出错的是角度约定:航向角是从北向东顺时针还是从东向北逆时针,俯仰角的正负指向,旋转顺序是ZYX还是XYZ。这些约定一旦在代码和图纸之间不一致,定位结果会出现系统性偏移。建议在项目开始时固定一套约定并写入文档,例如“ENU坐标,方位角自北向东为正,俯仰角向上为正,旋转顺序Z-Y-X”,全文和图示保持一致。切换参照系时必须显式说明,不能默认读者知道。
二、距离与角度测量模型
不同探测手段给出的约束形式不同,常见的有四类:
- 距离约束(测距):如到达时间(TOA)乘以波速、信号强度衰减反演距离。模型为 $d_i = \|\mathbf{p} - \mathbf{s}_i\| + \varepsilon_i$,其中 $\mathbf{s}_i$ 为第 $i$ 个传感器位置,$\varepsilon_i$ 为测距误差。
- 距离差约束(TDOA):两个传感器收到同一信号的时间差对应距离差,$d_{ij} = \|\mathbf{p}-\mathbf{s}_i\| - \|\mathbf{p}-\mathbf{s}_j\|$,消去了源发射时刻这一未知量,适合时钟不同步的场景。
- 角度约束(测向):阵列或准直器给出源相对传感器的方位角/俯仰角,约束为一条射线,$\tan\alpha_i = (p_y - s_{iy})/(p_x - s_{ix})$。
- 场强/浓度约束:点源在距离 $r$ 处的剂量率近似按 $1/r^2$ 衰减,$C_i = A/(4\pi r_i^2)\cdot k$,可反解出距离,但受本底、散射和源强未知的影响,误差通常较大。
测量模型的选择直接决定后续方程的形式。测距给出圆(二维)或球(三维)约束,测向给出直线约束,TDOA给出双曲线约束。实际反演中常把多类约束联合使用,例如“两站测向 + 一站测距”,以提高可观测性。
三、三角定位与多边定位
三角定位(triangulation)利用角度:两条测向射线在平面上交于一点,即为源位置。三维中两条射线一般不相交(异面),需要求最近点对或做最小二乘。两站测向在二维中恰好可解,但若两站与源接近共线,射线夹角很小,交点对角度误差极其敏感。
多边定位(multilateration)利用距离:二维中三个距离圆一般交于一点,三维中四个球面确定一点。实际测量有误差,圆不会严格共点,此时转为求解非线性最小二乘问题:
min_p Σ w_i · ( ||p - s_i|| - d_i )^2
// w_i 为按测量不确定度设定的权重,w_i = 1/σ_i^2
// 初值可用线性化方法(如将方程两两相减消去二次项)给出
线性化技巧值得掌握:将 $d_i^2 = \|\mathbf{p}-\mathbf{s}_i\|^2$ 展开,用第 $i$ 式减第 $1$ 式,可消去 $\|\mathbf{p}\|^2$ 项,得到关于 $\mathbf{p}$ 的线性方程组,先解出粗略初值,再用高斯-牛顿或列文伯格-马夸尔特迭代精化。这一步在工程实现中几乎是标配,能显著改善收敛性。
四、误差传播与几何精度因子
测量误差不可避免,关键是它如何传到位置估计上。设测量向量为 $\mathbf{z}$,位置解为 $\hat{\mathbf{p}}$,一阶近似下有 $\delta\mathbf{p} \approx (\mathbf{J}^T\mathbf{W}\mathbf{J})^{-1}\mathbf{J}^T\mathbf{W}\,\delta\mathbf{z}$,其中 $\mathbf{J}$ 为测量方程对位置的雅可比矩阵。位置协方差为:
Cov(p) = (J^T W J)^(-1)
GDOP = sqrt( trace(Cov(p)) ) // 几何精度因子
GDOP 只取决于几何构型,与测量误差大小无关:传感器包围源、分布均匀、张角大时 GDOP 小;传感器集中在源的一侧、或与源接近共线时 GDOP 急剧增大。工程上常把 GDOP 作为布站优化的目标函数,在候选点位中搜索使 GDOP 最小的组合。一个实用的经验是:二维定位至少需要三个站形成包围,且源到各站的方位角尽量分散;若只能单侧布站,应补充测向或先验约束。
五、实践中的几个要点
第一,先检查可观测性。在解方程之前,用雅可比矩阵的秩或条件数判断当前测量组合能否唯一确定三维位置;只有两个测距站在三维中是不够的,必须补充角度或先验信息。第二,权重设置要反映真实不确定度,测向和测距的量纲不同,不能简单等权相加,应按各自方差归一化。第三,结果要带不确定度输出,只给一个点坐标而不给置信椭圆或误差范围,在反演报告中是不完整的。第四,坐标系约定要贯穿全流程,从传感器标定、数据入库到结果展示,任何一次隐式转换都可能引入米级甚至更大的偏差。
把几何基础打牢之后,后续的源强反演、多源分辨和不确定度量化才有可靠的落脚点。位置估计的精度上限由几何构型决定,而实际达到的精度则由测量模型与误差处理的严谨程度决定。