一种基于迭代模型调整的完整CFD方法论，用于提高高密度森林区域的风力模拟精度 Edouard Leonard, Ru Li, Eric Tromeur, Marianne Dupont, Aurélien Gaussorgues, Gaetan Martellozzo, Stavros Koutsioumpas, Mustafa Akcakaya

《Energies》：A Complete CFD Methodology Based on Iterative Model Adjustment to Improve Wind Simulation Accuracy in Highly Dense Forest Area Edouard Leonard, Ru Li, Eric Tromeur, Marianne Dupont, Aurélien Gaussorgues, Gaetan Martellozzo, Stavros Koutsioumpas and Mustafa Akcakaya

　　**摘要**
在茂密森林和复杂地形中进行风资源评估（WRA）仍然具有挑战性，因为强冠层引起的湍流和增强的风切变显著影响风流特性，并增加了建模的不确定性。依赖于植物面积密度（PAD）或叶面积密度（LAD）估算的方法需要昂贵的航空调查和特定地点的校准，这限制了它们的工业应用性。基于Meteodyn和EDF Power之间的科学合作，本研究提出了一种完整且可复制的计算流体动力学（CFD）方法，该方法围绕在Meteodyn WT?中实现的迭代模型调整（IMA）程序构建，利用标准工业输入来提高高森林地区的风资源评估精度。IMA程序通过使用单个参考测风塔的风速剖面测量值，迭代校准冠层阻力系数和森林模型参数，直到模拟的风切变与观测结果相匹配。该方法在芬兰、法国和苏格兰的三个地点进行了评估，在不同森林和复杂地形条件下进行了六次校准和交叉预测实验。交叉预测的不确定性显著降低，水平平均速度误差从[1.0–9.5%]下降到[0.5–2.2%]，全局平均绝对误差约为1.1%。这项研究为基于RANS的森林模型中冠层阻力项的敏感性提供了新的物理见解，表明阻力系数和冠层高度对风切变模拟都有重要且共同必要的影响。这些发现表明，在不依赖遥感获得的冠层密度数据集的情况下，可以在复杂地形和森林地区实现稳健准确的风资源评估，提供了一种实用且可工业扩展的替代方案。

**1. 引言**
由于有利条件（较少的法规和社会阻力）以及空间限制，风力发电场越来越多地建在茂密森林地区。然而，由于森林对风流的重大影响，风电场开发商在实现高精度风速模拟方面遇到了新的困难。在大气边界层中，由高粗糙度表面（如森林）引起的扰动会产生高水平的湍流和强风切变，这显然会影响风电场开发的风险[1,2]。最近有许多研究考察了各种数值模型[3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]，例如线性模型、使用CFD的RANS方法或LES模型，用于模拟林地环境中的风流。一些研究表明，线性模型的假设不足以再现密集森林地形的垂直风切变和流动特性[3,9,19]。LES通过将冠层表示为分布式动量源，并在湍流方程中加入额外的生产和耗散项（由PAD或LAD剖面驱动），正确捕捉了森林边缘的流动动态[4,5,6]，并且经常被用作验证更简单模型的基准[9,20]。然而，LES在工业风资源评估中仍存在两个问题。首先，计算成本过高：Ivanell等人[9]报告称，在一个森林基准地点，LES模拟所需的计算资源是相同域内等效RANS模拟的50到100倍。相比之下，本研究中使用的RANS方法在标准四核工作站上每个扇区只需13到18分钟即可收敛，完成20个扇区的评估仅需不到六小时。其次，LES需要能够正确再现大气边界层三维结构的时变湍流入流边界条件——通常通过单独的前体模拟生成——这增加了计算开销和方法论复杂性，与常规工业实践不兼容。尽管最近有计算进展，如GPU并行化和实用模拟策略（包括[19]提出的基于PALM的框架，该框架试图通过从测风塔数据中选择代表性稳定时间段并从中尺度再分析数据中指定时变边界条件来降低计算成本——LES在工业WRA方面仍然比稳态RANS复杂且计算要求更高。鉴于这些限制，使用不同湍流封闭方式的RANS CFD方法已成为森林和复杂地形中风资源评估最广泛采用的方法，它在物理精度和计算效率之间提供了实际平衡[3,9,10,11,12,13,14,15,21,22,23,24,25]。与线性模型不同，RANS解决了完整的时均动量和湍流传输方程，能够通过分布式体积阻力项表示流动分离、再循环和冠层引起的动量吸收[21,22]。使用RANS CFD方法在林地和复杂地形中进行了多次成功的风资源评估[3,9,11,12,13,14,15,21,22,23,24,25]，Barber等人[19]在五个复杂地形地点的七个工作流程中确认，当森林和地形复杂性共存时，RANS CFD始终优于线性模型。然而，将RANS CFD应用于高森林地区仍具有挑战性。多项研究表明，当缺乏或对冠层密度输入的约束不足时，RANS的准确性会显著下降[9,21]。Ivanell等人[9]在Ryningsn?s森林基准地点证明，使用默认或基于文献的阻力系数的RANS模型产生的误差与线性模型相当。因此，大多数工业RANS应用依赖于默认或经验分配的阻力系数，这在风切变预测中引入了较大且依赖于地点的不确定性，尤其是在密集和异质森林条件下[13,14]。Katul等人[22]明确指出了在RANS森林模型精度中参数α（α）的中心作用。它表明α是模型不确定性的主要来源，它将森林密度、树种、冠层孔隙率和叶角分布的综合效应聚集为一个无法直接测量的单一体积参数。同一研究还表明，一阶混合长度封闭方法在再现冠层流中的平均风速、剪切应力和湍流动能剖面方面与二阶模型表现相当，而高阶封闭方法的理论优势被常规应用中α的较大不确定性所抵消。这一发现表明，准确确定α对于RANS森林模型的性能比选择湍流封闭方式更为关键——这直接激发了本研究提出的校准方法。目前没有实用的方法可以在没有遥感数据的情况下确定工业RANS森林建模的特定地点α值——本研究直接解决了这一空白。基于EDF Power和Meteodyn之间的科学合作，本研究旨在开发一种清晰且可复制的CFD方法，仅使用标准工业输入来提高高森林地区的风模拟精度。因此，本研究使用的材料仅限于通常在操作性风资源评估项目中可获得的材料：公开可用的地形和粗糙度地图，以及单个气象测风塔的风速剖面测量数据。本文提出了一种新型方法，使用基于单个参考测风塔的现场风速测量的迭代程序来确定CFD森林模型参数。以往的研究通过两种主要方法解决了α的确定问题：要么从文献或风洞实验中规定固定值[13,22,28]，要么从航空LiDAR调查中推导出空间变化的PAD产品[7,9,21]。相比之下，本研究直接从单个测风塔的多个高度处的风速剖面测量中推断α值，无需任何遥感数据。据我们所知，这是第一项仅通过单个气象测风塔的风速剖面测量来校准基于RANS的森林冠层模型，并在复杂地形中进行系统交叉预测的研究。Ivanell等人[9]指出，Sogachev等人[29]的湍流封闭常数在异质森林地点表现良好，同时也指出从ALS导出的PAD数据确定α是标准方法。Lopes等人[30]强调，标准参数化或通用系数值都不适用于任意地点条件——突显了α的根本地点依赖性，即固定的文献值无法捕捉这一点。本文提出的IMA程序通过从反映实际森林条件的测量中推断α值，直接解决了这一不确定性。该方法在芬兰、法国和苏格兰的三个地点进行了评估，每个地点都配备有两个气象测风塔、共同的测量周期和在不同地形复杂条件下的异质森林覆盖。这种配置产生了六个独立的校准和交叉预测案例，涵盖了广泛的森林密度、冠层高度和测风塔之间的距离（从2到6.6公里），为评估所提出方法的通用性和工业相关性提供了坚实的基础。本研究的主要贡献包括：一种完整且可复制的CFD方法（IMA），用于密集森林和复杂地形中的风资源评估，仅需要标准工业输入——公开可用的粗糙度和地形地图以及单个气象测风塔的风速剖面测量数据——而不依赖于遥感获得的冠层密度数据；对森林模型参数（冠层高度和阻力系数）的系统性敏感性分析，证明了这两个变量对密集森林条件下的风切变模拟有影响；在包括六次交叉预测案例在内的异质森林条件和复杂地形下验证了空间常数阻力系数的假设，发现森林引起的流动效应可延伸至周围开阔区域1公里范围内，表明IMA可以从森林边界外的测风塔应用。

**2. 地点描述和计算设置**
**2.1. 地点描述**
本研究中的气象测风塔数据由EDF Power在三个国家提供：芬兰、法国和苏格兰。这些地点因其具有挑战性的条件而被特别选中，它们结合了密集的森林覆盖和不同水平的地形复杂性。用于地形描述和从卫星数据中得出的公开可用的地形和土地覆盖地图分别列在表1中，并在图1、图2和图3中展示。为了尽可能匹配测量周期和现场观测，选择了两种不同的土地覆盖地图：对于站点2和3，使用了Corine土地覆盖（CLC）粗糙度地图[31]；对于站点1，使用了Copernicus土地覆盖（LC）粗糙度地图[32]。表1. 地点描述。图1. 站点1的地形地图（左）和粗糙度地图（右）。两个测风塔测得的主导风向均为230°。图2. 站点2的地形地图（左）和粗糙度地图（右）。两个测风塔测得的主导风向均为225°。图3. 站点3的地形地图（左）和粗糙度地图（右）。两个测风塔测得的主导风向均为210°。所有站点均使用NASADEM地形地图，分辨率为30米[33]，除了站点1使用EU-DEM 1.1地图[34]，因为NASADEM地形地图在芬兰不可用。地形条件从站点1的 moderately complex 地形到站点2和3的 complex 地形不等，站点1的最大高差约为100米，站点2为500米，站点3为520米。在风工程应用中，粗糙度长度值通常使用转换表从土地覆盖地图中推导出来[35,36,37,38,39]。这些转换表提供了有用的粗糙度信息，但被认为具有不确定性，并且强烈依赖于应用地区的地理区域，导致同一数据集的推荐值存在差异。在Meteodyn WT中，转换表是使用基于多个可用来源的全球平均方法确定的[23]。例如，图1、图2和图3中显示的对应于针叶林区域的0.7米粗糙度长度值是从八个不同的转换表来源确定的，这些转换表中的粗糙度长度值范围从0.5米到1.2米不等。三个站点展示了多种森林覆盖条件。在100米分辨率的粗糙度地图上，所有站点的粗糙度长度变化范围从0.003米到0.7米，显示出具有开阔地的异质森林覆盖。1号测风塔（M1）和2号测风塔（M2）附近的冠层高度是通过现场观察、照片以及Google Earth影像估算得出的。这些估算基于以每根测风塔位置为中心、半径约为400米的四分之一圆形区域，并且方向与主导风向一致。表1中呈现的观察结果与粗糙度地图一致，表明在主导风向大约200°到230°的范围内，1号地点的M1和M2之间的森林状况几乎相同。而在2号和3号地点，M1和M2之间的森林状况存在明显差异。第4.1节中的敏感性分析探讨了这些冠层高度估算的不确定性，结果表明，如果冠层高度估算误差为7米，则平均交叉预测误差约为0.8%。根据现场观察，这三个地点的主要树种都是针叶树。

每个地点都安装有两根气象测风塔，每根测风塔上在三个不同高度处安装了杯式风速计，具体信息见表1。所有测量数据覆盖了一年的时间段，数据质量高，覆盖范围超过98%。所有地点的风速测量不确定性假设为2.5%，这与风工程指南和类似测风塔的观测结果相符[40,41,42]。如表1所示，1号、2号和3号地点的测风塔之间的距离分别为6.6公里、4.4公里和3.1公里。

2.2. 计算设置
本研究使用的所有CFD参数详见表2。在网格生成方面，建模时采用了最小垂直分辨率为4米、水平分辨率为25米的网格。垂直扩展率为1.2，使得网格尺寸从地面逐渐向顶部扩展；水平扩展率为1.1，使得水平网格尺寸从中心向边缘扩展。此外，Meteodyn WT?（版本1.10.1）提供了表面网格细化选项，也称为“网格映射”，这确保了在任何感兴趣的区域都能获得细致的网格分辨率，而不会受到水平扩展率的影响。

2.3. 风向划分
风向被划分为20个风向扇区。例如，在第一个扇区内，风速从0°方向（正北方向）到18°方向（东北方向）进行分类；在第二个扇区内，风速从18°方向到36°方向进行分类。每个地点都针对20个风向扇区进行了CFD计算。所有地点和风向扇区的CFD计算都达到了完全数值收敛，残差值低于所需容差，从而得到了可靠且一致的流场解决方案。所有地点的每个扇区的平均CFD计算时间根据输入的土地覆盖图和达到适当收敛所需的计算时间不同，范围在13分钟到18分钟之间。所有计算都是使用标准性能计算机和多核并行计算（4个核心，64 GB RAM）来执行的。

3. 森林CFD方法与验证框架
本节介绍了本研究提出的用于密集森林区域风资源评估的完整CFD方法。该方法包括四个相互关联的组成部分，以下将依次进行介绍。第3.1节介绍了在Meteodyn WT?中实现的RANS控制方程，包括湍流闭合模型和边界条件。第3.2节描述了森林冠层模型，并确定了体积阻力系数作为主导森林地形中风切变模拟的参数。第3.3节介绍了迭代模型调整（IMA）程序——这是本研究的核心贡献之一，该程序能够在无需遥感数据的情况下，仅根据单个参考测风塔的风速剖面测量值来确定体积阻力系数。最后，第3.4节定义了交叉预测验证框架和误差指标，用于在第4节中评估IMA的性能。

这三个组成部分之间的逻辑关系在图4中进行了总结：控制方程和冠层模型构成了物理基础，IMA使用现场测量值校准关键参数，而交叉预测框架则在三个国际地点的六个独立案例中应用了验证策略。图4展示了森林CFD方法与验证框架的总体情况。标准工业输入数据（灰色）被输入到CFD森林模型中（蓝色，第3.1节和第3.2节）。IMA程序（橙色，第3.3节）通过迭代校准体积阻力系数，直到模拟的风切变指数在10%的容差范围内与观测值匹配。校准后的模型随后应用于交叉预测验证框架中（蓝绿色，第3.4节），其结果使用MAE和MAPE指标进行评估。

3.1. 森林建模
计算流体动力学（CFD）方法基于求解完整的雷诺平均纳维-斯托克斯方程（RANS），从而能够模拟复杂地形上的流动分离和再循环现象。Meteodyn WT?是一款商用场地评估软件，用于模拟大气边界层。湍流通过混合长度模型[43]用湍流粘性假设来表示。通过结合先进的数值方法——如耦合多网格求解器，Meteodyn WT?实现了高效的解决方案[44]，并且已在多种情况下得到了验证[3,13,14,25,45]。在Meteodyn WT?中，不可压缩的大气流动由雷诺平均纳维-斯托克斯方程描述：
[此处应插入方程式]

其中ρ是空气密度（假设为恒定值），u是平均速度分量，p是平均压力，η是有效粘度。μ项代表了包括冠层引起的阻力在内的外力。

雷诺应力张量是由非线性平流项的平均值产生的，反映了湍流对平均流动的影响。这一项引入了额外的未知数，因此需要一个湍流闭合模型。纳维-斯托克斯方程采用布辛涅斯克假设进行闭合：
[此处应插入方程式]

湍流粘度η是通过普朗特混合长度模型[44]获得的：
[此处应插入方程式]

湍流动能的产生项可以表示为：
[此处应插入方程式]

耗散项可以表示为：
[此处应插入方程式]

其中C是一个与热稳定性条件相关的常数系数。湍流长度λ由以下公式给出：
[此处应插入公式]

湍流通量通过湍流粘度的概念与平均变量的梯度相关联，后者可以表示为风速尺度（通常为湍流动能k的平方根）和湍流长度尺度的乘积。k通过下面的湍流动能传输方程求解：
[此处应插入方程式]

在入口处，对数风速剖面和湍流动能剖面取决于热稳定性条件。出口处采用自由出流条件。在底部边界，假设湍流动能的产生被耗散所平衡。在顶部边界，也采用自由出流条件；对于k也施加了对称边界条件。在侧边边界同样施加了对称边界条件。

3.2. 冠层模型
在Meteodyn WT?中，冠层区域被视为阻力作用的动量能量汇，湍流长度尺度在此处被修改。森林引起的扰动通过在动量守恒方程中加入阻力项来建模：
[此处应插入方程式]

其中v是主方向上的风速，c是与森林密度和垂直叶面积剖面成比例的体积阻力系数。该叶面积剖面用于描述树木的形状，考虑到树木顶部的叶片较少。因此，如图5所示，森林顶部25%区域的阻力系数被线性减小。这种线性减小反映了针叶林中观察到的LAD剖面的减少[7]。

3.3. 迭代模型调整（IMA）方法在森林区域的应用
3.3.1. 物理基础和反问题构建
准确确定体积阻力系数是基于RANS的森林风建模中的核心挑战。该系数同时控制了动量吸收、湍流动能的产生以及冠层区域内的湍流动能耗散——使其成为建模不确定性的主要来源，无论湍流闭合阶数如何[22]。然而，这一系数并非可以直接观测得到：它综合了森林密度、树种、冠层孔隙度以及叶面积分布等多种因素的影响，形成了一个无法通过实地直接测量的单一体积参数。

本研究提出的IMA方法通过将这一系数视为一个未知量来解决这一挑战，该系数需要根据单个参考测风塔的风速剖面测量值来推断。具体来说，IMA解决了一个反问题：给定参考测风塔观测到的风切变指数，确定特定地点的阻力系数，使得RANS模拟在规定的收敛容差范围内再现该指数。风切变指数描述了风速随高度的变化规律：
[此处应插入方程式]

在本研究中，from是每个测风塔最高和最低风速测量值计算得出的，并作为整个第3.4节中主要的垂直剖面验证指标。这种反问题构建在物理上是合理的，因为风切变指数与模拟的风切变指数之间存在单调关系：增大from会增强冠层区域的动量吸收能力，从而降低低处风速，同时保持毂高处的风速基本不变，进而增加from值。相反，减小from会减弱动量吸收，使得垂直风速剖面更加平坦，从而导致from值降低。这种单调关系在所有三个地点的敏感性分析中都得到了证实。该方法基于这样的假设：在整个森林环境区域，某个特定参数（此处未明确指出具体参数）在空间上是恒定的。尽管在米级尺度上森林密度明显具有异质性，但IMA假设认为，在风资源评估的尺度上（塔间距为2-6.6公里），森林阻力的空间变异性会趋于平均，并且可以用一个单一的、代表整个站点的平均值来表示。这一假设类似于在大型模型中使用单一的粗糙度长度来表示异质表面，但它并非显而易见，需要通过实验进行验证。

3.3.2 IMA程序的描述
IMA程序如图6所示，包括两个连续的步骤。完整的工作流程总结如下：
图6. 迭代模型调整（IMA）方法论的描述。
从初始森林模型参数开始，不断迭代调整阻力系数，直到模拟的风切变指数在参考塔处与观测值相差在10%的容差范围内。

步骤1——森林模型初始化：树冠高度和粗糙度长度的调整。
树冠高度由粗糙度长度和森林比例来确定。首先根据现场观测结果调整这些变量，初始的粗糙度地图值和森林比例是基于默认值。粗糙度长度值由土地覆盖图提供，而森林比例则是根据现场访问、照片或Google Earth影像估计的近似树冠高度来确定的[35,36,37,38]。
当整个区域的树冠高度被系统性地低估时，相应地重新调整森林比例。例如，对于一个以针叶树为主的站点，如果树冠高度约为30米，粗糙度长度为某个值，则相应的森林比例为另一个值。相反，如果观察到局部树冠高度不一致，则通过修改转换表或手动调整该区域的粗糙度长度值来纠正。

步骤2——阻力系数的迭代调整。
参考阻力系数直接控制树冠动量汇的强度，从而控制模拟的风切变指数。该参数被视为一个自由参数，并按以下方式迭代确定：
1. 初始化（Meteodyn WT? 默认值）。
2. 使用当前的阻力系数值对所有20个风向扇区进行完整的CFD模拟。提取参考塔位置的模拟平均风速剖面，并计算模拟的切变指数。
3. 计算切变收敛标准：（此处公式未给出）。
4. 如果满足收敛标准，则校准完成，进入交叉预测阶段；如果不满足，则增加阻力系数的值并返回步骤1。

公式（16）中的10%收敛标准被认为是一个实际可行的阈值，它比所有六个参考塔配置下观察到的默认模拟误差更为严格，同时由于阻力系数和某个参数之间的单调关系，这一标准在实践中是可以实现的，从而保证了收敛到唯一解。在三个地点的所有六个配置中都成功实现了收敛，迭代次数在2到7次之间（表3）。

表3. 用于IMA程序的森林模型参数。
阻力系数的增量值是校准精度和计算效率之间的折中选择。第4节中的敏感性分析表明，阻力系数的变化0.01会导致三个地点的平均切变指数变化约为0.02-0.05——远低于10%的收敛阈值——确保了迭代过程不会超出目标值。六个校准案例中所需的最大迭代次数为7次（站点2，塔1），平均每次迭代大约3到4次。

3.3.3 IMA程序的范围和假设
IMA方法基于两个假设，这些假设在结果和讨论中都有明确提及：
- **的空间均匀性**：假设该参数在整个站点范围内是恒定的，并不要求森林覆盖在空间上完全均匀。树冠的空间异质性——包括森林斑块、空地和不同植被密度区域之间的树冠高度变化——已经通过土地覆盖图中的空间变化粗糙度长度来表示，该粗糙度长度决定了每个计算单元的局部树冠高度。因此，森林模型的激活在区域内随空间变化，与粗糙度图保持一致。假设一个单一的参数可以代表整个站点的森林阻力特性。虽然该参数决定了阻力场的空间分布和垂直范围，但另一个参数则决定了其空气动力学的大小——具体来说是每单位体积树冠的阻力，这取决于树种特征，如针叶密度、树枝结构和树冠形状。这一假设在第4节的六个交叉预测案例中得到了验证，这些案例包括异质森林条件和塔间距从2到6.6公里的情况，总体平均绝对误差约为1.1%。
- **中性大气稳定性**：CFD模拟仅在中性稳定性条件下进行，基于全年平均值接近中性条件的假设[27]。虽然大气稳定性已被确定为森林地形建模中的一个重要不确定性来源[16]，但中性假设与标准工业WRA实践是一致的。热稳定性对IMA校准的影响仍然是一个待解决的问题，留待未来的研究。

最终使用调整后的森林模型，以参考塔作为入口边界条件来模拟整个站点的风流。IMA可以应用于方向平均值或扇区平均值的风速垂直剖面。在这项工作中，IMA应用于包括所有风向扇区的平均垂直剖面。为了减少CFD计算时间，建议在校准迭代阶段仅对主导风向扇区应用IMA，然后使用完整的多扇区模拟来验证结果。

3.4 WRA中的交叉预测方法
在这项工作中，使用交叉预测方法和基于测量点的风况比较来验证CFD模拟方法，这是风能行业指南[48]推荐的。具体来说，对于每个站点，使用一个塔的位置来估计区域内（以及其他塔位置在不同高度）的风速，或使用风切变指数。这一程序同时测试了模型的两个方面准确性：水平外推——CFD模型将风信息从参考塔正确传递到2-6.6公里异质森林地形中的第二个塔位置的能力；以及垂直外推——用于估计测量高度以上至风力涡轮机轮毂高度的风速的模拟风切变剖面的准确性。每个高度的风速值和风切变指数（在公式（15）中定义）都被用作验证指标。

这种交叉预测框架特别适合验证IMA方法，因为它直接测试了核心假设，即在一个参考塔处校准的参数可以用来准确预测站点内另一个位置的风况。对于每个站点，该过程应用两次——一次使用M1作为参考塔来预测M2的位置，另一次使用M2作为参考塔来预测M1的位置——从而在三个站点产生了六个独立的交叉预测案例。这种方法与标准工业WRA实践一致，通常是从附近的一个塔的测量数据来估计潜在风力涡轮机位置的风况[48]。
模型准确性使用两个互补的误差指标来量化。平均绝对误差（MAE）分别应用于最高塔高度的风速和风切变指数，而平均绝对百分比误差（MAPE）则表征了所有三个测量高度的整体垂直剖面准确性。

4. 结果与讨论
4.1 森林模型的敏感性
在本节中，我们提出研究CFD森林模型的敏感性。在三个站点上，分别使用第一个塔M1和第二个塔M2作为每个站点的唯一参考，分析了在其他塔处的交叉预测结果。例如，对于站点1，在选择了参考塔M1之后，将塔M2的模拟风速与测量值进行比较，并在图7中标记为“站点1-M2”。为此考虑了两种不同的情景：

图7. 输入CFD森林模型变量的敏感性影响：森林比例和阻力系数对风速交叉预测的影响，使用MAPE变异性来表示。对于每个站点，使用另一个塔作为参考来模拟风速值。这意味着M1被用作参考塔来估计M2，M2被用作参考塔来估计M1。
- **树冠高度对风流的影响**：在敏感性分析中，通过森林比例研究了树冠高度对风流的影响，使用了10、20、30、40和50这些值，同时应用了一个参考阻力系数。这意味着每个网格单元的树冠高度是通过将这些比例乘以粗糙度长度值来计算的。例如，如果考虑的粗糙度长度值为0.7米（对应于针叶树区域），则检查了对风流的影响：7米、14米、21米、28米和35米。
- **森林密度的影晌**：通过改变阻力系数的值（0.001、0.005、0.01、0.02和0.05）来研究森林密度的影晌，同时保持森林比例不变。图7使用平均绝对百分比误差（MAPE）（在公式（18）中描述）展示了输入CFD森林模型变量对风速交叉预测的影响。结果表明，这两个变量对风速交叉预测的影响相似。在森林比例[10-50]的范围内，结果的变异性介于1%到7%之间；而在阻力系数[0.001-0.05]的范围内，森林密度的变异性介于1%到10%之间。每个站点和每个塔的模型变异性非常大。这种站点依赖的变异性反映了三个站点之间的异质森林条件。值得注意的是，在站点2上，M1和M2分别观察到最大和最小的变异性，显示出该项目中最极端的森林条件。如表1和本研究的后续内容所描述，站点1位于非常高且密集的森林条件下，而M2位于靠近空地的位置，因此受森林影响较小。
- **树冠高度估计的敏感性**：这种分析也有助于估计每个塔位置周围树冠高度估计的敏感性影响。如表1所示，这些值是根据现场观察、照片和Google Earth影像估计的。如图7所示，由于7米树冠高度的变化，MAPE变异性在0.04%到2.3%之间，平均值为0.8%。值得注意的是，在考虑的较高阻力系数范围内也观察到相似的平均变异性约为1%。这意味着树冠高度估计的7米误差平均只会对交叉预测的MAPE产生影响约1%。因此，在塔周围的森林高度估计上观察到相对的容差。

正如第4.1节中的敏感性分析所示，单独调整阻力系数或树冠高度都不足以再现观测到的风切变指数——即使指定不切实际的树冠高度，当阻力系数保持在其默认值时，模拟的风切变也只能略微接近平观值。因此，这两个参数必须共同校准，这证实了两步IMA程序的必要性，并与之前的研究[14,20,22]一致。此分析强调了系统化方法的必要性，以共同校准树冠高度和阻力系数——即第3.3节中提出的IMA方法。从这种敏感性分析中可以得出几个具有物理意义的观察结果。这两个参数的默认值都位于测试范围的较低端，这表明默认配置系统性地低估了密集针叶林条件下的森林诱导风切变。这种预测通过IMA校准得到了后验验证，在所有三个站点，这两个参数都被一致地调整为更高值。与相比，产生了更大的MAPE敏感性范围（1-10%），尽管覆盖的相对范围更大。这两个参数都对风切变产生积极且不可补偿的影响——增加或都会增加模拟的风切变指数，并且没有一个参数可以替代另一个参数。这种正相关确保了IMA两步校准程序是合理的：首先根据站点观测数据固定，然后在第二步中调整，从而得到唯一解，因为没有相反参数变化的组合能够产生等效的风速剖面。

4.2. 在每个参考桅杆上应用的IMA方法
本节的目的是详细描述每个站点和参考桅杆上的IMA方法。在参考桅杆上评估IMA校准，因为参考桅杆是域内唯一有风速剖面测量的位置，这些测量数据用于确定校准完成的收敛标准——迭代循环需要观测数据来决定何时校准完成。在参考桅杆上评估校准后的模型确认IMA程序正确识别了能够再现局部风速剖面的森林模型参数。然而，仅在校准点再现条件并不足以构成验证——一个具有足够自由参数的模型总能拟合其训练数据。IMA的科学价值仅在第4.3节的第二个桅杆上进行交叉预测时得到证明，其中未经进一步调整地应用了校准后的参数来预测模型未进行校准的位置的风况。表3显示了每个站点上分别使用的模型参数。对于每个站点，已经应用了图6中描述的IMA方法（每个站点20个风CFD扇区的相同设置），以模拟森林对每个桅杆位置平均风速剖面的影响。从初始参数开始，然后使用IMA程序修改森林模型参数，以在参考桅杆的测量风速垂直剖面和CFD模拟之间达到合适的调整。如前所述，初始森林比率设置为，参考阻力系数设置为。然后修改这些森林模型参数以更好地适应测量的风速垂直剖面。在这项工作中，我们认为当模拟结果与测量值之间的剪切差异小于10%时，森林模型参数是可接受的。这个最终结果分别是在站点1、2和3进行了2次、7次和1次IMA程序迭代后获得的。以下是详细的迭代模型调整（IMA）程序：
在站点1和站点3，首先根据表1中呈现的树高推导出森林比率。例如，在站点3，森林比率为。
在站点2，森林比率设置为默认值，并根据现场观测对粗糙度图进行了局部修正，将针叶林区域的树高从14米增加到35 ± 5米。因此，在粗糙度图转换表中将针叶林区域对应的粗糙度长度从0.7米修改为1.7米。
在确定森林比率后，参考阻力系数进行了迭代调整，直到模拟的风切变指数在10%的收敛标准范围内与观测值匹配。

通过图8中的风切变指数比较和图9中的垂直剖面比较来检验IMA方法的效果。对于每个站点，选择最高桅杆高度（站点1和2为100米，站点3为70米）作为参考高度。图8显示了迭代模型调整（IMA）方法对每个参考桅杆剪切值的影响。图9显示了IMA方法对每个站点和参考桅杆的影响。相应的剪切值（）也显示在图例中。在图8中，使用默认的和得到的风切变被系统性地低估了。通过IMA程序校准的森林模型参数提供了每个站点局部树冠条件的更物理准确的表示，使模拟的风切变指数与所有站点参考桅杆配置的观测结果非常接近。站点2 M1的观测剪切值最高，而站点2 M2的剪切值最低。
如图9所示，在参考桅杆处，风速垂直剖面的改善在低海拔地区最为明显。这种改善在站点2 M1最为显著，这是本研究中最密集的森林配置，其中默认模拟低估了观测剪切指数0.3——这是所有六个案例中最大的差异。无论选择哪个桅杆作为参考，IMA程序都能产生一致的森林模型参数，表明校准后的代表了整个站点，而不仅仅是参考桅杆的局部森林条件。这一结果在站点2尤为显著，因为M1和M2位于明显不同的森林环境中——M1位于非常密集的针叶林中（米），而M2位于空地附近（米）。虽然M2位于空地附近，测得的剪切指数为0.25，但M1被密集的针叶林包围，估计的树冠高度为35米，显示出本研究中的最高剪切指数，接近0.6。这一点也通过图片和谷歌地球图像得到了证实（此处未展示）。IMA程序从两个桅杆得出了相同的校准结果。这表明，尽管M2位于森林之外且距离非常密集的森林超过1公里，森林诱导的动量效应仍延伸到森林边界以外1公里的风向。因此，IMA不仅适用于森林内，也适用于周边区域。尽管风能行业指南建议使用具有与风力涡轮机站点相似粗糙度条件的代表性桅杆[48]，但当前结果表明，只要位于森林效应范围之内，位于森林边界之外的桅杆也可以提供可靠的IMA校准数据。

4.3. 迭代模型调整（IMA）对交叉预测的影响
图10、图11和图12展示了所有站点和桅杆的交叉预测结果。对于每个站点，使用在参考桅杆上获得的校准森林模型来模拟第二个桅杆位置的风速垂直剖面，然后将其与所有三个风速计高度的测量值进行比较。图10显示了IMA方法对每个站点和桅杆的交叉预测风速剖面的影响。对于每个站点，使用另一个桅杆作为参考来模拟交叉预测的风速值。这意味着使用M1作为参考桅杆来估计M2，使用M2作为参考桅杆来估计M1。图11显示了IMA方法对交叉预测风切变值的影响。对于每个站点，使用另一个桅杆作为参考来模拟交叉预测的风速值。这意味着使用M1作为参考桅杆来估计M2，使用M2作为参考桅杆来估计M1。图12显示了IMA方法对交叉预测风速值的平均绝对百分比误差（MAPE）的影响。对于每个站点，使用另一个桅杆作为参考来模拟交叉预测的风速值。这意味着使用M1作为参考桅杆来估计M2，使用M2作为参考桅杆来估计M1。在所有情况下都观察到IMA对交叉预测垂直剖面的明显影响，这一点在图10的风速剖面形状和图11的剪切指数值中都非常明显。当增加时，低海拔处的风速降低，而高海拔处的风速增加。如图12所示，所有情况下的MAPE都减少了，证实了使用IMA方法可以减少CFD模型的不确定性。在适度森林化的条件下（站点1和3以及站点2 M1），MAPE的减少范围从0.8%到4.4%。对于站点2的极端森林条件M1，观察到了15.8%的改善。正如预期的那样，当默认森林参数显著低估实际树冠高度和密度时，IMA方法似乎最为相关。
测量值和模拟值之间的剪切差异从[0.08–0.3]显著减少到[0.01–0.07]，证实了IMA方法改善了垂直剖面模拟。图8和图11中观察到的相似模式证实了在参考桅杆校准的森林模型参数不仅控制了参考位置的风切变，还控制了整个站点域的风切变——验证了IMA校准的空间可转移性。IMA还提高了最高测量高度（站点1和2为100米，站点3为70米）处的水平外推精度。如图12所示，应用于水平外推风速的MAE在所有情况下都减少了，除了站点3 M1。在站点2处，100米处的风速误差从9.5%减少到0.7%。由于风力涡轮机轮毂高度通常超过80米，这一结果对于森林地形中的风电场能量产量估计特别相关。在低海拔以及轮毂高度附近，IMA方法的影响也很明显。所有六个案例中准确的交叉预测结果验证了IMA方法的假设——即参考阻力系数在整个站点范围内是空间恒定的。这一假设在文献中很常见，在物理上是有根据的，因为树木树冠引起的阻力形式占主导地位[30]。然而，交叉预测的范围仅限于[2–6.6公里]。IMA是否适用于更长的距离（>7公里）或针叶林和混交林类型仍是未来工作的开放性问题。
与之前的工作[7,9,12,15,21,22,26]相比，这些结果表明，即使在没有LAD数据的先进遥感的情况下，也可以使用单方程湍流闭合RANS CFD在高度森林化区域进行准确的速度交叉预测。这是通过使用IMA方法基于单个参考桅杆的全年风速剖面测量来仔细估计森林密度而实现的。虽然IMA-RANS方法提高了与测量风速剖面的吻合度，但其固有限制应当得到承认。RANS平均了湍流，无法捕捉树冠内的细小湍流结构或瞬态阵风，而且树冠由忽略详细树木几何形状和局部效应的总体阻力系数表示。从单个桅杆进行校准也限制了空间精度，对于树木异质性较强的站点，在桅杆间距较小的尺度上也是如此。尽管有这些限制，RANS在计算上仍然高效，并且适用于密集森林中的风能评估，而LES可以在需要高分辨率流动信息的区域作为补充工具。

5. 结论
本研究提出了一种基于迭代模型调整（IMA）程序的创新且可复制的CFD方法，以提高在密集森林和复杂地形中的风能资源评估精度，该方法仅需要标准工业输入——公开可用的粗糙度和/或地形图以及单个气象桅杆处的风速剖面测量数据。该方法在三个国际站点的六个校准和交叉预测案例中得到了验证，这些站点由EDF Power提供，使用了Meteodyn WT?在异质森林和复杂地形条件下进行了测试。关键发现总结如下：
基于IMA方法的准确风速模拟：使用IMA程序确定的适当森林高度和阻力系数值时，风速估计精度显著提高。模拟和测量值之间的剪切差异从[0.08–0.3]显著减少到[0.01–0.07]。水平平均速度误差已从 [1.0–9.5%] 降低到 [0.5–2.2%]，基于六个案例的计算，全球平均绝对速度误差为 1.1%。对森林建模的物理理解表明：阻力系数和树冠高度对风切变模拟都有类似且共同必要的影响——仅靠其中一个参数不足以再现密集森林条件下的观测到的风速分布。默认参数值系统性地低估了森林引起的气流剪切效应，而 IMA 在所有三个站点都始终校准出更高的值。假设树冠高度在空间上是恒定的（即整个站点的物种级空气动力学特性一致）在包括异质森林覆盖、以针叶林为主的森林、复杂地形以及预测距离为 2–6.6 公里在内的所有六个案例中都得到了验证。森林引起的气流效应可以观察到距离森林边缘 1 公里范围内的区域，这表明 IMA 可以应用于位于森林边界外的代表性测风塔——这扩展了使用具有相似粗糙度条件的测风塔的标准建议 [48]。虽然预测距离超过 7 公里超出了复杂地形下的既定指南 [48]，但 IMA 在更远距离和不同类型的森林（如落叶林、混交林）中的适用性仍然是一个未解决的问题。

使用基于卫星数据的土地利用图可以为森林风速建模提供有用的基础。然而，在所有三个站点都需要转换表和森林比例校正，以获得更符合物理实际情况的树冠高度。IMA 结果对中等范围的树冠高度估计不确定性具有鲁棒性——7 米的误差对交叉预测的平均绝对百分比误差 (MAPE) 的影响小于 1%，平均变异率为 0.8%。结果表明，仅使用标准工业输入数据（公共地图和单个测风塔的测量结果）就可以获得准确的计算流体动力学 (CFD) 结果。这种方法已在涵盖各种地形和森林条件的六个案例中得到验证。这里提出的观察结果被认为适用于其他林地研究。

与之前的研究 [9,21,22] 相比，当前结果表明，在高度森林化的地区，即使没有来自航空调查的 PAD 或 LAD 数据，只要通过 IMA 校准程序使用全年风速测量数据仔细估计森林密度，也可以实现准确的交叉预测。然而，也应承认几个局限性：RANS 方法无法捕捉到细尺度湍流结构、瞬态阵风或详细的树冠几何形状——LES 方法在这些特征至关重要的区域可以提供更高分辨率的时空流场信息。关于树冠高度恒定的假设仅针对以针叶林为主的森林得到了验证；对于落叶林或混交林，由于林下层（LAD）随季节变化且物种多样性导致树冠阻力特性在整个站点范围内存在差异，其适用性仍是一个未解决的问题。预测距离超过 7 公里超出了现有的 WRA 指南 [48]，并且尚未进行测试。此外，从单个测风塔进行校准也会限制在测风塔间距较小且森林异质性较强的站点上的空间精度。

在整个 CFD 模拟过程中假设了中性大气稳定条件，基于全年平均可以整合所有稳定性类别的假设，从而得到代表近中性条件的平均风速分布 [27]。这一假设是必要的，作为将森林引起的气流效应与热稳定性效应分开的首次近似，从而便于进行 IMA 校准。尽管所有三个站点使用的全年测量周期支持了中性假设，但现有方法并未明确考虑大气热稳定性对风场时间变化的影响。特别是，在非中性稳定条件下（即稳定分层或对流不稳定可能与树冠阻力效应相互作用的情况下），IMA 的适用性尚未确定，并被确定为未来工作的重点。