简介：面向无源被动雷达定位与椭圆法算法研究的一份MATLAB源码资源，解决多站观测下无源目标位置求解问题。工程中可利用信号到达时间差或频率差构造椭圆方程，通过解算多个椭圆的交点来确定目标坐标，适用于被动雷达、电子侦察等隐蔽探测场景。压缩包内仅含1个m文件，包体大小约1KB，代码精简，便于直接阅读和改造。已有345人学习下载。资源核心价值在于提供一个可运行的椭圆相交计算程序，帮助使用者快速验证TDOA/FDOA双曲面投影为椭圆后的几何求解流程；同时可作为定位精度评估、多站协作等后续研究的基础脚本，尤其适合算法初学者用于理解非线性方程迭代求解与椭圆参数估计的工程实现。
1. 无源定位里的椭圆法：为什么被动雷达盯上一条时延线就能算目标位置
做被动雷达和无源定位的工程师，大多数时候手里只有被动接收机，不能像主动雷达那样发射信号，所以所有的定位信息都来自目标辐射的信号或者外部辐射源被目标散射的信号。常见做法是用时差（TDOA）做双曲线定位、用测向角（AOA）做交叉定位，但外辐射源雷达里还有一种非常实用的路子——椭圆法。它利用的是双基地几何：发射站、接收站分置，目标把外部辐射源（广播、电视、手机基站等）的信号散射到接收站，通过直达波和散射波的时延差，可以算出目标到发射站与接收站的距离之和。距离之和固定的点，落在一个以发射站和接收站为焦点的椭圆上。于是被动定位就变成了“解椭圆或求椭圆交点”的问题，而 findEllIntersect 这个名字对应的就是这个几何问题。
对做无源目标定位的人来说，椭圆法的价值在于它只需要单部接收机加上已知位置的第三方辐射源就能形成一条目标位置约束；如果再多一个辐射源或一部接收机，就有了两条约束，目标位置就是椭圆交点。这套思路不依赖目标主动发射信号，也不依赖阵列测向精度，尤其适合城市环境里的低空小目标、无人机等“静默”目标。下面从几何推导、函数实现、参数设定到精度验证，把这个方向完整拆开。
2. findEllIntersect 的几何基础：双基地时延如何变成平面上的一副椭圆
2.1 从直达波和散射波的时延差推导椭圆方程
外辐射源无源雷达里，接收机同时收到两个信号：一个是发射站直达接收站的参考信号，另一个是发射站经目标散射到接收站的回波信号。设发射站位置为 T(x_t, y_t)，接收站位置为 R(x_r, y_r)，目标位置为 P(x, y)。两个信号的传播距离之差由时延差 Δt 乘光速 c 得到：
ΔR = c * Δt
这个 ΔR 与目标位置的关系是：
ΔR = (|TP| + |PR|) - |TR|
稍微挪一下：
|TP| + |PR| = |TR| + ΔR
右边是一个常数。所有满足“到两个固定点距离之和等于常数”的点，在平面上的轨迹就是椭圆，焦点分别是发射站和接收站，长轴长度等于 |TR| + ΔR。所以测到一个时延差，就得到一条椭圆约束线，目标一定落在椭圆上。
当只有这一个约束没有其他量测时，目标无法唯一确定，椭圆上任意一点都可能。但这并不代表没用——配合测向角（目标相对接收站的方位角）或者另一个独立时延差，就能把目标从整个椭圆上缩到一个点。
椭圆参数的标准形式是这样的：中心在焦点连线的中点，长半轴 a = (|TR| + ΔR)/2，焦距 c = |TR|/2，短半轴 b = sqrt(a² - c²)。推出这一步之后，就可以把“目标定位”转换为“椭圆与椭圆/直线的几何求交问题”，这正是 findEllIntersect 要做的事。
2.2 用发射机、接收机和目标位置建立椭圆的参数化形式
实际工程中，发射站和接收站不一定在同一坐标系里，所以写程序前先把所有位置统一到同一个平面坐标系。常见做法是选一个投影坐标系（比如高斯-克吕格投影后的局部平面坐标），把经纬度换成 x、y，再以米为单位跑椭圆求交。如果直接在经纬度下做，椭圆方程会因为地球曲率变成椭球面上的复杂曲线，工程上基本不这么干。
建立参数化椭圆时，推荐用旋转加平移的方式而不是直接用标准 x²/a² + y²/b² = 1，因为焦点连线并不总是沿坐标轴方向。设焦点连线的方向角为 θ，椭圆中心为 O(x_c, y_c)，目标坐标可以写成：
x = x_c + a*cos(t)cos(θ) - bsin(t)sin(θ) y = y_c + acos(t)sin(θ) + bsin(t)*cos(θ)
这里的参数 t 是偏心率角，不是目标的极角。代码里用这个参数形式做椭圆求交最稳定，因为它能自然处理旋转，也不会像隐式方程那样出现除以短半轴的数值风险。
参数说明如下：
- a（长半轴）：由时延差决定，(距离和)/2。
- b（短半轴）：sqrt(a² - c²)，其中 c = |TR|/2 是半焦距。
- θ：焦点连线方向角，atan2(y_t - y_r, x_t - x_r) 或用中心指向发射站的方向。
- O：椭圆中心，(x_t + x_r)/2, (y_t + y_r)/2。
工程实现时还有个细节：时延差测量值必须大于 0，并且 |TR| + ΔR 必须大于长轴下界，否则 a < c 就会得到虚短轴，说明时延测量异常或发射站与接收站位置不匹配，这个情况需要在调用求解前拦截。
3. findEllIntersect 的数学解法：解析求交与退化 case 处理
3.1 两个椭圆求交的本质：把几何问题变成多项式求根
findEllIntersect 绝大多数场景下要处理的是“两组椭圆约束求交点”。两个椭圆求交展开后是一个四次方程，理论上能解析求根，但实现复杂且容易遇到数值不稳定。工程上常见做法是用参数方程加数值迭代，或者把其中一个椭圆离散成点序列再求另一个椭圆的最近点，这种思路代码简单，但精度受离散密度限制，不适合高精度定位。
更稳妥的方案是化为求最小距离问题：在两个椭圆上分别采样参数 t1 和 t2，目标函数为 |P1(t1) - P2(t2)|²，这个函数的极小值点如果距离为零，就找到了交点；距离接近零就是近似交点。多维无约束优化可以用 scipy 的 least_squares 来做，因为这是二维参数空间上的最小二乘问题。
但直接优化也有问题：两个椭圆相切或非常接近时，目标函数非常平坦，优化器容易停在局部极小值。所以正规做法是解析求一个四次方程，或者用 Groebner 基消元后求解一元四次多项式。对工程应用，另一个做法是退而求其次——用等间隔参数采样初值，做多次局部优化，从最优初值出发收敛到全局交点。
3.2 findEllIntersect 核心函数设计：输入与输出
用 Python 实现时，我会把 findEllIntersect 设计成接收两个椭圆的参数，返回所有交点坐标。一个可用的函数签名如下：
import numpy as np from scipy.optimize import least_squares def ell_params_from_delay(tx_pos, rx_pos, delta_R): """ 根据时延差构建椭圆参数 tx_pos: 发射站坐标 [x, y] rx_pos: 接收站坐标 [x, y] delta_R: 双基地距离和减去基线长度，即 c*delta_t - |T-R| """ tx_pos = np.asarray(tx_pos, dtype=float) rx_pos = np.asarray(rx_pos, dtype=float) baseline = np.linalg.norm(tx_pos - rx_pos) sum_dist = baseline + delta_R a = 0.5 * sum_dist c_half = 0.5 * baseline b_sq = a*a - c_half*c_half if b_sq <= 0: raise ValueError("delta_R 太小或基线计算错误，无法形成椭圆") b = np.sqrt(b_sq) theta = np.arctan2(tx_pos[1] - rx_pos[1], tx_pos[0] - rx_pos[0]) center = 0.5 * (tx_pos + rx_pos) return dict(a=a, b=b, theta=theta, center=center) def _ell_point(params, t): """参数 t 对应的椭圆点坐标""" a = params['a'] b = params['b'] theta = params['theta'] cx, cy = params['center'] cos_t = np.cos(t) sin_t = np.sin(t) x = cx + a*cos_t*np.cos(theta) - b*sin_t*np.sin(theta) y = cy + a*cos_t*np.sin(theta) + b*sin_t*np.cos(theta) return np.array([x, y]) def find_ell_intersect(ell1, ell2, n_initial=36, tol=1e-3): """ 求两个椭圆的交点坐标 ell1, ell2: ell_params_from_delay 返回的椭圆参数字典 n_initial: 初始化采样数，越大越不易漏交点 tol: 残差阈值，小于该值认为两点重合 """ best_points = [] t1_scan = np.linspace(0, 2*np.pi, n_initial, endpoint=False) t2_scan = np.linspace(0, 2*np.pi, n_initial, endpoint=False) # 扫描网格，用每个网格点做初值 for t1 in t1_scan: for t2 in t2_scan: def residual(vars): p1 = _ell_point(ell1, vars[0]) p2 = _ell_point(ell2, vars[1]) return p1 - p2 res = least_squares(residual, [t1, t2], method='lm', xtol=1e-10, ftol=1e-10) if np.linalg.norm(res.fun) < tol: p = 0.5 * (_ell_point(ell1, res.x[0]) + _ell_point(ell2, res.x[1])) # 去重 if not any(np.linalg.norm(p - q) < 1.0 for q in best_points): best_points.append(p) return np.array(best_points)这段代码的逻辑并不复杂：先在两个椭圆的参数域 t1、t2 上做 36×36 的网格采样，把每个采样点组合作为初始猜测，用 Levenberg-Marquardt 最小化两个椭圆上的点差。收敛后如果残差足够小，就认为得到一个交点候选。最后按 1 米距离去重，避免同一交点被多个初值反复找到。
参数设置上有三个关键点：
- n_initial 不能太小，两个椭圆最多有 4 个交点，如果采样数过少，初值可能陷入错误局部极值。
- tol 的值一般取 0.1 到 1 米，具体和系统定位精度相关；如果定位误差要求到分米级，tol 应该设到 0.05 以下。
- least_squares 的 method 用 'lm'，也就是 Levenberg-Marquardt，适合只有两个参数的小规模问题，收敛快；不要用 'trf' 在这个场景，它处理边界约束更擅长，但在这里反而更慢。
3.3 解析求根与数值兜底的混合策略
前面这种网格扫描加优化的方法有个弱点：当两个椭圆非常接近但不相交时，优化器可能收敛到最小距离点，给出的“交点”其实是近焦点，造成假定位。因此在实战中我不会只用数值优化，而是叠加一个解析步骤：把两个椭圆的隐式方程写成二次锥形式，消元后得到一个单变量四次多项式，然后用 numpy.roots 求根。
解析消元的实现依赖于椭圆的标准表达：将每个椭圆写作二次型矩阵 Q_i，方程是 [x, y, 1] * Q_i * [x, y, 1]^T = 0。两个方程做线性组合，可以消去 x² 和 y² 中的一项，再用代入法降阶。这个推导过程冗长，但代码实现很固定。实际工程里我经常先把解析解跑出来的交点作为初值，再跑一次局部最小化修正，这样既不担心漏根，又能保证收敛精度。
混合策略有个好处：解析求根能识别“两椭圆分离”的情况，因为无实根意味着约束不可满足；而数值优化直接返回一个最小距离点，用户如果没有后续校验，很容易把“零距离”当作“目标就在那里”，然后整条定位链路就翻车了。所以在设计 findEllIntersect 时，返回值里我通常会额外带一个 status 字段，0 表示找到正常交点，1 表示两椭圆间最小距离超过阈值，2 表示解析求根无解但数值优化收敛到了距离极小值。调用方看到非 0 时，不应该把结果当作定位结果使用。
4. 把 findEllIntersect 用到被动雷达定位：从时延测量到目标坐标的完整流程
4.1 数据预处理：时延估计和量测对齐
椭圆定位对时延差的精度非常敏感，因为椭圆长轴由时延差决定。若时延误差为 Δτ，则距离和误差为 c*Δτ，叠加上基线距离后，长轴就会偏移。实际接收机系统里，直达波与散射波的时延用互相关函数估计，采样率 100 MHz 时理论时延分辨率为 10 ns，对应距离误差约 3 米。城市环境里的多径会让互相关峰出现多个峰值，挑错峰会让时延差相差几百纳秒甚至几微秒，这是定位错误的头号来源。
预处理阶段需要完成三件事：
- 时延粗估计：用滑动相关在搜索窗口内找峰。
- 时延精估计：对相关峰做抛物线插值或频域插值，把分辨力提到亚采样级。
- 量测一致性检查：计算 ΔR = c*Δt 后，验证 |TR| + ΔR 大于 |TR|，也就是长轴大于基线的一半对应的情况。
这里建议把量测送入 findEllIntersect 前先做一步“门限淘汰”：ΔR 必须大于 0，且小于最大探测距离与基线差。如果 ΔR 为负数或超过物理范围，直接丢弃。
4.2 单站单发条件下的椭圆与测向交叉
最常见的工程配置是：一部接收机、一个已知外部辐射源。此时只有一个椭圆，不足以定位，需要增加测向信息。接收机用阵列天线测目标方位角 φ，那么目标位于从接收站出发的一条射线上，该射线与椭圆求交，一般得到一个交点，特殊情况下有两个。这个项目标题只提到椭圆法，但在工程实施中常常是把椭圆和测向联合使用，因为单椭圆在数学上本来就欠定。
射线与椭圆的求交比椭圆求交更简单。把射线写作 P = R + u * d，其中 d = [cos φ, sin φ]，代入椭圆隐式方程，得到一个关于 u 的二次方程，解出 u 后即可得到位置。这个步骤用 findEllIntersect 的思想稍微改造一下就能实现：把射线视为退化的椭圆——也就是短半轴趋近于 0 的椭圆——这样求交变成射线与椭圆求交。
不过这里要提醒一个容易踩的坑：测向角有误差时，射线和椭圆可能不相交。这时如果直接解二次方程，判别式小于 0 就没结果。工程上的后悔药是把测向角视为约束而不是硬条件：在测向角 ±3σ 范围内搜索，求椭圆到射线的最小距离点。最小距离点虽然不是严格意义上的“交点”，但这个结果在定位精度上往往比硬相交更可靠。
4.3 多站联合定位：多椭圆交点的选择和残差筛选
当系统可用多于一条椭圆约束时，比如两个不同方向的辐射源，或者两部接收机同时看同一个目标，就会出现多个椭圆求交的问题。N 个椭圆两两相交会产生很多伪交点，需要一种方法选真点。我的做法是分两步：
- 第一步，列出所有椭圆对的交点。
- 第二步，计算每个交点到每个椭圆的“代数距离”或“几何距离”，累加作残差。真实目标位置在所有椭圆上残差都最小。
具体残差可以用距离差来衡量：对椭圆 i，把候选点 P 到该椭圆的距离定义为 |P 到焦点1 距离 + P 到焦点2 距离 - 该椭圆焦点距之和|。残差最小的候选点就是最终目标位置候选；残差偏大超过门限，说明该候选点可能由杂波造成。这个做法不需要解高维方程，只需要大量重复二维求交，工程上稳定性足够。
def solve_multi_ellipse(ell_list, threshold=10.0): from itertools import combinations candidates = [] for (ell_a, ell_b) in combinations(ell_list, 2): pts = find_ell_intersect(ell_a, ell_b) if pts.size: candidates.extend(pts) if not candidates: return None scores = [] for p in candidates: total = 0.0 for ell in ell_list: total += ellipse_point_distance(ell, p) scores.append(total) best_idx = int(np.argmin(scores)) if scores[best_idx] < threshold: return candidates[best_idx] return None # 超门限，所有候选都不可信参数 threshold 的取值要根据时延误差来定。如果时延精度是 10 ns，对应距离误差 3 米，那么点到椭圆的距离残差门限可以设到 6-10 米；如果系统精度到亚米级，可以收紧到 2 米。这里的思路是能在存在干扰点的场景下提供鲁棒性。交汇出的候选点数量在噪声较大会变多，残差筛选能压制大部分假点。
4.4 坐标转换和时间同步在流程中的位置
不要忘了，椭圆法本质上是几何解算，它的误差下限取决于坐标和时间的同步精度。工程上常见方案是接收站和发射站位置由 GPS/北斗定位获取，接收站时间由驯服时钟或网络时间同步保证。如果发射站和接收站之间时基偏差超过 10 ns，直接换算成 3 米的位置误差。做整套流程时，坐标转换建议集中在一个模块里：
- 原始输入：经纬度、海拔。
- 转换到局部直角坐标：使用 UTM 投影或自定义当地切平面。建议在几百公里范围内使用自定义切平面，让基线长度不会因投影变形而误差过大。
- 时延测量：在信号处理前端完成，通过数据帧给定位模块送 ΔR 值。
如果你在写代码时基站坐标锁定到 UTM，但时延测量时把光速用 2.99792458e8，结果单位全要统一成米，不然椭圆长短轴差之毫厘，定位结果失之千里。
5. findEllIntersect 实战避坑：5 个会让定位结果玄学化的真实问题
5.1 两个椭圆平行或同焦，导致无穷多交点假象
现象：代码跑出来交点坐标一大片，而且每个都离真实目标很远，看起来像随机分布。原因：当两个椭圆拥有同一个焦点的连线且长轴差很小时，两椭圆非常接近重合或平行，数值求解时参数 t1 和 t2 的映射关系极不稳定，局部极小值遍布整个参数域。解决：在调用 findEllIntersect 前增加几何条件检查：计算两个椭圆的中心距离、长轴差、焦距方向夹角；若中心距离小于 20 米且长轴差在 5 米以内，程序中直接告警，不进入求交流程，改成使用测向线约束或其他办法。
5.2 时延互相关峰选错：多径让椭圆长轴直接偏移几十米
现象：定位结果在某个方向系统性偏了 30-50 米。原因：城市环境里直达波经过建筑物反射形成多径，散射波路径也有多条；相关峰选择器选了错误的峰。解决：在时延估计环节改用高分辨率算法，比如 MUSIC 或 ESPRIT 估计多径时延，取最先到达的那一路作为直达波；散射波则需要在相关输出上保留前几个峰并分别送入椭圆解算。每次输出目标位置前，检查对应的时延差是否大于发射站到接收站直达波时延，如果小于，直接丢弃。
5.3 优化器收敛到局部极小值，交点个数不对
现象：明明两个椭圆应该交于两个点，但函数只返回一个点，或返回一个错点。原因：网格初值太粗，收敛落到了局部极小值；特别是两个椭圆相切和接近相切时会这样。解决：把 n_initial 提高，但代价是运行时间增加到 36×36 次优化。更快的做法是先用解析四次求根得到候选，再用候选点做精确化。如果你的代码没有解析求根部分，至少把初始网格依据椭圆偏心率自适应：偏心率大时在近焦点区域加密采样，偏心率小时均匀采样。
5.4 基线长度不准，整个椭圆系整体平移
现象：即使没有测量噪声，解算结果也在某个方向偏移。原因：发射站和接收站的位置没有经过校准，基线误差 10 米，直接导致椭圆焦点偏移，椭圆求交结果整体偏移。解决：工程上使用“自校准”流程：让已知位置的目标做几次定位，反推椭圆参数误差，然后在运行时做修正。这个校准过程本质上是最小化残差的位置修正量，建议定期在整条链路运行前自动执行。
5.5 目标在基线的延长线上，椭圆变成线段的退化场景
现象：目标靠近发射站与接收站连线的外侧或延长线附近时，求解结果剧烈跳变，前后帧的定位结果不连续。原因：椭圆短半轴趋于零，目标位置在焦线上时方程退化，交点计算病态。解决：对短半轴 b < 1 米的情况不直接用椭圆求交，而是改用焦点距离解算：在焦线上，目标到两个焦点的距离和等于长轴乘以 2，先算出目标距离接收站的距离，再沿焦线方向外推。这个特殊处理只需要在调用 findEllIntersect 前增加一个 if 分支。
6. 验证与进阶：用蒙特卡洛仿真评估椭圆定位的精度上限
椭圆定位到这一步，能跑通流程，但一个负责任的工程交付还得回答“精度到底多少”。 我习惯的做法是搭建一个蒙特卡洛仿真平台，把时延噪声、基站位置噪声和测向噪声分别注入系统，统计不同位置上的均方根误差（RMSE），并与克拉美罗下界（CRLB）对比。
仿真参数设定示例：
| 参数 | 值 | 说明 |
|---|---|---|
| 发射站位置 | (0, 0) km | 外部辐射源 |
| 接收站位置 | (20, 0) km | 被动接收机 |
| 目标位置 | (12, 15) km | 典型低空目标 |
| 时延噪声 σ_t | 10 ns / 20 ns / 50 ns | 对应约 3m / 6m / 15m 距离误差 |
| 测向噪声 σ_angle | 0.5° / 1° / 2° | 取决于阵列口径 |
| 蒙特卡洛次数 | 1000 | 每次独立随机注入噪声 |
仿真代码中把 findEllIntersect 作为黑匣子调用。每个采样点得到目标位置估计，经过 1000 次统计 RMSE。对比 CRLB 时，注意椭圆法的等效时延精度受双基地几何配置影响严重：目标在基线垂直方向时，椭圆曲率大，定位精度高；目标在基线延长方向时，精度急剧下降。所以工程部署时尽量让接收站和辐射源围绕重点监视区域呈三角形布局，而不是一字排开。
做完仿真后，还有两个进阶方向可以做：一是把椭圆法和测向法做混合解算，利用不同量测方程的互补性提高单帧定位精度；二是针对时延差量测引入加权最小二乘，让高精度量测在解算中占更大权重。用 matlab 的定位工具箱或 python 自带 numpy/scipy 就够了，最终交付时把仿真代码和 findEllIntersect 一起放进自动化测试套件，每次修改都跑一轮回归，防止后续调整把验证过的算法弄坏。
另一个值得做的验证是现场跑车实验：装载接收机，围绕辐射源转几圈，用 GPS 记录真实轨迹，再和椭圆定位输出轨迹对比。这个验证能发现仿真里没建模的系统误差，比如天线相位中心偏差、通道不一致、时钟漂移。跑车的血泪经验是，仿真再漂亮，不经过现场数据校验的算法上线后一定会出幺蛾子。希望这套从几何推导、函数实现到验证的路径能帮到你。
本文还有配套的精品资源，点击获取
