无源定位椭圆法：被动雷达多站融合定位算法解析

简介：面向无源雷达与被动定位研究场景，这份MATLAB源码实现椭圆法目标定位中的关键步骤——多站观测椭圆交点求解。它根据信号到达时间差/频率差信息构建椭圆模型，通过数值迭代计算目标平面位置，可避免手工解算非线性方程的繁琐并降低误差积累。代码接口清晰，使用者只需提供各观测站对应的椭圆参数或时差数据，即可获得候选交点坐标，便于直接嵌入定位流程；压缩包内含1个.m文件，体积约1KB，轻量紧凑，适合算法验证、教学演示或二次开发。目前累计已有345人学习浏览，主要面向雷达信号处理、无源定位领域的初学者与工程技术人员。借助该程序，可快速理解椭圆法几何原理，复现目标坐标估计过程，也可作为改进TDOA/FDOA定位精度或构建多站协同定位系统的参考起点。

1. 无源定位椭圆法：被动雷达不发声，也能用时间差画椭圆

在不发射任何能量的前提下定位一个未知目标，听起来违背直觉。被动雷达接收机旁边恰好有广播塔或通信基站作照射源，接收机同时收到直达波和目标反射的回波；两路信号做互相关测出时延差，乘光速后得到“距离和”——目标到照射源的距离加目标到接收机的距离。这个量把目标限制在以照射源和接收机为焦点的椭圆上，这就是无源定位椭圆法。单个椭圆定不了位，换一组“照射源-接收机”组合再测一次，两个椭圆的交点就是目标坐标，findEllIntersect 这类数值函数就是为算这个交点而生的。下面按双基测量建模、椭圆参数生成、交点求解、多站加权最小二乘、残差验证的顺序把链路走通，适合正在写无源定位算法或做被动雷达定位工程的人。

2. 椭圆法的几何基础：双基时延怎么变成椭圆参数

2.1 双基距离和测量：直达波、回波与时延差的换算

无源雷达系统里必须有一个机会照射源，典型如调频广播塔、模拟电视塔、4G/5G 基站。接收站有两个接收通道，一个指向照射源方向锁住直达波，一个指向监视空域锁住回波。设照射源位置为 T，接收机位置为 R，未知目标位置为 P。直达波到达时刻是||T - R|| / c，回波到达时刻是(||T - P|| + ||P - R||) / c，两者相减得到时延差 Δτ，于是双基距离和为：

L = ||T - P|| + ||P - R|| = ||T - R|| + c * Δτ

其中||T - R||是两个已知站点之间的距离，工程上可以精确计算，测量引入的未知量只有 Δτ。这里和 TDOA 双曲线定位有个容易混淆的点：TDOA 是同一目标辐射信号被两个接收站接收，测的是“距离差”，对应双曲线；椭圆法测的是“照射源-目标-接收机”的绕行路径和直达路径之间的时延差，得到的是“距离和”，对应椭圆。被动雷达定位天然具备直达波参考，所以椭圆法在无源雷达里比双曲线法更自然。实际系统中 Δτ 通过互模糊函数峰值搜索获得，同时还能给出多普勒频移，窄带调频信号的距离分辨率大约几十米量级，信号有效带宽越宽，测时延越准。

2.2 由焦点和距离和生成椭圆参数：半长轴、短轴与倾角

椭圆的标准定义是平面上到两个定点（焦点）距离之和等于常数 2a 的点的轨迹。在本场景里，两个焦点就是 T 和 R，“距离和”就是上一节算出的 L，因此半长轴a = L / 2。两个焦点之间的半焦距为c0 = ||T - R|| / 2，再通过直角三角形关系得到半短轴b = sqrt(a² - c0²)。焦点连线与全局坐标系的夹角就是椭圆长轴方向：

φ = atan2(R.y - T.y, R.x - T.x)

有了 center、a、b、φ 四个量，椭圆就完全确定了。下面的 Python 函数直接把这些量算出来，供后面的交点求解使用：

import numpy as np def ellipse_params(f1, f2, r_sum): """从焦点 f1, f2 和目标距离和 r_sum 生成椭圆参数 f1: 照射源坐标 (x, y)，单位 m f2: 接收机坐标 (x, y)，单位 m r_sum: 双基距离和 L，单位 m，由时延差乘光速推算 返回: center, a, b, angle """ center = (f1 + f2) / 2.0 c_half = np.linalg.norm(f2 - f1) / 2.0 a = r_sum / 2.0 b = np.sqrt(max(a * a - c_half * c_half, 1e-12)) angle = np.arctan2(f2[1] - f1[1], f2[0] - f1[0]) return center, a, b, angle

代码逻辑很直接：中心取两焦点中点，半焦距取焦点距离的一半，a 由距离和减半给出，b 用勾股关系从 a 和 c_half 求。max判断是为了防止测量误差导致a² - c_half²出现微小的负值，这是工程里常见的数值防御，随手写上比报错强。如果 r_sum 小于两焦点距离，说明时延差测量异常，这个椭圆不成立，应该在数据预处理阶段丢弃，而不是留到求交阶段。

2.3 椭圆带宽度与测量精度：解析求交何时失效

每个椭圆约束在几何上不是一条无宽度的曲线，时延测量误差σ_τ会直接转化为距离和误差σ_L = c * σ_τ，让椭圆变成一条“带”。两个椭圆在噪声下相交，实际是两个带相交成一个模糊区域。常用量级的换算关系如下：

当 σ_L 达到几百米而两焦点基线也只有几公里时，解析方法解出的“精确交点”在噪声意义下没有价值。判断标准可以这样定：当 σ_L 超过椭圆短轴的 5% 到 10% 时，就不要再用严格求交的思路去报点，直接改用加权最小二乘，解析求交只负责产生初值和数据关联的候选点。

3. 两个椭圆求交：findEllIntersect 的数值实现与初值选取

3.1 把求交转成一维找根问题

两个任意位置、任意倾角的椭圆求交，解析上可以展开成两个二元二次方程联立，消元后得到一个一元四次方程，理论上最多 4 个实交点。但工程代码里很少有人真的去解四次方程，原因有两个：椭圆方程在全局坐标下的系数包含大量三角函数，推导和写代码都容易出错；噪声场景下两个椭圆画出来有交点，数值上却不严格满足等式，四次方程求根的数值稳定性也差。常见做法是把其中一个椭圆写成参数方程，代入另一个椭圆的距离和约束，把求交变成单变量方程求根。

具体来说，椭圆 1 的参数方程是p(t) = center1 + R(φ1) * [a1 * cos t, b1 * sin t]，t 是偏近点角，取值范围[0, 2π)。把这个 p(t) 代入椭圆 2 的约束，定义函数：

g(t) = ||p(t) - T2|| + ||p(t) - R2|| - L2

找椭圆交点就是找 g(t) 的零点。g(t) 在[0, 2π)上连续，零点最多 4 个，所以先粗网格扫描找出符号变化的区间，再用 Brent 方法精确求根。这个方案实现简单、对初值不敏感，还能一次性找回所有交点而不是只找一个。

3.2 findEllIntersect 核心实现：参数扫描加精确求根

下面的代码沿用无源雷达处理软件里常见的函数名 findEllIntersect，方便和手里的老工程代码对照。输入是两个椭圆的焦点对和距离和，输出是所有交点坐标：

from scipy.optimize import brentq def ellipse_point(t, center, a, b, angle): """椭圆参数方程：t 为偏近点角，返回全局坐标""" ca, sa = np.cos(angle), np.sin(angle) x_loc, y_loc = a * np.cos(t), b * np.sin(t) return np.array([ center[0] + ca * x_loc - sa * y_loc, center[1] + sa * x_loc + ca * y_loc ]) def findEllIntersect(ell1, ell2, n_scan=90): """求两个椭圆的全部交点 ell1, ell2: (f1, f2, r_sum) 三元组，f 为焦点坐标 n_scan: 粗扫描点数，越大越不容易漏根，默认 90 """ center1, a1, b1, ang1 = ellipse_params(*ell1) ts = np.linspace(0.0, 2.0 * np.pi, n_scan + 1) roots = [] def g(t): p = ellipse_point(t, center1, a1, b1, ang1) d1 = np.linalg.norm(p - ell2[0]) d2 = np.linalg.norm(p - ell2[1]) return d1 + d2 - ell2[2] for i in range(n_scan): g0, g1 = g(ts[i]), g(ts[i + 1]) if g0 * g1 < 0.0: t_root = brentq(g, ts[i], ts[i + 1]) p = ellipse_point(t_root, center1, a1, b1, ang1) if all(np.linalg.norm(p - q) > 1e-6 for q in roots): roots.append(p) return np.array(roots) if roots else np.empty((0, 2))

几个容易被忽略的参数值得单独说明。n_scan决定粗扫描密度，两个椭圆在长轴方向拉得很长时，交点区域可能只落在很窄的 t 区间内，n_scan 过小会漏根，一般取 90 到 180 足够。brentq要求区间端点处函数值异号，而 g(t) 在端点处恰好为零时（相切）不会被捕获，所以近似相切场景下建议加密扫描后重试。去重阈值1e-6是按坐标单位米设置的，如果坐标系换成经纬度，要改成1e-11量级，否则同一个交点会被重复报出来。

提示：若发现漏根，先不要怀疑求根算法，先把 n_scan 从 90 提高到 180。多花的代价只是几百次距离计算，但常能避开那些肉眼可见、代码却扫不到的窄交点区间。

3.3 运行示例与交点选择逻辑

给一组可复现的输入验证算法行为。照射源 T1=(0,0)，接收机 R1=(20,0)，距离和 L1=26；第二组 T2=(10,5)，R2=(30,5)，距离和 L2=21。运行 findEllIntersect 后，返回的可能会有 2 个或 4 个交点。

ell1 = (np.array([0.0, 0.0]), np.array([20.0, 0.0]), 26.0) ell2 = (np.array([10.0, 5.0]), np.array([30.0, 5.0]), 21.0) pts = findEllIntersect(ell1, ell2) print(pts)

两个椭圆最多有 4 个交点，物理上有意义的点只有一个，多出来的点来自几何多解。工程上的筛选依据有三个：目标高度、多普勒频移和轨迹连续性。目标高度把二维交点投影回三维需要先验高度或者用第三组椭圆去卡；多普勒频移与目标相对观测几何的径向速度有关；轨迹连续性用于跟踪滤波，比如卡尔曼滤波的预测门限。多站无源定位里一般不强制在几何求交阶段选出唯一点，而是把所有候选点全送进下一级的数据关联模块，让多帧观测来裁决哪个是真实目标。

4. 多站融合定位：被动雷达定位的加权最小二乘与 GDOP 控制

4.1 从严格交点到残差最小化：目标函数与雅可比

实际多站无源定位里，N 组观测对应 N 个椭圆，这些椭圆通常不会交于同一点。原因除了时延测量误差，还有站点坐标误差、多径导致的额外传播路径、以及目标高度被忽略带来的系统偏差。继续依赖几何求交会陷入“选哪个交点”的循环，更稳妥的做法是放弃几何交点概念，定义加权残差：

r_i(p) = (||p - T_i|| + ||p - R_i|| - L_i) / σ_i

目标函数写成加权平方和J(p) = Σ r_i(p)²。σ_i 是第 i 站的双基距离和标准差，由时延估计方差和站址误差分量合成。这个目标函数对 p 的雅可比很容易推导：

∂r_i/∂p = (p - T_i) / ||p - T_i|| + (p - R_i) / ||p - R_i||

每一项都是一个指向焦点方向的单位向量之和，几何含义很直观——残差梯度由目标指向照射源和接收机的单位矢量合成。这个解析形式可以直接交给优化器用，也可以忽略改用数值差分，站点数少时数值差分足够稳定，代码上也省事。

4.2 多站无源定位的 least_squares 实现

用 scipy.optimize.least_squares 实现多椭圆融合定位的代码很短：

from scipy.optimize import least_squares def locate_ellipse(obs, x0, sigma_L=None): """多站无源定位椭圆法融合 obs: 数组，每行 (Tx_x, Tx_y, Rx_x, Rx_y, L) x0: 初值，可用 findEllIntersect 的交点或站点几何中心 sigma_L: 各站距离和标准差，None 时等权 """ if sigma_L is None: sigma_L = np.ones(len(obs)) def resid(p): out = [] for (Tx, Rx, L), s in zip(obs, sigma_L): d = np.linalg.norm(p - Tx) + np.linalg.norm(p - Rx) out.append((d - L) / s) return np.array(out) res = least_squares(resid, x0, method='lm') return res.x, res.cost, res.jac

调用时 x0 用第 3 节算出的所有候选交点中残差最小的那个，或者直接取所有接收站位置的平均值；sigma_L 的单位和距离一致。method='lm'适合 m 个残差、2 个未知数的中小规模问题，不需要显式提供雅可比。least_squares 内部按残差向量做优化，这里的加权已经让不同精度的站在同一尺度下参与拟合，后续协方差计算可以直接使用返回的雅可比。locate_ellipse返回平均意义下的定位点、残差平方和的一半以及雅可比三项，后面两项在精度评估时要用。

sigma_L 的整定没有统一标准，常见做法是按三项合成：时延估计标准差给出的c * σ_τ、站址坐标误差在两个焦点方向的投影（通常取两个站点各自位置不确定度的半数）、以及直达波通道多径造成的固定偏移。前两项是随机量，第三项在城市环境里往往是系统性偏置，最好用已知位置的静态目标做一次标定，把偏置残差均值减掉后再进入后续实时处理。

提示：least_squares 返回的 cost 是0.5 * Σr²，计算误差方差时系数 2 不要漏掉，否则 GDOP 会系统性偏小。

4.3 协方差估计与布站建议

最小二乘收敛后，用残差雅可比估算定位协方差矩阵：

def gdop_from_jac(res_jac, res_cost, m, n_free=2): """由 least_squares 结果估算 GDOP res_jac: 加权残差雅可比 res_cost: least_squares 的 cost，数值上等于 0.5 * sum(r^2) """ sigma2 = 2.0 * res_cost / (m - n_free) cov = sigma2 * np.linalg.inv(res_jac.T @ res_jac) return np.sqrt(np.trace(cov)), cov

m 是站数，n_free 是待估参数个数 2。sigma2 的估计假设加权后的残差是零均值白噪声，如果系统里有未消除的多径误差，这个估计会偏大，物理上反而是好事，如实反映定位结果不可信。GDOP 的单位与坐标单位一致，表示定位误差的均方根半径。布站经验可以总结成下表：

一个常见误操作是所有站设成等权。正确的做法是把时延测量方差、站址误差、甚至直达波多径不确定性都折算进 sigma_L 再代入 least_squares。加权和不加权，在代码上只差一个参数，在最终定位误差上经常是几百米和一公里的差别。

5. 实战技巧：后验残差剔除坏观测，把无源目标钉在地图上

5.1 用卡方检验判断定位结果是否可信

多站定位算出一个坐标后，只报坐标不报可信度没有意义。把估计点回代每个站的残差，构造统计量χ² = Σ(r_i/σ_i)²。在二维定位、N 站观测的场景下，自由度为 N - 2；给定置信度 0.99，查卡方分布临界值，低于临界值说明残差幅度和噪声假设一致，定位结果可信，否则就要怀疑有坏观测：

from scipy.stats import chi2 def verify_loc_chi2(p_est, obs, sigma_L, alpha=0.01): """卡方检验：判断定位结果是否与噪声假设一致 obs: 每行 (Tx_x, Tx_y, Rx_x, Rx_y, L) """ r = [] for (Tx, Rx, L), s in zip(obs, sigma_L): d = np.linalg.norm(p_est - Tx) + np.linalg.norm(p_est - Rx) r.append((d - L) / s) r = np.array(r) chi2_val = np.sum(r ** 2) df = len(obs) - 2 return chi2_val < chi2.ppf(1 - alpha, df), chi2_val

自由度减 2 是因为二维坐标耗掉了两个自由度。若只有两个站，自由度为 0，卡方检验退化为要求两个残差严格同号且幅度一致，此时只能靠第 3 节的交点筛选逻辑做补充判断。

5.2 坏值剔除与重定位：一次只剔一站

多径是椭圆法最大的实际威胁。城市环境里目标回波被建筑物反射后，等效路径变长，L 偏大，直接污染距离和约束。剔除策略不应该是“残差最大就删”，因为最小二乘会把坏值影响分摊到多个残差上；正确做法是一次剔除残差最大的那一站，剩余站重新定位，再重新做卡方检验，直到检验通过或站数少于 3 为止。每轮重定位后残差会重新分配，上一轮看似正常的站可能在剔除后暴露问题，所以必须迭代，不能一轮定案。相关经验是：剔除后目标位置移动超过 3 倍 GDOP，说明被剔除的站确实在拉偏结果。

5.3 坐标投影与工程落地细节

椭圆定位计算全部在局部平面坐标系内进行。站点经纬度要先做投影转换，常见的是 UTM 或高斯-克吕格投影，定位结果再反投影回经纬度。不要在经纬度上直接算欧氏距离和椭圆参数，纬度 60 度处经度方向 1 度的实际距离不到纬度方向的一半，直接算会把椭圆焦点距离和距离和全部算错。另一个容易忽略的细节是目标高度：二维椭圆假设目标与站点在同一平面上，目标飞过接收机上方时，距离和里多出的高度分量会被误判为水平距离，导致定位向站点方向收缩。处理办法是给每个 L 加高度修正，用目标先验高度 h 近似减去h² / (2 * R_target)量级的修正项，或者把目标高度也放进未知数里扩展为三维定位。把时延环宽度、站址误差和高度不确定性全部折算进 sigma_L，再按 5.1 的卡方门限筛点，误报率会显著下降，这也是无源雷达数据链路上比“算得准”更要紧的——报得准。

本文还有配套的精品资源，点击获取