尧图网站建设 
1. 项目概述从“光污染”到量化分析去年带队做美赛E题“光污染”的经历现在回想起来依然觉得是一次非常典型的建模实战。题目给了一个看似宏大的社会与环境议题但内核却要求我们将其转化为一系列可量化、可计算、可预测的数学问题。很多队伍一开始就懵了要么陷在“光污染危害”的文献综述里出不来要么对着“评价指标”和“控制模型”这些词无从下手。其实这道题的核心思路非常清晰定义问题 - 建立评价体系 - 构建预测模型 - 提出干预策略。整个过程MATLAB作为核心计算与可视化工具贯穿始终。简单来说这道题要求我们做三件事第一为光污染建立一个综合性的评价指标体系不能只看亮度还得考虑光谱、时间、空间和对生态/社区的影响。第二基于这个体系建立一个能预测光污染水平如何随城市发展如人口、GDP、土地类型变化而演变的数学模型。第三利用这个模型设计一套成本效益最优的干预策略告诉决策者“在哪里、用什么方法、花多少钱”能最有效地控制光污染。这不仅仅是套个公式那么简单。你需要处理多源异构数据卫星遥感夜间灯光数据、城市统计数据、地理信息数据需要融合环境科学、经济学和公共政策的知识更需要用MATLAB将这一切想法高效地实现出来从数据清洗、空间分析、模型拟合到策略模拟。接下来我就结合我们当时的解题思路和具体的MATLAB代码实现拆解每一个关键环节希望能给未来参赛的你一些实实在在的参考。2. 解题核心思路与模型框架设计面对“光污染”这种开放性问题最忌讳的就是一上来就埋头写代码。我们花了将近半天的时间来搭建整个解题的逻辑框架确保每一步都有的放矢。我们的整体思路可以概括为“一个核心两个维度三层模型”。2.1 问题拆解与核心定义题目中的“光污染”是一个模糊概念。我们首先要将其操作化为可测量的指标。我们将其分解为两个核心维度物理强度维度即光本身的性质。包括亮度最直接的指标来自卫星夜间灯光数据如VIIRS DNB数据的辐射值。光谱构成不同波长的光对生态和人类的影响不同。例如短波蓝光对天文观测干扰大且更易影响人体节律。我们通过查阅文献为常见光源高压钠灯、LED白光、LED暖光赋予不同的光谱权重系数。时间动态通宵照明与午夜后关闭的照明其污染程度截然不同。我们引入了“夜间光照持续时间”作为时间维度指标。影响效应维度即光对受体产生的影响。包括生态影响重点考虑对夜间迁徙鸟类、海龟、昆虫的干扰。我们根据土地利用类型如湿地、森林、海岸线设定生态敏感区光辐射到这些区域则产生“生态影响分值”。社区影响考虑对人类健康睡眠干扰、交通安全眩光、夜空观赏性星空可见度的影响。这部分更多依赖于人口密度数据和问卷调查数据的模拟。2.2 综合评价指数CEI的构建将上述多维指标合成一个综合性的“光污染评价指数”Comprehensive Evaluation Index, CEI是关键一步。我们采用了加权线性综合法但权重的确定不是主观臆断。我们使用了层次分析法AHP结合熵权法来确定权重。AHP用于整合我们通过文献和专家咨询模拟得到的各指标相对重要性判断熵权法则基于实际数据如各区域亮度值的差异度计算客观权重。最后将主客观权重结合得到最终权重。在MATLAB中AHP的实现需要构建判断矩阵并进行一致性检验。% 示例AHP权重计算核心片段 function weights calculateAHP(matrix) % matrix: 判断矩阵 [n, ~] size(matrix); [V, D] eig(matrix); % 求特征值和特征向量 [~, idx] max(diag(D)); % 找到最大特征值 w V(:, idx); % 对应的特征向量 weights w / sum(w); % 归一化得到权重 % 一致性检验 CI, RI, CR lambda_max max(diag(D)); CI (lambda_max - n) / (n - 1); RI [0, 0, 0.58, 0.90, 1.12, 1.24, 1.32, 1.41, 1.45]; % 平均随机一致性指标 CR CI / RI(n); if CR 0.1 error(判断矩阵一致性不满足要求请调整); end end2.3 预测模型的选择系统动力学SD与多元回归的融合题目要求预测未来光污染趋势。我们放弃了单一的回归模型因为光污染系统内部存在复杂的反馈关系例如经济发展推动照明需求光污染加剧可能引发政策调控。因此我们引入了系统动力学System Dynamics, SD模型作为核心预测框架。我们将系统划分为几个关键存量光污染综合指数CEI社会经济驱动因子如人均GDP、城市化率政策干预强度然后定义流量和反馈回路。例如“经济发展”正反馈于“照明基础设施投资”从而增加“光辐射源强度”而“公众环保意识”和“已实施政策”会负反馈于“新照明项目的审批尺度”。实操心得直接用Vensim等SD专用软件建模当然可以但在美赛有限时间内为了与MATLAB数据处理和后续优化无缝衔接我们在MATLAB中用差分方程组来近似实现这个SD模型的核心逻辑。这要求我们对模型进行合理的简化抓住主要矛盾。同时SD模型中的一些参数关系如GDP增长率对灯光亮度增长率的影响系数需要通过历史数据进行标定。这里我们用到了多元非线性回归。MATLAB的fitnlm函数在这里大显身手。% 示例使用多元非线性回归拟合SD模型参数 % 假设我们有历史数据年份Year GDP 灯光亮度Brightness % 建立关系Brightness a * GDP^b c % 准备数据 tbl table(GDP, Brightness, VariableNames, {GDP,Brightness}); % 定义非线性模型公式 modelfun (b,x) b(1) * x.GDP.^b(2) b(3); % 初始参数猜测 beta0 [0.01, 0.5, 10]; % 拟合模型 nlmodel fitnlm(tbl, modelfun, beta0); % 查看拟合结果和参数 disp(nlmodel.Coefficients.Estimate); plot(nlmodel); % 可视化拟合效果这个融合框架的优势在于既有系统思考的宏观视角SD又有数据驱动的微观精度回归使得我们的预测模型既有说服力又具备一定的稳健性。3. 数据获取、处理与空间分析实战模型框架搭好了接下来就是喂数据。数据决定了模型的上限。E题通常不提供现成数据需要自己寻找。我们主要使用了以下几类数据并全部在MATLAB中完成处理。3.1 核心数据源夜间灯光数据NASA的VIIRS DNB可见光红外成像辐射计套件-日/夜波段月度数据是研究光污染的黄金标准。可以从NOAA官网下载https://eogdata.mines.edu/products/vnl/。我们下载了目标区域我们自选了一个大都市圈最近5年的月度平均辐射合成数据.tif格式。3.2 数据处理全流程数据读取与裁剪使用geotiffread读取TIFF文件获取图像数据和地理参考信息。然后用目标区域的矢量边界Shapefile通过mapcrop或地理坐标匹配进行裁剪。[DNB_data, R] readgeoraster(viirs_202201.tif); % 假设有边界坐标向量latlim, lonlim [DNB_cropped, R_cropped] mapcrop(DNB_data, R, latlim, lonlim);异常值与背景值处理夜间灯光数据包含背景噪声气辉、月光等和异常高值火光、油气燃烧。我们采用分位数阈值法进行过滤。lower_prc prctile(DNB_cropped(DNB_cropped0), 1); % 取正值数据的1%分位数 upper_prc prctile(DNB_cropped(:), 99.9); % 取全体数据的99.9%分位数 DNB_cleaned DNB_cropped; DNB_cleaned(DNB_cleaned lower_prc) 0; % 低于下限的视为背景噪声 DNB_cleaned(DNB_cleaned upper_prc) upper_prc; % 高于上限的进行截断避免极端值影响时间序列整合与趋势分析将月度数据按年度平均得到年际序列。使用detrend函数或线性拟合来观察亮度总体增长趋势。yearly_mean zeros(5,1); for i 1:5 % 读取第i年的12个月数据求平均 yearly_mean(i) mean(annual_data_i(:)); end % 线性趋势拟合 p polyfit((1:5), yearly_mean, 1); trend polyval(p, (1:5));空间网格化与指标计算将研究区域划分为1km x 1km的网格。对每个网格计算其平均亮度、亮度年增长率、以及其所属的土地利用类型来自全球土地覆盖数据如FROM-GLC。% 假设已将区域网格化每个网格有ID和坐标范围 % 将灯光数据聚合到网格 grid_brightness zeros(num_grids, 1); for g 1:num_grids grid_mask (lon_grid grid_lon_min(g) lon_grid grid_lon_max(g) ... lat_grid grid_lat_min(g) lat_grid grid_lat_max(g)); grid_brightness(g) mean(DNB_cleaned(grid_mask), omitnan); end3.3 多源数据融合将处理好的灯光网格数据与人口密度网格数据WorldPop、GDP空间化数据夜间灯光常被用作代理变量、生态保护区矢量数据、道路网络数据等进行空间关联。这里大量使用了MATLAB的Mapping Toolbox和地理信息处理函数如shaperead读取矢量数据geoshow进行叠加展示。注意事项不同数据源的空间参考坐标系、分辨率必须统一。我们使用geotiffwrite输出统一坐标系和分辨率的栅格数据或使用vec2mtx将矢量数据转换为与灯光数据匹配的栅格这是后续所有空间分析的基础。坐标转换错误是导致整个空间分析失败的最常见原因。4. 核心模型实现与MATLAB代码详解有了干净的数据和清晰的框架就可以开始实现核心模型了。这里重点讲两个部分综合评价指数CEI的计算和系统动力学SD模型的模拟。4.1 综合评价指数CEI计算模块这个模块的输入是每个网格的多维指标数据输出是一个综合得分。function [CEI_scores, weights] calculateCEI(grid_data, params) % grid_data: 表格或矩阵每一行是一个网格列是各指标值 % 例如[亮度, 蓝光比例, 照明时长, 距生态区距离, 人口密度...] % params: 结构体包含权重计算方法选择、AHP判断矩阵等参数 [num_grids, num_indicators] size(grid_data); % 步骤1数据标准化消除量纲 % 采用极差标准化对于成本型指标如距生态区距离需取倒数或负向处理 normalized_data zeros(size(grid_data)); for i 1:num_indicators if params.is_benefit(i) % 效益型指标越大越好 min_val min(grid_data(:, i)); max_val max(grid_data(:, i)); normalized_data(:, i) (grid_data(:, i) - min_val) / (max_val - min_val eps); else % 成本型指标越小越好 min_val min(grid_data(:, i)); max_val max(grid_data(:, i)); normalized_data(:, i) (max_val - grid_data(:, i)) / (max_val - min_val eps); end end % 步骤2计算权重结合AHP与熵权法 % 计算熵权 p normalized_data ./ sum(normalized_data, 1); % 计算比重 e -sum(p .* log(p eps), 1) / log(num_grids); % 计算信息熵 d 1 - e; % 计算信息效用值 w_entropy d / sum(d); % 熵权 % 计算AHP权重调用前面定义的函数 w_ahp calculateAHP(params.ahp_matrix); % 综合权重这里采用简单线性组合也可用乘法合成 alpha 0.6; % 主观权重占比可根据实际情况调整 weights alpha * w_ahp (1-alpha) * w_entropy; weights weights / sum(weights); % 归一化 % 步骤3计算综合得分 CEI_scores normalized_data * weights; end4.2 系统动力学SD模型模拟我们在MATLAB中用ODE求解器来模拟这个包含反馈的系统。将SD的流图转化为一阶微分方程组。function dYdt lightPollutionSD(t, Y, params, intervention_func) % t: 时间 % Y: 状态变量向量 [CEI; GDP_per_capita; Policy_Strength; Public_Awareness] % params: 模型参数结构体增长率、影响系数等 % intervention_func: 外部干预策略函数句柄接受t和Y返回干预力度 CEI Y(1); GDP Y(2); Policy Y(3); Awareness Y(4); % 计算当前时刻的外部干预例如从某年开始实施强力的灯光管理 intervention intervention_func(t, Y); % 定义微分方程 % 1. CEI的变化率 经济驱动增长 - 政策抑制 - 自然衰减 干预影响 growth_driven params.k1 * GDP; % 经济驱动项 policy_suppression params.k2 * Policy * CEI; % 政策抑制项与现有污染水平和政策强度成正比 natural_decay params.k3 * CEI; % 自然衰减如技术自然更替 intervention_effect params.k4 * intervention; dCEIdt growth_driven - policy_suppression - natural_decay intervention_effect; % 2. GDP增长率简化假设与CEI有微弱正相关但受资源约束 dGDPdt params.gdp_growth_rate * GDP * (1 - GDP/params.gdp_capacity) params.k5 * CEI; % 3. 政策强度变化率随公众意识提高和CEI升高而增强但有延迟和饱和 dPolicydt params.k6 * Awareness params.k7 * max(CEI - params.CEI_threshold, 0) - params.k8 * Policy; % 4. 公众意识变化率受媒体报道与CEI相关和教育投入影响 dAwarenessdt params.k9 * CEI params.k10 * intervention - params.k11 * Awareness; dYdt [dCEIdt; dGDPdt; dPolicydt; dAwarenessdt]; end % 主模拟程序 params.k1 0.05; % 经济驱动系数 params.k2 0.1; % 政策抑制系数 % ... 其他参数通过历史数据拟合或文献设定 % 定义干预策略例如从第10年t10开始实施强度为0.5的干预 intervention_strategy (t, Y) 0.5 * (t 10); % 初始状态 Y0 [10; 10000; 0.1; 0.2]; % 初始CEI, GDP, 政策强度公众意识 tspan [0, 30]; % 模拟30年 % 使用ode45求解 [t, Y] ode45((t,Y) lightPollutionSD(t, Y, params, intervention_strategy), tspan, Y0); % 可视化结果 figure; subplot(2,2,1); plot(t, Y(:,1)); xlabel(年份); ylabel(综合污染指数 CEI); title(光污染趋势预测); grid on; % ... 绘制其他状态变量这个SD模型虽然经过了大量简化但它清晰地展示了系统中各要素的动态相互作用并能通过调整参数和干预函数来模拟不同发展情景下的污染走势为策略分析提供了“数字沙盘”。5. 干预策略的成本效益分析与优化模型预测出未来光污染会加剧那么该怎么办题目要求提出干预策略。我们将其构建为一个带约束的优化问题在有限的预算下选择一组干预措施如更换路灯类型、设定照明宵禁、划定暗天空保护区使得在规划期如20年末的CEI降低量最大或总成本效益比最高。5.1 策略量化与成本函数我们将空间网格作为决策单元。对每个网格i定义决策变量x_i表示施加的干预力度0到1之间。干预措施分为几类每类有单位成本c_j和单位效果系数e_j即降低CEI的能力。措施1将路灯更换为低色温LED暖黄光。成本较高对生态影响光谱改善大。措施2实施午夜后关闭半数路灯。成本低主要是管理成本对亮度改善直接。措施3在生态敏感区周边设置遮光罩。成本中等效果局限于局部。每个网格可以施加多种措施其总效果是叠加的总成本也是叠加的。5.2 优化模型建立目标函数最大化总效益CEI降低总量。 约束条件总预算不超过B。每个网格的干预力度x_i在[0,1]之间。某些敏感区域必须实施最低限度的干预x_i x_min。这是一个线性规划Linear Programming问题因为目标函数和约束条件关于决策变量x_i都是线性的。我们使用MATLAB的linprog函数求解。但需要注意的是linprog默认是求解最小化问题我们的最大化问题需要转换。% 假设有N个网格M种措施 N 1000; M 3; % 定义参数 cost_per_unit [50, 5, 20]; % 每种措施的单位成本万元/平方公里 effect_per_unit [0.3, 0.15, 0.25]; % 每种措施的单位CEI降低效果 budget 10000; % 总预算万元 sensitive_grids [10, 20, 30]; % 生态敏感网格编号 min_intervention 0.7; % 敏感网格最低干预力度 % 构建优化问题 % 决策变量 x 是一个 N*M 维的列向量按网格顺序排列[网格1_措施1, 网格1_措施2, ..., 网格N_措施M] f -repmat(effect_per_unit, 1, N); % 目标函数系数求最大化故取负 f f(:); % 不等式约束 A*x b % 1. 预算约束总成本 budget A_budget zeros(1, N*M); for i 1:N A_budget(1, (i-1)*M1 : i*M) cost_per_unit; end b_budget budget; % 2. 每个决策变量上限 1 A_ub eye(N*M); b_ub ones(N*M, 1); % 3. 敏感网格下限约束 -x_i -min_intervention A_lb_sensitive zeros(length(sensitive_grids), N*M); for idx 1:length(sensitive_grids) grid_id sensitive_grids(idx); for m 1:M A_lb_sensitive(idx, (grid_id-1)*M m) -1; % 构造 -x -min 即 x min end end b_lb_sensitive -min_intervention * ones(length(sensitive_grids), 1); % 合并不等式约束 A [A_budget; A_ub; A_lb_sensitive]; b [b_budget; b_ub; b_lb_sensitive]; % 变量下界非负 lb zeros(N*M, 1); % 求解线性规划 options optimoptions(linprog, Display, iter, Algorithm, dual-simplex); [x_opt, fval, exitflag] linprog(f, A, b, [], [], lb, [], options); if exitflag 0 disp(优化求解成功); total_cost A_budget * x_opt; total_effect -fval; % 恢复为正的效益值 fprintf(最优方案总成本%.2f 万元预计降低CEI%.4f\n, total_cost, total_effect); % 将解向量 x_opt 重塑为 N x M 矩阵便于分析每个网格的措施组合 decision_matrix reshape(x_opt, M, N); else error(优化求解失败); end5.3 结果可视化与策略解读求解后我们得到了一张“干预地图”。通过imagesc或geoshow将每个网格的最优干预力度或首选措施类型进行空间可视化。% 将决策矩阵中每个网格的最大干预措施类型找出 [~, optimal_measure] max(decision_matrix, [], 2); % 假设有网格的经纬度坐标网格矩阵 lat_grid, lon_grid % 创建一个与地理坐标匹配的干预力度图层 intervention_layer zeros(size(lat_grid)); for i 1:N % 找到该网格对应的所有像素位置 mask (grid_id_map i); % grid_id_map是每个像素所属网格的编号图 intervention_layer(mask) sum(decision_matrix(i, :)); % 或用 optimal_measure(i) 显示措施类型 end figure; worldmap([min(lat_lim) max(lat_lim)], [min(lon_lim) max(lon_lim)]); geoshow(lat_grid, lon_grid, intervention_layer, DisplayType, texturemap); colorbar; title(最优光污染干预策略空间分布图);这张图直观地告诉决策者钱应该优先花在哪里。我们的分析结果显示在市中心商业区推广智能调光措施2性价比最高而在城市边缘的湿地公园周边则必须投入较高成本更换光源并加装遮光设施措施1和3。我们还可以进行敏感性分析模拟预算增减10%、20%对总效果的影响为决策提供弹性空间。实操心得在论文中我们不仅给出了优化结果还讨论了模型的局限性。例如我们的成本函数是线性的但现实中可能存在规模效应批量采购折扣或阈值效应低于某个投资额度无效。我们在附录中简要提及了可以考虑非线性规划或整数规划是否对某个区域实施干预作为模型扩展这体现了思考的深度。6. 论文写作要点与可视化技巧美赛论文模型和结果再漂亮也需要通过论文清晰有力地表达出来。MATLAB不仅在计算上强大在生成论文所需的图表方面也是一把利器。6.1 图表可视化要点时空动态展示对于灯光数据不要只放一张静态图。我们制作了时间序列动画展示研究区域2018-2022年光污染的蔓延过程。使用for循环更新imagesc的数据并用getframe和VideoWriter保存为视频或GIF插入论文附录或提供链接。writerObj VideoWriter(light_pollution_evolution.mp4, MPEG-4); open(writerObj); fig figure; for year 1:5 % 绘制第year年的灯光数据 imagesc(lon_range, lat_range, annual_data(:,:,year)); axis equal tight; colormap(jet); colorbar; clim([0 max_light]); % 统一颜色刻度 title(sprintf(Annual Mean Light Pollution - Year %d, 2017year)); xlabel(Longitude); ylabel(Latitude); frame getframe(fig); writeVideo(writerObj, frame); end close(writerObj);多图对比与组合使用subplot将CEI空间分布图、GDP分布图、干预策略图放在一起对比清晰展示其空间相关性。figure(Position, [100, 100, 1200, 400]); subplot(1,3,1); geoshow(CEI_map, R_CEI, DisplayType, texturemap); title((a) 光污染综合指数 (CEI)); colorbar; subplot(1,3,2); geoshow(GDP_density_map, R_GDP, DisplayType, texturemap); title((b) GDP空间密度); colorbar; subplot(1,3,3); geoshow(intervention_layer, R, DisplayType, texturemap); title((c) 推荐干预优先级); colormap(jet); colorbar;统计图表专业化预测趋势图一定要加上置信区间。使用polyfit和polyval进行拟合后用predint函数计算预测区间。[p, S] polyfit(year, CEI_trend, 1); [y_fit, delta] polyval(p, year_future, S); % year_future是未来年份向量 plot(year, CEI_trend, ko, MarkerFaceColor, b); hold on; plot(year_future, y_fit, r-, LineWidth, 2); fill([year_future, fliplr(year_future)], ... [y_fitdelta, fliplr(y_fit-delta)], ... r, FaceAlpha, 0.2, EdgeColor, none); xlabel(Year); ylabel(CEI); legend(Historical Data, Fitted Trend, 95% Prediction Interval, Location, best); grid on;6.2 论文行文建议摘要用“我们建立了...模型M1该模型包含...指标采用...方法合成。基于此我们构建了...预测模型M2并模拟了三种情景。最后我们提出了一个...优化模型M3以...为目标...为约束给出了分区的策略建议。敏感性分析表明...。我们的模型优点在于...”模型假设清晰列出并说明其合理性。例如“假设不同光源的光谱权重在短期内保持不变”、“假设政策执行力度与政策文本强度成正比”。灵敏度分析必须做。改变AHP判断矩阵、调整SD模型中的关键参数如k1, k2观察CEI预测结果和优化方案的变化幅度。这能极大增强模型的鲁棒性。优缺点与推广客观评价。优点可以是“模型融合了主观与客观赋权结合了宏观系统动态与微观数据驱动”。缺点要诚实如“模型未考虑大气传输对光污染扩散的影响”、“成本数据基于估算可能与实际有出入”。推广可以提到“该框架可应用于噪声污染、水污染等其它环境问题的评估与管理”。7. 常见问题与代码调试心得在实战中肯定会遇到各种坑。这里分享几个我们当时遇到并解决了的问题。7.1 数据内存不足夜间灯光数据分辨率高多年多期数据直接读入内存可能导致MATLAB崩溃。解决方案使用blockproc函数进行分块处理或者使用imread时只读取需要的波段和区域。对于时间序列分析可以先将每年的数据聚合为年度均值一个较低分辨率的数据集再进行后续分析。7.2 空间匹配错误多源数据叠加时出现错位这是坐标系CRS不统一导致的。解决方案务必使用geotiffinfo查看每个.tif文件的GeoTIFFTags信息确认其Projection和PCS平面坐标系或GCS地理坐标系。使用projfwd和projinv进行坐标转换或利用QGIS等专业软件将所有数据统一转换为同一坐标系如WGS84 Web Mercator后再导入MATLAB。7.3 模型预测结果不合理SD模型跑出来的CEI指数爆炸增长或快速归零。排查步骤检查参数量纲确保所有参数增长率、系数的量级是合理的。例如一个0.5的年增长率已经很高了。检查反馈回路符号确认是正反馈还是负反馈。政策抑制项应该是-k * Policy * CEI如果符号错了效果就反了。缩小时间步长ode45是变步长算法通常很稳健。但如果方程刚性很强可以尝试ode15s。进行单因素测试将其他变量固定只改变一个变量看输出是否符合预期逻辑。7.4 优化问题无解linprog提示“No feasible solution found”。原因与解决约束过严预算太低而敏感区域的最低干预要求又太高导致没有可行域。尝试放松约束如降低min_intervention或增加预算。变量上下界冲突检查lb和ub是否合理。将问题可视化对于简单的二维或三维问题可以用plot画出约束区域和目标函数等高线直观查看是否存在可行域。7.5 绘图美观性问题论文图表看起来不专业。技巧统一风格所有图表使用相同的配色方案如parula,viridis字体大小一致。去除白边保存图片时使用exportgraphics(gca, figure.png, Resolution, 300)或print -dpng -r300 figure.png并调整Position和PaperPosition属性。使用紧凑布局subplot之前使用tiledlayout能更好地控制子图间距。最后想说的是美赛E题这类开放性问题没有标准答案。评委看重的是你将现实问题转化为数学问题的能力、模型构建的逻辑严谨性、解决方案的创造性与可行性以及结果呈现的清晰度。MATLAB是你实现这一切想法的强大工具但工具背后的思考才是灵魂。从理解“光污染”这个现象开始一步步拆解、量化、建模、优化直到给出一个有理有据、图文并茂的解决方案这个过程本身就是一次完整的科研训练。希望这份结合了思路和代码的复盘能让你在未来的比赛中少走弯路更高效地将你的创意转化为实实在在的论文和代码。 MATLAB建模实战：从光污染量化到干预策略优化
1. 项目概述从“光污染”到量化分析去年带队做美赛E题“光污染”的经历现在回想起来依然觉得是一次非常典型的建模实战。题目给了一个看似宏大的社会与环境议题但内核却要求我们将其转化为一系列可量化、可计算、可预测的数学问题。很多队伍一开始就懵了要么陷在“光污染危害”的文献综述里出不来要么对着“评价指标”和“控制模型”这些词无从下手。其实这道题的核心思路非常清晰定义问题 - 建立评价体系 - 构建预测模型 - 提出干预策略。整个过程MATLAB作为核心计算与可视化工具贯穿始终。简单来说这道题要求我们做三件事第一为光污染建立一个综合性的评价指标体系不能只看亮度还得考虑光谱、时间、空间和对生态/社区的影响。第二基于这个体系建立一个能预测光污染水平如何随城市发展如人口、GDP、土地类型变化而演变的数学模型。第三利用这个模型设计一套成本效益最优的干预策略告诉决策者“在哪里、用什么方法、花多少钱”能最有效地控制光污染。这不仅仅是套个公式那么简单。你需要处理多源异构数据卫星遥感夜间灯光数据、城市统计数据、地理信息数据需要融合环境科学、经济学和公共政策的知识更需要用MATLAB将这一切想法高效地实现出来从数据清洗、空间分析、模型拟合到策略模拟。接下来我就结合我们当时的解题思路和具体的MATLAB代码实现拆解每一个关键环节希望能给未来参赛的你一些实实在在的参考。2. 解题核心思路与模型框架设计面对“光污染”这种开放性问题最忌讳的就是一上来就埋头写代码。我们花了将近半天的时间来搭建整个解题的逻辑框架确保每一步都有的放矢。我们的整体思路可以概括为“一个核心两个维度三层模型”。2.1 问题拆解与核心定义题目中的“光污染”是一个模糊概念。我们首先要将其操作化为可测量的指标。我们将其分解为两个核心维度物理强度维度即光本身的性质。包括亮度最直接的指标来自卫星夜间灯光数据如VIIRS DNB数据的辐射值。光谱构成不同波长的光对生态和人类的影响不同。例如短波蓝光对天文观测干扰大且更易影响人体节律。我们通过查阅文献为常见光源高压钠灯、LED白光、LED暖光赋予不同的光谱权重系数。时间动态通宵照明与午夜后关闭的照明其污染程度截然不同。我们引入了“夜间光照持续时间”作为时间维度指标。影响效应维度即光对受体产生的影响。包括生态影响重点考虑对夜间迁徙鸟类、海龟、昆虫的干扰。我们根据土地利用类型如湿地、森林、海岸线设定生态敏感区光辐射到这些区域则产生“生态影响分值”。社区影响考虑对人类健康睡眠干扰、交通安全眩光、夜空观赏性星空可见度的影响。这部分更多依赖于人口密度数据和问卷调查数据的模拟。2.2 综合评价指数CEI的构建将上述多维指标合成一个综合性的“光污染评价指数”Comprehensive Evaluation Index, CEI是关键一步。我们采用了加权线性综合法但权重的确定不是主观臆断。我们使用了层次分析法AHP结合熵权法来确定权重。AHP用于整合我们通过文献和专家咨询模拟得到的各指标相对重要性判断熵权法则基于实际数据如各区域亮度值的差异度计算客观权重。最后将主客观权重结合得到最终权重。在MATLAB中AHP的实现需要构建判断矩阵并进行一致性检验。% 示例AHP权重计算核心片段 function weights calculateAHP(matrix) % matrix: 判断矩阵 [n, ~] size(matrix); [V, D] eig(matrix); % 求特征值和特征向量 [~, idx] max(diag(D)); % 找到最大特征值 w V(:, idx); % 对应的特征向量 weights w / sum(w); % 归一化得到权重 % 一致性检验 CI, RI, CR lambda_max max(diag(D)); CI (lambda_max - n) / (n - 1); RI [0, 0, 0.58, 0.90, 1.12, 1.24, 1.32, 1.41, 1.45]; % 平均随机一致性指标 CR CI / RI(n); if CR 0.1 error(判断矩阵一致性不满足要求请调整); end end2.3 预测模型的选择系统动力学SD与多元回归的融合题目要求预测未来光污染趋势。我们放弃了单一的回归模型因为光污染系统内部存在复杂的反馈关系例如经济发展推动照明需求光污染加剧可能引发政策调控。因此我们引入了系统动力学System Dynamics, SD模型作为核心预测框架。我们将系统划分为几个关键存量光污染综合指数CEI社会经济驱动因子如人均GDP、城市化率政策干预强度然后定义流量和反馈回路。例如“经济发展”正反馈于“照明基础设施投资”从而增加“光辐射源强度”而“公众环保意识”和“已实施政策”会负反馈于“新照明项目的审批尺度”。实操心得直接用Vensim等SD专用软件建模当然可以但在美赛有限时间内为了与MATLAB数据处理和后续优化无缝衔接我们在MATLAB中用差分方程组来近似实现这个SD模型的核心逻辑。这要求我们对模型进行合理的简化抓住主要矛盾。同时SD模型中的一些参数关系如GDP增长率对灯光亮度增长率的影响系数需要通过历史数据进行标定。这里我们用到了多元非线性回归。MATLAB的fitnlm函数在这里大显身手。% 示例使用多元非线性回归拟合SD模型参数 % 假设我们有历史数据年份Year GDP 灯光亮度Brightness % 建立关系Brightness a * GDP^b c % 准备数据 tbl table(GDP, Brightness, VariableNames, {GDP,Brightness}); % 定义非线性模型公式 modelfun (b,x) b(1) * x.GDP.^b(2) b(3); % 初始参数猜测 beta0 [0.01, 0.5, 10]; % 拟合模型 nlmodel fitnlm(tbl, modelfun, beta0); % 查看拟合结果和参数 disp(nlmodel.Coefficients.Estimate); plot(nlmodel); % 可视化拟合效果这个融合框架的优势在于既有系统思考的宏观视角SD又有数据驱动的微观精度回归使得我们的预测模型既有说服力又具备一定的稳健性。3. 数据获取、处理与空间分析实战模型框架搭好了接下来就是喂数据。数据决定了模型的上限。E题通常不提供现成数据需要自己寻找。我们主要使用了以下几类数据并全部在MATLAB中完成处理。3.1 核心数据源夜间灯光数据NASA的VIIRS DNB可见光红外成像辐射计套件-日/夜波段月度数据是研究光污染的黄金标准。可以从NOAA官网下载https://eogdata.mines.edu/products/vnl/。我们下载了目标区域我们自选了一个大都市圈最近5年的月度平均辐射合成数据.tif格式。3.2 数据处理全流程数据读取与裁剪使用geotiffread读取TIFF文件获取图像数据和地理参考信息。然后用目标区域的矢量边界Shapefile通过mapcrop或地理坐标匹配进行裁剪。[DNB_data, R] readgeoraster(viirs_202201.tif); % 假设有边界坐标向量latlim, lonlim [DNB_cropped, R_cropped] mapcrop(DNB_data, R, latlim, lonlim);异常值与背景值处理夜间灯光数据包含背景噪声气辉、月光等和异常高值火光、油气燃烧。我们采用分位数阈值法进行过滤。lower_prc prctile(DNB_cropped(DNB_cropped0), 1); % 取正值数据的1%分位数 upper_prc prctile(DNB_cropped(:), 99.9); % 取全体数据的99.9%分位数 DNB_cleaned DNB_cropped; DNB_cleaned(DNB_cleaned lower_prc) 0; % 低于下限的视为背景噪声 DNB_cleaned(DNB_cleaned upper_prc) upper_prc; % 高于上限的进行截断避免极端值影响时间序列整合与趋势分析将月度数据按年度平均得到年际序列。使用detrend函数或线性拟合来观察亮度总体增长趋势。yearly_mean zeros(5,1); for i 1:5 % 读取第i年的12个月数据求平均 yearly_mean(i) mean(annual_data_i(:)); end % 线性趋势拟合 p polyfit((1:5), yearly_mean, 1); trend polyval(p, (1:5));空间网格化与指标计算将研究区域划分为1km x 1km的网格。对每个网格计算其平均亮度、亮度年增长率、以及其所属的土地利用类型来自全球土地覆盖数据如FROM-GLC。% 假设已将区域网格化每个网格有ID和坐标范围 % 将灯光数据聚合到网格 grid_brightness zeros(num_grids, 1); for g 1:num_grids grid_mask (lon_grid grid_lon_min(g) lon_grid grid_lon_max(g) ... lat_grid grid_lat_min(g) lat_grid grid_lat_max(g)); grid_brightness(g) mean(DNB_cleaned(grid_mask), omitnan); end3.3 多源数据融合将处理好的灯光网格数据与人口密度网格数据WorldPop、GDP空间化数据夜间灯光常被用作代理变量、生态保护区矢量数据、道路网络数据等进行空间关联。这里大量使用了MATLAB的Mapping Toolbox和地理信息处理函数如shaperead读取矢量数据geoshow进行叠加展示。注意事项不同数据源的空间参考坐标系、分辨率必须统一。我们使用geotiffwrite输出统一坐标系和分辨率的栅格数据或使用vec2mtx将矢量数据转换为与灯光数据匹配的栅格这是后续所有空间分析的基础。坐标转换错误是导致整个空间分析失败的最常见原因。4. 核心模型实现与MATLAB代码详解有了干净的数据和清晰的框架就可以开始实现核心模型了。这里重点讲两个部分综合评价指数CEI的计算和系统动力学SD模型的模拟。4.1 综合评价指数CEI计算模块这个模块的输入是每个网格的多维指标数据输出是一个综合得分。function [CEI_scores, weights] calculateCEI(grid_data, params) % grid_data: 表格或矩阵每一行是一个网格列是各指标值 % 例如[亮度, 蓝光比例, 照明时长, 距生态区距离, 人口密度...] % params: 结构体包含权重计算方法选择、AHP判断矩阵等参数 [num_grids, num_indicators] size(grid_data); % 步骤1数据标准化消除量纲 % 采用极差标准化对于成本型指标如距生态区距离需取倒数或负向处理 normalized_data zeros(size(grid_data)); for i 1:num_indicators if params.is_benefit(i) % 效益型指标越大越好 min_val min(grid_data(:, i)); max_val max(grid_data(:, i)); normalized_data(:, i) (grid_data(:, i) - min_val) / (max_val - min_val eps); else % 成本型指标越小越好 min_val min(grid_data(:, i)); max_val max(grid_data(:, i)); normalized_data(:, i) (max_val - grid_data(:, i)) / (max_val - min_val eps); end end % 步骤2计算权重结合AHP与熵权法 % 计算熵权 p normalized_data ./ sum(normalized_data, 1); % 计算比重 e -sum(p .* log(p eps), 1) / log(num_grids); % 计算信息熵 d 1 - e; % 计算信息效用值 w_entropy d / sum(d); % 熵权 % 计算AHP权重调用前面定义的函数 w_ahp calculateAHP(params.ahp_matrix); % 综合权重这里采用简单线性组合也可用乘法合成 alpha 0.6; % 主观权重占比可根据实际情况调整 weights alpha * w_ahp (1-alpha) * w_entropy; weights weights / sum(weights); % 归一化 % 步骤3计算综合得分 CEI_scores normalized_data * weights; end4.2 系统动力学SD模型模拟我们在MATLAB中用ODE求解器来模拟这个包含反馈的系统。将SD的流图转化为一阶微分方程组。function dYdt lightPollutionSD(t, Y, params, intervention_func) % t: 时间 % Y: 状态变量向量 [CEI; GDP_per_capita; Policy_Strength; Public_Awareness] % params: 模型参数结构体增长率、影响系数等 % intervention_func: 外部干预策略函数句柄接受t和Y返回干预力度 CEI Y(1); GDP Y(2); Policy Y(3); Awareness Y(4); % 计算当前时刻的外部干预例如从某年开始实施强力的灯光管理 intervention intervention_func(t, Y); % 定义微分方程 % 1. CEI的变化率 经济驱动增长 - 政策抑制 - 自然衰减 干预影响 growth_driven params.k1 * GDP; % 经济驱动项 policy_suppression params.k2 * Policy * CEI; % 政策抑制项与现有污染水平和政策强度成正比 natural_decay params.k3 * CEI; % 自然衰减如技术自然更替 intervention_effect params.k4 * intervention; dCEIdt growth_driven - policy_suppression - natural_decay intervention_effect; % 2. GDP增长率简化假设与CEI有微弱正相关但受资源约束 dGDPdt params.gdp_growth_rate * GDP * (1 - GDP/params.gdp_capacity) params.k5 * CEI; % 3. 政策强度变化率随公众意识提高和CEI升高而增强但有延迟和饱和 dPolicydt params.k6 * Awareness params.k7 * max(CEI - params.CEI_threshold, 0) - params.k8 * Policy; % 4. 公众意识变化率受媒体报道与CEI相关和教育投入影响 dAwarenessdt params.k9 * CEI params.k10 * intervention - params.k11 * Awareness; dYdt [dCEIdt; dGDPdt; dPolicydt; dAwarenessdt]; end % 主模拟程序 params.k1 0.05; % 经济驱动系数 params.k2 0.1; % 政策抑制系数 % ... 其他参数通过历史数据拟合或文献设定 % 定义干预策略例如从第10年t10开始实施强度为0.5的干预 intervention_strategy (t, Y) 0.5 * (t 10); % 初始状态 Y0 [10; 10000; 0.1; 0.2]; % 初始CEI, GDP, 政策强度公众意识 tspan [0, 30]; % 模拟30年 % 使用ode45求解 [t, Y] ode45((t,Y) lightPollutionSD(t, Y, params, intervention_strategy), tspan, Y0); % 可视化结果 figure; subplot(2,2,1); plot(t, Y(:,1)); xlabel(年份); ylabel(综合污染指数 CEI); title(光污染趋势预测); grid on; % ... 绘制其他状态变量这个SD模型虽然经过了大量简化但它清晰地展示了系统中各要素的动态相互作用并能通过调整参数和干预函数来模拟不同发展情景下的污染走势为策略分析提供了“数字沙盘”。5. 干预策略的成本效益分析与优化模型预测出未来光污染会加剧那么该怎么办题目要求提出干预策略。我们将其构建为一个带约束的优化问题在有限的预算下选择一组干预措施如更换路灯类型、设定照明宵禁、划定暗天空保护区使得在规划期如20年末的CEI降低量最大或总成本效益比最高。5.1 策略量化与成本函数我们将空间网格作为决策单元。对每个网格i定义决策变量x_i表示施加的干预力度0到1之间。干预措施分为几类每类有单位成本c_j和单位效果系数e_j即降低CEI的能力。措施1将路灯更换为低色温LED暖黄光。成本较高对生态影响光谱改善大。措施2实施午夜后关闭半数路灯。成本低主要是管理成本对亮度改善直接。措施3在生态敏感区周边设置遮光罩。成本中等效果局限于局部。每个网格可以施加多种措施其总效果是叠加的总成本也是叠加的。5.2 优化模型建立目标函数最大化总效益CEI降低总量。 约束条件总预算不超过B。每个网格的干预力度x_i在[0,1]之间。某些敏感区域必须实施最低限度的干预x_i x_min。这是一个线性规划Linear Programming问题因为目标函数和约束条件关于决策变量x_i都是线性的。我们使用MATLAB的linprog函数求解。但需要注意的是linprog默认是求解最小化问题我们的最大化问题需要转换。% 假设有N个网格M种措施 N 1000; M 3; % 定义参数 cost_per_unit [50, 5, 20]; % 每种措施的单位成本万元/平方公里 effect_per_unit [0.3, 0.15, 0.25]; % 每种措施的单位CEI降低效果 budget 10000; % 总预算万元 sensitive_grids [10, 20, 30]; % 生态敏感网格编号 min_intervention 0.7; % 敏感网格最低干预力度 % 构建优化问题 % 决策变量 x 是一个 N*M 维的列向量按网格顺序排列[网格1_措施1, 网格1_措施2, ..., 网格N_措施M] f -repmat(effect_per_unit, 1, N); % 目标函数系数求最大化故取负 f f(:); % 不等式约束 A*x b % 1. 预算约束总成本 budget A_budget zeros(1, N*M); for i 1:N A_budget(1, (i-1)*M1 : i*M) cost_per_unit; end b_budget budget; % 2. 每个决策变量上限 1 A_ub eye(N*M); b_ub ones(N*M, 1); % 3. 敏感网格下限约束 -x_i -min_intervention A_lb_sensitive zeros(length(sensitive_grids), N*M); for idx 1:length(sensitive_grids) grid_id sensitive_grids(idx); for m 1:M A_lb_sensitive(idx, (grid_id-1)*M m) -1; % 构造 -x -min 即 x min end end b_lb_sensitive -min_intervention * ones(length(sensitive_grids), 1); % 合并不等式约束 A [A_budget; A_ub; A_lb_sensitive]; b [b_budget; b_ub; b_lb_sensitive]; % 变量下界非负 lb zeros(N*M, 1); % 求解线性规划 options optimoptions(linprog, Display, iter, Algorithm, dual-simplex); [x_opt, fval, exitflag] linprog(f, A, b, [], [], lb, [], options); if exitflag 0 disp(优化求解成功); total_cost A_budget * x_opt; total_effect -fval; % 恢复为正的效益值 fprintf(最优方案总成本%.2f 万元预计降低CEI%.4f\n, total_cost, total_effect); % 将解向量 x_opt 重塑为 N x M 矩阵便于分析每个网格的措施组合 decision_matrix reshape(x_opt, M, N); else error(优化求解失败); end5.3 结果可视化与策略解读求解后我们得到了一张“干预地图”。通过imagesc或geoshow将每个网格的最优干预力度或首选措施类型进行空间可视化。% 将决策矩阵中每个网格的最大干预措施类型找出 [~, optimal_measure] max(decision_matrix, [], 2); % 假设有网格的经纬度坐标网格矩阵 lat_grid, lon_grid % 创建一个与地理坐标匹配的干预力度图层 intervention_layer zeros(size(lat_grid)); for i 1:N % 找到该网格对应的所有像素位置 mask (grid_id_map i); % grid_id_map是每个像素所属网格的编号图 intervention_layer(mask) sum(decision_matrix(i, :)); % 或用 optimal_measure(i) 显示措施类型 end figure; worldmap([min(lat_lim) max(lat_lim)], [min(lon_lim) max(lon_lim)]); geoshow(lat_grid, lon_grid, intervention_layer, DisplayType, texturemap); colorbar; title(最优光污染干预策略空间分布图);这张图直观地告诉决策者钱应该优先花在哪里。我们的分析结果显示在市中心商业区推广智能调光措施2性价比最高而在城市边缘的湿地公园周边则必须投入较高成本更换光源并加装遮光设施措施1和3。我们还可以进行敏感性分析模拟预算增减10%、20%对总效果的影响为决策提供弹性空间。实操心得在论文中我们不仅给出了优化结果还讨论了模型的局限性。例如我们的成本函数是线性的但现实中可能存在规模效应批量采购折扣或阈值效应低于某个投资额度无效。我们在附录中简要提及了可以考虑非线性规划或整数规划是否对某个区域实施干预作为模型扩展这体现了思考的深度。6. 论文写作要点与可视化技巧美赛论文模型和结果再漂亮也需要通过论文清晰有力地表达出来。MATLAB不仅在计算上强大在生成论文所需的图表方面也是一把利器。6.1 图表可视化要点时空动态展示对于灯光数据不要只放一张静态图。我们制作了时间序列动画展示研究区域2018-2022年光污染的蔓延过程。使用for循环更新imagesc的数据并用getframe和VideoWriter保存为视频或GIF插入论文附录或提供链接。writerObj VideoWriter(light_pollution_evolution.mp4, MPEG-4); open(writerObj); fig figure; for year 1:5 % 绘制第year年的灯光数据 imagesc(lon_range, lat_range, annual_data(:,:,year)); axis equal tight; colormap(jet); colorbar; clim([0 max_light]); % 统一颜色刻度 title(sprintf(Annual Mean Light Pollution - Year %d, 2017year)); xlabel(Longitude); ylabel(Latitude); frame getframe(fig); writeVideo(writerObj, frame); end close(writerObj);多图对比与组合使用subplot将CEI空间分布图、GDP分布图、干预策略图放在一起对比清晰展示其空间相关性。figure(Position, [100, 100, 1200, 400]); subplot(1,3,1); geoshow(CEI_map, R_CEI, DisplayType, texturemap); title((a) 光污染综合指数 (CEI)); colorbar; subplot(1,3,2); geoshow(GDP_density_map, R_GDP, DisplayType, texturemap); title((b) GDP空间密度); colorbar; subplot(1,3,3); geoshow(intervention_layer, R, DisplayType, texturemap); title((c) 推荐干预优先级); colormap(jet); colorbar;统计图表专业化预测趋势图一定要加上置信区间。使用polyfit和polyval进行拟合后用predint函数计算预测区间。[p, S] polyfit(year, CEI_trend, 1); [y_fit, delta] polyval(p, year_future, S); % year_future是未来年份向量 plot(year, CEI_trend, ko, MarkerFaceColor, b); hold on; plot(year_future, y_fit, r-, LineWidth, 2); fill([year_future, fliplr(year_future)], ... [y_fitdelta, fliplr(y_fit-delta)], ... r, FaceAlpha, 0.2, EdgeColor, none); xlabel(Year); ylabel(CEI); legend(Historical Data, Fitted Trend, 95% Prediction Interval, Location, best); grid on;6.2 论文行文建议摘要用“我们建立了...模型M1该模型包含...指标采用...方法合成。基于此我们构建了...预测模型M2并模拟了三种情景。最后我们提出了一个...优化模型M3以...为目标...为约束给出了分区的策略建议。敏感性分析表明...。我们的模型优点在于...”模型假设清晰列出并说明其合理性。例如“假设不同光源的光谱权重在短期内保持不变”、“假设政策执行力度与政策文本强度成正比”。灵敏度分析必须做。改变AHP判断矩阵、调整SD模型中的关键参数如k1, k2观察CEI预测结果和优化方案的变化幅度。这能极大增强模型的鲁棒性。优缺点与推广客观评价。优点可以是“模型融合了主观与客观赋权结合了宏观系统动态与微观数据驱动”。缺点要诚实如“模型未考虑大气传输对光污染扩散的影响”、“成本数据基于估算可能与实际有出入”。推广可以提到“该框架可应用于噪声污染、水污染等其它环境问题的评估与管理”。7. 常见问题与代码调试心得在实战中肯定会遇到各种坑。这里分享几个我们当时遇到并解决了的问题。7.1 数据内存不足夜间灯光数据分辨率高多年多期数据直接读入内存可能导致MATLAB崩溃。解决方案使用blockproc函数进行分块处理或者使用imread时只读取需要的波段和区域。对于时间序列分析可以先将每年的数据聚合为年度均值一个较低分辨率的数据集再进行后续分析。7.2 空间匹配错误多源数据叠加时出现错位这是坐标系CRS不统一导致的。解决方案务必使用geotiffinfo查看每个.tif文件的GeoTIFFTags信息确认其Projection和PCS平面坐标系或GCS地理坐标系。使用projfwd和projinv进行坐标转换或利用QGIS等专业软件将所有数据统一转换为同一坐标系如WGS84 Web Mercator后再导入MATLAB。7.3 模型预测结果不合理SD模型跑出来的CEI指数爆炸增长或快速归零。排查步骤检查参数量纲确保所有参数增长率、系数的量级是合理的。例如一个0.5的年增长率已经很高了。检查反馈回路符号确认是正反馈还是负反馈。政策抑制项应该是-k * Policy * CEI如果符号错了效果就反了。缩小时间步长ode45是变步长算法通常很稳健。但如果方程刚性很强可以尝试ode15s。进行单因素测试将其他变量固定只改变一个变量看输出是否符合预期逻辑。7.4 优化问题无解linprog提示“No feasible solution found”。原因与解决约束过严预算太低而敏感区域的最低干预要求又太高导致没有可行域。尝试放松约束如降低min_intervention或增加预算。变量上下界冲突检查lb和ub是否合理。将问题可视化对于简单的二维或三维问题可以用plot画出约束区域和目标函数等高线直观查看是否存在可行域。7.5 绘图美观性问题论文图表看起来不专业。技巧统一风格所有图表使用相同的配色方案如parula,viridis字体大小一致。去除白边保存图片时使用exportgraphics(gca, figure.png, Resolution, 300)或print -dpng -r300 figure.png并调整Position和PaperPosition属性。使用紧凑布局subplot之前使用tiledlayout能更好地控制子图间距。最后想说的是美赛E题这类开放性问题没有标准答案。评委看重的是你将现实问题转化为数学问题的能力、模型构建的逻辑严谨性、解决方案的创造性与可行性以及结果呈现的清晰度。MATLAB是你实现这一切想法的强大工具但工具背后的思考才是灵魂。从理解“光污染”这个现象开始一步步拆解、量化、建模、优化直到给出一个有理有据、图文并茂的解决方案这个过程本身就是一次完整的科研训练。希望这份结合了思路和代码的复盘能让你在未来的比赛中少走弯路更高效地将你的创意转化为实实在在的论文和代码。