简介面向无源雷达与被动定位研究场景这份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典型场景工程判断10 ns3 m宽带扩频信号、高精度同步站可做严格几何求交0.1 µs30 m通用通信信号、GPS 驯钟同步求交后需最小二乘复核1 µs300 m调频广播带宽受限场景直接走统计定位流程10 µs3000 m窄带信号、模糊函数主峰很宽仅能维持航迹级定位当 σ_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_scan90): 求两个椭圆的全部交点 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)距离和 L126第二组 T2(10,5)R2(30,5)距离和 L221。运行 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_LNone): 多站无源定位椭圆法融合 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, methodlm) return res.x, res.cost, res.jac调用时 x0 用第 3 节算出的所有候选交点中残差最小的那个或者直接取所有接收站位置的平均值sigma_L 的单位和距离一致。methodlm适合 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_free2): 由 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)), covm 是站数n_free 是待估参数个数 2。sigma2 的估计假设加权后的残差是零均值白噪声如果系统里有未消除的多径误差这个估计会偏大物理上反而是好事如实反映定位结果不可信。GDOP 的单位与坐标单位一致表示定位误差的均方根半径。布站经验可以总结成下表布站特征GDOP 表现工程建议目标落在两焦点连线或其延长线上椭圆退化GDOP 极大布站时避开目标主航路与基线共线基线长度接近目标距离椭圆交叉角大误差椭圆较圆优先保证基线长度足够各站时延精度差异大高精度站被低精度站拖累必须加权不能等权目标高度未建模残差带系统偏差协方差偏小用 DEM 或高度先验修正 L一个常见误操作是所有站设成等权。正确的做法是把时延测量方差、站址误差、甚至直达波多径不确定性都折算进 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, alpha0.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 的卡方门限筛点误报率会显著下降这也是无源雷达数据链路上比“算得准”更要紧的——报得准。本文还有配套的精品资源点击获取