简介：本资源是一套面向海洋工程、环境科学及气候研究领域的SWAN波浪模型Python实现方案，聚焦风浪边界条件生成与二维能谱构建，适用于具备Python基础与数据分析经验的科研人员和工程师。资源包共13个文件，含6个Jupyter Notebook（如swan_stats.ipynb用于谱分析与可视化、write_GFS_wind.ipynb和write_TPAR_era5.ipynb等用于多源风场与再分析数据格式转换）、3个Python脚本（manuscript_functions.py与spec2d.py提供核心算法封装）、2个备份文件及1个README说明文档，整体压缩后仅279KB，轻量易部署。已有84人学习下载。读者可直接复用完整的数据预处理—谱计算—模型输入生成全流程代码，获得CSIRO历史海况、ERA5/GFS风场适配、TPAR格式输出等关键能力，并通过Notebook交互式调试理解波浪分割逻辑与二维能谱统计建模原理。
1. 项目概述：从“拍脑袋”到“算出来”的风浪模拟革新
在海岸工程、海洋预报和船舶设计领域，准确模拟海浪是确保一切计算和决策可靠性的基石。过去，很多项目在设置波浪边界条件时，往往依赖于经验公式、历史统计资料，甚至是“拍脑袋”定一个参数。这种方法在面对复杂多变的海况，特别是风浪与涌浪混合的实际情况时，常常力不从心，导致模拟结果与实测数据偏差巨大，轻则影响设计冗余度造成浪费，重则可能埋下安全隐患。我参与过的一个近海风电项目就曾吃过亏，初期用简化谱模型算出来的基础受力偏小，后期补强成本飙升。这个痛点，正是“SWAN代码实现：基于波浪分割与数据分析的风浪边界条件生成及二维谱构建方法”要解决的核心问题。
简单来说，这不是一个简单的软件使用教程，而是一套从原始观测数据出发，通过严谨的数据处理和分析，逆向工程出符合物理规律且能被第三代海浪模型SWAN直接使用的、高精度二维波浪谱边界条件的完整方法论。SWAN模型本身功能强大，但其输入边界——二维波浪谱的构建，一直是门槛所在。本方法将“数据分析”作为桥梁，把杂乱无章的浮标、测波仪数据，转化为蕴含风浪、涌浪各自能量、方向、频率特征的“结构化信息”，进而生成SWAN所需的谱文件。其价值在于，它将波浪边界条件的设定从一门“艺术”变成了可重复、可验证的“科学”，特别适合需要精细化模拟的科研与工程项目，比如台风浪过程再现、复杂地形下的波浪折射绕射分析、以及海洋可再生能源的场址评估。
2. 核心思路拆解：为何要分割风浪与涌浪？
要理解这个方法，首先得打破一个常见误区：不是所有海浪都是一回事。一片海面上的波浪，其能量来源主要有两种：当地风直接吹拂产生的风浪，以及从遥远风暴区传播过来的涌浪。这两者在特性上截然不同：风浪周期短、波陡大、方向谱较宽；涌浪周期长、波面平滑、方向性集中。在SWAN这类谱模型中，如果将它们混为一谈，用一个单一的谱形（如JONSWAP谱）来概括，就相当于用一首歌的旋律去描述整个交响乐团，必然会丢失大量关键信息，导致模型无法准确模拟波浪在传播过程中的演变，尤其是非线性相互作用和方向分布的变化。
因此，本方法的核心思路可以概括为“先分家，再建模，后合成”。
第一步：先分家（波浪分割）。利用时间序列分析（如波面高程数据）或谱分析方法，从观测到的混合波浪场中，将风浪和涌浪的能量分离出来。这不仅仅是能量大小的分离，更是对各自主导频率范围和主波向的识别。常用的分割方法有频谱分割法（在频率域设定分割阈值）和参数化方法（如利用风速、风时等气象参数辅助判断）。这一步是后续所有工作的基础，分割的准确性直接决定最终边界条件的质量。
第二步：再建模（单成分谱构建）。对于分割出来的风浪部分，我们通常采用与风场参数相关的谱模型，如JONSWAP谱或TMA谱，并根据分割结果调整其尖度因子、峰频等参数。对于涌浪部分，则可能采用更简单的谱形（如Pierson-Moskowitz谱的变形），或直接对分割后的涌浪频谱进行参数化拟合。关键在于，要为风浪和涌浪分别构建一个一维频率谱。
第三步：后合成（二维谱构建）。一维频率谱只告诉我们能量在不同频率上的分布，但波浪是有方向的。因此，需要将上一步得到的一维谱与一个方向分布函数结合，扩展成二维谱。方向分布函数通常采用余弦幂次方模型，其分布宽度参数需要根据数据分析结果（如测向浮标数据）或经验关系来确定。最终，风浪二维谱和涌浪二维谱会进行线性叠加（能量相加），生成一个完整的、包含双峰甚至多峰特征的二维方向频率谱，这个谱文件就是SWAN模型可以直接读取的边界条件。
注意：这里说的“线性叠加”是指能量谱的相加，在物理上基于线性波理论的假设，即不同波成分之间独立传播。在实际操作中，对于强非线性的情况需要谨慎评估。
3. 数据处理与波浪分割实战
理论清晰后，我们进入实战环节。假设我们拥有某站点一段时间的逐时观测数据，包括有效波高、谱峰周期、平均波向、以及风速和风向。我们的目标是为一个持续三天的风暴过程生成SWAN边界条件。
3.1 数据准备与质控
原始数据往往存在缺失、异常或格式不统一的问题。第一步永远是数据清洗。
格式标准化：将来自不同仪器（如波浪浮标、ADCP、气象站）的数据，统一到相同的时间戳（例如UTC时间）和数据结构中。我习惯使用Python的Pandas库来处理，创建一个以时间为索引的DataFrame。
import pandas as pd # 假设有波浪数据和风数据两个CSV文件 wave_df = pd.read_csv('wave_buoy.csv', parse_dates=['time'], index_col='time') wind_df = pd.read_csv('wind_station.csv', parse_dates=['time'], index_col='time') # 按小时重采样并对齐 wave_hourly = wave_df.resample('1H').mean() wind_hourly = wind_df.resample('1H').mean() # 合并数据框 merged_df = pd.merge(wave_hourly, wind_hourly, left_index=True, right_index=True, how='inner', suffixes=('_wave', '_wind'))异常值处理：识别并处理物理上不可能的值。例如，有效波高为负值或极大值，周期超出合理范围（如小于3秒或大于25秒）。可以采用阈值法或统计方法（如3σ原则）进行筛选和插补。
# 示例：剔除有效波高异常大的数据 reasonable_max_hs = 20.0 # 根据海域设定合理上限 merged_df = merged_df[merged_df['significant_wave_height'] <= reasonable_max_hs] # 对于缺失值，可采用前后时刻线性插值，但需谨慎，连续长时间缺失应考虑分段处理 merged_df.interpolate(method='linear', inplace=True, limit=4) # 最多连续插值4个点
3.2 风浪与涌浪的能量分割
这是最具技术含量的一步。这里介绍一种基于风速和波龄关系的常用参数化分割方法。
核心原理：波龄是波速与风速的比值。风浪处于成长或充分成长状态，其波速与风速有较强的相关性，波龄通常在一个特定范围内（例如，充分成长风浪的波龄约0.83）。而涌浪的波速远大于当地风速，波龄很大。我们可以设定一个波龄阈值来区分。
实操步骤：
计算波速与波龄：首先需要由谱峰周期估算谱峰频率fp，然后利用线性波理论计算谱峰对应的波速Cp。
Cp = g / (2 * π * fp)，其中g为重力加速度。 波龄β = Cp / U10，其中U10为海面10米高处的风速。设定分割阈值：根据大量观测研究和海域特性，确定一个临界波龄β_crit。常见取值范围在1.2到1.5之间。当β < β_crit时，认为波浪以风浪为主；当β > β_crit时，认为涌浪能量占主导。更精细的方法会计算风浪和涌浪的有效波高：
H_s,wind = H_s * F_wind和H_s,swell = H_s * (1 - F_wind)。 其中，风浪能量比例F_wind可以通过经验公式计算，例如一个常用的公式是：F_wind = 1 / (1 + (β / β_crit)^n )，参数n通常取4或5。代码实现示例：
import numpy as np # 假设merged_df中已有‘peak_period’(Tp, 秒)和‘wind_speed’(U10, m/s)列 g = 9.81 merged_df['peak_freq'] = 1.0 / merged_df['peak_period'] # 谱峰频率 merged_df['peak_phase_speed'] = g / (2 * np.pi * merged_df['peak_freq']) # 谱峰波速 merged_df['wave_age'] = merged_df['peak_phase_speed'] / merged_df['wind_speed'] # 定义分割参数 beta_crit = 1.3 n = 4 # 计算风浪能量比例 merged_df['wind_energy_fraction'] = 1.0 / (1.0 + (merged_df['wave_age'] / beta_crit) ** n) # 计算风浪和涌浪的有效波高 (能量与波高平方成正比) merged_df['Hs_wind'] = merged_df['significant_wave_height'] * np.sqrt(merged_df['wind_energy_fraction']) merged_df['Hs_swell'] = merged_df['significant_wave_height'] * np.sqrt(1 - merged_df['wind_energy_fraction'])
实操心得：波龄阈值β_crit不是金科玉律，需要根据具体海域进行率定。最好能有一批同时包含分离风浪、涌浪波高的人工分析数据作为基准，来反推和验证最适合本海域的β_crit和n值。在没有这类数据时，可以参考邻近海域的文献值，并通过对比模拟结果与独立观测数据来间接评估分割效果。
4. 二维波浪谱的构建与参数化
得到每小时的风浪和涌浪波高后，我们需要为它们分别构建二维谱。SWAN模型通常接受两种形式的边界谱：参数化谱（如JONSWAP）或离散化的谱表。这里以更灵活的参数化谱为例。
4.1 风浪谱模型（JONSWAP谱）的参数确定
对于风浪部分，我们采用JONSWAP谱，其形式为：S(f) = α * g^2 * (2π)^{-4} * f^{-5} * exp[-1.25*(fp/f)^4] * γ^{exp[-(f-fp)^2/(2*σ^2*fp^2)]}其中，α是能量尺度参数，fp是谱峰频率，γ是峰升因子，σ是峰形参数。
我们的任务是从数据中确定这些参数：
谱峰频率fp_wind：不能直接使用混合波的Tp。一个实用的方法是，假设风浪的谱峰周期与风速相关。可采用经验公式估算，例如
Tp_wind ≈ 0.729 * U10（适用于充分成长风浪）。更准确的做法是，利用分割后的风浪波高Hs_wind和风速U10，通过风浪成长关系（如SMB公式或其变体）反推风浪的谱峰周期或频率。峰升因子γ：描述谱的尖锐程度。在北海JONSWAP实验中，平均值为3.3。在实际应用中，可根据风速或波龄给出范围，通常介于1到7之间。对于成长中的风浪，γ值较大；对于充分成长的风浪，γ值接近3.3。可以建立一个查找表或简单关系式，例如
γ = 3.3 + 0.1*(U10/Cp - 0.83)，并限制在合理范围内。能量尺度参数α：这个参数决定了谱的总能量（从而决定Hs）。它可以通过将JONSWAP谱的零阶矩m0与风浪波高关联来反算：
Hs_wind = 4.0 * sqrt(m0)。因此，我们可以先预设fp和γ，然后数值求解α，使得由该谱计算出的Hs等于我们分割得到的Hs_wind。这是一个简单的单变量求根问题。
4.2 涌浪谱模型的简化处理
涌浪谱的形态通常比风浪谱更窄、更规则。在没有详细的涌浪频谱观测时，常采用修改的PM谱或简单的单峰谱形。
谱峰频率fp_swell：可以直接使用原始谱峰频率，或者根据涌浪波高Hs_swell和假设的涌浪陡度进行估算。更可靠的方法是，如果有原始的频谱数据，在对频谱进行分割后，直接读取涌浪部分的峰值频率。
谱形：可以采用JONSWAP谱但设置较低的γ值（如1~2），或者采用形式更简单的谱，如
S(f) = A * f^{-5} * exp[-B * f^{-4}]，通过调整A、B两个参数来匹配Hs_swell和fp_swell。
4.3 方向分布函数的确定
无论是风浪谱还是涌浪谱，都需要与方向分布函数D(θ)结合，形成二维谱S(f, θ) = S(f) * D(θ)。
方向分布函数通常采用余弦幂次模型：D(θ) = G * cos^s((θ-θ_m)/2)， 当 |θ-θ_m| <= π；否则为0。 其中，θ_m是主波向，s是控制分布宽度的参数，G是归一化系数使得∫D(θ)dθ = 1。
主波向θ_m：
- 对于风浪，主波向通常与风向高度一致，但存在一个偏转角（因科氏力等因素），通常可取风向偏右（北半球）10-30度。可以直接使用观测的平均波向，但需注意这是混合波向。更严谨的做法是，使用分割后的数据或基于风浪成长模型估算。
- 对于涌浪，主波向应使用观测的波向，并假设在短时间内（如几小时）变化缓慢。
方向分布宽度参数s：这个参数最难确定。s值越大，方向分布越集中。
- 风浪：s值与风区、风速有关。成长中的风浪方向较集中（s大，如10-20），充分成长后变宽（s小，如2-5）。可参考WAM模型等第三代海浪模型中的经验关系。
- 涌浪：涌浪在传播中由于弥散和散射，方向性会变得更集中，s值通常比风浪大，可能达到20以上。
- 实操技巧：在没有明确观测数据时，这是一个重要的调参项。一个稳妥的初始设置是：风浪s=5-10，涌浪s=15-25。然后通过对比SWAN模拟结果与站点观测的方向分布（如果有）来进行校准。
5. SWAN边界条件文件生成与集成
构建好每个时刻、每种波成分的二维谱参数后，我们需要将其转化为SWAN能识别的边界条件文件。SWAN支持多种边界条件输入方式，对于非均匀、随时间变化的边界，最常用的是BOUNDSPEC命令配合谱数据文件。
5.1 生成谱数据文件
SWAN期望的谱数据文件格式通常是表格形式的。对于每个边界点（本例假设为单点），每个输出时刻，需要提供谱的离散化数据。但更常用的方式是提供参数化谱的参数，让SWAN内部生成谱。
我们可以为每个时刻生成一个包含边界条件的指令块。例如，在SWAN的输入文件（.swn）中：
BOUND SHAPESPEC JONSWAP 3.3 PEAK DSPR POWER ... BOUNDSPEC SEGMENT IJ ... & CONSTANT PAR Hs_wind Per_wind Dir_wind s_wind & CONSTANT PAR Hs_swell Per_swell Dir_swell s_swell然而，对于复杂的时变边界，更高效的方式是使用BOUNDSPEC命令读取外部的“谱参数文件”。
我们需要创建一个文本文件（如boundary_conditions.bnd），其格式大致如下：
TIME 时间点1 QUANTITY 参数1 参数2 ... ! 对应第一个边界谱 参数1‘ 参数2’ ... ! 对应第二个边界谱（如涌浪） 时间点2 QUANTITY ...具体参数取决于选择的谱形。例如，对于JONSWAP谱，可能需要提供：Hs, Tp/ fp, 主波向，γ，s (方向分布宽度)，有时还有σ_a和σ_b。
自动化脚本示例： 我们可以用Python将之前计算好的时间序列参数写入特定格式的文件。
# 假设我们有一个DataFrame `boundary_params`，索引是时间，列包括： # 'Hs_wind', 'Tp_wind', 'dir_wind', 'gamma_wind', 's_wind' # 'Hs_swell', 'Tp_swell', 'dir_swell', 's_swell' (涌浪假设gamma=1.0) with open('swan_boundary.bnd', 'w') as f: f.write('TIME\n') for time, row in boundary_params.iterrows(): # 转换时间为SWAN接受的格式，例如从1970年1月1日开始的秒数 swan_time = (time - pd.Timestamp('1970-01-01')).total_seconds() f.write(f' {swan_time:.0f}\n') f.write(' QUANTITY\n') # 写入风浪谱参数：假设顺序为 Hs, Per, Dir, Gamma, s f.write(f' {row["Hs_wind"]:.3f} {row["Tp_wind"]:.2f} {row["dir_wind"]:.1f} {row["gamma_wind"]:.2f} {row["s_wind"]:.1f}\n') # 写入涌浪谱参数：假设涌浪用JONSWAP但gamma固定为1.0 f.write(f' {row["Hs_swell"]:.3f} {row["Tp_swell"]:.2f} {row["dir_swell"]:.1f} 1.0 {row["s_swell"]:.1f}\n')5.2 在SWAN输入文件中配置
在SWAN的.swn主输入文件中，需要正确引用这个边界文件并设置边界线段。
PROJECT 'My_Wave_Simulation' ... BOUND SHAPESPEC JONSWAP 3.3 PEAK DSPR POWER BOUNDSPEC SEGMENT IJ 1 1 100 1 VARIABLE FILE 0 'swan_boundary.bnd'这里SEGMENT IJ定义了边界线段的位置（网格索引）。VARIABLE FILE 0表示从文件swan_boundary.bnd中读取随时间变化的边界条件，文件中的每个“QUANTITY”块对应两个谱（因为我们在每个时间点写了两行参数）。
5.3 关键配置解析与经验
边界嵌套：如果计算域很大，通常采用嵌套网格。大区域模型（如WW3）的输出可以作为本区域SWAN模型的边界。此时，本方法生成的边界条件可能用于更内层的小尺度精细网格，或者用于驱动一个单站点的谱模型来验证边界条件。
谱的分辨率：在
BOUNDSPEC或GENPAR命令中，需要设置边界谱的频率和方向分辩率（FRANGE和DIRANGE）。这应与计算域内的谱离散设置协调。通常边界谱的分辨率可以略低于计算域内部，以节省I/O，但不能太低以致丢失关键谱峰信息。时间插值：SWAN在读取时间序列边界文件时，会在给定的时间点之间进行线性插值。因此，输入文件的时间间隔需要足够小，以捕捉风浪和涌浪的快速变化，特别是在风暴过境期间。通常1小时间隔是常见的，对于快速变化的台风，可能需要缩短到15-30分钟。
踩坑记录：最大的一个坑是单位一致性。SWAN默认使用国际单位制（米、秒、度），但风向波向的定义（相对于正北顺时针还是逆时针？数学角度还是地理角度？）必须仔细核对。我遇到过因为波向定义搞反（数学角度0度是向东，地理角度0度是向北），导致波浪完全从错误的方向传入，模拟结果完全错误。务必查阅所用SWAN版本的说明书，并在生成边界文件时进行正确的转换。一个检查方法是，用一个小测试用例，输入一个简单的单向谱，看模拟的波浪传播方向是否符合预期。
6. 验证、常见问题与调参策略
生成边界条件并运行SWAN后，工作只完成了一半。验证和调参是确保模拟结果可信的关键。
6.1 验证方法
单点验证：在计算域内，选取与观测站点位置对应的网格点，输出该点的波浪参数（Hs, Tp, 波向）时间序列，与实测数据进行对比。绘制时间序列对比图，计算偏差统计量（如偏差Bias、均方根误差RMSE、相关系数R）。
- 重点看什么：是否捕捉到了主要的波动过程？风浪和涌浪的峰值是否对齐？波向的变化趋势是否一致？如果分割正确，模型应能更好地再现混合浪中不同成分的演变。
谱验证：如果有可能，输出该点的二维谱或一维频率谱、方向谱，与观测谱进行对比。这是最严格的检验，能直接看出构建的边界谱是否真实反映了能量的频率-方向分布。
6.2 常见问题排查表
| 问题现象 | 可能原因 | 排查与解决思路 |
|---|---|---|
| 模拟的Hs整体偏大/偏小 | 1. 边界谱总能量（Hs）设置错误。 2. 风场输入同时存在，能量双重输入。 | 1. 检查分割计算和谱参数反推过程，确认Hs计算无误。 2. 如果SWAN同时使用了风场驱动，检查 WIND命令和BOUNDSPEC命令是否造成了能量重复。通常，使用详细边界谱时，计算域内不应用风场再生成风浪，除非是模拟边界内的风区。 |
| 模拟的波向与观测持续存在固定偏差 | 1. 输入的主波向定义错误（地理角/数学角）。 2. 方向分布函数过于集中或分散。 | 1.彻底检查角度转换。SWAN通常使用地理角度（顺时针从正北起算）。确认观测数据、内部计算和边界文件的角度约定。 2. 调整方向分布宽度参数 s。偏角系统偏差可能需修正风向-波向偏转角参数。 |
| 无法模拟出涌浪事件 | 1. 波浪分割失败，涌浪能量被归入风浪。 2. 涌浪谱参数（如Tp）设置不合理，能量在传播中耗散过快。 3. 边界线段方向设置错误，涌浪无法传入。 | 1. 回顾分割算法，检查波龄阈值β_crit是否过大。对比分割前后的Hs_wind和Hs_swell时间序列，看涌浪事件是否被识别出来。 2. 检查涌浪的底摩擦、白冠耗散等源项设置是否过于激进。可尝试在模型中对涌浪成分调低耗散系数。 3. 检查 BOUNDSPEC SEGMENT的定位和方向，确保其法线方向与涌浪来向大致垂直。 |
| 模拟谱形与观测谱形差异大 | 1. 选择的参数化谱形（JONSWAP）不适合该海域或该波成分。 2. 峰升因子γ、方向分布宽度s等形状参数取值不当。 | 1. 尝试其他谱形，如TMA谱（适用于有限水深）。对于涌浪，可尝试双峰谱或直接输入离散谱数据。 2. 进行参数敏感性分析。固定Hs和Tp，微调γ和s，观察对模拟结果的影响，寻找最优组合。 |
| 模型运行不稳定或报错 | 1. 边界谱参数存在非物理值（如负的Hs，极端的Tp）。 2. 边界文件格式错误，时间顺序错乱。 | 1. 在生成边界文件的脚本中加入严格的合理性检查，过滤掉非物理值。 2. 仔细对照SWAN手册检查边界文件格式，确保时间严格递增，每个时间块内的数据行数正确。用文本编辑器查看文件开头和结尾部分。 |
6.3 系统性调参策略
当模拟结果与观测存在偏差时，应遵循由外及内、由主到次的顺序进行调参：
- 第一步：校准边界总能量。确保输入的Hs（风浪+涌浪）时间序列与观测的总体Hs在量级和变化趋势上基本吻合。这是最重要的基础。
- 第二步：校准谱峰周期。调整风浪和涌浪的Tp参数化方案，使模拟的Tp与观测的Tp（或分割后的Tp）匹配。
- 第三步：校准波向。修正系统性的波向偏差，并微调方向分布宽度参数s，以改善波向的分布和集中程度。
- 第四步：校准谱形。最后调整峰升因子γ等精细参数，以改善模拟的谱峰尖锐度、谱宽等特征。
整个调参过程应是一个迭代的过程，并且最好使用独立的验证数据集（即未参与边界条件生成的数据时段）来评估最终模型的性能，以避免过拟合。
7. 从脚本到流程：构建自动化工具链
手动处理数据、分割、计算参数、写边界文件是非常繁琐且易出错的。在实际项目中，我将上述所有步骤整合成了一个半自动化的Python工具链，核心模块包括：
- 数据预处理模块：统一读取多种格式的原始数据，进行质量控制和插补。
- 波浪分割引擎：实现多种分割算法（参数化法、频谱分割法），并提供可视化界面用于对比和选择分割结果。
- 谱参数计算模块：根据选定的谱模型，从分割后的波高、周期、风速等数据反算谱参数（α, fp, γ, s等）。
- SWAN边界文件生成器：将时间序列的谱参数按照指定格式写入SWAN边界文件，并自动生成对应的SWAN命令片段。
- 验证与后处理模块：自动运行SWAN（可调用其命令行接口），提取指定点的模拟结果，与观测数据对比，生成标准化的验证图表和报告。
这个工具链不仅将过去需要数天的工作缩短到几小时，更重要的是保证了处理过程的一致性和可重复性，减少了人为错误。它允许我快速进行敏感性试验，比如：“如果我把波龄阈值从1.3改成1.4，对最终的风电场疲劳载荷估算会有多大影响？”
最后，分享一个深刻的体会：基于数据分析的风浪边界条件生成，其精度上限取决于观测数据的质量和代表性。再精巧的算法也无法弥补糟糕数据带来的问题。因此，在项目初期，投入精力进行数据评估、了解观测仪器的局限性和数据的时空代表性，与后续的模型构建同等重要。这套方法的价值，在于它为我们提供了一套严谨的框架，将有限的数据信息最大化地利用起来，让数值模拟的边界不再是模糊的“灰色地带”，而是建立在坚实数据分析基础上的“透明窗口”。
本文还有配套的精品资源，点击获取
