1. 赛题核心从“云中的海盐”到数学模型的构建最近在准备“认证杯”数学建模网络挑战赛第二阶段C题“云中的海盐”这个题目挺有意思它把大气科学、环境物理和数学建模结合在了一起。题目背景大概是研究海盐气溶胶简单理解就是海水蒸发后被风带到空中的微小盐粒如何通过云的形成和降水过程最终被清除出大气。这个过程听起来很自然但要用数学模型去量化它里面涉及到的变量和机制就复杂了。很多同学拿到这种题第一反应可能是去搜现成的代码或者论文但更关键的是先理解题目到底在问什么以及背后的物理图景是什么。这直接决定了你模型的方向和精度。“云中的海盐”本质上是一个大气颗粒物气溶胶的源-输送-转化-沉降过程模拟问题。海盐是源大气运动和云物理过程是输送与转化的载体湿沉降随降水落下和干沉降直接沉降是最终的汇。题目通常会提供或暗示一些观测数据比如不同高度、不同时间的海盐浓度、气象数据等要求我们建立一个模型来描述海盐浓度随时间、空间的变化并可能预测其在特定条件下的归宿。所以我们的核心任务不是去复现一个复杂的气候模式而是抓住主要矛盾用相对简洁但物理意义明确的数学模型来刻画这个过程。理解这一点后我们的建模思路就不能停留在简单的数据拟合上。我们需要考虑几个核心环节海盐的排放源强如何计算、在大气中的扩散用什么方程描述、云内过程云如何捕获这些颗粒云水含量、降水强度如何影响、以及最终的沉降。每一个环节都可以对应一个或一组数学方程。接下来的内容我会结合常见的建模方法和Matlab工具一步步拆解如何构建这个模型并分享一些在数据处理和数值求解中容易踩的坑。2. 模型框架选择与关键物理过程数学描述面对这样一个多过程耦合的问题直接上手编代码很容易迷失。我建议先搭建一个清晰的模型框架。对于时间尺度在几天到一周、空间尺度在区域范围内的海盐传输问题一个基于欧拉框架的箱模型或一维/二维平流-扩散方程是比较务实且有效的起点。我们不需要像WRF-Chem那样做完整的三维模拟但模型必须包含核心物理过程。2.1 控制方程平流-扩散-沉降方程这是描述大气中污染物浓度演化的基石。对于海盐气溶胶的质量浓度 ( C(x, y, z, t) )单位通常是 (\mu g/m^3)其控制方程可以写为[\frac{\partial C}{\partial t} -\nabla \cdot (\mathbf{u} C) \nabla \cdot (K \nabla C) S - R]这个方程看起来复杂我们拆开看(\frac{\partial C}{\partial t})浓度随时间的变化率。这就是我们要求解的东西。(-\nabla \cdot (\mathbf{u} C))平流项。表示风风速矢量 (\mathbf{u})把海盐从一个地方带到另一个地方。这是浓度空间分布变化的主要驱动力。在Matlab中处理这项时要特别注意数值格式的选择比如迎风格式否则容易产生数值震荡非物理的浓度波动。(\nabla \cdot (K \nabla C))扩散项。表示由于大气湍流造成的海盐从高浓度区向低浓度区的扩散。(K) 是湍流扩散系数它通常不是常数而是随高度、大气稳定度变化。简化模型中我们有时会把它设为一个经验常数。(S)源项。这就是海盐从海面的排放。源强 (S) 的计算本身就是一个子模型通常与风速特别是10米高度风速 (U_{10})的某次方成正比因为风越大海浪破碎产生海盐气泡越多。一个经典的公式是 (S A \cdot U_{10}^B)其中 (A) 和 (B) 是经验参数需要根据题目给出的数据或参考文献进行率定。(R)清除项。这是本题的关键主要指湿沉降。干沉降直接撞到地面或海面在短时间模拟中有时可以忽略或者合并为一个简单的沉降速度参数。湿沉降 (R_{wet}) 通常表示为(R_{wet} \Lambda \cdot C \cdot q)。这里 (\Lambda) 是清除系数或冲刷比(q) 是降水强度单位时间内的降水量。这个公式的物理意义是降水强度越大空气中能被雨滴捕获带走的颗粒物比例就越高。在编程实现时我们往往需要对上述方程进行离散化。例如如果我们简化为一维垂直柱模型只考虑高度z方向的变化方程可以简化为 [\frac{\partial C}{\partial t} -w \frac{\partial C}{\partial z} \frac{\partial}{\partial z}(K_z \frac{\partial C}{\partial z}) S(z) - \Lambda q C] 这里 (w) 是垂直速度通常很小可假设为0(K_z) 是垂直扩散系数。这个一维模型足以研究海盐从海面排放向上输送然后被云和降水清除的垂直廓线演变对于回答一些宏观问题已经足够。2.2 云过程的参数化连接气溶胶与降水的桥梁题目叫“云中的海盐”云的作用至关重要。但直接模拟云微物理云滴凝结、碰并等对于数模竞赛来说过于沉重。因此我们需要参数化。也就是说用一些相对简单的公式来表征云对海盐的净影响。一个常见的思路是引入“云水含量” (L)单位g/m³和“降水效率” (E) 的概念。我们可以认为当相对湿度达到一定阈值如100%时水汽凝结形成云海盐颗粒作为凝结核被包裹进云滴。此时部分海盐质量会从“气溶胶相”转移到“云水相”。随后部分云水通过自动转化等过程形成降水落下。这个过程可以粗略地用两个方程描述云激活/捕获假设在云区内由相对湿度或液态水含量判断气溶胶浓度 (C) 以一定速率 (k_{in}) 转化为云中物质或认为云水中的海盐浓度与气溶胶浓度成正比。湿沉降清除云中物质携带海盐以降水强度 (q) 和清除效率 (\eta) 被移除。即 (R_{wet} \eta \cdot q \cdot C_{cloud})其中 (C_{cloud}) 是云中海盐相关物质的浓度。在实际编程中我们往往将这两个过程合并直接使用前面提到的 (R_{wet} \Lambda \cdot C \cdot q) 公式。这里的 (\Lambda) 就是一个综合了云捕获效率和降水效率的参数。它的取值需要参考文献典型值可能在 (10^{-5}) 到 (10^{-4}) /s/mm/h量级。这里一个重要的经验是(\Lambda) 的值对模拟结果极其敏感。在参数率定或敏感性分析部分必须测试 (\Lambda) 在不同量级下对海盐寿命和垂直分布的影响。3. 基于Matlab的数值求解策略与代码实现要点有了数学模型接下来就是用Matlab把它实现出来。这里我们以一维垂直模型为例展示核心的求解流程和代码片段。选择一维模型是因为它概念清晰计算量小适合在竞赛有限时间内调试并得到有物理意义的结果同时其数值方法可以很容易推广到二维。3.1 空间离散与时间推进有限差分法我们采用有限差分法来离散化控制方程。将垂直高度从海面z0到模型顶层zH划分为N层每层厚度为 (\Delta z)。浓度 (C) 在网格点 (j) 上取值 (C_j)代表第j层中心的浓度。以简化的垂直扩散-沉降方程为例假设无风有恒定源和湿沉降 [\frac{\partial C}{\partial t} K_z \frac{\partial^2 C}{\partial z^2} S - \Lambda q C]其显式差分格式为 [\frac{C_j^{n1} - C_j^n}{\Delta t} K_z \frac{C_{j1}^n - 2C_j^n C_{j-1}^n}{(\Delta z)^2} S_j - \Lambda q C_j^n]这里上标 (n) 代表时间步。整理后得到时间推进公式 [C_j^{n1} C_j^n \Delta t \left[ K_z \frac{C_{j1}^n - 2C_j^n C_{j-1}^n}{(\Delta z)^2} S_j - \Lambda q C_j^n \right]]注意显式格式有条件稳定。稳定性要求 (\Delta t \leq \frac{(\Delta z)^2}{2K_z})。如果扩散系数 (K_z) 很大或网格很细时间步长必须非常小否则计算会发散。这是第一个大坑。% 参数设置 H 2000; % 模型顶高度 (m) N 100; % 网格数 dz H/N; % 网格间距 dt 10; % 时间步长 (s)需满足稳定性条件 total_time 24*3600; % 模拟总时间 24小时 nt round(total_time / dt); Kz 10.0; % 垂直扩散系数 (m^2/s)这是一个简化假设实际随高度变化 Lambda 5e-5; % 湿沉降清除系数 (1/s) q 2.0 / 3600 / 1000; % 假设降水强度 2 mm/h转换为 m/s % 初始化 z linspace(dz/2, H-dz/2, N); % 网格中心高度 C zeros(N, 1); % 初始浓度为零 S zeros(N, 1); S(1:10) 1e-6; % 假设源集中在近地面10层单位: kg/m^3/s (示例值) % 时间循环 for it 1:nt C_new C; % 内部网格点计算扩散项使用循环清晰但较慢实际可用矩阵运算优化 for j 2:N-1 diff_term Kz * (C(j1) - 2*C(j) C(j-1)) / (dz^2); C_new(j) C(j) dt * (diff_term S(j) - Lambda * q * C(j)); end % 边界条件处理顶部通量为零底部为海面可设为定浓度或通量边界 j 1; % 底层 % 假设底层浓度梯度由海面通量平衡简化处理下边界采用虚拟点法或直接设定 % 这里简化底层扩散项采用单边差分并加上源 diff_term_bottom Kz * (C(2) - C(1)) / (dz^2); % 近似处理 C_new(1) C(1) dt * (diff_term_bottom S(1) - Lambda * q * C(1)); j N; % 顶层 % 顶层通量为零边界条件dC/dz 0意味着 C(N1) C(N-1) 在虚拟点 % 因此扩散项为Kz * (C(N-1) - 2*C(N) C(N-1)) / dz^2 2*Kz*(C(N-1)-C(N))/dz^2 diff_term_top 2 * Kz * (C(N-1) - C(N)) / (dz^2); C_new(N) C(N) dt * (diff_term_top S(N) - Lambda * q * C(N)); C C_new; % 更新浓度场 % 可选每模拟一段时间输出一次结果 if mod(it, round(3600/dt)) 0 fprintf(模拟时间: %d 小时\n, it*dt/3600); end end % 绘制最终浓度垂直廓线 figure; plot(C, z/1000, b-o, LineWidth, 1.5); xlabel(海盐浓度 (kg/m^3)); ylabel(高度 (km)); title(24小时模拟后海盐浓度垂直分布); grid on;这段代码提供了一个最基础的骨架。在实际竞赛中你需要根据题目给出的具体条件进行大幅修改和增强比如源项S需要根据题目给出的风速公式计算随时间或空间变化的源。扩散系数Kz很少是常数。一个常见的参数化是随高度增加而增大近地面湍流强达到最大值后再减小。可以参考“K理论”或使用边界层高度相关的经验公式。清除项降水强度 (q) 很可能不是常数而是随时间变化的例如有降水事件。你需要根据题目提供的气象数据如果有来定义 (q(t))。清除系数 (\Lambda) 也可能与降水类型对流性降水/层状云降水有关。边界条件底边界海面的处理非常关键。通常有两种一是给定浓度Dirichlet条件如海面浓度恒定二是给定通量Neumann条件即排放通量已知。需要根据题目表述选择。顶边界一般设为通量为零。3.2 提升稳定性与效率隐式格式与矩阵运算显式格式的稳定性限制是个麻烦。对于扩散问题采用隐式格式如Crank-Nicolson格式可以无条件稳定允许使用更大的时间步长 (\Delta t)虽然每步计算量稍大但总体效率可能更高。隐式格式将方程写为 [\frac{C_j^{n1} - C_j^n}{\Delta t} \theta \cdot [K_z \frac{\partial^2 C}{\partial z^2} S - \Lambda q C]^{n1} (1-\theta) \cdot [K_z \frac{\partial^2 C}{\partial z^2} S - \Lambda q C]^{n}] 当 (\theta 0.5) 时即为Crank-Nicolson格式。这将导致一个线性方程组 [A \mathbf{C}^{n1} B \mathbf{C}^{n} \mathbf{b}] 其中 (A) 和 (B) 是三对角矩阵因为只涉及相邻网格点(\mathbf{b}) 是源项相关的向量。在Matlab中可以使用稀疏矩阵和反斜杠运算符高效求解。% 使用Crank-Nicolson格式θ0.5求解一维扩散-沉降方程 % 假设S和q为常数仅演示格式 % 构造系数矩阵A和B忽略源项和沉降项中的隐式部分简化完整形式更复杂 r Kz * dt / (2 * dz^2); main_diag_A (1 2*r 0.5*dt*Lambda*q) * ones(N,1); % 主对角线元素 sub_diag -r * ones(N-1, 1); % 次对角线元素 % 构建三对角矩阵A (I - θ*L)其中L是扩散算子矩阵 A spdiags([sub_diag, main_diag_A, sub_diag], [-1, 0, 1], N, N); % 处理边界条件以顶层通量零为例会影响A矩阵的最后一个元素 A(N, N-1) -2*r; % 顶层通量零边界修正 A(N, N) 1 2*r 0.5*dt*Lambda*q; % 主对角线不变 % 矩阵B (I (1-θ)*L) main_diag_B (1 - 2*r - 0.5*dt*Lambda*q) * ones(N,1); B spdiags([-sub_diag, main_diag_B, -sub_diag], [-1, 0, 1], N, N); B(N, N-1) 2*r; % 顶层边界修正 B(N, N) 1 - 2*r - 0.5*dt*Lambda*q; % 时间推进 for it 1:nt b dt * S; % 源项贡献这里S是向量 C_new A \ (B * C b); % 求解线性方程组 C C_new; end这里的关键心得是在竞赛有限时间内如果模型维度不高如一维显式格式编程简单易于调试只要时间步长设得足够小是可以接受的。但如果模型复杂或需要长时间模拟花点时间实现隐式格式是值得的它能避免很多因稳定性问题导致的诡异结果。另外一定要善用Matlab的稀疏矩阵spdiags功能来存储和运算A、B矩阵这对于二维甚至三维问题能极大节省内存和计算时间。4. 数据驱动、参数率定与模型验证的实战思路数学模型建好了代码也跑起来了但你怎么知道你的模型是对的呢在数学建模竞赛中模型验证和参数率定是区分优秀论文和普通论文的关键。对于“云中的海盐”这类问题通常题目会提供一部分数据可能是某个站点的观测浓度随时间变化或者不同高度的平均浓度廓线。我们的目标就是让模型输出尽可能贴近这些观测数据。4.1 参数敏感性分析与率定模型中有很多参数是不确定的比如源强公式中的系数 (A) 和 (B)扩散系数 (K_z)湿沉降清除系数 (\Lambda) 等。我们不能随意给它们赋值。一个系统的方法是敏感性分析在合理范围内变动某个参数其他参数固定观察模型输出如地面浓度、柱总量、垂直分布形状的变化程度。这能告诉我们哪个参数对结果影响最大需要重点率定。在Matlab中可以写一个循环来实现。Lambda_range logspace(-6, -4, 20); % 测试Λ从1e-6到1e-4 ground_concentration zeros(size(Lambda_range)); for i 1:length(Lambda_range) Lambda_test Lambda_range(i); % 运行一次模拟使用Lambda_test % ... [模拟代码记录模拟结束时的地面浓度C(1)] ground_concentration(i) C(1); end figure; loglog(Lambda_range, ground_concentration, s-); xlabel(\Lambda (湿沉降清除系数)); ylabel(模拟地面浓度); title(地面浓度对Λ的敏感性); grid on;如果曲线很陡说明敏感性高如果平缓说明模型对该参数不敏感甚至可以取一个文献中的典型值。参数率定找到一组参数使得模型输出与观测数据的误差最小。这本质上是一个优化问题。最常用的方法是最小二乘法即最小化模拟值与观测值之差的平方和。% 假设我们有观测数据 obs_data (时间序列或垂直廓线) % 和对应的模拟输出 sim_data(param) % 定义误差函数 function error myErrorFunction(param) % param是一个向量包含要率定的参数如 [A, B, Lambda] A param(1); B param(2); Lambda param(3); % 调用你的模型函数得到模拟结果 sim_result sim_result run_my_model(A, B, Lambda, ...); % 计算与观测数据 obs_data 的均方根误差 (RMSE) error sqrt(mean((sim_result - obs_data).^2)); end % 使用Matlab优化工具箱进行参数寻优 initial_guess [1.3e-5, 3.4, 5e-5]; % 初始猜测值 lb [1e-6, 2.0, 1e-6]; % 参数下界 ub [1e-4, 4.0, 1e-4]; % 参数上界 options optimoptions(fmincon, Display, iter); [best_param, min_error] fmincon(myErrorFunction, initial_guess, [], [], [], [], lb, ub, [], options);常用的优化函数有fmincon有约束优化、lsqnonlin非线性最小二乘等。这里有个大坑优化结果可能陷入局部最优。因此初始猜测值initial_guess最好基于文献或物理意义来设定并且可以尝试多组不同的初始值来验证结果的稳健性。4.2 模型验证与不确定性讨论率定好参数后需要用另一部分未参与率定的观测数据来验证模型。这叫“交叉验证”。如果模型在验证数据上也表现良好说明模型具有一定的泛化能力。在论文中这部分结果通常用图表展示时间序列对比图将模拟的海盐浓度时间序列与观测数据画在同一张图上直观对比趋势和峰值。散点图与拟合线横坐标为观测值纵坐标为模拟值。如果点都分布在1:1线附近说明模拟效果好。可以计算相关系数 (R^2)、均方根误差RMSE等定量指标。垂直廓线对比对比模拟和观测的浓度随高度变化曲线。必须讨论模型的不确定性。在结论部分要坦诚地指出模型的局限性比如我们忽略水平输送一维模型的局限、简化了复杂的云微物理过程、参数存在不确定性等。并可以做一个简单的情景分析如果未来风速增加20%或者降水强度加倍根据你的模型海盐的沉降通量会如何变化这能体现模型的应用价值。5. 从一维到二维/三维的扩展思路与竞赛实战建议虽然一维模型足以应对很多问题但如果题目明确要求考虑水平分布例如研究海盐从海岸向内陆的输送就需要扩展到二维甚至三维。思路是相通的但计算复杂度会指数级增加。5.1 二维模型构建要点对于二维x-z剖面模型控制方程变为 [\frac{\partial C}{\partial t} -u \frac{\partial C}{\partial x} \frac{\partial}{\partial z}(K_z \frac{\partial C}{\partial z}) S(x,z) - \Lambda q C] 这里增加了水平平流项 (-u \frac{\partial C}{\partial x})(u) 是水平风速可能随高度变化。数值求解时扩散项的处理和一维类似但平流项的数值处理需要格外小心。前面提到的迎风格式几乎是必须的它能保证数值稳定性。在Matlab中实现二维模型浓度场 (C) 变成一个矩阵代码中会用到双重循环或者更高效的矩阵运算。% 二维模型伪代码结构示意 nx 100; nz 50; C zeros(nz, nx); U ... % 水平风速场 (nz, nx) Kz ... % 垂直扩散系数场 S ... % 源项场在近海区域设置源 for it 1:nt C_new C; for ix 2:nx-1 for iz 2:nz-1 % 水平平流项迎风差分 if U(iz, ix) 0 adv_x U(iz, ix) * (C(iz, ix) - C(iz, ix-1)) / dx; else adv_x U(iz, ix) * (C(iz, ix1) - C(iz, ix)) / dx; end % 垂直扩散项 diff_z Kz(iz) * (C(iz1, ix) - 2*C(iz, ix) C(iz-1, ix)) / (dz^2); % 源和汇 source S(iz, ix); sink Lambda * q * C(iz, ix); C_new(iz, ix) C(iz, ix) dt * (-adv_x diff_z source - sink); end end % 处理边界条件海岸、地面、顶部、远场 C C_new; end二维模拟的输出可以用pcolor或imagesc函数绘制浓度空间分布图非常直观。5.2 给参赛者的几点核心建议先简后繁确保核心流程跑通不要一开始就追求复杂的二维、三维模型。先用一维模型把排放、扩散、沉降的整个流程用代码实现并能在简单条件下如恒定风、恒定降水跑出一个合理的结果比如浓度随高度递减有降水时浓度降低。这是你的“基本盘”。数据预处理是重中之重题目给的数据风速、降水、初始浓度往往不是直接可用的。可能需要单位换算、插值到你的模型网格上、处理缺失值。在Matlab中interp1一维插值、interp2二维插值、meshgrid等函数会非常有用。花在数据清洗上的时间往往比写模型本身还多。可视化贯穿始终每完成一个步骤就画图看看。初始化后的浓度场、风场、源项分布、每模拟一段时间的浓度廓线……图形能帮你快速发现代码错误比如出现负浓度、不合理的空间分布和物理上的不合理之处。论文写作与代码并重数学建模竞赛最终提交的是论文。在论文中你的模型思路、公式推导、参数选取理由、结果分析和图表展示比代码本身更重要。代码是工具论文才是作品。确保你的论文逻辑清晰图表美观对每一个结果都有合理解释。团队分工要明确通常三人一组建议一人主攻模型构建与公式推导理论一人主攻编程实现与调试编程一人主攻论文写作与数据可视化写作。但三者需要紧密沟通理论者要理解编程的可行性编程者要理解模型的物理意义写作者要准确表达两者的工作。最后关于“云中的海盐”这个具体赛题由于没有看到具体的题目描述和数据以上内容是一个通用的、高可行性的建模框架。你需要将上述通用原理与题目给出的具体条件、数据相结合。比如如果题目给出了具体的风速时间序列你的源项 (S) 就要随之变化如果给出了不同天气案例有云无云、有雨无雨你的清除项 (\Lambda) 和 (q) 就要分情况讨论。记住所有模型都是对现实的简化关键在于你的简化是否抓住了主要矛盾并且能用数学语言自洽地描述出来。