三维路面不平度生成与RoadRunner导入：谐波叠加法实践指南

前两年做整车平顺性仿真，我一开始只用两条独立车辙剖面来应付路面输入，左右轮各一条一维序列，跑起来倒是能算，但真到了要把路面数据放进场景级驾驶仿真工具的时候，问题全冒出来了——车辆并不是只在两条轮辙上运动，悬架几何、侧倾响应、转向输入、路面渲染，全都需要一块有横向起伏的三维路面。后来我把生成方法从一维扩展到三维，核心算法用的还是大家熟悉的谐波叠加法，但把中间的数据链路彻底打通了：Matlab生成三维路面不平度，分别落成TXT和RDF格式，最后导入RoadRunner这类场景工具继续做仿真验证。

这篇文章就把整条链路展开讲。内容包含谐波叠加法的物理背景与公式系数来源、三维化展开的思路、Matlab完整可运行代码、TXT和RDF两种格式的细节，以及我实际使用中踩到的几个参数坑。适合正在做车辆平顺性、悬架载荷、耐久性分析，或者想把路面不平度数据弄进RoadRunner做场景级仿真的人参考。

1. 为什么非要做成“三维路面”：从两道车辙到一张地面网格

1.1 一维路面到底缺了什么

传统整车动力学里最常用的路面激励模型，是直接给左轮和右轮分别施加一条纵向不平度剖面。这种方法有个很实际的优势：计算量小、参数少、实现快。用谐波叠加法随手就能生成两条随机序列，中间甚至可以用相干函数控制左右轮的相关程度。做平顺性或疲劳分析时，这种输入确实够用。

但如果你想做的是轨迹级仿真，或者要把路面数据丢进可视化仿真场景，一维剖面的短板就非常明显。侧向有坡度的车道、路面横坡、轮迹之外的高程信息，这些都是两条一维曲线完全表达不了的。更关键的是，车辆在变道、转弯时，轮胎并不总是压在固定的两条轮迹上。路面必须是一张连续的三维网格，任何位置都能查得到高度。RoadRunner这类工具加载路面时，要的也是一个带顶点的曲面网格，而不是两条点列。

所以三维路面不平度，本质上解决的问题是：把一条沿纵向分布的随机高程，扩展成一张覆盖道路上所有位置的二维随机场，并且横向不同位置之间存在合理相关性。

1.2 三维路面不平度需要的目标：不平度等级与空间分辨率

要生成随机路面，第一件事是确定目标等级。国内车辆动力学领域最常用的还是国标GB/T 7031里给出的路面功率谱密度分级标准。核心定义是：

[ G_q(n)=G_q(n_0)\left(\frac{n}{n_0}\right)^{-2} ]

其中 (n_0=0.1 \text{ m}^{-1}) 为参考空间频率，(G_q(n_0)) 是路面不平度系数，单位是 (\text{m}^3)。A到H级路面对应的 (G_q(n_0)) 从低到高排列，其中车辆平顺性仿真最常用的几个等级我在下面给出。

做仿真之前先明确你要模拟的是什么路面。沥青路面选B级通常比较保守；如果研究的是耐久性和冲击载荷，C级甚至D级更容易激发悬架限位工况。有了这个系数，后面的幅值计算才有依据。

空间分辨率同样重要。路面纵向采样间隔我一般取0.05m，横向取0.1m到0.2m，这样既能保留足够的高频不平度信息，又不会让数据量爆炸。具体怎么和最高空间频率匹配，我在最后一部分专门讲。

2. 谐波叠加法核心公式与二维场构造思路

2.1 从国标功率谱密度出发：幅值系数为什么是 ( \sqrt{2G(n)\Delta n} )

谐波叠加法的思路是：把随机路面看成无数个正弦波的叠加，每个正弦波有自己的空间频率、振幅和随机相位。离散化之后，纵向剖面写成：

[ q(x)=\sum_{k=1}^{m} \sqrt{2G_q(n_k)\Delta n}\sin(2\pi n_k x+\varphi_k) ]

这里的 (n_k) 是第k个空间频率，(\varphi_k) 是均匀分布在 ([0,2\pi)) 的随机相位，(\Delta n) 是频率间隔。

幅值系数为什么是 (\sqrt{2G_q(n)\Delta n})，这一点很多代码直接抄，但没有理解。功率谱密度表示单位频率内的平均功率贡献。单边功率谱下，第k个频率区间 ([n_k-\Delta n/2, n_k+\Delta n/2]) 对路面方差的贡献是 (G_q(n_k)\Delta n)。写成单频正弦 (A_k\sin(2\pi n_k x+\varphi_k)) 后，它的均方值是 (A_k^2/2)。让这两者相等：

[ \frac{A_k^2}{2}=G_q(n_k)\Delta n ]

所以 (A_k=\sqrt{2G_q(n_k)\Delta n})。网上有些代码直接写 (\sqrt{G_q(n)\Delta n})，那是把单边谱和双边谱搞混了，生成的剖面整体幅值偏小，频谱积分后等效路面等级至少低了一档，这个问题在结果比对时很容易被忽略。

随机相位的作用是让不同频段之间没有确定关系。每次运行都是一张不同路面，但统计特性一样。这就带来一个工程问题：如果希望两次仿真之间只有路面种子不同，其他条件完全一致，那就必须在所有随机相位之前固定随机数发生器种子。我习惯把种子作为外部输入参数而不是函数内部随机生成，这样后续做正交试验、多组工况轮次时，每一轮路面都能精确复现。

2.2 怎样把一维叠加结果铺成二维面：横向剖线与横向平滑

把一维谐波叠加扩展成三维，最直接的方法并不是真的去计算二维功率谱。二维功率谱的测量数据不像一维国标谱那么规范，强行套二维模型，理论推导会很复杂，且很多参数没有实测支撑。

我采用的工程化做法是：先把道路在横向上分成若干条纵向剖线，每条剖线单独用一维谐波叠加生成，然后再沿横向做一次低通平滑。这样做有三个好处。

第一，每条剖线在纵向上严格保持目标功率谱，因为你没有改动任何沿x方向的数据处理过程。第二，横向剖线之间采用不同的随机相位序列，天然产生横向变化，左右轮迹不会完全一样。第三，横向平滑的作用相当于给道路表面加了一个空间低通滤波器，避免相邻剖线之间出现“刀切”感。这里有一个很容易误解的点：横向平滑确实会让相邻剖线变得相关，但并不会改变某一条固定剖线沿纵向的频谱，因为平滑是在每个x位置沿y方向做的，不是沿x方向做的。这点保证了最终路面在纵向上的统计特征和国标谱一致。

每条剖线完全独立时，路面在横向上会太碎，人眼看过去像砂纸，轮胎压过时左右轮高频激励完全不相关，和真实路面不符。所以剖线间距不宜太宽，两剖线之间最好通过插值或平滑建立过渡。代码里我用了横向移动平均进行平滑，核宽度大约0.5m，实测效果比较自然。

2.3 左右轮迹相关性：一个不花大代价的实用处理

有经验的读者会追问一个问题：左右轮迹在高频和低频段的相干程度本来就不一样。空间频率很低的长波，左右轮基本是刚性同步的，相当于整个车身在同一个起伏上；空间频率很高的短波，左右轮受到的地面激励几乎没有关系。

严格处理这个问题，需要引入相干函数模型。工程近似上，可以按空间频率分段处理：低于某个截止频率的成分，所有横向剖线共用同一组随机相位；高于截止频率的成分，每条剖线用独立的随机相位。这个截止频率通常和轮距有关，轮距越大，左右轮完全相干的频率上限越低。如果你只是做一般的垂向平顺性仿真，可以用比较折中的办法：让所有剖线频率下限附近的低频相位共享，高频相位自由，再通过横向平滑过渡。这样生成的路面既不会在数字上违反频谱特性，视觉上也比完全独立剖线自然得多。

下面的Matlab代码用了相对保守的方案，每一条剖线相位独立，但加入横向低通平滑来补偿相关性。如果想做更精细的相干控制，可以把随机相位部分按频段拆开处理。

3. Matlab实操：生成三维路面并写出TXT

3.1 总体参数设计

先定义参数块。路面长度、宽度、纵向和横向采样间隔、路面等级都是外部控制的。把路面等级作为字符串输入，函数内部用一个switch结构把等级映射到国标系数。

这里有个参数设计的关键点：空间频率下限和上限不能乱选。最低空间频率建议不低于 (1/L)，也就是说道路长度至少要能容纳最低频成分的一个完整波长。如果你用0.011 m⁻¹做下限，但路长只有20m，那这个超长波在20m长度内根本展不开，生成的轮廓其实是截断的。理论上谐波叠加并不是周期函数强制约束，但最低频波长远大于路长，对实际路段影响有限，只是路面整体会偏移一个低频趋势，往往表现为两端高度差明显。更合理的是让最低频对应的波长接近道路长度，或者干脆把 (n_{\min}=\max(0.011, 1/L)) 作为下限使用。

最高空间频率则由采样间隔决定。纵向采样间隔 (dx) 能表示的最高空间频率不超过 (1/(2dx))。当dx取0.05m时，可表示的最高空间频率是10 cycle/m，对应最短波长为0.1m。把谐波上限设成10而不是5，能保留更多路面细节；但如果上限设得太高而采样间隔不够密，就会产生混叠。

频率数量m也很关键。谐波叠加法的频率数量越多，相位分布越连续，生成路面越接近真实随机过程。m取500是能用，取1000以上更从容，生成时间也不会明显拉长。

3.2 主函数代码

下面这段函数完成了从路面等级参数到三维高程矩阵的全部计算。

function [x, y, Z, V, F] = gen3DRoadSurface(grade, L, B, dx, dy, m, seed) % gen3DRoadSurface 利用谐波叠加法生成三维路面不平度 % 输入: % grade - 路面等级字符, 如 'A','B','C','D' % L - 道路纵向长度, m % B - 道路横向宽度, m % dx - 纵向采样间隔, m % dy - 横向采样间隔, m % m - 谐波叠加频率数量 % seed - 随机数种子 % 输出: % x, y - 纵向、横向坐标向量 % Z - 高程矩阵, size = [length(x), length(y)] % V, F - 用于网格重建的顶点与三角面索引, 写TXT/RDF时用 n0 = 0.1; % 参考空间频率 switch upper(grade) case 'A', G0 = 16e-6; case 'B', G0 = 64e-6; case 'C', G0 = 256e-6; case 'D', G0 = 1024e-6; case 'E', G0 = 4096e-6; otherwise, error('unsupported road grade'); end rng(seed); % --- 空间频率参数 --- n_high = 1/(2*dx); % 采样最高可表示频率 n_low = max(0.011, 1/L); % 最低频率与路长匹配 n_list = linspace(n_low, n_high, m)'; dn = n_list(2) - n_list(1); % 幅值系数: sqrt(2 * G_q(n) * dn) G_list = G0 * (n_list / n0).^(-2); A_list = sqrt(2 * G_list * dn); % --- 坐标网格 --- x = (0:dx:L)'; y = (-B/2:dy:B/2); nx = length(x); ny = length(y); % 逐条纵剖面生成 Z = zeros(nx, ny); for j = 1:ny phases = 2 * pi * rand(m, 1); % 利用矩阵广播: 频率列向量 x 位移行向量 Z(:, j) = A_list' * sin(2 * pi * n_list * x' + phases); end % 横向平滑, 降低相邻剖线的“切刀感” win = max(3, round(0.5 / dy)); if mod(win, 2) == 0 win = win + 1; end for i = 1:nx Z(i, :) = smoothdata(Z(i, :), 'movmean', win); end % --- 转换成顶点与三角面 --- [Xmesh, Ymesh] = ndgrid(x, y); ids = reshape(1:nx*ny, nx, ny); numV = nx * ny; V = zeros(numV, 4); V(:, 1) = ids(:); V(:, 2) = Xmesh(:); V(:, 3) = Ymesh(:); V(:, 4) = Z(:); numF = 2 * (nx - 1) * (ny - 1); F = zeros(numF, 3); cnt = 0; for i = 1:nx-1 for j = 1:ny-1 a = ids(i, j); b = ids(i+1, j); c = ids(i+1, j+1); d = ids(i, j+1); cnt = cnt + 1; F(cnt, :) = [a, b, c]; cnt = cnt + 1; F(cnt, :) = [a, c, d]; end end end

这段代码里有一个容易被忽略的细节：(G_q(n)) 随频率按 (n^{-2}) 衰减，低频段能量占主导。所以空间频率列表用等间隔linspace时，低频段每个频段间隔内的能量变化很剧烈，需要足够密的频率点才能准确覆盖谱密度曲线。如果把m设得太小，比如几十个点，生成的路面会丢失低频包络细节，看起来像简单的几个大波形叠加。通常m取500以上，代码中取1000也很快。

3.3 验证与可视化

生成完三维矩阵，别急着写文件，先做一个快速验证。最有效的验证是把某一条纵向剖线提取出来，做功率谱估计，和目标功率谱密度画在同一张对数坐标图上。

zCenter = Z(:, floor(ny/2)); zCenter = zCenter - mean(zCenter); fs = 1/dx; [pxx, f] = periodogram(zCenter, hann(length(zCenter)), 2^nextpow2(length(zCenter)), fs); targetG = G0 * (f(n_list range?) / n0).^(-2); loglog(f, pxx); hold on; % 在地图画目标谱曲线

目标谱 (G_q(n)) 一般是单调递减的直线（在双对数坐标下斜率为-2）。如果生成剖面的功率谱曲线围绕这条直线上下波动，且波动幅度不夸张，说明路面频谱能量分布是对的。如果整体偏低，就要检查幅值系数是否少了根号2；如果高频区域离直线太远，多半是频率数量太少或采样间隔不够小。

也可以直接算标准差。B级路面理论标准差通常在4~10mm级别，如果算出来偏差过大，优先检查代码里谱密度单位和频率范围的设置。

三维路面生成后，还可以用surf函数快速扫一眼：

surf(y, x, Z, 'EdgeColor', 'none'); xlabel('横向位置/m'); ylabel('纵向位置/m'); zlabel('高程/m'); axis equal;

如果视觉上横向能看到明显条带，说明横向平滑窗口设得不够，或者剖线间距太远。二维路面图里出现明显横向条纹，会在仿真时给轮胎施加不真实的高频横向激励。

3.4 落盘TXT：节点表和面索引表分开写

TXT文件不需要搞得很复杂。路面三维数据给下游工具时，最通用的形式是“顶点表+三角面索引表”。顶点表保存每个网格点的ID，x坐标，y坐标，z高程。面索引表保存每个三角形由哪三个顶点构成。

function writeRoadTXT(nodeFile, faceFile, V, F) % 写出节点文件和面索引文件 fid = fopen(nodeFile, 'w'); fprintf(fid, 'ID x y z\n'); for i = 1:size(V,1) fprintf(fid, '%d %.6f %.6f %.6f\n', V(i,1), V(i,2), V(i,3), V(i,4)); end fclose(fid); fid = fopen(faceFile, 'w'); fprintf(fid, 'a b c\n'); for i = 1:size(F,1) fprintf(fid, '%d %d %d\n', F(i,1), F(i,2), F(i,3)); end fclose(fid); end

这里我特意把节点和面分开保存，而不是混在一个大文件里。原因很实际：后续转RDF格式时，顶点表可以直接映射成节点列表，面索引表直接映射成三角形网格；如果中间想做抽稀、裁剪或网格重建，分文件处理也更方便。

4. TXT如何转RDF：RoadRunner导入前的数据换装

4.1 RDF需要哪些信息

很多做动力学的人一听“RDF格式”会觉得这是某个标准规范里的专用格式。实际工程中，RDF往往指目标仿真工具定义的道路描述文件，不同工具的字段不完全一样。RoadRunner生态里经常接触到的RDF，作用是把路面的几何、等级、网格等描述信息，以规范化文本形式组织起来供场景工具读取。

所以TXT转RDF，绝对不是在文件名上把.txt改成.rdp或.rdf就算完成，而是要做一次结构性映射。TXT里只有纯数据，RDF里还需要有元信息：路面等级、网格步长、单位、顶点数量、三角面数量。这些信息对仿真工具来说并不是冗余，RoadRunner或其他工具导入时要用它们来配置场景和碰撞物理属性。

我这里给出一个可扩展的RDF组织方式，字段设计为节点+索引+元数据。如果你使用的RoadRunner版本或第三方工具对字段名有特殊要求，按我下面这个映射关系改标签名即可，核心内容不变。

4.2 从TXT到RDF的转换代码

在Matlab里直接由V和F生成RDF，代码结构如下：

function writeRoadRDF(rdfFile, grade, V, F, dx, dy, L, B) % 把三维路面网格写成 RoadRunner 可导入的道路描述文件 fid = fopen(rdfFile, 'w'); fprintf(fid, '<?xml version="1.0" encoding="UTF-8"?>\n'); fprintf(fid, '<RoadDescription FileType="RoadSurfaceMesh" Version="1.0">\n'); fprintf(fid, ' <Attributes>\n'); fprintf(fid, ' <Grade>%s</Grade>\n', grade); fprintf(fid, ' <Length>%.3f</Length>\n', L); fprintf(fid, ' <Width>%.3f</Width>\n', B); fprintf(fid, ' <Unit>m</Unit>\n'); fprintf(fid, ' <Dx>%.4f</Dx>\n', dx); fprintf(fid, ' <Dy>%.4f</Dy>\n', dy); fprintf(fid, ' </Attributes>\n'); fprintf(fid, ' <Vertices Count="%d">\n', size(V, 1)); for i = 1:size(V, 1) fprintf(fid, ' <Vertex ID="%d" X="%.6f" Y="%.6f" Z="%.6f"/>\n', ... V(i,1), V(i,2), V(i,3), V(i,4)); end fprintf(fid, ' </Vertices>\n'); fprintf(fid, ' <Faces Count="%d">\n', size(F, 1)); for i = 1:size(F, 1) fprintf(fid, ' <Face ID="%d" V1="%d" V2="%d" V3="%d"/>\n', ... i, F(i,1), F(i,2), F(i,3)); end fprintf(fid, ' </Faces>\n'); fprintf(fid, '</RoadDescription>\n'); fclose(fid); end

这个文件结构完全可以从TXT节点表和面索引表直接转换，不需要重新采样。RoadRunner导入网格时，真正关心的是三件事：几何坐标、三角连接关系、单位。你在TXT里已经包含了前两项，最后把单位、等级、网格间距塞进Attributes段，RDF就完整了。

4.3 坐标系和单位检查

这一节是实际导入时最容易翻车的部分。Matlab里三维数据默认是右手坐标系，x代表道路延伸方向，y代表横向，z代表高程。但RoadRunner的地面或道路网格系不一定默认就是x纵y横z上。不同版本、不同导入通道，轴系可能不同。不要等到RoadRunner里路面整个竖起来或者倒过来，再回去排查。

一个务实的方法：先构造一个非常小的路面，比如只有3x3个顶点的平面，输出成RDF，导入目标工具。如果显示正常，再把完整路面写进去。如果发现路面被翻转，多半是y轴或z轴需要交换；如果道路方向不对，则可能是x方向定义和工具默认前进方向不同。

单位问题同样隐蔽。我在代码里强制输出坐标单位为米，并且把单位写入RDF的Attributes段。但有些工具导入时不读这个字段，而是直接按工具默认单位解释，也就是“数字是多少就按多少米/厘米处理”。如果导入后路面尺度明显不对，先确认工具默认单位是否和你输出单位一致，而不是急着改Matlab代码。

5. 几个真正影响仿真结果的参数“坑”

5.1 路面长度与空间频率下限的匹配

谐波叠加法的频率范围选择，直接影响道路两端是否会产生明显的整体趋势。如果空间频率下限过低，比如取0.005m⁻¹，而这个值的半波长已经大于路长，那这条长波在道路长度内只会显示出一个单调上升或下降的斜坡，让一段明明应该统计平稳的路面在视觉上变成“上坡”或者“下坡”。车辆在这种路面上行驶会产生额外重力分量。

处理方式是让 (n_{\min}) 和路长相匹配。最省事的做法是 (n_{\min}=1/L)，意思是道路长度刚好对应最低频成分的一个波长。这样最低频长波在一个路段内可以完整展开，不会出现强制的单调趋势。如果L