辐射换热仿真-主题049_城市热岛效应辐射分析-主题049_城市热岛效应辐射分析

kkchenjj

144人浏览 · 2026-04-27 21:31:32

kkchenjj · 2026-04-27 21:31:32 发布

主题049：城市热岛效应辐射分析

摘要

城市热岛效应是城市化进程中产生的重要环境问题，其本质是城市地表辐射能量平衡的改变。本主题从辐射换热的角度深入分析城市热岛效应的形成机制，研究不同地表材料（沥青、混凝土、植被、水体等）的辐射特性对城市热环境的影响。通过建立城市冠层辐射模型，分析太阳短波辐射的吸收、地表长波辐射的发射与吸收、以及建筑之间的多次反射过程。结合Python仿真程序，模拟不同城市规划方案下的热岛强度变化，为城市微气候优化和热环境改善提供科学依据。

关键词： 城市热岛、地表辐射、反照率、热惯性、微气候、辐射平衡、冠层模型

在这里插入图片描述

1. 城市热岛效应概述

1.1 热岛效应的定义与分类

城市热岛（Urban Heat Island, UHI）是指城市区域的气温明显高于周边郊区的现象。这一概念最早由英国气象学家Luke Howard于1818年在研究伦敦城市气候时提出。经过近两个世纪的研究，科学家们发现城市热岛是一个复杂的多尺度现象，可以从不同角度进行分类：

按空间尺度分类：

城市尺度热岛（Urban Scale UHI）：覆盖整个城市区域，水平范围从几公里到几十公里，热岛强度（城乡温差）通常在1-5°C之间，极端情况下可达10°C以上。这种热岛主要影响城市边界层气候，与天气系统相互作用。
街区尺度热岛（Neighborhood Scale UHI）：覆盖几个街区，水平范围几百米到几公里，热岛强度可达2-8°C。这种热岛与局地建筑布局、街道走向、绿化分布密切相关。
微尺度热岛（Micro Scale UHI）：发生在建筑尺度，水平范围几十米到几百米，热岛强度可达5-15°C。这种热岛直接影响人体热舒适，与行人高度处的热环境密切相关。

按形成机制分类：

地表热岛（Surface UHI）：指城市地表温度高于郊区的现象，通过卫星遥感或航空热像仪测量。地表热岛强度通常大于气温热岛，白天可达10-20°C，夜间5-10°C。
冠层热岛（Canopy Layer UHI）：指城市冠层（地面到建筑平均高度）内的气温高于郊区的现象，通过气象站观测获得。这是人们日常生活中直接感受到的热岛。
边界层热岛（Boundary Layer UHI）：指城市边界层（建筑平均高度到混合层顶）内的气温高于郊区的现象，通过气象气球或飞机观测获得。

1.2 热岛效应的时空特征

城市热岛具有明显的日变化和季节变化特征：

日变化特征：

城市热岛的日变化呈现典型的"夜间强、白天弱"特征。白天，城市地表吸收大量太阳辐射而升温，但由于城市下垫面热容量大，升温速度反而比郊区慢，此时热岛强度较弱甚至可能出现"冷岛"现象。夜间，城市储存的热量缓慢释放，而郊区地表降温迅速，城乡温差逐渐增大，通常在夜间22:00至凌晨04:00达到最大值。

这种日变化特征与辐射过程密切相关：

白天：太阳短波辐射主导能量平衡。城市反照率通常低于郊区（沥青路面0.05-0.15 vs 草地0.20-0.30），吸收更多太阳辐射。但城市蒸散发量小，感热通量比例高，同时城市粗糙度大，热量向大气扩散的效率相对较低。
夜间：长波辐射冷却主导能量平衡。城市天空视角因子小（建筑遮挡），地表向天空的长波辐射损失减少。同时，城市人为热排放（空调、交通、工业）在夜间持续，进一步维持城市高温。

季节变化特征：

城市热岛的季节变化因气候类型而异：

温带地区：夏季热岛最强，冬季最弱。夏季太阳辐射强，城市蓄热多，同时空调使用频繁，人为热排放大。
热带地区：热岛全年存在，雨季强度略低。雨季云量多，太阳辐射减弱，同时降水蒸发带走热量。
干旱地区：夜间热岛强，白天可能出现冷岛。白天城市蒸散发几乎为零，感热通量极大，但城市反照率差异和粗糙度效应可能导致城市升温较慢。

空间分布特征：

城市热岛的空间分布与土地利用类型密切相关：

高密度商业区：热岛强度最大，建筑密集、绿地稀少、人为热排放集中。
工业区：热岛强度大，工厂设备散热、高蓄热性地表、植被覆盖率低。
居住区：热岛强度中等，与建筑密度、绿化率、街道走向有关。
城市公园：热岛强度最小，甚至出现"冷岛"效应，植被蒸散发和遮阴作用显著。
水体：城市中的河流、湖泊具有显著的降温效应，是缓解热岛的重要资源。

1.3 热岛效应的影响

城市热岛对城市环境、居民健康和能源消耗产生深远影响：

对人体健康的影响：

热相关死亡率增加：高温加剧热应激，导致中暑、心脑血管疾病发作。研究表明，当夜间温度超过25°C时，死亡率显著上升，而城市热岛使夜间降温困难。
空气质量恶化：高温加速光化学反应，增加地面臭氧浓度。同时，热岛形成的城市热低压阻碍污染物扩散，加剧雾霾。
传染病传播：高温延长蚊虫活动季节，增加登革热、疟疾等虫媒传染病风险。

对能源消耗的影响：

空调负荷增加：城市温度每升高1°C，空调能耗增加约5-10%。夏季高峰期，热岛效应可使城市总电力负荷增加数吉瓦。
峰值负荷压力：热岛效应延长了高温时段，使电力峰值负荷持续时间更长，对电网造成压力。

对城市生态的影响：

物候改变：城市热岛使春季提前、秋季延后，改变植物开花、落叶时间，影响城市生态系统。
生物多样性下降：高温环境不适合某些物种生存，导致城市生物多样性降低。
水体生态恶化：城市径流温度升高，影响受纳水体的水生生态系统。

2. 城市地表辐射能量平衡

2.1 辐射能量平衡方程

城市地表的能量平衡是理解热岛效应的基础。对于任意地表单元，能量平衡方程为：

$Rn+QF=H+LE+G+ΔSR_n + Q_F = H + LE + G + \Delta S$

其中：

$R_n$ ：净辐射通量（W/m²）
$Q_F$ ：人为热通量（W/m²）
$H$ ：感热通量（W/m²）
$L E$ ：潜热通量（W/m²）， $L$ 为汽化潜热， $E$ 为蒸散发率
$G$ ：土壤/建筑储热通量（W/m²）
$ΔS\Delta S$ ：能量储存变化率（W/m²）

在稳态条件下， $ΔS≈0\Delta S \approx 0$ ，方程简化为：

$R_n + Q_F = H + LE + G$

净辐射通量 $R_n$ 由短波辐射和长波辐射组成：

$Rn=(1−α)K↓+L↓−L↑=(1−α)K↓+L∗R_n = (1-\alpha)K_\downarrow + L_\downarrow - L_\uparrow = (1-\alpha)K_\downarrow + L^*$

其中：

$K↓K_\downarrow$ ：到达地表的太阳短波辐射（W/m²）
$α\alpha$ ：地表反照率（0-1）
$L↓L_\downarrow$ ：大气逆辐射（W/m²）
$L↑L_\uparrow$ ：地表长波辐射发射（W/m²）
$L∗=L↓−L↑L^* = L_\downarrow - L_\uparrow$ ：净长波辐射

地表长波辐射发射遵循斯蒂芬-玻尔兹曼定律：

$L↑=εσTs4+(1−ε)L↓L_\uparrow = \varepsilon \sigma T_s^4 + (1-\varepsilon)L_\downarrow$

其中：

$ε\varepsilon$ ：地表发射率（0-1）
$σ\sigma$ ：斯蒂芬-玻尔兹曼常数（5.67×10⁻⁸ W/m²·K⁴）
$T_s$ ：地表温度（K）
$(1−ε)L↓(1-\varepsilon)L_\downarrow$ ：地表反射的长波辐射

因此，净长波辐射为：

$L∗=L↓−εσTs4−(1−ε)L↓=ε(L↓−σTs4)L^* = L_\downarrow - \varepsilon \sigma T_s^4 - (1-\varepsilon)L_\downarrow = \varepsilon(L_\downarrow - \sigma T_s^4)$

2.2 不同地表材料的辐射特性

城市热岛效应的根本原因在于城市下垫面材料与自然地表的辐射特性差异。以下是主要城市地表材料的辐射特性：

（1）沥青路面

沥青是城市道路的主要材料，其辐射特性表现为：

反照率：0.05-0.15（新沥青0.05-0.10，老化沥青0.10-0.15）
发射率：0.90-0.95
热容量：约1.9 MJ/m³·K
导热系数：约0.75 W/m·K
热扩散率：约0.4×10⁻⁶ m²/s

沥青路面反照率低，吸收大量太阳辐射；热容量大，储存热量多；导热性适中，热量向深层传递。这些特性使沥青路面成为城市热岛的重要贡献者。

（2）混凝土

混凝土广泛用于人行道、广场和建筑外墙：

反照率：0.20-0.35（普通混凝土0.20-0.30，白色混凝土0.50-0.70）
发射率：0.85-0.95
热容量：约1.7 MJ/m³·K
导热系数：约1.4 W/m·K
热扩散率：约0.8×10⁻⁶ m²/s

混凝土反照率高于沥青，但仍低于自然地表。其导热系数较高，热量容易向深层传递并储存。

（3）植被（草地、树木）

植被是缓解城市热岛的关键要素：

反照率：0.20-0.30（草地），0.15-0.25（树叶）
发射率：0.95-0.99
蒸散发率：高（白天可达70-90%的净辐射转化为潜热）
表面温度：比周围铺装低5-15°C

植被通过以下机制降温：

遮阴：树冠遮挡太阳辐射，减少地表得热
蒸散发：水分蒸发带走大量热量（潜热通量）
高发射率：有效发射长波辐射散热
低粗糙度：减少热量向大气的湍流输送

（4）水体

城市水体（湖泊、河流、喷泉）具有显著的降温效应：

反照率：0.05-0.15（取决于太阳高度角和水体浑浊度）
发射率：0.95-0.98
热容量：4.18 MJ/m³·K（远高于固体材料）
蒸散发率：高

水体的巨大热容量使其温度变化缓慢，白天吸收大量热量而升温有限，夜间缓慢释放热量，起到"热缓冲"作用。

（5）建筑屋顶和墙面

建筑表面是城市下垫面的重要组成部分：

屋顶反照率：0.10-0.35（深色屋顶0.10-0.20，浅色屋顶0.30-0.50）
墙面反照率：0.20-0.50
发射率：0.85-0.95

建筑表面之间多次反射（峡谷效应）使城市冠层内辐射 trapping 增强，加剧热岛效应。

2.3 城市冠层辐射模型

城市冠层是指地面到建筑平均高度的空间层。与平坦地表不同，城市冠层内辐射传输受建筑几何形状显著影响。城市冠层辐射模型需要考虑以下特殊过程：

（1）天空视角因子（Sky View Factor, SVF）

天空视角因子定义为地表某点可见天空的立体角与半球立体角（2π）之比：

$\frac{1}{\pi} \int_{0}^{2\pi} \int_{0}^{\theta_{max}(\phi)} \cos\theta \sin\theta \, d\theta \, d\phi$

对于无限长街道峡谷，天空视角因子可解析计算：

$\cos\beta - \sin\beta \cdot \arctan\left(\frac{W/H}{\sin\beta}\right) + \frac{H}{W} \ln\left(\frac{W/H + \sqrt{1+(W/H)^2}}{1+\sin\beta}\right)$

其中：

$W$ ：街道宽度
$H$ ：建筑高度
$β\beta$ ：街道走向与太阳方位角的夹角

简化情况下，对于东西走向的街道：

$\cos\left(\arctan\frac{H}{W/2}\right) = \frac{W/2}{\sqrt{H^2+(W/2)^2}}$

天空视角因子对热岛效应的影响：

白天：SVF小意味着建筑遮挡太阳辐射，减少地表得热，有利于降温
夜间：SVF小意味着地表向天空的长波辐射损失减少，不利于散热，加剧热岛

（2）建筑表面多次反射

在城市峡谷中，太阳辐射在墙面和地面之间多次反射，增加辐射 trapping。考虑一个简单的两墙面峡谷模型：

设峡谷高宽比为 $H / W$ ，墙面反照率为 $αw\alpha_w$ ，地面反照率为 $αg\alpha_g$ 。入射太阳辐射 $K↓K_\downarrow$ 在峡谷内的多次反射过程可用几何级数描述。

墙面接收的太阳辐射包括直接辐射和地面反射辐射：

$K_w = K_{w,direct} + K_{w,reflected}$

地面接收的太阳辐射包括直接辐射、两侧墙面反射辐射：

$Kg=Kg,direct+2⋅Kg,reflectedK_g = K_{g,direct} + 2 \cdot K_{g,reflected}$

经过多次反射后，峡谷有效反照率 $αeff\alpha_{eff}$ 低于材料本身反照率：

$αeff=αbase−Δαtrapping\alpha_{eff} = \alpha_{base} - \Delta\alpha_{trapping}$

其中 $Δαtrapping\Delta\alpha_{trapping}$ 随高宽比增大而增大。

（3）长波辐射交换

夜间，城市地表通过长波辐射冷却。在城市峡谷中，地表不仅向天空辐射，还向周围建筑墙面辐射。长波辐射能量平衡为：

$Lnet=εs(Lsky−σTs4)⋅SVF+∑iεsεiFsi(σTi4−σTs4)L_{net} = \varepsilon_s (L_{sky} - \sigma T_s^4) \cdot SVF + \sum_{i} \varepsilon_s \varepsilon_i F_{si} (\sigma T_i^4 - \sigma T_s^4)$

其中：

$L_{sky}$ ：大气逆辐射
$F_{si}$ ：地表与第 $i$ 个建筑表面的视角因子
$T_i$ ：第 $i$ 个建筑表面的温度

由于建筑表面温度通常高于天空等效温度，峡谷内地表的长波辐射净损失小于开阔地表，降温较慢。

3. 城市热岛辐射机制详解

3.1 白天热岛形成的辐射机制

白天，太阳短波辐射主导城市地表能量平衡。城市热岛的形成与以下辐射过程密切相关：

（1）反照率差异

城市下垫面平均反照率（0.10-0.20）明显低于郊区自然地表（草地0.20-0.30，森林0.15-0.25）。这意味着城市吸收更多的太阳辐射：

$ΔQabs=(1−αurban)K↓−(1−αrural)K↓=(αrural−αurban)K↓\Delta Q_{abs} = (1-\alpha_{urban})K_\downarrow - (1-\alpha_{rural})K_\downarrow = (\alpha_{rural} - \alpha_{urban})K_\downarrow$

以正午太阳辐射800 W/m²为例：

城市（反照率0.12）吸收：800 × (1-0.12) = 704 W/m²
郊区（反照率0.25）吸收：800 × (1-0.25) = 600 W/m²
差值：104 W/m²

这部分额外吸收的能量如果全部转化为感热，可使城市气温显著升高。

（2）蒸散发差异

郊区植被覆盖率高，蒸散发强烈，将大量净辐射转化为潜热：

$LErural=βrural⋅RnLE_{rural} = \beta_{rural} \cdot R_n$

其中 Bowen 比 $βrural=H/LE\beta_{rural} = H/LE$ 通常在0.2-0.5之间，意味着60-80%的净辐射用于蒸散发。

城市蒸散发受限（铺装不透水、植被稀少），Bowen 比高达2-5，意味着70-80%的净辐射转化为感热：

$Hurban=11+βurban−1Rn≈0.7−0.8RnH_{urban} = \frac{1}{1+\beta_{urban}^{-1}} R_n \approx 0.7-0.8 R_n$

（3）建筑几何效应

城市建筑几何形状对白天辐射过程的影响是双刃剑：

降温效应：

建筑遮挡直接太阳辐射，减少地表得热
墙面阴影降低地表温度
高宽比越大，遮阳效果越明显

增温效应：

峡谷效应增加辐射 trapping，降低有效反照率
墙面吸收辐射后通过对流加热空气
粗糙度增大降低通风效率

综合效应取决于高宽比和太阳高度角。夏季正午，高宽比大于1的街道可获得显著遮阳效果；冬季太阳高度角低，遮阳效应减弱。

（4）热容量效应

城市材料（混凝土、沥青）热容量大，白天吸收热量后升温较慢：

$ΔT=Qabs⋅Δtρcpd\Delta T = \frac{Q_{abs} \cdot \Delta t}{\rho c_p d}$

其中 $d$ 为热扩散深度。虽然城市吸收更多辐射，但由于热容量大，地表温度上升反而可能比郊区慢，导致白天热岛强度较弱甚至出现"冷岛"。

3.2 夜间热岛形成的辐射机制

夜间，长波辐射冷却主导地表能量平衡。城市热岛在夜间最为显著，其辐射机制包括：

（1）长波辐射冷却差异

郊区开阔地表直接向天空辐射冷却：

$Lnet,rural=εσ(Tsky4−Ts4)≈εσ(Ts4−Tsky4)L_{net,rural} = \varepsilon \sigma (T_{sky}^4 - T_s^4) \approx \varepsilon \sigma (T_s^4 - T_{sky}^4)$

城市峡谷内地表向天空的辐射受建筑遮挡：

$Lnet,urban=εσ(Tsky4−Ts4)⋅SVF+εσ∑Fsi(Ti4−Ts4)L_{net,urban} = \varepsilon \sigma (T_{sky}^4 - T_s^4) \cdot SVF + \varepsilon \sigma \sum F_{si}(T_i^4 - T_s^4)$

由于 $T_i > T_{sky}$ 且 $S V F < 1$ ，城市地表净长波辐射损失小于郊区：

$L_{net,urban}| < |L_{net,rural}|$

以典型值为例：

郊区SVF ≈ 1.0， $T_s$ = 25°C， $T_{sky}$ = 10°C， $L_{net}$ ≈ -85 W/m²
城市SVF = 0.3， $T_s$ = 28°C， $T_{wall}$ = 26°C， $L_{net}$ ≈ -45 W/m²

城市辐射冷却速率仅为郊区的50-60%，温度下降缓慢。

（2）储热释放

白天城市储存的大量热量在夜间缓慢释放：

$Grelease=−λ∂T∂zG_{release} = -\lambda \frac{\partial T}{\partial z}$

沥青路面白天储存的热量可在夜间持续释放数小时，维持地表高温。储热量估算：

$Qstored=ρcp∫0δ(T(z)−Tinitial)dzQ_{stored} = \rho c_p \int_0^{\delta} (T(z) - T_{initial}) dz$

其中 $δ\delta$ 为热扩散深度（沥青约0.1-0.3m，混凝土约0.2-0.5m）。

（3）人为热排放

夜间人为热排放是维持城市高温的重要因素：

$Q_F = Q_{traffic} + Q_{buildings} + Q_{industry}$

交通散热：车辆发动机、刹车、轮胎与路面摩擦
建筑散热：空调排热、照明、电器使用
工业散热：工厂设备、数据中心

典型城市夜间人为热通量：

商业区：50-100 W/m²
高密度居住区：20-50 W/m²
工业区：30-80 W/m²
郊区：< 5 W/m²

（4）大气逆辐射增强

城市大气中温室气体（CO₂、水汽）和气溶胶浓度高，大气逆辐射增强：

$L↓=εaσTa4L_\downarrow = \varepsilon_a \sigma T_a^4$

其中大气发射率 $εa\varepsilon_a$ 与温室气体浓度、云量、气溶胶有关。城市大气发射率比郊区高5-15%，进一步减少地表净长波辐射损失。

3.3 热岛强度的定量分析

热岛强度 $ΔTUHI\Delta T_{UHI}$ 定义为城市与郊区的气温差：

$ΔTUHI=Turban−Trural\Delta T_{UHI} = T_{urban} - T_{rural}$

基于能量平衡，可以建立热岛强度的简化模型。假设城乡净辐射差异主要由反照率和蒸散发差异引起，感热通量与温差成正比：

$H = h_c (T_s - T_a)$

其中 $h_c$ 为对流换热系数。稳态能量平衡给出：

$ΔTUHI=(1−αu)K↓−(1−αr)K↓+QF−(LEu−LEr)hc,u+hc,r\Delta T_{UHI} = \frac{(1-\alpha_u)K_\downarrow - (1-\alpha_r)K_\downarrow + Q_F - (LE_u - LE_r)}{h_{c,u} + h_{c,r}}$

简化形式：

$ΔTUHI∝Δα⋅K↓+QF−ΔLEhc\Delta T_{UHI} \propto \frac{\Delta \alpha \cdot K_\downarrow + Q_F - \Delta LE}{h_c}$

这表明热岛强度与反照率差异、人为热排放正相关，与蒸散发差异、对流效率负相关。

Oke（1982）基于观测数据提出了热岛强度的经验公式：

$ΔTUHI,max=k⋅ln⁡(P)−c\Delta T_{UHI,max} = k \cdot \ln(P) - c$

其中 $P$ 为城市人口， $k$ 和 $c$ 为经验常数。对于北美城市， $\approx 2.5$ ， $\approx 6.6$ 。

Park（1986）考虑了城市几何形状的影响：

$ΔTUHI=a⋅ln⁡(HW)+b⋅(1−SVF)+c\Delta T_{UHI} = a \cdot \ln(\frac{H}{W}) + b \cdot (1-SVF) + c$

其中 $H / W$ 为高宽比， $a$ 、 $b$ 、 $c$ 为与气候相关的常数。

4. 缓解城市热岛的辐射策略

4.1 高反照率材料（冷屋顶、冷路面）

提高城市地表反照率是缓解热岛最直接有效的辐射策略。

冷屋顶（Cool Roofs）：

冷屋顶使用高反照率材料（白色涂料、金属屋面、反射瓦片），反照率可达0.50-0.80，远高于传统深色屋顶（0.10-0.20）。

降温效果估算：

太阳辐射800 W/m²时，反照率从0.20提高到0.60，减少吸收320 W/m²
假设50%转化为感热，可减少160 W/m²的感热通量
对流换热系数10 W/m²·K时，可降低气温约16°C（屋顶表面）
建筑内部降温约2-5°C，空调能耗减少20-40%

冷屋顶的辐射特性要求：

高太阳反照率：> 0.50（ENERGY STAR标准）
高发射率：> 0.80（促进夜间散热）
耐久性：反照率保持率> 0.70（3年后）

冷路面（Cool Pavements）：

冷路面技术包括：

浅色骨料：使用浅色花岗岩、石英替代深色石灰岩
反射涂层：在现有路面上涂覆高反照率涂层
透水铺装：增加水分蒸发，降低表面温度
相变材料：利用材料相变潜热吸收热量

冷路面反照率可从0.05-0.10提高到0.30-0.50，表面温度降低10-20°C。但需注意：

高反照率可能增加行人眩光
冬季可能增加采暖能耗（但通常夏季收益大于冬季损失）
路面磨损后反照率下降，需要定期维护

4.2 绿化策略（绿色屋顶、城市森林）

植被通过辐射和非辐射机制缓解热岛：

绿色屋顶（Green Roofs）：

绿色屋顶在屋顶铺设土壤和植被，具有以下辐射特性：

反照率：0.15-0.30（低于冷屋顶，但通过蒸散热降温）
蒸散发：白天将60-80%净辐射转化为潜热
隔热：土壤层减少热量向建筑内部传递
发射率：0.95-0.99（有效长波辐射散热）

绿色屋顶表面温度比传统屋顶低15-30°C，建筑内部降温2-8°C。此外，绿色屋顶还具有雨水管理、空气净化、生物多样性保护等生态效益。

城市森林（Urban Forests）：

城市森林的降温机制：

遮阴：树冠遮挡太阳辐射，减少地表得热50-90%
蒸散发：单棵大树日蒸散发量可达100-400升，带走大量热量
长波辐射：树冠温度低于周围建筑，减少环境长波辐射
风场调节：树木引导气流，改善通风

单棵树的降温效果：

遮阴区域地表温度降低10-20°C
周围空气温度降低1-3°C（距离树冠10m范围内）
综合降温效应可达2-5°C（街区尺度）

4.3 水体与喷泉设计

水体通过高蒸发率和热容量缓解热岛：

城市湖泊与河流：

水体表面温度比周围铺装低5-15°C
蒸发带走大量热量（潜热通量200-500 W/m²）
巨大热容量缓冲温度波动
降温效应可延伸至下风向数百米

喷泉与喷雾系统：

水雾蒸发直接冷却空气（蒸发1kg水带走约2.4MJ热量）
增加空气湿度，提高人体舒适度
可用于广场、步行街等开放空间

4.4 城市形态优化

通过优化城市几何形状改善辐射环境：

街道走向与宽度：

东西走向街道夏季遮阳效果好，但冬季也减少太阳得热
南北走向街道冬季阳光充足，但夏季东西墙面受热严重
增加街道宽度提高天空视角因子，促进夜间散热
优化高宽比（H/W = 0.5-1.5）平衡遮阳与通风

建筑布局：

错列布局增加通风，带走热量
建筑高度渐变，引导气流
保留通风廊道，连接城市绿地和水体

表面材料组合：

人行道使用高反照率材料
车行道使用透水沥青或冷路面
建筑立面使用反射涂料或绿化墙

5. Python仿真案例分析

5.1 案例1：不同地表材料的辐射能量平衡对比

本案例模拟夏季典型日（夏至前后）不同地表材料的辐射能量平衡和温度变化。

模拟条件：

日期：6月21日（夏至）
纬度：30°N（中国华东地区）
天气：晴朗无云
大气条件：标准大气

地表材料：

沥青路面（反照率0.10，发射率0.95）
混凝土（反照率0.30，发射率0.90）
草地（反照率0.25，发射率0.98，蒸散发活跃）
高反照率路面（反照率0.50，发射率0.90）

模拟结果分析：

正午时刻（12:00）各材料表面温度：

沥青：58-62°C
混凝土：48-52°C
草地：32-36°C
高反照率路面：38-42°C

沥青与草地的温差可达25-30°C，这是城市热岛的重要来源。

能量平衡分量对比（正午）：

沥青：净辐射600 W/m²，感热500 W/m²，储热100 W/m²，潜热≈0
草地：净辐射450 W/m²，感热100 W/m²，储热50 W/m²，潜热300 W/m²

草地通过蒸散发将大部分净辐射转化为潜热，显著降低感热通量和表面温度。

夜间（02:00）各材料表面温度：

沥青：28-30°C
混凝土：25-27°C
草地：20-22°C
高反照率路面：22-24°C

沥青夜间降温缓慢，与草地的温差仍有8-10°C，维持夜间热岛。

5.2 案例2：城市峡谷辐射传输模拟

本案例模拟不同高宽比（H/W）的城市峡谷内的辐射传输过程。

模拟条件：

街道走向：东西向
高宽比：0.5, 1.0, 2.0, 3.0
墙面反照率：0.30
地面反照率：0.15（沥青）
太阳位置：夏至正午，太阳高度角83°

模拟结果分析：

天空视角因子（SVF）：

H/W = 0.5：SVF = 0.71
H/W = 1.0：SVF = 0.45
H/W = 2.0：SVF = 0.24
H/W = 3.0：SVF = 0.16

随着高宽比增加，天空视角因子迅速减小。

峡谷内辐射分布（正午）：

H/W = 0.5：地面接收太阳辐射约85%的开阔地表值
H/W = 1.0：地面接收太阳辐射约60%的开阔地表值
H/W = 2.0：地面接收太阳辐射约35%的开阔地表值
H/W = 3.0：地面接收太阳辐射约20%的开阔地表值

高宽比大于2的街道峡谷在正午可获得显著的遮阳效果，地面温度比开阔沥青路面低15-25°C。

峡谷有效反照率：

考虑多次反射后，峡谷有效反照率低于材料本身反照率：

H/W = 0.5：有效反照率0.13（材料0.15）
H/W = 1.0：有效反照率0.11
H/W = 2.0：有效反照率0.09
H/W = 3.0：有效反照率0.08

峡谷效应使城市冠层吸收更多太阳辐射。

夜间长波辐射冷却：

夜间峡谷内地表长波辐射净损失：

H/W = 0.5：约为开阔地表的70%
H/W = 1.0：约为开阔地表的50%
H/W = 2.0：约为开阔地表的35%
H/W = 3.0：约为开阔地表的25%

深峡谷严重阻碍夜间辐射散热，加剧夜间热岛。

5.3 案例3：城市热岛时空演变模拟

本案例模拟一个简化城市区域在典型夏季日的热岛时空演变。

城市模型：

5km × 5km 城市区域
中心商业区：高密度建筑，低绿化率
周边居住区：中密度建筑，中等绿化率
郊区：农田和林地

模拟结果分析：

热岛强度日变化：

06:00：热岛强度1-2°C（夜间热岛持续）
09:00：热岛强度降至0-1°C（郊区升温快）
12:00：热岛强度0°C或负值（郊区表面温度超过城市）
15:00：热岛强度1-2°C（城市储热释放）
18:00：热岛强度2-3°C
00:00：热岛强度达到峰值4-5°C

空间分布特征：

中心商业区热岛最强（4-6°C）
工业区次之（3-5°C）
居住区中等（2-4°C）
城市公园出现"冷岛"（-2至-4°C）
河流沿岸降温明显（-1至-3°C）

不同缓解措施效果：

冷屋顶普及50%：热岛强度降低0.8-1.2°C
城市绿化率提高20%：热岛强度降低1.0-1.5°C
冷路面普及30%：热岛强度降低0.5-0.8°C
综合措施：热岛强度降低2.5-3.5°C

5.4 案例4：缓解策略对比分析

本案例对比不同热岛缓解策略的辐射机制和降温效果。

对比策略：

基准情景：传统城市（沥青路面、深色屋顶、低绿化率）
高反照率策略：冷屋顶（反照率0.60）+ 冷路面（反照率0.40）
绿化策略：绿色屋顶（覆盖率50%）+ 城市森林（覆盖率30%）
综合策略：高反照率 + 绿化 + 水体

模拟结果分析：

白天表面温度（正午）：

基准情景：沥青58°C，屋顶55°C，平均气温35°C
高反照率策略：路面42°C，屋顶38°C，平均气温32°C
绿化策略：路面52°C（遮阴区38°C），屋顶40°C，平均气温30°C
综合策略：路面40°C，屋顶35°C，平均气温28°C

夜间气温（02:00）：

基准情景：平均气温28°C
高反照率策略：平均气温26°C（夜间长波散热改善）
绿化策略：平均气温24°C（储热少，蒸散发降温）
综合策略：平均气温22°C

能量平衡对比（日平均）：

策略	净辐射	感热	潜热	储热
基准	180	140	20	20
高反照率	120	90	15	15
绿化	160	60	80	20
综合	100	40	50	10

高反照率策略主要减少净辐射吸收；绿化策略主要增加潜热通量；综合策略两者兼顾，效果最佳。

成本效益分析：

高反照率策略：初期成本低，维护简单，但仅对白天有效
绿化策略：生态效益多，降温效果好，但需要维护和水资源
综合策略：初期成本高，但长期效益最大，适应性强

附录：Python仿真程序说明

本主题配套的Python仿真程序包含以下功能模块：

A.1 程序结构

run_simulation.py
├── 案例1：不同地表材料辐射能量平衡对比
│   ├── 太阳辐射计算（考虑太阳高度角、大气透过率）
│   ├── 地表能量平衡求解
│   └── 温度变化模拟（24小时）
├── 案例2：城市峡谷辐射传输模拟
│   ├── 天空视角因子计算
│   ├── 峡谷内太阳辐射分布
│   └── 多次反射模型
├── 案例3：城市热岛时空演变模拟
│   ├── 城市区域网格划分
│   ├── 不同土地利用类型参数
│   └── 热岛强度计算
└── 案例4：缓解策略对比分析
    ├── 多种缓解情景设置
    ├── 效果对比分析
    └── 成本效益评估

A.2 关键算法

太阳位置计算：

def solar_position(day_of_year, hour, latitude, longitude):
    """计算太阳高度角和方位角"""
    # 太阳赤纬角
    declination = 23.45 * np.sin(np.radians(360 * (284 + day_of_year) / 365))
    
    # 时角
    hour_angle = 15 * (hour - 12)
    
    # 太阳高度角
    sin_altitude = (np.sin(np.radians(latitude)) * np.sin(np.radians(declination)) +
                    np.cos(np.radians(latitude)) * np.cos(np.radians(declination)) * 
                    np.cos(np.radians(hour_angle)))
    altitude = np.degrees(np.arcsin(sin_altitude))
    
    return altitude

地表能量平衡求解：

def surface_energy_balance(R_n, Q_f, T_a, wind_speed, surface_type):
    """求解地表能量平衡方程"""
    # 根据地表类型获取参数
    params = surface_parameters[surface_type]
    
    # 迭代求解表面温度
    T_s = T_a  # 初始猜测
    for iteration in range(100):
        # 计算各能量分量
        H = sensible_heat(T_s, T_a, wind_speed, params)
        LE = latent_heat(T_s, T_a, wind_speed, params)
        G = ground_heat(T_s, params)
        
        # 能量平衡残差
        residual = R_n + Q_f - H - LE - G
        
        if abs(residual) < 0.1:
            break
            
        # 更新表面温度
        dT = residual / (dH_dT + dLE_dT + dG_dT)
        T_s += dT
    
    return T_s, H, LE, G

A.3 可视化输出

程序生成以下可视化结果：

能量平衡分量图：展示各材料24小时能量平衡分量变化
温度对比图：不同材料表面温度和气温对比
峡谷辐射分布图：不同高宽比峡谷内的辐射分布
热岛时空演变动画：城市热岛强度的时空演变
缓解策略对比图：不同策略的降温效果对比

A.4 使用说明

运行程序前请确保安装以下依赖：

pip install numpy matplotlib scipy pillow

运行程序：

python run_simulation.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
主题049：城市热岛效应辐射分析
城市热岛效应辐射换热仿真程序

本程序模拟城市热岛效应的辐射机制，包括：
1. 不同地表材料的辐射能量平衡对比
2. 城市峡谷辐射传输模拟
3. 城市热岛时空演变模拟
4. 缓解策略对比分析

作者：仿真模拟专家
日期：2026-03-06
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')  # 使用非交互式后端
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle, FancyBboxPatch
from matplotlib.collections import PatchCollection
import matplotlib.animation as animation
from PIL import Image
import os
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 创建输出目录
output_dir = os.path.dirname(os.path.abspath(__file__))

# ============================================================================
# 物理常数和参数
# ============================================================================
SIGMA = 5.67e-8  # 斯蒂芬-玻尔兹曼常数 W/m^2/K^4
CP_AIR = 1005  # 空气比热 J/kg/K
RHO_AIR = 1.225  # 空气密度 kg/m^3
L_V = 2.45e6  # 水的汽化潜热 J/kg

# 地表材料参数定义
SURFACE_MATERIALS = {
    'asphalt': {  # 沥青路面
        'albedo': 0.10,
        'emissivity': 0.95,
        'thermal_conductivity': 0.75,  # W/m/K
        'heat_capacity': 1.9e6,  # J/m^3/K
        'evaporation_rate': 0.0,  # 无蒸发
        'name': '沥青路面',
        'color': '#2C2C2C'
    },
    'concrete': {  # 混凝土
        'albedo': 0.30,
        'emissivity': 0.90,
        'thermal_conductivity': 1.4,
        'heat_capacity': 1.7e6,
        'evaporation_rate': 0.0,
        'name': '混凝土',
        'color': '#A0A0A0'
    },
    'grass': {  # 草地
        'albedo': 0.25,
        'emissivity': 0.98,
        'thermal_conductivity': 0.3,
        'heat_capacity': 2.0e6,
        'evaporation_rate': 0.7,  # 高蒸发率
        'name': '草地',
        'color': '#228B22'
    },
    'cool_pavement': {  # 高反照率路面
        'albedo': 0.50,
        'emissivity': 0.90,
        'thermal_conductivity': 1.0,
        'heat_capacity': 1.8e6,
        'evaporation_rate': 0.0,
        'name': '高反照率路面',
        'color': '#F5F5DC'
    },
    'water': {  # 水体
        'albedo': 0.08,
        'emissivity': 0.97,
        'thermal_conductivity': 0.6,
        'heat_capacity': 4.18e6,
        'evaporation_rate': 0.8,
        'name': '水体',
        'color': '#4169E1'
    }
}

# ============================================================================
# 太阳辐射计算
# ============================================================================
def solar_declination(day_of_year):
    """计算太阳赤纬角（度）"""
    return 23.45 * np.sin(np.radians(360 * (284 + day_of_year) / 365))

def solar_hour_angle(hour):
    """计算时角（度）"""
    return 15 * (hour - 12)

def solar_altitude(day_of_year, hour, latitude):
    """计算太阳高度角（度）"""
    dec = solar_declination(day_of_year)
    ha = solar_hour_angle(hour)
    
    sin_alt = (np.sin(np.radians(latitude)) * np.sin(np.radians(dec)) +
               np.cos(np.radians(latitude)) * np.cos(np.radians(dec)) * 
               np.cos(np.radians(ha)))
    
    sin_alt = np.clip(sin_alt, -1, 1)
    altitude = np.degrees(np.arcsin(sin_alt))
    
    return np.maximum(altitude, 0)  # 负值表示夜晚

def solar_azimuth(day_of_year, hour, latitude):
    """计算太阳方位角（度，从正北顺时针）"""
    dec = solar_declination(day_of_year)
    ha = solar_hour_angle(hour)
    alt = solar_altitude(day_of_year, hour, latitude)
    
    if alt <= 0:
        return 0
    
    cos_az = ((np.sin(np.radians(dec)) * np.cos(np.radians(latitude)) -
               np.cos(np.radians(dec)) * np.sin(np.radians(latitude)) * 
               np.cos(np.radians(ha))) / np.cos(np.radians(alt)))
    
    cos_az = np.clip(cos_az, -1, 1)
    azimuth = np.degrees(np.arccos(cos_az))
    
    if ha > 0:
        azimuth = 360 - azimuth
    
    return azimuth

def solar_radiation(day_of_year, hour, latitude, cloud_cover=0.0):
    """
    计算太阳辐射通量（W/m^2）
    
    参数:
    - day_of_year: 一年中的第几天（1-365）
    - hour: 小时（0-24）
    - latitude: 纬度（度）
    - cloud_cover: 云量（0-1）
    
    返回:
    - 直接辐射、散射辐射、总辐射（W/m^2）
    """
    altitude = solar_altitude(day_of_year, hour, latitude)
    
    if altitude <= 0:
        return 0, 0, 0
    
    # 太阳常数
    S0 = 1367  # W/m^2
    
    # 日地距离修正
    dr = 1 + 0.033 * np.cos(np.radians(360 * day_of_year / 365))
    
    # 大气透过率（简化模型）
    tau = 0.75  # 晴朗天气
    
    # 直接辐射
    Ib = S0 * dr * tau ** (1 / np.sin(np.radians(altitude)))
    
    # 散射辐射
    Id = 0.3 * (1 - tau) * S0 * dr * np.sin(np.radians(altitude))
    
    # 云量修正
    cloud_factor = 1 - 0.75 * cloud_cover
    Ib *= cloud_factor
    Id *= (0.5 + 0.5 * cloud_factor)
    
    # 总辐射
    Ig = Ib + Id
    
    return Ib, Id, Ig

def atmospheric_longwave(T_air, cloud_cover=0.0, humidity=0.6):
    """
    计算大气逆辐射（W/m^2）
    
    使用Brunt公式
    """
    # 饱和水汽压（hPa）
    es = 6.11 * np.exp(17.27 * T_air / (T_air + 237.3))
    
    # 实际水汽压
    ea = humidity * es
    
    # 大气发射率
    eps_a = 0.65 * (ea / (T_air + 273.15)) ** 0.14
    
    # 云量修正
    eps_a = eps_a * (1 + 0.26 * cloud_cover)
    
    # 大气逆辐射
    L_down = eps_a * SIGMA * (T_air + 273.15) ** 4
    
    return L_down

# ============================================================================
# 地表能量平衡计算
# ============================================================================
def convective_heat_transfer(T_surface, T_air, wind_speed):
    """
    计算对流换热系数和感热通量
    
    使用简化模型
    """
    # 对流换热系数（W/m^2/K）
    h_c = 5.7 + 3.8 * wind_speed
    
    # 感热通量
    H = h_c * (T_surface - T_air)
    
    return H, h_c

def latent_heat_flux(T_surface, T_air, wind_speed, evaporation_rate, RH=0.6):
    """
    计算潜热通量（W/m^2）
    """
    if evaporation_rate <= 0:
        return 0
    
    # 表面饱和水汽压
    es_surface = 6.11 * np.exp(17.27 * T_surface / (T_surface + 237.3))
    
    # 空气水汽压
    es_air = 6.11 * np.exp(17.27 * T_air / (T_air + 237.3))
    ea_air = RH * es_air
    
    # 水汽压差
    delta_e = es_surface - ea_air
    
    # 蒸散发系数
    h_e = 0.01 * (5.7 + 3.8 * wind_speed)  # s/m
    
    # 潜热通量
    LE = evaporation_rate * h_e * delta_e * 100 * L_V / 101.3
    
    return np.maximum(LE, 0)

def ground_heat_flux(T_surface, T_deep=20, thermal_conductivity=1.0, depth=0.5):
    """
    计算土壤/建筑储热通量（W/m^2）
    
    使用简化的一维热传导模型
    """
    G = thermal_conductivity * (T_surface - T_deep) / depth
    
    return G

def solve_surface_temperature(material, solar_rad, L_down, T_air, wind_speed, 
                               Q_f=0, T_deep=20, initial_T=None):
    """
    求解地表能量平衡，计算表面温度
    
    能量平衡方程: R_n + Q_f = H + LE + G
    
    参数:
    - material: 材料参数字典
    - solar_rad: 太阳总辐射（W/m^2）
    - L_down: 大气逆辐射（W/m^2）
    - T_air: 气温（°C）
    - wind_speed: 风速（m/s）
    - Q_f: 人为热通量（W/m^2）
    - T_deep: 深层温度（°C）
    - initial_T: 初始表面温度猜测（°C）
    
    返回:
    - T_surface: 表面温度（°C）
    - energy_components: 各能量分量字典
    """
    albedo = material['albedo']
    emissivity = material['emissivity']
    lambda_s = material['thermal_conductivity']
    evap_rate = material['evaporation_rate']
    
    # 初始猜测
    if initial_T is None:
        T_s = T_air + 5
    else:
        T_s = initial_T
    
    # 迭代求解
    for iteration in range(100):
        # 净短波辐射
        K_net = (1 - albedo) * solar_rad
        
        # 长波辐射
        L_up = emissivity * SIGMA * (T_s + 273.15) ** 4 + (1 - emissivity) * L_down
        L_net = L_down - L_up
        
        # 净辐射
        R_n = K_net + L_net
        
        # 感热通量
        H, h_c = convective_heat_transfer(T_s, T_air, wind_speed)
        
        # 潜热通量
        LE = latent_heat_flux(T_s, T_air, wind_speed, evap_rate)
        
        # 储热通量
        G = ground_heat_flux(T_s, T_deep, lambda_s)
        
        # 能量平衡残差
        residual = R_n + Q_f - H - LE - G
        
        # 收敛判断
        if abs(residual) < 0.5:
            break
        
        # 更新表面温度（牛顿迭代法）
        dR_n_dT = -4 * emissivity * SIGMA * (T_s + 273.15) ** 3
        dH_dT = h_c
        dLE_dT = evap_rate * 0.01 * h_c * L_V * 100 * 17.27 * 237.3 / (T_s + 237.3) ** 2 / 101.3
        dG_dT = lambda_s / 0.5
        
        dT = residual / (dR_n_dT - dH_dT - dLE_dT - dG_dT)
        T_s -= dT * 0.5  # 欠松弛
        
        # 限制温度范围
        T_s = np.clip(T_s, T_air - 10, T_air + 60)
    
    energy_components = {
        'R_n': R_n,
        'K_net': K_net,
        'L_net': L_net,
        'H': H,
        'LE': LE,
        'G': G,
        'Q_f': Q_f
    }
    
    return T_s, energy_components

# ============================================================================
# 案例1：不同地表材料的辐射能量平衡对比
# ============================================================================
def case1_surface_materials():
    """
    案例1：不同地表材料的辐射能量平衡对比
    模拟夏季典型日（夏至）不同地表材料的温度变化和能量平衡
    """
    print("\n" + "=" * 70)
    print("案例1：不同地表材料的辐射能量平衡对比")
    print("=" * 70)
    
    # 模拟参数
    day_of_year = 172  # 夏至（6月21日）
    latitude = 30  # 北纬30度（中国华东地区）
    hours = np.linspace(0, 24, 49)  # 每30分钟一个时间点
    
    # 气象条件
    T_air_day = 35  # 白天最高气温
    T_air_night = 25  # 夜间最低气温
    wind_speed = 2.0  # m/s
    cloud_cover = 0.0  # 晴朗
    humidity = 0.6  # 相对湿度60%
    
    # 材料列表
    materials = ['asphalt', 'concrete', 'grass', 'cool_pavement']
    
    # 存储结果
    results = {}
    
    for mat_name in materials:
        print(f"\n计算 {SURFACE_MATERIALS[mat_name]['name']}...")
        
        material = SURFACE_MATERIALS[mat_name]
        T_surface = np.zeros(len(hours))
        T_air = np.zeros(len(hours))
        energy_data = {'R_n': [], 'H': [], 'LE': [], 'G': [], 'K_net': [], 'L_net': []}
        
        T_s_prev = T_air_night
        
        for i, hour in enumerate(hours):
            # 气温日变化（正弦曲线）
            T_air[i] = T_air_night + (T_air_day - T_air_night) * \
                       np.maximum(0, np.sin(np.pi * (hour - 6) / 12))
            
            # 太阳辐射
            Ib, Id, Ig = solar_radiation(day_of_year, hour, latitude, cloud_cover)
            
            # 大气逆辐射
            L_down = atmospheric_longwave(T_air[i], cloud_cover, humidity)
            
            # 求解表面温度
            T_s, energy = solve_surface_temperature(
                material, Ig, L_down, T_air[i], wind_speed, 
                Q_f=0, T_deep=20, initial_T=T_s_prev
            )
            
            T_surface[i] = T_s
            T_s_prev = T_s
            
            for key in energy_data:
                energy_data[key].append(energy[key])
        
        results[mat_name] = {
            'T_surface': T_surface,
            'T_air': T_air,
            'energy': energy_data
        }
    
    # 绘制结果
    print("\n生成可视化图表...")
    
    # 图1：表面温度对比
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 表面温度时间序列
    ax1 = axes[0, 0]
    for mat_name in materials:
        mat_info = SURFACE_MATERIALS[mat_name]
        ax1.plot(hours, results[mat_name]['T_surface'], 
                label=mat_info['name'], color=mat_info['color'], linewidth=2)
    ax1.plot(hours, results['asphalt']['T_air'], 'k--', label='气温', linewidth=1.5, alpha=0.7)
    ax1.set_xlabel('时间 (小时)', fontsize=11)
    ax1.set_ylabel('温度 (°C)', fontsize=11)
    ax1.set_title('不同地表材料表面温度日变化', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper left', fontsize=9)
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(0, 24)
    
    # 正午能量平衡对比（12:00）
    ax2 = axes[0, 1]
    noon_idx = 24  # 12:00
    
    energy_types = ['K_net', 'L_net', 'R_n', 'H', 'LE', 'G']
    energy_labels = ['净短波', '净长波', '净辐射', '感热', '潜热', '储热']
    x_pos = np.arange(len(materials))
    width = 0.12
    
    colors_bar = ['#FFD700', '#87CEEB', '#FF6B6B', '#FF8C00', '#4169E1', '#8B4513']
    
    for j, (e_type, e_label, color) in enumerate(zip(energy_types, energy_labels, colors_bar)):
        values = [results[mat]['energy'][e_type][noon_idx] for mat in materials]
        ax2.bar(x_pos + j * width, values, width, label=e_label, color=color, alpha=0.8)
    
    ax2.set_xlabel('地表材料', fontsize=11)
    ax2.set_ylabel('热流密度 (W/m²)', fontsize=11)
    ax2.set_title('正午(12:00)能量平衡分量对比', fontsize=12, fontweight='bold')
    ax2.set_xticks(x_pos + width * 2.5)
    ax2.set_xticklabels([SURFACE_MATERIALS[m]['name'] for m in materials], fontsize=9)
    ax2.legend(loc='upper right', fontsize=8, ncol=2)
    ax2.grid(True, alpha=0.3, axis='y')
    ax2.axhline(y=0, color='k', linewidth=0.5)
    
    # 净辐射日变化
    ax3 = axes[1, 0]
    for mat_name in materials:
        mat_info = SURFACE_MATERIALS[mat_name]
        ax3.plot(hours, results[mat_name]['energy']['R_n'], 
                label=mat_info['name'], color=mat_info['color'], linewidth=2)
    ax3.set_xlabel('时间 (小时)', fontsize=11)
    ax3.set_ylabel('净辐射 (W/m²)', fontsize=11)
    ax3.set_title('净辐射日变化', fontsize=12, fontweight='bold')
    ax3.legend(loc='upper left', fontsize=9)
    ax3.grid(True, alpha=0.3)
    ax3.set_xlim(0, 24)
    ax3.axhline(y=0, color='k', linewidth=0.5)
    
    # 感热与潜热对比
    ax4 = axes[1, 1]
    
    # 计算日平均
    H_avg = [np.mean(results[mat]['energy']['H']) for mat in materials]
    LE_avg = [np.mean(results[mat]['energy']['LE']) for mat in materials]
    
    x_pos = np.arange(len(materials))
    width = 0.35
    
    bars1 = ax4.bar(x_pos - width/2, H_avg, width, label='感热通量', color='#FF8C00', alpha=0.8)
    bars2 = ax4.bar(x_pos + width/2, LE_avg, width, label='潜热通量', color='#4169E1', alpha=0.8)
    
    ax4.set_xlabel('地表材料', fontsize=11)
    ax4.set_ylabel('日平均热流密度 (W/m²)', fontsize=11)
    ax4.set_title('感热与潜热通量对比（日平均）', fontsize=12, fontweight='bold')
    ax4.set_xticks(x_pos)
    ax4.set_xticklabels([SURFACE_MATERIALS[m]['name'] for m in materials], fontsize=9)
    ax4.legend(loc='upper right', fontsize=10)
    ax4.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case1_surface_materials.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"图表已保存: {output_dir}/case1_surface_materials.png")
    
    # 打印关键结果
    print("\n关键结果汇总:")
    print("-" * 60)
    print(f"{'材料':<15} {'最高温度(°C)':<12} {'日平均感热':<12} {'日平均潜热':<12}")
    print("-" * 60)
    for mat_name in materials:
        T_max = np.max(results[mat_name]['T_surface'])
        H_avg = np.mean(results[mat_name]['energy']['H'])
        LE_avg = np.mean(results[mat_name]['energy']['LE'])
        print(f"{SURFACE_MATERIALS[mat_name]['name']:<12} {T_max:>10.1f} {H_avg:>12.1f} {LE_avg:>12.1f}")
    print("-" * 60)
    
    return results

# ============================================================================
# 案例2：城市峡谷辐射传输模拟
# ============================================================================
def sky_view_factor(H, W):
    """
    计算街道峡谷的天空视角因子
    
    参数:
    - H: 建筑高度 (m)
    - W: 街道宽度 (m)
    
    返回:
    - SVF: 天空视角因子 (0-1)
    """
    # 简化公式：无限长街道峡谷
    aspect_ratio = H / W
    
    # 使用近似公式
    SVF = np.cos(np.arctan(2 * aspect_ratio))
    
    return SVF

def canyon_solar_radiation(H, W, solar_alt, solar_az, street_az=90):
    """
    计算峡谷内太阳辐射分布
    
    参数:
    - H: 建筑高度 (m)
    - W: 街道宽度 (m)
    - solar_alt: 太阳高度角 (度)
    - solar_az: 太阳方位角 (度)
    - street_az: 街道走向方位角 (度，默认90度即东西向)
    
    返回:
    - 地面接收的太阳辐射比例
    - 墙面接收的太阳辐射比例
    """
    if solar_alt <= 0:
        return 0, 0, 0
    
    # 太阳相对于街道的入射角
    relative_az = np.abs(solar_az - street_az)
    if relative_az > 180:
        relative_az = 360 - relative_az
    
    # 计算墙面投影
    # 简化模型：考虑两侧墙面遮挡
    tan_alt = np.tan(np.radians(solar_alt))
    
    # 地面接收的直接辐射比例（考虑墙面阴影）
    if relative_az < 90:  # 太阳大致垂直于街道
        shadow_length = H / tan_alt * np.cos(np.radians(relative_az))
        if shadow_length >= W:
            ground_ratio = 0  # 完全阴影
        else:
            ground_ratio = 1 - shadow_length / W
    else:  # 太阳平行于街道
        ground_ratio = 1.0
    
    # 考虑散射辐射（天空视角因子）
    SVF = sky_view_factor(H, W)
    diffuse_ratio = SVF
    
    # 总辐射
    total_ground = ground_ratio * 0.7 + diffuse_ratio * 0.3
    
    # 墙面接收的辐射
    wall_ratio = (1 - SVF) * 0.5 * np.sin(np.radians(solar_alt))
    
    return total_ground, wall_ratio, SVF

def effective_albedo(H, W, wall_albedo, ground_albedo):
    """
    计算峡谷有效反照率（考虑多次反射）
    
    使用简化模型
    """
    SVF = sky_view_factor(H, W)
    
    # 峡谷内表面平均反照率
    # 地面占1份，两侧墙面各占H/W份
    total_area = 1 + 2 * H / W
    avg_albedo = (ground_albedo + 2 * (H / W) * wall_albedo) / total_area
    
    # 多次反射因子
    trapping_factor = (1 - SVF) * avg_albedo
    
    # 有效反照率
    alpha_eff = SVF * ground_albedo + (1 - SVF) * avg_albedo * trapping_factor
    
    return alpha_eff

def case2_urban_canyon():
    """
    案例2：城市峡谷辐射传输模拟
    分析不同高宽比峡谷的辐射特性
    """
    print("\n" + "=" * 70)
    print("案例2：城市峡谷辐射传输模拟")
    print("=" * 70)
    
    # 模拟参数
    day_of_year = 172  # 夏至
    latitude = 30
    street_az = 90  # 东西走向
    wall_albedo = 0.30
    ground_albedo = 0.15  # 沥青
    
    # 不同高宽比
    aspect_ratios = [0.5, 1.0, 2.0, 3.0]
    W = 20  # 街道宽度固定20m
    
    hours = np.linspace(6, 18, 25)  # 白天6:00-18:00
    
    results = {}
    
    for aspect_ratio in aspect_ratios:
        H = aspect_ratio * W
        print(f"\n计算高宽比 H/W = {aspect_ratio} (H={H}m, W={W}m)...")
        
        SVF = sky_view_factor(H, W)
        alpha_eff = effective_albedo(H, W, wall_albedo, ground_albedo)
        
        ground_rad = []
        wall_rad = []
        total_solar = []
        
        for hour in hours:
            alt = solar_altitude(day_of_year, hour, latitude)
            az = solar_azimuth(day_of_year, hour, latitude)
            
            Ib, Id, Ig = solar_radiation(day_of_year, hour, latitude)
            
            g_ratio, w_ratio, _ = canyon_solar_radiation(H, W, alt, az, street_az)
            
            ground_rad.append(g_ratio * Ig)
            wall_rad.append(w_ratio * Ig)
            total_solar.append(Ig)
        
        results[aspect_ratio] = {
            'H': H,
            'SVF': SVF,
            'alpha_eff': alpha_eff,
            'ground_rad': np.array(ground_rad),
            'wall_rad': np.array(wall_rad),
            'total_solar': np.array(total_solar)
        }
    
    # 绘制结果
    print("\n生成可视化图表...")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：峡谷内地面接收辐射
    ax1 = axes[0, 0]
    colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
    for i, aspect_ratio in enumerate(aspect_ratios):
        ax1.plot(hours, results[aspect_ratio]['ground_rad'], 
                label=f'H/W = {aspect_ratio}', color=colors[i], linewidth=2)
    ax1.plot(hours, results[0.5]['total_solar'], 'k--', label='开阔地表', 
            linewidth=1.5, alpha=0.5)
    ax1.set_xlabel('时间 (小时)', fontsize=11)
    ax1.set_ylabel('太阳辐射 (W/m²)', fontsize=11)
    ax1.set_title('峡谷地面接收太阳辐射', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right', fontsize=9)
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(6, 18)
    
    # 图2：天空视角因子和有效反照率
    ax2 = axes[0, 1]
    x_pos = np.arange(len(aspect_ratios))
    width = 0.35
    
    SVFs = [results[ar]['SVF'] for ar in aspect_ratios]
    alphas = [results[ar]['alpha_eff'] for ar in aspect_ratios]
    
    ax2_twin = ax2.twinx()
    
    bars1 = ax2.bar(x_pos - width/2, SVFs, width, label='天空视角因子', 
                    color='#4169E1', alpha=0.8)
    bars2 = ax2_twin.bar(x_pos + width/2, alphas, width, label='有效反照率', 
                         color='#FF8C00', alpha=0.8)
    
    ax2.set_xlabel('高宽比 (H/W)', fontsize=11)
    ax2.set_ylabel('天空视角因子', fontsize=11, color='#4169E1')
    ax2_twin.set_ylabel('有效反照率', fontsize=11, color='#FF8C00')
    ax2.set_title('天空视角因子与有效反照率', fontsize=12, fontweight='bold')
    ax2.set_xticks(x_pos)
    ax2.set_xticklabels([f'{ar}' for ar in aspect_ratios])
    ax2.tick_params(axis='y', labelcolor='#4169E1')
    ax2_twin.tick_params(axis='y', labelcolor='#FF8C00')
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 添加图例
    lines1, labels1 = ax2.get_legend_handles_labels()
    lines2, labels2 = ax2_twin.get_legend_handles_labels()
    ax2.legend(lines1 + lines2, labels1 + labels2, loc='upper right', fontsize=9)
    
    # 图3：不同高宽比峡谷示意图
    ax3 = axes[1, 0]
    ax3.set_xlim(0, 100)
    ax3.set_ylim(0, 80)
    ax3.set_aspect('equal')
    ax3.axis('off')
    ax3.set_title('不同高宽比峡谷示意图', fontsize=12, fontweight='bold')
    
    y_offset = 0
    for i, aspect_ratio in enumerate(aspect_ratios):
        H = results[aspect_ratio]['H']
        W = 20
        
        # 绘制峡谷
        # 地面
        rect_ground = Rectangle((10, y_offset), W, 2, 
                                facecolor='#555555', edgecolor='black')
        ax3.add_patch(rect_ground)
        
        # 左侧建筑
        rect_left = Rectangle((10-W/4, y_offset+2), W/4, H/3, 
                              facecolor='#8B4513', edgecolor='black')
        ax3.add_patch(rect_left)
        
        # 右侧建筑
        rect_right = Rectangle((10+W, y_offset+2), W/4, H/3, 
                               facecolor='#8B4513', edgecolor='black')
        ax3.add_patch(rect_right)
        
        # 标注
        ax3.text(5, y_offset + H/6, f'H/W={aspect_ratio}', 
                fontsize=10, va='center', fontweight='bold')
        ax3.text(20+W/2, y_offset + 1, f'W={W}m', 
                fontsize=8, ha='center', va='center', color='white')
        ax3.text(10-W/8, y_offset + H/6 + 2, f'H={H:.0f}m', 
                fontsize=8, ha='center', va='bottom', rotation=90)
        
        y_offset += 18
    
    # 图4：长波辐射冷却对比（夜间）
    ax4 = axes[1, 1]
    
    # 夜间长波辐射冷却
    T_surface = 25  # 地表温度
    T_sky = 10  # 天空等效温度
    T_wall = 23  # 墙面温度
    
    cooling_rates = []
    for aspect_ratio in aspect_ratios:
        SVF = results[aspect_ratio]['SVF']
        
        # 开阔地表
        L_net_open = 0.95 * SIGMA * ((T_sky + 273.15)**4 - (T_surface + 273.15)**4)
        
        # 峡谷内
        L_net_canyon = SVF * 0.95 * SIGMA * ((T_sky + 273.15)**4 - (T_surface + 273.15)**4) + \
                       (1 - SVF) * 0.95 * 0.90 * SIGMA * ((T_wall + 273.15)**4 - (T_surface + 273.15)**4)
        
        cooling_rates.append(abs(L_net_canyon) / abs(L_net_open) * 100)
    
    bars = ax4.bar(range(len(aspect_ratios)), cooling_rates, 
                   color=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728'], alpha=0.8)
    ax4.axhline(y=100, color='k', linestyle='--', linewidth=1.5, label='开阔地表基准')
    ax4.set_xlabel('高宽比 (H/W)', fontsize=11)
    ax4.set_ylabel('相对辐射冷却率 (%)', fontsize=11)
    ax4.set_title('夜间长波辐射冷却对比', fontsize=12, fontweight='bold')
    ax4.set_xticks(range(len(aspect_ratios)))
    ax4.set_xticklabels([f'{ar}' for ar in aspect_ratios])
    ax4.legend(loc='upper right', fontsize=9)
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for i, (bar, val) in enumerate(zip(bars, cooling_rates)):
        ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2, 
                f'{val:.0f}%', ha='center', va='bottom', fontsize=10, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case2_urban_canyon.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"图表已保存: {output_dir}/case2_urban_canyon.png")
    
    # 打印关键结果
    print("\n关键结果汇总:")
    print("-" * 70)
    print(f"{'高宽比':<10} {'建筑高度(m)':<12} {'天空视角因子':<12} {'有效反照率':<12}")
    print("-" * 70)
    for aspect_ratio in aspect_ratios:
        H = results[aspect_ratio]['H']
        SVF = results[aspect_ratio]['SVF']
        alpha = results[aspect_ratio]['alpha_eff']
        print(f"{aspect_ratio:<10.1f} {H:<12.1f} {SVF:<12.3f} {alpha:<12.3f}")
    print("-" * 70)
    
    return results

# ============================================================================
# 案例3：城市热岛时空演变模拟
# ============================================================================
def case3_urban_heat_island():
    """
    案例3：城市热岛时空演变模拟
    模拟简化城市区域的热岛时空演变
    """
    print("\n" + "=" * 70)
    print("案例3：城市热岛时空演变模拟")
    print("=" * 70)
    
    # 城市网格设置
    nx, ny = 50, 50
    dx = 100  # 网格间距 100m
    
    # 创建城市土地利用类型网格
    # 0: 郊区, 1: 居住区, 2: 商业区, 3: 工业区, 4: 公园, 5: 水体
    land_use = np.zeros((ny, nx), dtype=int)
    
    # 设置城市区域
    center_x, center_y = nx // 2, ny // 2
    
    for j in range(ny):
        for i in range(nx):
            dist = np.sqrt((i - center_x)**2 + (j - center_y)**2)
            
            if dist < 5:  # 中心商业区
                land_use[j, i] = 2
            elif dist < 12:  # 居住区
                land_use[j, i] = 1
            elif dist < 18:  # 工业区
                land_use[j, i] = 3
            elif i > nx - 8:  # 东侧水体（河流）
                land_use[j, i] = 5
            elif dist > 25 and np.random.random() < 0.3:  # 郊区公园
                land_use[j, i] = 4
    
    # 添加一些公园
    parks = [(15, 15), (35, 35), (10, 40)]
    for px, py in parks:
        for j in range(max(0, py-3), min(ny, py+4)):
            for i in range(max(0, px-3), min(nx, px+4)):
                if np.sqrt((i-px)**2 + (j-py)**2) < 3:
                    land_use[j, i] = 4
    
    # 土地利用参数
    land_params = {
        0: {'name': '郊区', 'albedo': 0.25, 'emissivity': 0.95, 'evap': 0.6, 'Q_f': 5, 'color': '#90EE90'},
        1: {'name': '居住区', 'albedo': 0.20, 'emissivity': 0.92, 'evap': 0.3, 'Q_f': 25, 'color': '#FFD700'},
        2: {'name': '商业区', 'albedo': 0.15, 'emissivity': 0.90, 'evap': 0.1, 'Q_f': 80, 'color': '#FF6347'},
        3: {'name': '工业区', 'albedo': 0.18, 'emissivity': 0.88, 'evap': 0.05, 'Q_f': 60, 'color': '#A9A9A9'},
        4: {'name': '公园', 'albedo': 0.22, 'emissivity': 0.97, 'evap': 0.8, 'Q_f': 0, 'color': '#228B22'},
        5: {'name': '水体', 'albedo': 0.08, 'emissivity': 0.97, 'evap': 0.9, 'Q_f': 0, 'color': '#4169E1'}
    }
    
    # 模拟参数
    day_of_year = 172
    latitude = 30
    hours = np.linspace(0, 24, 25)
    
    # 存储结果
    T_surface = np.zeros((len(hours), ny, nx))
    T_air = np.zeros((len(hours), ny, nx))
    UHI_intensity = np.zeros(len(hours))
    
    print("\n计算城市热岛时空演变...")
    
    for t, hour in enumerate(hours):
        print(f"  计算时间: {hour:.1f}h")
        
        # 气象条件
        T_air_base = 25 + 10 * np.maximum(0, np.sin(np.pi * (hour - 6) / 12))
        wind_speed = 2.0
        
        # 太阳辐射
        Ib, Id, Ig = solar_radiation(day_of_year, hour, latitude)
        L_down = atmospheric_longwave(T_air_base)
        
        # 计算每个网格的表面温度
        for j in range(ny):
            for i in range(nx):
                lu = land_use[j, i]
                params = {
                    'albedo': land_params[lu]['albedo'],
                    'emissivity': land_params[lu]['emissivity'],
                    'thermal_conductivity': 1.0,
                    'heat_capacity': 2.0e6,
                    'evaporation_rate': land_params[lu]['evap']
                }
                
                Q_f = land_params[lu]['Q_f']
                
                T_s, _ = solve_surface_temperature(
                    params, Ig, L_down, T_air_base, wind_speed, Q_f=Q_f
                )
                
                T_surface[t, j, i] = T_s
                T_air[t, j, i] = T_air_base + 0.3 * (T_s - T_air_base)  # 简化气温计算
        
        # 计算热岛强度（城市与郊区平均温差）
        urban_mask = (land_use == 1) | (land_use == 2) | (land_use == 3)
        rural_mask = land_use == 0
        
        if np.any(urban_mask) and np.any(rural_mask):
            UHI_intensity[t] = np.mean(T_air[t][urban_mask]) - np.mean(T_air[t][rural_mask])
    
    # 绘制结果
    print("\n生成可视化图表和动画...")
    
    # 创建静态图
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1：土地利用分布
    ax1 = axes[0, 0]
    colors_land = [land_params[i]['color'] for i in range(6)]
    cmap_land = matplotlib.colors.ListedColormap(colors_land)
    im1 = ax1.imshow(land_use, cmap=cmap_land, origin='lower', vmin=-0.5, vmax=5.5)
    ax1.set_title('城市土地利用分布', fontsize=12, fontweight='bold')
    ax1.set_xlabel('X (100m)', fontsize=10)
    ax1.set_ylabel('Y (100m)', fontsize=10)
    
    # 添加图例
    from matplotlib.patches import Patch
    legend_elements = [Patch(facecolor=land_params[i]['color'], 
                            label=land_params[i]['name']) for i in range(6)]
    ax1.legend(handles=legend_elements, loc='upper right', fontsize=8)
    
    # 图2：正午表面温度分布
    ax2 = axes[0, 1]
    noon_idx = 12  # 12:00
    im2 = ax2.imshow(T_surface[noon_idx], cmap='hot', origin='lower', 
                     vmin=20, vmax=70)
    ax2.set_title('正午(12:00)表面温度分布', fontsize=12, fontweight='bold')
    ax2.set_xlabel('X (100m)', fontsize=10)
    ax2.set_ylabel('Y (100m)', fontsize=10)
    plt.colorbar(im2, ax=ax2, label='温度 (°C)')
    
    # 图3：夜间表面温度分布
    ax3 = axes[1, 0]
    night_idx = 22  # 22:00
    im3 = ax3.imshow(T_surface[night_idx], cmap='hot', origin='lower', 
                     vmin=15, vmax=35)
    ax3.set_title('夜间(22:00)表面温度分布', fontsize=12, fontweight='bold')
    ax3.set_xlabel('X (100m)', fontsize=10)
    ax3.set_ylabel('Y (100m)', fontsize=10)
    plt.colorbar(im3, ax=ax3, label='温度 (°C)')
    
    # 图4：热岛强度日变化
    ax4 = axes[1, 1]
    ax4.plot(hours, UHI_intensity, 'b-', linewidth=2, marker='o', markersize=4)
    ax4.axhline(y=0, color='k', linestyle='--', linewidth=1, alpha=0.5)
    ax4.fill_between(hours, 0, UHI_intensity, where=(UHI_intensity > 0), 
                     alpha=0.3, color='red', label='热岛')
    ax4.fill_between(hours, 0, UHI_intensity, where=(UHI_intensity < 0), 
                     alpha=0.3, color='blue', label='冷岛')
    ax4.set_xlabel('时间 (小时)', fontsize=11)
    ax4.set_ylabel('热岛强度 (°C)', fontsize=11)
    ax4.set_title('城市热岛强度日变化', fontsize=12, fontweight='bold')
    ax4.legend(loc='upper left', fontsize=9)
    ax4.grid(True, alpha=0.3)
    ax4.set_xlim(0, 24)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case3_urban_heat_island.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"静态图表已保存: {output_dir}/case3_urban_heat_island.png")
    
    # 创建动画
    print("  创建热岛演变动画...")
    
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
    
    def update(frame):
        ax1.clear()
        ax2.clear()
        
        # 表面温度图
        im = ax1.imshow(T_surface[frame], cmap='hot', origin='lower', 
                        vmin=15, vmax=65)
        ax1.set_title(f'表面温度分布 ({hours[frame]:.0f}:00)', 
                     fontsize=12, fontweight='bold')
        ax1.set_xlabel('X (100m)', fontsize=10)
        ax1.set_ylabel('Y (100m)', fontsize=10)
        
        # 热岛强度曲线
        ax2.plot(hours[:frame+1], UHI_intensity[:frame+1], 'b-', linewidth=2)
        ax2.axhline(y=0, color='k', linestyle='--', linewidth=1, alpha=0.5)
        ax2.set_xlabel('时间 (小时)', fontsize=11)
        ax2.set_ylabel('热岛强度 (°C)', fontsize=11)
        ax2.set_title('热岛强度日变化', fontsize=12, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        ax2.set_xlim(0, 24)
        ax2.set_ylim(-2, 6)
        
        # 添加当前时间标记
        ax2.axvline(x=hours[frame], color='r', linestyle='--', alpha=0.5)
        
        return im,
    
    ani = animation.FuncAnimation(fig, update, frames=len(hours), 
                                  interval=200, blit=False)
    
    # 保存为GIF
    ani.save(f'{output_dir}/case3_uhi_animation.gif', writer='pillow', fps=5)
    plt.close()
    
    print(f"动画已保存: {output_dir}/case3_uhi_animation.gif")
    
    # 打印关键结果
    print("\n关键结果汇总:")
    print("-" * 60)
    print(f"最大热岛强度: {np.max(UHI_intensity):.2f}°C (发生在 {hours[np.argmax(UHI_intensity)]:.0f}:00)")
    print(f"最小热岛强度: {np.min(UHI_intensity):.2f}°C (发生在 {hours[np.argmin(UHI_intensity)]:.0f}:00)")
    print(f"日平均热岛强度: {np.mean(UHI_intensity):.2f}°C")
    print("-" * 60)
    
    return {
        'land_use': land_use,
        'T_surface': T_surface,
        'T_air': T_air,
        'UHI_intensity': UHI_intensity,
        'hours': hours
    }

# ============================================================================
# 案例4：缓解策略对比分析
# ============================================================================
def case4_mitigation_strategies():
    """
    案例4：缓解策略对比分析
    对比不同热岛缓解策略的效果
    """
    print("\n" + "=" * 70)
    print("案例4：缓解策略对比分析")
    print("=" * 70)
    
    # 基准城市参数
    base_params = {
        'albedo': 0.15,  # 平均反照率
        'emissivity': 0.92,
        'evap_rate': 0.2,
        'Q_f': 40,  # 平均人为热
        'green_coverage': 0.15  # 绿化覆盖率
    }
    
    # 缓解策略
    strategies = {
        'baseline': {
            'name': '基准情景',
            'albedo': 0.15,
            'evap_rate': 0.2,
            'Q_f': 40,
            'green_coverage': 0.15
        },
        'high_albedo': {
            'name': '高反照率策略',
            'albedo': 0.35,  # 冷屋顶+冷路面
            'evap_rate': 0.2,
            'Q_f': 40,
            'green_coverage': 0.15
        },
        'greening': {
            'name': '绿化策略',
            'albedo': 0.18,
            'evap_rate': 0.5,  # 增加绿化
            'Q_f': 35,  # 绿化降低部分人为热
            'green_coverage': 0.40
        },
        'combined': {
            'name': '综合策略',
            'albedo': 0.30,
            'evap_rate': 0.45,
            'Q_f': 30,
            'green_coverage': 0.35
        }
    }
    
    # 模拟参数
    day_of_year = 172
    latitude = 30
    hours = np.linspace(0, 24, 49)
    
    results = {}
    
    print("\n计算不同缓解策略...")
    
    for strat_name, strat in strategies.items():
        print(f"  计算: {strat['name']}")
        
        T_surface = np.zeros(len(hours))
        T_air = np.zeros(len(hours))
        energy_balance = {'R_n': [], 'H': [], 'LE': [], 'G': []}
        
        params = {
            'albedo': strat['albedo'],
            'emissivity': 0.92,
            'thermal_conductivity': 1.0,
            'heat_capacity': 2.0e6,
            'evaporation_rate': strat['evap_rate']
        }
        
        T_s_prev = 25
        
        for i, hour in enumerate(hours):
            T_air_base = 25 + 10 * np.maximum(0, np.sin(np.pi * (hour - 6) / 12))
            wind_speed = 2.0
            
            Ib, Id, Ig = solar_radiation(day_of_year, hour, latitude)
            L_down = atmospheric_longwave(T_air_base)
            
            T_s, energy = solve_surface_temperature(
                params, Ig, L_down, T_air_base, wind_speed, 
                Q_f=strat['Q_f'], initial_T=T_s_prev
            )
            
            T_surface[i] = T_s
            T_air[i] = T_air_base + 0.3 * (T_s - T_air_base)
            T_s_prev = T_s
            
            for key in energy_balance:
                energy_balance[key].append(energy[key])
        
        results[strat_name] = {
            'T_surface': T_surface,
            'T_air': T_air,
            'energy': energy_balance,
            'params': strat
        }
    
    # 绘制结果
    print("\n生成可视化图表...")
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：表面温度对比
    ax1 = axes[0, 0]
    colors = ['#d62728', '#ff7f0e', '#2ca02c', '#1f77b4']
    for i, (strat_name, strat) in enumerate(strategies.items()):
        ax1.plot(hours, results[strat_name]['T_surface'], 
                label=strat['name'], color=colors[i], linewidth=2)
    ax1.set_xlabel('时间 (小时)', fontsize=11)
    ax1.set_ylabel('表面温度 (°C)', fontsize=11)
    ax1.set_title('不同策略表面温度对比', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper left', fontsize=9)
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(0, 24)
    
    # 图2：气温对比
    ax2 = axes[0, 1]
    for i, (strat_name, strat) in enumerate(strategies.items()):
        ax2.plot(hours, results[strat_name]['T_air'], 
                label=strat['name'], color=colors[i], linewidth=2)
    ax2.set_xlabel('时间 (小时)', fontsize=11)
    ax2.set_ylabel('气温 (°C)', fontsize=11)
    ax2.set_title('不同策略气温对比', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper left', fontsize=9)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 24)
    
    # 图3：能量平衡对比（日平均）
    ax3 = axes[1, 0]
    
    strat_names = list(strategies.keys())
    strat_labels = [strategies[s]['name'] for s in strat_names]
    
    R_n_avg = [np.mean(results[s]['energy']['R_n']) for s in strat_names]
    H_avg = [np.mean(results[s]['energy']['H']) for s in strat_names]
    LE_avg = [np.mean(results[s]['energy']['LE']) for s in strat_names]
    G_avg = [np.mean(results[s]['energy']['G']) for s in strat_names]
    
    x = np.arange(len(strat_names))
    width = 0.2
    
    ax3.bar(x - 1.5*width, R_n_avg, width, label='净辐射', color='#FFD700', alpha=0.8)
    ax3.bar(x - 0.5*width, H_avg, width, label='感热', color='#FF8C00', alpha=0.8)
    ax3.bar(x + 0.5*width, LE_avg, width, label='潜热', color='#4169E1', alpha=0.8)
    ax3.bar(x + 1.5*width, G_avg, width, label='储热', color='#8B4513', alpha=0.8)
    
    ax3.set_xlabel('缓解策略', fontsize=11)
    ax3.set_ylabel('日平均热流密度 (W/m²)', fontsize=11)
    ax3.set_title('能量平衡分量对比（日平均）', fontsize=12, fontweight='bold')
    ax3.set_xticks(x)
    ax3.set_xticklabels(strat_labels, fontsize=9, rotation=15, ha='right')
    ax3.legend(loc='upper right', fontsize=9)
    ax3.grid(True, alpha=0.3, axis='y')
    ax3.axhline(y=0, color='k', linewidth=0.5)
    
    # 图4：降温效果对比
    ax4 = axes[1, 1]
    
    # 计算相对于基准的降温效果
    T_max_baseline = np.max(results['baseline']['T_air'])
    T_avg_baseline = np.mean(results['baseline']['T_air'])
    
    cooling_effects_max = []
    cooling_effects_avg = []
    
    for strat_name in strat_names:
        T_max = np.max(results[strat_name]['T_air'])
        T_avg = np.mean(results[strat_name]['T_air'])
        cooling_effects_max.append(T_max_baseline - T_max)
        cooling_effects_avg.append(T_avg_baseline - T_avg)
    
    x = np.arange(len(strat_names))
    width = 0.35
    
    bars1 = ax4.bar(x - width/2, cooling_effects_max, width, 
                    label='最高温度降低', color='#FF6347', alpha=0.8)
    bars2 = ax4.bar(x + width/2, cooling_effects_avg, width, 
                    label='平均温度降低', color='#4169E1', alpha=0.8)
    
    ax4.set_xlabel('缓解策略', fontsize=11)
    ax4.set_ylabel('降温幅度 (°C)', fontsize=11)
    ax4.set_title('不同策略降温效果对比', fontsize=12, fontweight='bold')
    ax4.set_xticks(x)
    ax4.set_xticklabels(strat_labels, fontsize=9, rotation=15, ha='right')
    ax4.legend(loc='upper right', fontsize=9)
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar in bars1:
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                f'{height:.1f}°C', ha='center', va='bottom', fontsize=9)
    for bar in bars2:
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                f'{height:.1f}°C', ha='center', va='bottom', fontsize=9)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case4_mitigation_strategies.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"图表已保存: {output_dir}/case4_mitigation_strategies.png")
    
    # 创建缓解策略效果动画
    print("  创建缓解策略对比动画...")
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    lines = []
    for i, strat_name in enumerate(strat_names):
        line, = ax.plot([], [], label=strategies[strat_name]['name'], 
                       color=colors[i], linewidth=2)
        lines.append(line)
    
    ax.set_xlim(0, 24)
    ax.set_ylim(20, 45)
    ax.set_xlabel('时间 (小时)', fontsize=11)
    ax.set_ylabel('气温 (°C)', fontsize=11)
    ax.set_title('不同缓解策略气温对比动画', fontsize=12, fontweight='bold')
    ax.legend(loc='upper left', fontsize=10)
    ax.grid(True, alpha=0.3)
    
    def init():
        for line in lines:
            line.set_data([], [])
        return lines
    
    def update(frame):
        for i, strat_name in enumerate(strat_names):
            lines[i].set_data(hours[:frame+1], results[strat_name]['T_air'][:frame+1])
        return lines
    
    ani = animation.FuncAnimation(fig, update, init_func=init,
                                  frames=len(hours), interval=100, blit=True)
    
    ani.save(f'{output_dir}/case4_mitigation_animation.gif', writer='pillow', fps=10)
    plt.close()
    
    print(f"动画已保存: {output_dir}/case4_mitigation_animation.gif")
    
    # 打印关键结果
    print("\n关键结果汇总:")
    print("-" * 80)
    print(f"{'策略':<15} {'最高温度(°C)':<12} {'平均温度(°C)':<12} {'降温幅度(°C)':<12}")
    print("-" * 80)
    for strat_name in strat_names:
        T_max = np.max(results[strat_name]['T_air'])
        T_avg = np.mean(results[strat_name]['T_air'])
        cooling = T_avg_baseline - T_avg
        print(f"{strategies[strat_name]['name']:<12} {T_max:>10.1f} {T_avg:>12.1f} {cooling:>12.1f}")
    print("-" * 80)
    
    return results

# ============================================================================
# 主程序
# ============================================================================
def main():
    """主程序入口"""
    print("\n" + "=" * 70)
    print("城市热岛效应辐射分析仿真程序")
    print("主题049：城市热岛效应辐射分析")
    print("=" * 70)
    print(f"\n输出目录: {output_dir}")
    print("\n开始运行仿真案例...")
    
    # 运行案例1
    try:
        results1 = case1_surface_materials()
        print("\n案例1完成!")
    except Exception as e:
        print(f"案例1出错: {e}")
    
    # 运行案例2
    try:
        results2 = case2_urban_canyon()
        print("\n案例2完成!")
    except Exception as e:
        print(f"案例2出错: {e}")
    
    # 运行案例3
    try:
        results3 = case3_urban_heat_island()
        print("\n案例3完成!")
    except Exception as e:
        print(f"案例3出错: {e}")
    
    # 运行案例4
    try:
        results4 = case4_mitigation_strategies()
        print("\n案例4完成!")
    except Exception as e:
        print(f"案例4出错: {e}")
    
    print("\n" + "=" * 70)
    print("所有仿真案例运行完成!")
    print(f"结果文件保存在: {output_dir}")
    print("=" * 70)

if __name__ == "__main__":
    main()

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

当工程软件走进智能体时代，MATLAB昂首开启Agent之旅

AtomGit开源社区

基于 HarmonyOS 6.0 的跨端页面构建实战：停车管理系统开发解析

AtomGit开源社区

山东大学创新实训——群面智伴5（优化报告导出、新增工具调用、解决当前项目的部分问题）

本轮工作把群面系统从“AI 候选人可以自动讨论”推进到“用户可以真实参与、报告可以稳定生成、Agent 行为可以被观察和约束”的阶段。报告生成 Agent 的优化解决的是工程能力上的瓶颈：上下文如何压缩（InferencePack）、事实如何保留、工具如何稳定执行（ReAct 约束 + lenient 解析）、外部知识如何接入（Tavily）、长任务如何让用户可感知（SSE），以及失败时如何降级（