一种选取数值模型最优方案组合的方法与流程

本发明涉及地球科学技术领域，尤其涉及一种选取数值模型最优方案组合的方法。

背景技术：

在地球科学领域，尤其是大气科学领域，使用数值模型进行模拟是一种常用的实验手段。科学家们对自然界中的物理现象进行抽象，用数学方程将该现象进行模拟，并用计算机程序予以实现，从而形成了数值模型。通过数值模型可以预测未来气候变化的影响，可以进行气候变化的归因分析等。数值模型的建立大大增进了人类对大自然的认识，对人类科学的发展起到了极其重要的推动作用。

在大气科学领域，实际模拟中，模型需要将二维的研究区域划分为等分的网格，每个网格内看作同样的情况来进行模型。由于模拟的研究区域一般比较大且模拟花费的时间较长，模拟的网格通常不会太小，一般情况下，网格的间隔为0.5°，即50km，这样每个网格的大小是50km×50km。此时新的问题出现了，有些自然现象的影响范围是小于50km的，即数值模型完全没有能力将这些现象进行模拟，这些过程称为“次网格过程”，即小于模拟网格的过程。典型的次网格过程包括：辐射过程，湍流过程，对流过程等。为解决上述问题，科学家们开发了参数化方案来进行弥补。这样，模型最终的表达方式是，在原始数值模型格点的模拟上加入了相应参数化方案的模拟，从而更真实的反应自然现象，得到更精准的模拟结果。

通常一个数值模型包含很多个需要参数化方案来描述的物理过程，以大气科学领域的标志性模型wrf(weatherresearchandforecasting)为例，共有7个物理过程需要使用到参数化方案。这7个物理过程分别为：微物理过程，陆面模型过程，长波辐射过程，短波辐射过程，边界层过程，近地层过程以及积云过程。对于上述每一个物理过程，都有一些可供选择的方案。我们最终的目标是为每一个物理过程选择一个或多个最适合的方案，从而得到七个物理过程最优的参数化方案组合，达到更好的模拟效果。

然而，随着人类对大自然认识的不断加深，模型的结构越来越复杂，每个物理过程里可供选择的方案数也越来越多。每个物理过程中包括的方案从简单到复杂，从运行时间短到需要耗费较长时间，从经历过严谨多次测试的当前业务广泛使用的方案到最新开发鲜有人使用的方案。以wrf模型3.7.1版本为例，该模型共包括24个微物理过程方案，7个陆面模型过程方案，8个长波辐射过程方案，8个短波辐射过程方案，13个边界层过程方案，8个近地层过程方案，和13个积云过程方案。理论上，每个物理过程中的方案都可以和其他物理过程中的方案组合，基于此，总共有24×7×8×8×13×8×13＝14,536,704种排列组合。巨大的方案排列组合数实际上是一个复杂的组合优化问题，而且将所有方案组合全部运行从中选取最好的一个是几乎不可能实现的。如何从数量巨大的排列组合中选出最合适的组合是一个亟待解决的问题。

模型模拟本身含有一定的模拟误差，模型模拟的误差来源主要有以下几方面：第一，输入数据误差；第二，模型结构误差；第三，模型参数误差。上述问题即带来了模型模拟过程中的结构误差。大量前期研究表明选择了不正确的参数化方案带来的模型结构误差会对模型模拟结果造成很大影响，选取合理的参数化方案组合将大大提高模型的模拟效果。因此，为特定的模拟过程选择合适的参数化方案组合是本领域的一个重要的研究内容，也是一个亟待解决的难题。

传统的研究主要是技术人员根据经验人工选择最合适的方案组合，常用的方法是试错法。即通过人工经验尝试不同的参数化方案组合，在多次尝试中选择模拟效果较好的方案组合来完成全部模拟过程。然而，根据经验人工选择最合适的方案组合存在以下缺陷：

1)传统研究通过人工选择合适的方案组合，需要相关专业人员拥有充足的相关经验。一般来说，完成方案手动选择的专业技术人员需要对每个物理过程的每个方案有充分的了解，包括方案的物理机制，该方案适合在哪些特定的区域内使用，以及了解模型模拟值不佳的情况下需要调整哪些方案能够改进模拟值。这个过程需要多年的经验积累，培养专业技术人员需要花费大量时间和金钱。

2)一个模型可能有数十万种方案组合，全部采用手调的方式费时费力。随着模型结构越来越复杂，每个物理过程中包含的参数化方案不断增多，以上述wrf模型为例，共有14,536,704种组合方式，全部采用人工手动选取方案的话，具有很大的工作量，需要耗费很长的时间来完成这个过程。

3)实际研究中，我们可能不止关注模型的一个输出变量。如果同时关注模型的多个输出变量，需要每个变量都进行手调参数人工选出最适合的方案组合，但是对输出变量x综合效果好的参数值不一定对输出变量y的综合模拟效果也好。这种情况下，人工选出的参数化方案组合很难同时顾忌多个输出变量的综合输出效果，从而降低了数值模型的模拟效果。

技术实现要素：

本发明的目的在于解决上述现有技术存在的缺陷，提供一种相比于传统人工手调的方法，能够更加准确及高效的选取数值模型最优方案组合的方法。

一种选取数值模型最优方案组合的方法，包括以下步骤：

步骤1：利用先验知识对方案进行初筛；

步骤2：利用均一化实验设计完成对模型所有方案的抽样；

步骤3：采用多变量方差分析的方法计算每个物理过程对于输出变量的主效应和交互效应，从而得到物理过程的敏感性；所述物理过程的敏感性包括非敏感物理过程和敏感物理过程；

其中，所述主效应，即每个物理过程对模型输出变量的影响；

所述交互效应，即两个或多个物理过程的交互作用对模型输出变量的影响；

步骤4：对于非敏感物理过程，根据经验或模拟效果确定每个组的最优方案；

步骤5：对于敏感物理过程，使用图西检验来检测该过程里每两个方案之间是否有显著差异，如果两个方案有显著差异，则将两个方案分类到两个不同的组里；反之，则将两个方案归为同一个组内；

步骤6：针对敏感物理过程，分类过后，每个组里有一个或多个方案，根据泰勒评分计算每个组内所有方案评价分数的均值，平均值最高的一个组内的方案即为这个物理过程的最优方案；

步骤7：重复步骤5-6，直到每个敏感物理过程都得到它的最优方案；

步骤8：将非敏感物理过程和敏感物理过程中的最优方案排列组合，从而得到模型的最优方案组合。

进一步地，如上所述的选取数值模型最优方案组合的方法，所述步骤2包括：

21)确定好实验所需的因子数和每个因子中的水平数；其中，对于同一个因子，里面每一个水平的重复次数是一致的；

22)根据因子数和水平数选取合适的实验设计。

进一步地，如上所述的选取数值模型最优方案组合的方法，所述方差分析具体计算过程如下：

步骤31：对每个观测值，我们定义误差函数为：

deviation＝individualvalue–overallmean

步骤32：我们定义变量tss，它的计算方式为每个观测值相对于均值的误差的平方和，tss代表了单个观测相对于均值的离差程度；

步骤33：定义自由度df这个变量，主效应中的自由度等于每个因子中的水平数减1，交互效应中的自由度为两个因子自由度的乘积。

步骤34：使用ms来表征每个物理过程对因变量的影响效应，ms的定义为tss除以相应的自由度；

假设两个因子分别为a和b，因子a有a个水平，因子b有b个水平，那么，因子a和因子b的主效应以及他们的交互效应定义如下：

进一步地，如上所述的选取数值模型最优方案组合的方法，所述图西检测具体的工作方式为：

假设一个物理过程a共含有a个参数化方案，那么是每个对应方案的平均泰勒评分值，图西检验检测每一对两两的方案是否有显著差异，图西检测用来进行检测的公式如下所示：

其中是用于比较的两个方案中较大的目标函数均值，是两个方案中较小的目标函数均值，mse指的是两个方案中所有样本点相对于均值的平均误差，n指的是每个方案中的样本数，公式中qs的值需要和学生氏分布中的q值相比较，如果qs值比特定的qα值大，那么两个方案在水平α(0≤α≤1)下是有显著差异的。

进一步地，如上所述的选取数值模型最优方案组合的方法，使用泰勒评分t来表征模型模拟值和真实观测值之间的差异，泰勒评分的表达式如下所示：

其中r＝corr(sim，obs)，泰勒评分的范围是0到1，评分越高，代表模型模拟的越准确。

有益效果：

1.使用多变量方差分析的方法对数值模型中的物理过程进行敏感性分析，筛选出对结果影响大的物理过程。将非敏感的物理过程凭借经验确定为一个固定的最优方案，从而大大减少了总体的方案组合数；

2.采用图西检验将敏感物理过程中的方案划分为不同的组别，从而通过选出综合评价函数值最高的组别选取出对于某个物理过程的最优参数化方案或方案的集合，为最终选取参数化方案的组合提供必要的支持；

3.通过改变评价函数的表达方式可以方便的对多变量同时进行方案优选，使得选取参数化方案组合的时候可以同时考虑多个输出变量，从而达到获得更准确的模型模拟值的目的。

附图说明

图1为本发明选取数值模型最优方案组合的方法流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供的选取数值模型最优方案组合的方法，具体包括以下步骤：

步骤1：利用先验知识对方案进行初筛，以筛选掉不适合的方案。

在实际的模型运行中，科学家可能遇到的问题是某些方案组合搭配在一起不能工作或者某个方案不适用于特定的研究区域。所以，需要利用先验知识对方案进行初筛。对于前者情况，主要通过模型的用户手册进行筛选，对于后者情况，可以通过特定研究区域的特征进行筛选。

步骤2：利用均一化实验设计完成对模型所有方案的抽样

为了满足后续的多变量方差分析，需要进行合理的方案实验设计作为前提。实验设计可以理解为在所有可能的方案组合中进行抽样。抽样的目标是对于同一个因子(物理过程)，里面每一个水平(方案)的重复次数是一致的。这样，满足条件的实验设计才能用于后续方差分析的计算。在本申请中，我们使用均一化实验设计来完成此步骤。

具体地，假设一共有s个因子(即wrf模型中的物理过程数)，且每个因子中的水平数为[q1，q2，…，qs](即每个物理过程中的方案数)，那么方案组合的全部实验设计总数为q1×q2×…×qs。当s和q的值都比较大时，这个组合数是非常庞大的。部分实验设计是另一类实验设计方法，它巧妙的选取全部实验设计的一个子集，使得这个子集能用最少的实验数目达到最大限度表征问题特征的目的。均一化实验设计是部分实验设计中的一种，它的目标是对于同一个因子，里面每一个水平的重复次数是一致的。在实际应用中，大部分均一化实验设计已经被开发者设计好了，只需登录官方网站查找即可。完成均一化实验设计需要两个步骤：1)确定好实验所需的因子数和每个因子中的水平数；2)在官方网站上根据因子数和水平数选取合适的实验设计。均一化实验设计可以实现用相对较少的实验数目来完成方差分析的主效应和交互效应的计算。

步骤3：采用多变量方差分析的方法计算每个物理过程对于输出变量的主效应和交互效应，从而得到物理过程的敏感性；所述物理过程的敏感性包括非敏感物理过程和敏感物理过程；

采用多变量方差分析的方法计算每个物理过程对于输出变量的主效应和交互效应。一般来说，并不是所有的物理过程对于结果影响都很大，这个步骤相当于是对物理过程进行敏感性分析的过程。

具体地，方差分析是一种分析离散的自变量如何影响连续的因变量的统计方法。具体来说，方差分析用来分析离散的因子数(离散变量)对相应因变量(连续变量)的影响。其中，每一个因子里含有特定数量的水平数，不同因子里包含的水平数可以不同。在本申请中，我们关注的wrf模型的输出变量为因变量，物理过程(例如：微物理过程，长波辐射过程，短波辐射过程，陆面过程，边界层过程，近地层过程以及积云过程)为方差分析中的因子，每个物理过程中的方案称为方差分析里面的水平。因为大部分的模型有多于两个物理过程，所以我们使用多变量方差分析而不是单变量或双变量方差分析。在多变量方差分析中，我们关注两个结果：1)主效应，即每个物理过程对模型输出变量的影响；2)交互效应，即两个或多个物理过程的交互作用对模型输出变量的影响。

方差分析的具体计算过程如下：

对每个观测值，我们定义误差函数为：

deviation＝individualvalue–overallmean

进而，我们定义变量tss(totalsumofsquares)，它的计算方式为每个观测值相对于均值的误差的平方和，tss代表了单个观测相对于均值的离差程度。接下来定义自由度(thedegreesoffreedom，df)这个变量。主效应中的自由度等于每个因子中的水平数减1。交互效应中的自由度为两个因子自由度的乘积。接下来，我们使用ms(meansquare)来表征每个物理过程对因变量的影响效应，ms的定义为tss除以相应的自由度。以两个因子的计算过程为例，两个因子分别为a和b，因子a有a个水平，因子b有b个水平。那么，因子a和因子b的主效应以及他们的交互效应定义如下：

步骤4：对于非敏感物理过程，根据经验或模拟效果确定每个组的最优方案；

步骤5：对于敏感物理过程，使用图西检验来检测该过程里每两个方案之间是否有显著差异，如果两个方案有显著差异，则将两个方案分类到两个不同的组里；反之，则将两个方案归为同一个组内；

具体地，对于敏感的物理过程，使用图西检验来检测该过程里每两个方案之间是否有显著差异。如果两个方案有显著差异，则将两个方案分类到两个不同的组里；反之，则将两个方案归为同一个组内。分类过后，每个组里有一个或多个方案。

图西检测用来检测两个方案间是否有显著差异。如果两个方案间没有显著差异，那么他们可以归为同一个类别，反之，则需要归为不同的组别。因此，同一个物理过程中的方案最终被分类为不同的组别。图西检测具体的工作方式为：假设一个物理过程a共含有a个参数化方案，那么是每个对应方案的平均泰勒评分值。图西检验检测每一对两两的方案是否有显著差异，图西检测用来进行检测的公式如下所示：

其中是用于比较的两个方案中较大的目标函数均值，是两个方案中较小的目标函数均值，mse(meansquareerror)指的是两个方案中所有样本点相对于均值的平均误差，n指的是每个方案中的样本数。公式中qs的值需要和学生氏分布中的q值相比较。如果qs值比特定的qα值大，那么我们可以说两个方案在水平α(0≤α≤1)下是有显著差异的。在本申请中，我们把α值设定为0.05，即在图西检验95％的置信区间下，两个方案有/无显著差异。

步骤6：针对敏感物理过程，分类过后，每个组里有一个或多个方案，根据泰勒评分计算每个组内所有方案评价分数的均值，平均值最高的一个组内的方案即为这个物理过程的最优方案；

具体地，根据泰勒评分计算每个组内方案评价分数的均值。平均值最高的一个组内的方案即为这个物理过程的最佳方案。在这里，我们可以通过改变泰勒评分的表达方式方便的考虑关注多个模型输出变量的情形，即将多变量的综合评价融合到泰勒评分中，最终挑选出对多个模型输出变量综合最优的方案组合。

具体地，在这里，我们使用泰勒评分(t)来表征模型模拟值和真实观测值之间的差异，泰勒评分的表达式如下所示：

其中r＝corr(sim，obs)，泰勒评分的范围是0到1，评分越高，代表模型模拟的越准确。在这里，我们可以通过改变泰勒评分的表达方式方便的考虑关注多个模型输出变量的情形，即将多变量的综合评价融合到泰勒评分中，最终挑选出对多个模型输出变量综合最优的方案组合。

步骤7：重复步骤5-6，直到每个敏感物理过程都得到它的最优方案；

步骤8：将非敏感物理过程和敏感物理过程中的最优方案排列组合，从而得到模型的最优方案组合。

具体的，将非敏感物理过程中固定下来的最优方案以及敏感物理过程中筛选出的最优方案进行排列组合，最终得到最优的物理方案组合或物理方案组合的集合。此处，如果最终得到的结果是一个排列组合，则最适合的物理方案为一个排列组合构成的最终结果；如果最终得到的结果是多个排列组合，那么得到的最优物理方案是多个排列组合构成的集合。

实施例

假设一个模型共有五个物理过程a,b,c,d,e，每个物理过程里面有10个可供选择的方案，那么总共有10⁵种方案组合数。

此时，我们先对全部方案进行实验设计，然后用多目标方差分析的方法对物理过程进行敏感性分析，假设得出b和e为非敏感物理过程，a,c,d为敏感物理过程。

对于两个非敏感物理过程确定最优方案的方法：里面的方案无论选择哪个对结果的影响都是很小的，所以本申请通过经验或者比较模拟效果人为地为这个物理过程选择一个最优方案，这样就把总的组合数减少了。此时，由于b和e已经确定出一个最优方案，总组合数变成了10³。

对于敏感物理过程确定最优方案的方法：因为敏感物理过程中方案的选择对结果的影响很大，而且我们很难从10个备选方案里通过某个指标就选择出某个比其他都优的方案，有可能里面的2-3个方案表现差不多，都比剩余的方案要好。所以我们使用图西检验对10个方案进行了分组，选出表现最好的那个组，组里的全部方案都是最优方案。假设敏感物理过程a中10个方案分别为(a1,a2,…,a10)，通过图西检验根据10个方案模拟的相似性进行分组，得到四个组(a1,a3,a5)(a2,a6,a8),(a4),(a7,a9,a10)。随后使用泰勒评分，发现这四个组的评价分数均值为第三组(a4)最高。那么a4即为物理过程a的最优方案。同理，敏感物理过程c分组后得到五个组(c1,c5)，(c2,c8),(c4,c7),(c3,c6),(c9,c10)，其中第四组(c3,c6)的综合评分最高，那么c3和c6即为物理过程c的最优方案。敏感物理过程d分组后得到三个组(d1,d3,d4,d6),(d2,d8,d9,d10),(d5,d7)，其中第三组(d5,d7)的综合评分最高，那么d5和d7即为物理过程d的最优方案。将每个物理过程得到的最优方案排列组合，最终我们得到了四组(1╳1╳2╳2╳1)最优方案排列组合结果。

由于在模型模拟领域，大家普遍认为多个方案组合模拟结果的平均比单个方案组合模拟的效果要好。对于敏感物理过程，它对模拟结果很重要，既然我们不一定能找出里面的某个方案远远好于其他方案，那么就选择最好的2-3个方案都进行模拟。将敏感物理过程中每个组选出的最优方案和非敏感物理选出的每个组里的最优方案进行排列组合，最终得到的即为模型的最优方案组合。

综上，本发明使用严谨的统计学方法完成了数值模型的最优参数化方案组合筛选过程。相比于传统人工手调的方法，本申请所使用的方法更加准确及高效。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。