简介：交通规划中的交通量预测是城市路网设计的关键环节。压缩包面向交通规划课程学习者与相关从业者，以4个MATLAB脚本分别实现弗兰特法、底特律法、平均系数法与佛尼斯法，覆盖基于历史增长、设施饱和度、平均增长率和路网迭代分配等多种预测思路，适合作为课后习题编程的参考模板。弗兰特法适合快速估算，底特律法更关注道路饱和度的影响，平均系数法简单易行，佛尼斯法则可细致刻画路网拓扑与流量分配。包内共4个m文件，整体仅2KB，结构精简，便于直接读取和修改参数，四个脚本均包含从数据读取到结果输出的完整流程。已有931人学习下载。通过运行这些脚本，可快速掌握各方法的计算流程、适用场景及参数调整方式，对比不同模型的预测差异，为实际交通规划中的模型选择提供直观依据。

1. 交通规划里的四种增长系数法：一张OD表怎么跑出未来分布

做交通规划课后习题或者课程设计时，最容易卡住的不是四阶段法本身的概念，而是“手头只有一张现状OD表，外加未来各交通小区的发生量和吸引量，怎么用弗兰特法、底特律法、平均系数法、佛尼斯法跑出未来OD矩阵，还要保证行和、列和都对得上”。这套资源把四种增长系数法串在同一个数据集上，从公式拆解、手算核对到可运行的Python代码全部打包，适合刚学完交通分布预测、正在刷课后编程题的学生，也适合要快速出预测结果却不想从零调公式的从业者。直接换数据、调收敛参数就能复现，省掉大部分试错时间。

2. 四种方法公式拆解：从平均系数到佛尼斯法的迭代逻辑

2.1 平均系数法：双系数取平均，入门最简单但收敛最慢

平均系数法也叫均匀增长法，它的核心逻辑是假设未来OD量同时受出发小区发生量增长和到达小区吸引量增长的影响，取两者的平均值作为增长因子。设现状OD矩阵为T，未来发生量为U，未来吸引量为V，那么每一行的小区i有一个发生量增长系数α_i = U_i / O_i，每一列的小区j有一个吸引量增长系数β_j = V_j / D_j，其中O_i是现状第i行的行和，D_j是现状第j列的列和。

单次迭代的更新公式是：

t_ij = T_ij × (α_i + β_j) / 2

这看起来很简单，但注意两点：第一，α_i和β_j是用上一轮的矩阵重新算的，而不是用初始OD矩阵算一次就结束；第二，平均系数法算完一轮之后，行和与列和不会自动等于U和V，必须继续迭代，直到行列和的误差小于设定的阈值。

实操中我一般把收敛条件设成“行和与U的相对偏差最大值小于1e-4”或者“不平衡率小于1%”，而不是追求绝对相等。平均系数法对初始数据不敏感，即使OD矩阵里有零元素也能跑，但因为增长因子只做了简单平均，收敛速度偏慢，通常需要十几轮甚至几十轮才能稳住。它适合快速验证量级，或者作为其他方法的结果对照基准，不适合精度要求高的正式预测。

2.2 底特律法与弗兰特法：全域校准因子与局部修正因子的取舍

底特律法是对平均系数法的直接修正。它不再取平均值，而是把发生量增长系数和吸引量增长系数相乘，同时引入一个全域校准因子γ来保证总量守恒。单步公式为：

t_ij = T_ij × α_i × β_j / γ

其中γ = ΣΣ(T_ij × α_i) / ΣΣT_ij，也就是以发生量为权重算出的全区平均增长率。这里有个容易混淆的点：不同教材对γ的写法不太一样，有的是用发生量加权，有的是用吸引量加权，有的甚至直接用ΣU/ΣO。实际迭代时，因为每一轮都会重新计算γ，初始定义略有不同并不会改变最终收敛结果，只会影响中间步的数值轨迹。

底特律法的优点是每轮更新都带着全区总量约束，整体收敛比平均系数法稳定。缺点是它用同一个γ修正全矩阵，无法反映某个特定小区因为新建开发区、人口导入等原因产生的局部增长差异，如果区域内部发展不均衡，底特律法的预测结果会偏“钝”。

弗兰特法正是针对这个问题提出的。它在底特律法乘积项的基础上，增加了两个局部校准因子L_i和L_j：

t_ij = T_ij × α_i × β_j × (L_i + L_j) / 2

其中L_i = O_i / Σ_j(T_ij × β_j)，L_j = D_j / Σ_i(T_ij × α_i)。L_i衡量的是第i行在考虑吸引量增长后的残余不平衡，L_j衡量的是第j列在考虑发生量增长后的残余不平衡。这两个因子逐行、逐列地修正，相当于在全局总量基础上叠加了局部特征。

弗兰特法比底特律法更容易收敛，对非均衡发展的区域适用性更强，但代价是计算量增加，每轮迭代都需要分别算行方向和列方向的校准因子。如果数据里出现某一行或某一列全为零，L_i或L_j会直接除零，这是用弗兰特法最常见的翻车点之一，后面避坑章节会专门讲。

推荐的使用场景是：区域间发展差异明显、个别小区的增长率是其他区域两倍以上时，优先用弗兰特法；如果各小区增长率相差不大、总量平衡更重要，底特律法就够了。

2.3 佛尼斯法：行、列交替调整，收敛快但顺序有讲究

佛尼斯法走的是另一条路。它不直接构造一个统一的增长公式，而是把行约束和列约束分开，交替施加。每一轮分两步：先用发生量U_i调整各行，使行和等于U_i；再用吸引量V_j调整各列，使列和等于V_j。反复交替，行列误差同时缩到阈值以内。

具体步骤是：

第一步，行调整：t'_ij = T_ij × α_i，其中α_i = U_i / Σ_j T_ij。算完后每一行的行和等于对应的U_i，但列和已经被破坏。
第二步，列调整：t''_ij = t'_ij × β_j，其中β_j = V_j / Σ_i t'_ij。算完后每一列的列和等于对应的V_j，但行和又偏了。

之后反复以t''_ij为基础重复第一步和第二步。理论上两步交替可以同时逼近行列约束，实际迭代次数一般在5到20次之间，明显少于平均系数法。

先调行还是先调列，会影响中间步骤的数值，但对最终收敛结果影响很小，只要最终行列不平衡率都低于阈值即可。不过建议统一按“先调行再调列”的流程写，方便别人复核代码。佛尼斯法有一个特性：它天然保证每一步矩阵元素非负，只要初始OD非负，中间过程和最终结果都不会出现负值，这在处理包含零流量的OD对时非常友好。

选择建议：如果给到的数据已经是总量守恒的（ΣU = ΣV），佛尼斯法是最稳的选择，收敛快、无负数、不除零。如果总量本身不平衡，先做总量平衡再套佛尼斯法；硬跑的话列约束永远满足不了，会陷入死循环。

3. 用Python实现四种增长系数法：代码、参数与输出对照

3.1 基础数据准备：OD矩阵、发生量、吸引量的数据结构

实现之前先把数据格式固定下来。我常用的做法是用NumPy二维数组存OD矩阵，用两个一维数组存未来发生量和吸引量，顺序必须与矩阵的行列索引一致。这里以4个小区的对称OD矩阵为例给出完整代码。

import numpy as np # 现状OD矩阵，4个交通小区，T[i][j]表示从i区到j区的现状出行量 T = np.array([ [0, 150, 180, 100], [150, 0, 120, 90], [180, 120, 0, 110], [100, 90, 110, 0] ], dtype=float) # 未来各小区发生量（行和的目标值） U = np.array([420, 460, 520, 380], dtype=float) # 未来各小区吸引量（列和的目标值） V = np.array([450, 480, 500, 350], dtype=float) # 校验基础条件：现状总量与未来总量是否守恒 print("现状总量:", T.sum()) print("未来总量(U和):", U.sum()) print("未来总量(V和):", V.sum())

这段代码里，T必须用float类型，否则后续迭代中的除法会直接把小数截断。U和V的顺序不能错，U[0]对应第一行的小区1，V[0]对应第一列的小区1。打印总量是为了确认ΣU等于ΣV，这是所有增长系数法能够正常收敛的前提；如果不相等，先做总量平衡，不能直接往下跑。

3.2 四种方法的函数实现：迭代主体与收敛条件

下面把四种方法各写成一个函数。我统一用“当前矩阵T_curr的某一行和或列和与目标值之间的最大相对误差”作为收敛判断，阈值默认1e-4。这样比直接用绝对差更合理，因为不同量级的数据对绝对差的敏感度完全不同。

def check_convergence(T_curr, U, V, tol=1e-4): """计算行和与列和的最大相对误差，返回是否收敛""" row_err = np.max(np.abs(T_curr.sum(axis=1) - U) / U) col_err = np.max(np.abs(T_curr.sum(axis=0) - V) / V) return max(row_err, col_err) < tol, row_err, col_err def average_growth(T, U, V, max_iter=1000, tol=1e-4): """平均系数法：t_ij = T_ij * (alpha_i + beta_j) / 2""" T_curr = T.copy() for _ in range(max_iter): O = T_curr.sum(axis=1) D = T_curr.sum(axis=0) alpha = U / O beta = V / D T_next = T_curr * ((alpha[:, None] + beta[None, :]) / 2) converged, r_err, c_err = check_convergence(T_next, U, V, tol) T_curr = T_next if converged: break return T_curr, r_err, c_err def detroit_method(T, U, V, max_iter=1000, tol=1e-4): """底特律法：t_ij = T_ij * alpha_i * beta_j / gamma""" T_curr = T.copy() for _ in range(max_iter): O = T_curr.sum(axis=1) D = T_curr.sum(axis=0) alpha = U / O beta = V / D # 全域校准因子：按发生量加权的平均增长率 gamma = (T_curr * alpha[:, None]).sum() / T_curr.sum() T_next = T_curr * alpha[:, None] * beta[None, :] / gamma converged, r_err, c_err = check_convergence(T_next, U, V, tol) T_curr = T_next if converged: break return T_curr, r_err, c_err def fratar_method(T, U, V, max_iter=1000, tol=1e-4): """弗兰特法：t_ij = T_ij * alpha_i * beta_j * (L_i + L_j) / 2""" T_curr = T.copy() for _ in range(max_iter): O = T_curr.sum(axis=1) D = T_curr.sum(axis=0) alpha = U / O beta = V / D # 行校准因子和列校准因子 L_i = O / (T_curr * beta[None, :]).sum(axis=1) L_j = D / (T_curr * alpha[:, None]).sum(axis=0) T_next = T_curr * alpha[:, None] * beta[None, :] * (L_i[:, None] + L_j[None, :]) / 2 converged, r_err, c_err = check_convergence(T_next, U, V, tol) T_curr = T_next if converged: break return T_curr, r_err, c_err def furness_method(T, U, V, max_iter=1000, tol=1e-4): """佛尼斯法：先调行再调列，交替迭代""" T_curr = T.copy() for _ in range(max_iter): # 行调整：让行和等于U alpha = U / T_curr.sum(axis=1) T_curr = T_curr * alpha[:, None] # 列调整：让列和等于V beta = V / T_curr.sum(axis=0) T_curr = T_curr * beta[None, :] converged, r_err, c_err = check_convergence(T_curr, U, V, tol) if converged: break return T_curr, r_err, c_err

逐段解释一下逻辑。average_growth函数中，alpha[:, None]把一维数组扩成列向量，beta[None, :]扩成行向量，两者相加再除以2，就得到了一个与T同形状的增长因子矩阵，与T逐元素相乘即可。check_convergence里的(U)分母不能省略，它把绝对误差归一化成相对误差；如果U中某个值为0，需要先做数据清洗，否则这里会除零。detroit_method里的gamma每次迭代都在变，这是正常的，它负责把总量拉回平衡。fratar_method中L_i和L_j的计算要特别注意括号，乘beta[None, :]是让beta广播到每一行，然后按行求和；如果这里错写成(beta)就会得到标量，整个方法就退化了。furness_method是四个函数里最简单的，因为不需要显式计算增长系数，直接对矩阵做缩放就行。

3.3 输出对比与收敛判断：行列和不平衡率的检验逻辑

四个函数跑完后，直接把结果并排输出，同时打印各自的迭代轮数和最终行列不平衡率。

# 统一执行并输出对比 methods = { "平均系数法": average_growth, "底特律法": detroit_method, "弗兰特法": fratar_method, "佛尼斯法": furness_method, } for name, func in methods.items(): T_result, r_err, c_err = func(T, U, V) print(f"{name}:") print(np.round(T_result, 2)) print(f"行和误差: {r_err:.2e}, 列和误差: {c_err:.2e}") print(f"总量: {T_result.sum():.2f}\n")

运行这段代码后会发现，四种方法在满足收敛条件时输出的总量都接近未来总量1780，但矩阵分布差别明显。平均系数法的分布最“均匀”，因为增长因子被取平均后拉平了差异；佛尼斯法对局部增长的反应最敏感，某个增长率特别高的小区，它的对外出行量会被快速放大。这里我建议把收敛阈值从1e-4放宽到1e-3再对比一次，观察矩阵变化是否被压在一个合理范围内，如果1e-4和1e-3算出的矩阵差异小于1%，说明结果稳定，阈值再小也没有实际意义。

注意这个输出结果只是验证代码正确性的中间件，不是最终成果。要判断四种方法谁更合理，需要结合第5章的守恒校验和实际业务背景，单纯比数字没有意义。

4. 避坑指南：迭代不收敛、除零错误与结果对不上的五个排查点

4.1 现象：迭代震荡，行列和来回跳

跑avg_growth或detroit_method时，行和误差不是单调递减，而是在某个区间来回震荡，甚至越跑误差越大。原因通常是收敛阈值设得太严，比如1e-8，远超过了浮点计算能达到的精度范围；另外当初始矩阵某些元素极小（接近0）时，单轮更新的增益过大，也会产生过冲。解决方法是先把阈值放宽到1e-3跑一遍，观察误差变化曲线；如果震荡，取最近20轮的平均值作为最终结果，而不是用最后一轮的瞬时值。我在实际刷题时会打印每轮的r_err，看到误差在1e-4附近反复横跳就直接收手，取上一轮的值。

4.2 现象：除零报错或出现NaN

弗兰特法最容易踩这个坑。L_i = O_i / Σ_j(T_ij × β_j)，如果某一行现状OD全是0，分母就是0，得到inf，矩阵直接变成NaN。另外，当某个小区未来吸引量V_j恰好等于0时，β_j也是0，分母依然有问题。解决方法是跑之前先做一次数据完整性检查：凡是行和或列和为0的小区，要么从矩阵中剔除，要么给它补一个很小的基准值（比如1e-6），并在结果输出后注明该小区数据为占位值。千万不要直接给分母加个随机小量，那样会把矩阵数值污染掉，后面算不平衡率时误差全被掩盖了。

4.3 现象：底特律法结果和教材示例差很多

同一道课后题，用底特律法算出来和教材答案差了5%甚至更多，但平均系数法和佛尼斯法都能对上。常见原因是底特律法的γ取法不一致。有的教材γ = ΣΣ(T_ij × α_i) / ΣΣT_ij，有的教材γ = ΣU / ΣO，有的直接取(ΣU + ΣV) / (ΣO + ΣD)。这三者在迭代初期给出的矩阵完全不同，虽然理论上都能收敛到同一结果，但如果你只迭代了有限轮次（比如教材手算只做了3轮），差异就会残留下来。我的做法是：教材只有3轮手算过程的习题，就按教材用的γ定义手动对齐；代码实现则统一用“以发生量加权”的版本并迭代到收敛，不跟手算结果做逐轮比对，只比对最终行列和。

4.4 现象：佛尼斯法先调行还是先调列，结果不一样

同一个初始矩阵，先调行得到的结果和先调列得到的结果在中间轮次里明显不同，尤其在迭代次数少的情况下。这个现象是正常的，因为行调整先满足了行约束，会破坏列约束；列调整再修正方向，两者此消彼长。当迭代充分收敛后，结果趋向同一组解，但收敛速度不同。我的建议是固定成“先调行再调列”，并且在代码注释里写清楚，这样别人复核时能复现你的数值轨迹。如果你发现两种顺序下最终结果仍然差很多，说明收敛条件没有真正满足，应把阈值的判断同时放到行和列上，而不是只看其中一个。

4.5 现象：代码算出来和课后习题答案对不上

这个是刷题场景下最高频的问题。原因通常是三个：一是教材手算时每一步都四舍五入了，代码里没有做舍入，误差在小数点后累积；二是教材只迭代了3到5轮，代码自动迭代到1e-4才停，数值自然不同；三是教材答案用平均系数法但全程没有迭代，只做了单轮计算。排查时先看习题题干的迭代要求，如果写了“迭代直至收敛”，就用代码默认行为；如果答案只给了三步手算结果，就把max_iter固定成3或5跑一遍做对照，别让代码自动迭代。另外比较时不要逐元素比，而是比较行和、列和是否与原题给出的U、V匹配，匹配就是做对了。

5. 一个值得养成的习惯：用三重守恒校验锁定计算误差

四种方法跑完，得到四张未来OD矩阵，哪张能用？我每次都会强制过一遍三重守恒校验，不通过直接弃用。第一重是总量守恒，未来矩阵的总和必须等于ΣU也等于ΣV；第二重是行约束，每一行的行和对应的U_i偏差不能超过收敛阈值；第三重是列约束，每一列的列和对应的V_j偏差同样要达标。三重校验用一段代码就能完成。

def verify_conservation(T_pred, U, V, tol=1e-3): """返回校验结果：总量、行最大误差、列最大误差""" total = T_pred.sum() row_err = np.max(np.abs(T_pred.sum(axis=1) - U) / U) col_err = np.max(np.abs(T_pred.sum(axis=0) - V) / V) total_pass = abs(total - U.sum()) / U.sum() < tol row_pass = row_err < tol col_pass = col_err < tol print(f"总量: {total:.2f} vs 目标: {U.sum():.2f}, 通过: {total_pass}") print(f"行最大相对误差: {row_err:.2e}, 通过: {row_pass}") print(f"列最大相对误差: {col_err:.2e}, 通过: {col_pass}") return total_pass and row_pass and col_pass # 以弗兰特法结果为例做校验 T_fratar, _, _ = fratar_method(T, U, V) verify_conservation(T_fratar, U, V)

还有一个容易被忽略的进阶技巧：确认结果不仅守恒，还要符合“非对称增长”的业务直觉。把预测矩阵和现状矩阵逐元素求增长倍率，如果某个OD对的倍率明显异常（比如从1.2突然跳到0.7），大概率是数据录入错误或初始某对OD值过小导致的数值噪声。我会把增长倍率矩阵打印出来，人工扫一遍，比写一百行校验代码都管用。

从那以后，我每次跑增长系数法都强制走一遍“总量守恒、行约束、列约束、增长倍率检查”四步，缺一不可。希望这套思路也能帮你少走弯路。

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