CN103019704B - 一种用于地球系统模式耦合过程的高阶插值方法 - Google Patents

Abstract

本发明涉及一种用于地球系统模式耦合过程的高阶插值方法，属于地球系统模式耦合技术领域，该方法分为两个部分，第一个部分根据网格A、网格B信息生成给定高阶插值算法的一阶形式的插值系数W′₁，并形成一阶形式的插值系数文件；第二个部分利用网格A上的物理量数据通过插值系数W′₁插值计算得到网格B相应的物理量数据；本发明按照高阶插值算法的计算表达式，结合偏导数系数和高阶插值系数，合成高阶插值算法对应的一阶形式的插值系数，应用该发明，已有的高阶插值算法可以应用于地球系统模式的耦合过程，进而提高插值结果的数值精度，潜在地对整个地球系统模式的模拟结果产生有利的影响。

Description

一种用于地球系统模式耦合过程的高阶插值方法

技术领域

本发明属于地球系统模式耦合技术领域,涉及一种用于地球系统模式耦合过程的高阶插值方法，应用该方法，当前无法应用于地球系统模式耦合过程的高阶插值算法可以得以应用，特别是在维持耦合器插值模块使用方法不变的前提下应用已有的高阶插值算法。

背景技术

地球系统模式是一个在模拟全球气候变化过程中用于全球变化研究的重要工具，主要由耦合器和分量模式构成。其中，分量模式是模拟地球各子系统的数值模拟模块，包括了大气、海洋、海冰、大陆、生物圈等等分量模式，而耦合器是一个用于连接分量模式最终构成地球系统模式的重要组件。在整个变化过程中，能量会在地球各子系统之间传递，分量模式之间需要通过耦合器传递数据，用于各自的数值模拟计算。然而，由于不同的分量模式根据自身的物理模型采用的网格是不同的，在耦合过程中需要通过插值计算对不同网格下的物理量数据进行转换，即将物理量数据从对方分量模式所采用的网格通过插值拟合到本地分量模式所采用的网格。

具体地，分量模式网格由若干互不相交的网格单元组成，每个网格单元有一个中心点坐标，对于给定的物理量数据，每个网格单元包含一个浮点数据，代表网格单元中心点上的物理量数据。在不同网格下，对物理量数据进行插值就是利用网格A的各网格单元中心点上的物理量数据计算网格B的各网格单元中心点上的物理量数据。对于固定的网格A和网格B，需要插值的物理量一般有很多，对每个物理量进行单独处理的开销巨大；因而，用于地球系统模式耦合过程的插值过程被分为两个部分：第一个部分为插值系数生成，即通过插值算法，从源网格以及目标网格生成一组插值系数，插值系数给出利用源网格物理量数据计算目标网格相应物理量数据的与物理量数据无关的系数。该部分方法通常为离线操作，即不用于地球系统模式的运行过程中，生成的插值系数被保存在数据文件里面供耦合器使用；第二个部分为插值计算，在地球系统模式运行过程中，当分量模式之间发生耦合时，耦合器根据相应的插值系数将对方分量模式发送来的物理量数据通过插值计算转化成本地分量模式的物理量数据，该部分主要通过稀疏矩阵乘实现，可以对不同的物理量进行统一处理；两部分之间的关系为：第一部分离线生成的插值系数提供给第二部分在耦合过程进行插值计算时使用。

目前应用于地球系统模式的耦合器均利用一阶形式的插值系数，即插值结果仅线性依赖于网格A上的物理量数据，其特征为：可以用一组与物理量数据无关的插值系数描述给定网格之间的插值，把仅有一组的插值系数称为一阶形式的插值系数。对于给定的网格A和网格B，如果将网格A上的物理量数据记为列向量f，插值到网格B后的物理量数据记为F,那么一阶形式的插值系数可以表示为稀疏矩阵W（矩阵的行数等于网格B的网格单元个数，列数等于网格A的网格单元个数）且满足矩阵运算表达式F＝Wf。由于插值系数与物理量数据无关，因而一阶形式的插值系数的好处在于可以使用同样的插值系数完成网格A和网格B之间不同物理量的插值，另外，不同网格之间、不同物理量的各种一阶插值算法可以使用一致的处理流程，使得耦合器中的插值计算变得简单、一致。

一阶插值算法生成的插值系数都是一阶形式的，可以直接应用于地球系统模式的耦合过程；所谓一阶插值算法，是指其计算表达式中只涉及物理量数据而不涉及物理量偏导数数据。但一阶插值算法也有其不足之处，由于计算插值系数时只使用物理量数据，忽略偏导数数据，使用一阶插值算法得到的插值结果的数值精度比较低，可能无法达到预期目标；特别是对于守恒插值算法，忽略物理量的偏导数数据会明显减低插值结果的数值精度。

设计高精度插值算法时需要考虑物理量的偏导数数据，以便更准确地反映物理量的插值过程，这样的插值算法一般被称为高阶插值算法。采用高阶插值算法进行计算能够增加插值结果的数值精度，进而提高相关模式计算（例如通量计算）结果的数值精度，并可能促进整体模式的预测。然而高阶插值算法一般会生成多组插值系数，每一组插值系数形式上都是一阶形式的插值系数，稀疏矩阵表示为W₁,...，W_i,...，W_N,其中W₁对应于物理量数据f，W_i，i＝2,...，N对应N-1个不同类型的物理量偏导数数据p_i，i＝2,...，N，N为高阶插值系数的个数（取值与所采用的高阶插值算法有关）；高阶插值算法插值结果的计算表达式为相比一阶插值算法，完成高阶插值额外需要物理量偏导数数据p_i，i＝2,...，N。

一阶插值算法已经应用于地球系统模式的耦合过程，如图1所示，按插值的两个部分表示如下：

第一部分：生成一阶形式的插值系数文件

1）根据网格配置文件，调用相应的API，读取网格A和网格B文件中的网格单元信息，包括网格单元个数、维度，网格单元坐标、遮盖数据，储存在相应的数据结构中；

2）利用已读取的网格A和网格B的坐标、遮盖数据，调用选用的一阶插值算法，生成一阶形式的插值系数W₁；调用相应的API将W₁中的数据写入数据文件,作为一阶形式的插值系数文件；

第二部分：利用第一部分得到的一阶形式的插值系数文件和网格A的物理量数据计算网格B的物理量数据

3）根据物理量配置信息，读入网格A上的物理量数据，利用步骤3中得到插值系数文件中的W₁计算得到网格B上的物理量数据，即计算F＝W₁f。

如果不应用于地球系统模式的耦合过程（即无需按插值的两个阶段进行），一般高阶插值算法按照计算表达式，整体流程如图2所示，包括以下步骤：

1）根据网格配置文件，调用相应的API，读取网格A和网格B文件中的网格单元信息，包括网格单元个数、维度，网格单元坐标、遮盖数据，储存在相应的数据结构中；

2）利用已读取的网格A和网格B的坐标、遮盖数据，调用选用的高阶插值算法，生成高阶插值系数W₁，...，W_i,...，W_N并储存在任何一种选定的稀疏矩阵数据结构中；

3）根据物理量配置信息，读入网格A上的物理量数据f；

4）利用网格A的网格单元坐标和物理量数据，使用选用的偏导数计算方法，计算选用的高阶插值算法所需的物理量偏导数数据p_i，i＝2,...，N；

5）按照计算表达式计算插值结果，即

高阶插值算法目前没有广泛地应用于当前的耦合器中，主要难点在于当前的耦合器只能处理一阶形式的插值系数，另外如何高效地计算同一网格上不同物理量的偏导数也给高阶插值算法在耦合过程中的应用带来挑战。直接使用图2中的高阶插值流程，必须修改耦合器的插值模块代码，使其可以计算或以参数形式传入物理量偏导数数据；但这种方法代码难以维护且影响了耦合器插值模块的结构，一般不予考虑。另外，由于同一网格上有不同的物理量，如果步骤4）中的偏导数算法不能高效地计算偏导数，偏导数计算将引入较大的计算开销。

综上，一阶插值算法生成一组插值系数W₁，W₁作为一阶形式的插值系数可以直接应用于地球系统模式的耦合过程；但一阶插值算法的缺点是插值结果的数值精度可能较低。高阶插值算法利用物理量偏导数数据，可以提高插值结果的数值精度，可能对地球系统模式计算产生有利影响；但高阶插值算法的缺点是无法直接给出一阶形式的插值系数，另外如何高效计算同一网格上不同物理量的偏导数数据也是一个难点。耦合器插值模块只能处理一阶形式的插值系数，而高阶插值算法没有直接给出一阶形式的插值系数，因而难以广泛、方便地应用于地球系统模式的耦合过程。

发明内容

本发明的目的是为克服已有技术的不足之处，提出一种用于地球系统模式耦合过程的高阶插值方法，对于给定的高阶插值算法，本方法可以给出一阶形式的插值系数，使得给定的高阶插值算法可以应用于地球系统模式的耦合过程，进而提高插值结果的数值精度。

本发明的特点及有益效果:

本发明采用最小二乘法计算偏导数系数的方法，该方法得到的偏导数系数与物理量数据无关，因而可以高效地计算同一网格内不同物理量的偏导数数据；本发明按照高阶插值算法的计算表达式，结合偏导数系数和高阶插值系数，合成高阶插值算法对应的一阶形式的插值系数。

本发明先计算与物理量数据无关的偏导数系数，使得同一网格内不同物理量的偏导数计算开销大大降低。本发明中的系数合并的方法结合高阶插值的计算表达式，得到一阶形式的高阶插值系数，应用该发明，高阶插值算法可以应用于地球系统模式的耦合过程，进而提高插值结果的数值精度，潜在地对整个地球系统模式的模拟结果产生有利的影响。

附图说明

图1为已有的一阶插值算法在地球系统模式耦合过程中的应用流程图。

图2为已有的高阶插值算法的使用流程图。

图3为本发明的方法流程图。

图4为实施例中的经纬网格。

图5为实施例中网格A的邻居单元。

具体实施方式

本发明提出的一种用于地球系统模式耦合过程的高阶插值方法,基于最小二乘法的偏导数系数计算和系数合并方法，生成已有高阶插值算法（未应用于地球系统模式的耦合过程）的一阶形式的插值系数，使其可以应用于地球系统模式的耦合过程。下面通过附图并结合实施例来阐述本发明的具体实施方式。

本发明提出的一种用于地球系统模式耦合过程的高阶插值方法分为两个部分，第一个部分根据网格A（源网格）、网格B（目标网格）信息生成给定高阶插值算法的一阶形式的插值系数W′₁，并形成一阶形式的插值系数文件，第二个部分利用网格A上的物理量数据通过插值系数文件中的W′₁插值计算得到网格B相应的物理量数据；

第一部分具体包括以下步骤：

1）根据网格配置文件，调用相应的API，读取网格A和网格B文件中的网格单元信息，包括网格单元个数、维度，网格单元坐标、遮盖数据，储存在相应的数据结构中；

2）根据读取的网格A文件中的网格信息对网格A中的每个网格单元进行搜索得到每个网格单元的多个邻居单元；具体包括：

2-1）读取已储存在数据结构中的网格A文件中的网格信息，包括网格单元个数、维度、顶点坐标、中心点坐标、遮盖；

2-2）若网格单元维度为2（即网格A是局部矩形网格），根据网格单元的二维索引，在二维索引空间中找到网格单元的网格邻居；

2-3）若网格单元维度为1（即网格A是无结构网格），根据网格单元中心点坐标将网格A划分为若干区域，然后对于每个网格单元，计算其与指定区域内网格单元的距离，按距离从小到大排序，选取排序位置处于前面的多个（预先指定的个数）相应网格单元，作为网格邻居；

2-4）利用步骤2-2）或2-3）的结果，依次遍历网格单元的网格邻居及网格邻居的网格邻居组成邻居单元，使得到的邻居单元数目超过预定个数（预定个数是与选用的高阶插值算法相应数目）；

3）利用已读取的网格A和网格B的坐标、遮盖数据，调用选用的高阶插值算法，生成高阶插值系数W₁,...，W_i,...，W_N并储存在任何一种选定的稀疏矩阵数据结构中，N为高阶插值系数的个数（取值与所采用的高阶插值算法有关）；

4）利用网格A中网格单元的邻居单元及其坐标数据，基于最小二乘法，计算网格A的每个网格单元的各个偏导数系数,具体包括:

4-1）根据步骤2）得到的网格A中的网格单元的若干个邻居单元，组成邻居集合Neighbor(a)，对邻居集合Neighbor(a)中的每个邻居单元j∈Neighbor(a)，将网格单元j的物理量数据在网格单元a的中心点进行Taylor展开如式（1）：

$f_j=f_a+{\mathrm{\Σ}}_{i=2}^N{p}_{i,a}{\mathrm{\λ}}_i({\mathrm{\θ}}_j-{\mathrm{\θ}}_a,{\mathrm{\φ}}_j-{\mathrm{\φ}}_a)+o_{j,a,N}---\left(1\right),$

展开式（1）中p_i，a为网格单元a上选定高阶插值算法所需的所有偏导数，θ_j，φ_j为邻居单元j的中心点坐标，θ_a，φ_a为网格单元a的中心点坐标,λ_i(θ_j-θ_a，φ_j-φ_a)为相应偏导数对应的表达式，o_j，a，N为展开式中的高阶小量。展开式中所有偏导数p_i,a,i＝2,...，N为待求解未知数，按照步骤2），展开式的个数大于待求解未知数的个数N-1；

4-2）基于最小二乘法，将高阶小量的平方和最小作为最优化目标；Δ是待求解未知数p_i,a,i＝2,...,N的二次多项式，多项式系数为常数（网格单元坐标是常数）且二次项系数为正，求解多元一次方程组i＝2,...，N即可得到所有偏导数p_i,a,i＝2,...，N；最终表达式形如式（2）：

p_i，a＝ω_a,i,af_a+∑_{j∈Neighbor(a)}ω_a，i，jf_j（2）；

4-3）从步骤4-2）的计算表达式p_i,a,a∈A，i＝2,...，N中提取偏导数系数ω_a，i，a和ω_a，i，j，对于每个偏导数p_i，将ω_a，i，a，ω_a，i，j，a∈A，j∈Neighbor(a)构成的稀疏矩阵G_i储存在给定的稀疏矩阵数据结构中,使得步骤4-2）的表达式转化为：p_i＝G_if；

5）将高阶插值系数和偏导数系数合并生成一阶形式的插值系数文件

5-1）结合高阶插值算法的计算表达式和4-3）中的计算表达式，利用矩阵运算给出一阶形式的插值系数矩阵如式（3）：即

5-2）按照插值系数文件的标准格式，调用相应的API将W′₁中的数据写入数据文件,作为一阶形式的插值系数文件；

第二部分：利用第一部分得到的一阶形式的插值系数文件和网格A的物理量数据计算网格B的物理量数据，具体为：

6）根据物理量配置信息，读入网格A上的物理量数据f，利用步骤5）中得到的插值系数文件中的W′₁计算得到网格B上的物理量数据F，即计算F＝W′₁f。

实施例

下面以SCRIP二阶守恒插值算法为例，分为两个部分给出本发明方法的实施例，第一个部分根据网格A、网格B信息生成SCRIP二阶守恒插值算法的一阶形式的插值系数W′₁，第二个部分利用网格A上的物理量数据通过插值系数W′₁插值计算得到网格B相应的物理量数据。在本实施例中，网格A和网格B是经纬网格，经纬网格的特征为球面被经线和纬线划分为若干网格单元，如图4所示。具体实施步骤如下：

第一部分：生成一阶形式的插值系数文件

1）根据网格配置文件，调用netCDF文件API，读取网格A和网格B文件中的网格单元信息，包括网格单元个数、维度，网格单元坐标、遮盖数据，储存在相应的数据结构中；

2）根据读取的网格A文件中的网格信息对网格A中的每个单元进行搜索得到每个单元的多个邻居单元；具体包括：

2-1）读取已储存在数据结构中的网格A文件中的网格信息，包括网格单元个数、维度、顶点坐标、中心点坐标、遮盖；

2-2）根据网格单元的二维索引，在二维索引空间中找到网格单元的8个或5个网格邻居；如图5所示，对于不靠近极区的网格单元a，找到上下左右及斜线方向共8个网格邻居，即a₁,a₂,…a₈,；对于靠近北极点的网格单元，找到共远离极区的5个网格邻居，对应图5中的a₁,a₂,a₃，a₄和a₈；对于靠近南极点的网格单元，找到远离极区的5个网格邻居，对应图5中的a₄,a₅,a₆，a₇和a₈；

2-3）步骤2-2）中每个网格单元的网格邻居个数超过该算法涉及的偏导数类型个数，步骤2-2）中得到的网格邻居即为所需邻居单元；

3）利用已读取的网格A和网格B的坐标、遮盖数据，调用SCRIP二阶守恒插值算法，生成高阶插值系数W₁，W₂，W₃并储存在任何一种选定的稀疏矩阵数据结构中；

4）利用网格A中网格单元的邻居单元及其坐标数据，基于最小二乘法，计算网格A的每个网格单元的各个偏导数系数,具体包括:

4-1）对网格A中的网格单元a，步骤2）得到a的8或5个邻居单元，组成集合Neighbor(a)，对Neighbor(a)中的每个网格单元j∈Neighbor(a)，将网格单元j的物理量数据在网格单元a的中心点进行Taylor展开，

f_j＝f_a+p_1，a(θ_j-θ_a)+p_1，a(φ_j-φ_a)+o_j，a，3,

展开式中p_i，a为网格单元a上SCRIP二阶守恒插值算法所需的两个偏导数，θ_j，φ_j为网格单元j的中心点坐标，θ_a，φ_a为网格单元a的中心点坐标，o_j，a，3为展开式中的高阶小量。展开式中p_1,a,p_2，a为待求解未知数，按照步骤2），展开式的个数（8或5）大于待求解未知数的个数2；

4-2）基于最小二乘法，将高阶小量的平方和最小作为最优化目标；Δ是待求解未知数p_1,a,p_2，a的二次多项式，多项式系数为常数（网格单元坐标是常数）且二次项系数为正，求解多元一次方程组i＝2，3即可得到p_1,a,p_2，a；最终表达式形如p_i，a＝ω_a，i，af_a+∑_{j∈Neighbor(a)}ω_a，i，jf_j；

4-3）从步骤4-2）的计算表达式p_i，a，a∈A，i＝2，3中提取偏导数系数ω_a，i，a和ω_a，i，j，对于每个i，将ω_a，i,a,ω_a，i，j，a∈A，j∈Neighbor(a)构成的稀疏矩阵G_i储存在给定的稀疏矩阵数据结构中,使得步骤4-2）的表达式转化为p_i＝G_if；

5）将高阶插值系数和偏导数系数合并生成一阶形式的插值系数文件

5-1）结合高阶插值算法的计算表达式和4-3）中的计算表达式，利用矩阵运算给出一阶形式的插值系数矩阵，即

5-2）按照插值系数文件的标准格式，调用相应的netCDFAPI将W′₁中的数据写入数据文件,作为一阶形式的插值系数文件；

第二部分：利用第一部分得到的一阶形式的插值系数文件和网格A的物理量数据计算网格B的物理量数据

6）根据物理量配置信息，读入网格A上的物理量数据f，利用步骤5）中得到的插值系数文件中的W′₁计算得到网格B上的物理量数据F，即计算F＝W′₁f。

本发明给出基于最小二乘法和系数合并方法的高阶插值方法，区别技术特征为：首先，给出基于最小二乘法的偏导数系数计算方法，生成的偏导数系数与物理量数据无关，使得同一网格内不同物理量的偏导数计算开销大大降低；其次，通过合并偏导数系数与高阶插值系数，得到一阶形式的高阶插值系数，使得高阶插值算法在经过本方法的流程后可以应用于地球系统模式的耦合过程，其中无须对耦合器插值模块代码进行修改。

本发明给出了一种针对地球系统模式耦合过程的高阶插值方法，该方法基于最小二乘法和系数合并方法，使得当前无法广泛应用于耦合过程的高阶插值算法得以应用，可以大大提高耦合过程中插值结果的数值精度，对整个模式模拟结果具有潜在的有利影响。任何熟悉该技术的人在本发明所揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护之内。

Claims (3)

1.一种用于地球系统模式耦合过程的高阶插值方法，其特征在于，该方法分为两个部分，第一个部分根据网格A、网格B信息生成给定高阶插值算法的一阶形式的插值系数W′₁，并形成一阶形式的插值系数文件；第二个部分利用网格A上的物理量数据通过插值系数W′₁插值计算得到网格B相应的物理量数据；

第一部分具体包括以下步骤：

1)根据网格配置文件，调用相应的API，读取网格A和网格B文件中的网格单元信息，包括网格单元个数、维度，网格单元坐标、遮盖数据，储存在相应的数据结构中；

2)根据读取的网格A文件中的网格单元信息，对网格A中的每个网格单元进行搜索得到每个网格单元的多个邻居单元；

3)利用已读取的网格A和网格B的网格单元坐标、遮盖数据，调用选用的高阶插值算法，生成高阶插值系数W₁,…,W_i,…,W_N并储存在任何一种选定的稀疏矩阵数据结构中，N为高阶插值系数的个数；

4)利用网格A中网格单元的邻居单元及其坐标数据，基于最小二乘法，计算网格A的每个网格单元的各个偏导数系数；

5)将高阶插值系数和偏导数系数合并生成一阶形式的插值系数文件；

第二部分具体为：

6)根据物理量配置信息，读入网格A上的物理量数据f，利用步骤5)中得到的插值系数文件中的W′₁计算得到网格B上的物理量数据F，即计算F＝W′₁f；

所述步骤4)具体包括:

4-1)根据步骤2)得到的网格A中的网格单元a的若干个邻居单元，组成邻居集合Neighbor(a)，对邻居集合Neighbor(a)中的每个邻居单元j∈Neighbor(a)，将网格单元j的物理量数据在网格单元a的中心点进行Taylor展开如式(1)：

展开式(1)中p_i,a为网格单元a上选定高阶插值算法所需的所有偏导数，θ_j,φ_j为邻居单元j的中心点坐标，θ_a,φ_a为网格单元a的中心点坐标,λ_i(θ_j-θ_a,φ_j-φ_a)为相应偏导数对应的表达式，o_j，a，N为展开式中的高阶小量；展开式中所有偏导数p_i，a,i＝2,…,N为待求解未知数，按照步骤2)，展开式的个数大于待求解未知数的个数N-1；

4-2)基于最小二乘法，将高阶小量的平方和最小作为最优化目标；Δ是待求解未知数p_i，a,i＝2,…,N的二次多项式，多项式系数为常数且二次项系数为正，求解多元一次方程组即可得到所有偏导数p_i，a,i＝2,…,N；最终表达式形如式(2)：

p_i，a＝ω_a，i,af_a+∑_{j∈Neighbor(a)}ω_a,i,jf_j(2)；

4-3)从步骤4-2)的计算表达式p_i,a,a∈A,i＝2,…,N中提取偏导数系数ω_a,i,a和ω_a,i,j，对于每个偏导数p_i，将ω_a,i,a,ω_a,i,j,a∈A,j∈Neighbor(a)构成的稀疏矩阵G_i储存在给定的稀疏矩阵数据结构中,使得步骤4-2)的表达式(2)转化为：p_i＝G_if。

2.如权利要求1所述方法，其特征在于，所述步骤2)具体包括：

2-1)读取已储存在数据结构中的网格A文件中的网格单元信息，包括网格单元个数、维度、顶点坐标、中心点坐标、遮盖；

2-2)若网格单元维度为2，根据网格单元的二维索引，在二维索引空间中找到网格单元的网格邻居；

2-3)若网格单元维度为1，根据网格单元中心点坐标将网格A划分为若干区域，然后对于每个网格单元，计算其与指定区域内网格单元的距离，按距离从小到大排序，选取排序位置处于前面的多个相应网格单元，作为网格邻居；

2-4)利用步骤2-2)或2-3)的结果，依次遍历网格单元的网格邻居及网格邻居的网格邻居组成邻居单元，使得到的邻居单元数目超过预定个数。

3.如权利要求2所述方法，其特征在于，所述步骤5)具体包括:

5-1)结合高阶插值算法的计算表达式和4-3)中的计算表达式，利用矩阵运算给出一阶形式的插值系数矩阵如式(3)：即

5-2)按照插值系数文件的标准格式，调用相应的API将W′₁中的数据写入数据文件,作为一阶形式的插值系数文件。