基于KD树的非结构网格地球系统模式观测稀疏化方法

技术领域

本发明属于计算机和气象、海洋的交叉技术领域，尤其涉及基于KD树的非结构网格地球系统模式观测稀疏化方法。

背景技术

地球系统模式采用数值模拟的手段预报地球未来的大气、海洋、海冰、陆面等演变状态，对气候变化、防灾减灾、经济活动等具有重要意义。地球系统模式的预报效果依赖于初始场的质量，资料同化能够将模式背景场与观测资料相融合，从而改进初始场。为确保资料同化的有效性，当观测的空间分布相较于模式格点更为密集时，需要采用观测稀疏化方法去除观测冗余信息和观测误差相关性。

观测稀疏化方法在每个模式网格中保留至多一个观测，排除冗余观测和位于模式格点之外的观测。对于规则结构网格，模式网格采用二维数组方式存储，模式网格在东西方向和南北方向以等间距分布，能够高效实现观测稀疏化。对于非结构网格，模式网格采用一维数组方式存储，模式网格是无序排列的，目前的高效观测稀疏化方法还很少。然而，非结构网格地球系统模式能够灵活地贴合复杂不规则几何形状的岸线和陡峭的底部地形，且能够针对重点区域进行加密而无需嵌套，这些是结构网格不具备的优势。随着中小尺度过程研究的不断发展与进步，以及超级计算机计算能力不断提高带来的模式分辨率的飞速提升，工程实践对非结构网格地球系统模式观测稀疏化方法的需求也越来越迫切。

现有的观测稀疏化方法对于非结构网格通常采用穷举法。假设共有M个观测、N个模式网格和K个模式网格顶点，每个模式网格为P边形，即具有P个顶点。穷举法首先读取观测

可见传统的观测稀疏化方法需要对观测数M和模式网格数N进行两层循环，计算时间较长，计算效率较低。而目前典型的全球地球系统模式网格点数和观测数分别可达10

发明内容

如何提供高效、易行、稳定的观测稀疏化方法是本领域急需解决的技术问题。有鉴于此，本发明提出了基于KD树的非结构网格地球系统模式观测稀疏化方法。

本发明公开的基于KD树的非结构网格地球系统模式观测稀疏化方法，包括以下步骤：

S1：获取M个观测位置和N个模式网格位置数据数据，并对数据进行预处理；

S2：初始化模式网格保留观测序号的存储数组NE

S3：根据模式网格中心点E

S4：针对观测

S5：根据模式网格中心点E

S6：根据模式网格顶点的位置

S7：循环到i+1执行S4步，直到i＝1,…,M循环完成，执行第S8步；

S8：遍历模式网格保留观测序号的存储数组NE

进一步的，获取观测位置和模式网格位置数据数据包括：读取观测

进一步的，所述对数据进行预处理，包括：

将观测和模式网格的经度范围调整一致，统一为[-180°W,180°E]；

若观测或模式网格为循环边界，则复制经度范围在[-180°W,(-180+h)°W]之间的观测并将它们的经度值加360，复制经度范围在[(180-h)°E,180°E]之间的观测并将它们的经度值减360，其中h为指定范围。

进一步的，步骤S3包括：

S301：判断待划分的模式网格中心点集

S302：设置分割维度D，若KD树为空则设置分割维度D为经度，否则设置当前步分割维度与上一步分割维度不同，即将分割维度D交替设置为经度、纬度；

S303：选用快速排序法对模式网格点的D维度进行排序，选取D维度上位于中位数的模式网格中心点，构造KD树的当前节点；

S304：根据当前节点，将剩余的待划分点放入左子树点集和右子树的待划分点集中，循环执行步骤S301。

进一步的，所述将剩余的待划分点放入左子树点集和右子树的待划分点集中，具体规则为将D维度值小于当前节点D维度值的待划分点放入左子树中，其他的则放入右子树中。

进一步的，步骤S4包括：

S401：以KD树T为初始节点，设置当前节ND＝T；

S402：判断当前节点ND是否为空，为空则结束搜索，否则继续搜索；

S403：计算当前节点ND与观测位置

S404：判断当前节点ND的分割维度D，计算当前节点ND与观测位置

S405：当搜索到叶子节点时，搜索过程开始回溯；判断观测位置

本发明的有益效果如下：

本发明与传统的穷举法相比，加速效果显著，能够剔除模式网格外的观测，且能够有效去除观测冗余信息。

附图说明

图1本发明的整体流程图；

图2本发明的详细流程图；

图3本发明的构建KD树流程图；

图4本发明查找KD树中与观测距离在指定范围内的模式网格中心点的流程图。

具体实施方式

下面结合附图对本发明作进一步的说明，但不以任何方式对本发明加以限制，基于本发明教导所作的任何变换或替换，均属于本发明的保护范围。

本发明要解决的技术问题是针对现有非结构网格观测稀疏化方法存在的计算效率较低问题，提供一种高效、易行、稳定的观测稀疏化方法。当模式网格数N远大于log

本发明的技术方案如图1所示：首先读取模式网格和观测的位置信息并处理循环边界，其次针对非结构模式网格构建KD树，然后根据KD树查找观测所在的模式网格并保存其相对于模式网格中心点的距离，最后保留距模式网格中心点最近的观测。

参考图2，本发明具体技术方案包括以下步骤：

1.读取观测

2.将观测和模式网格的经度范围调整一致，统一为[-180°W,180°E]；

3.若观测或模式网格为循环边界，则复制经度范围在[-180°W,(-180+h)°W]之间的观测并将它们的经度值加360，复制经度范围在[(180-h)°E,180°E]之间的观测并将它们的经度值减360，其中h为指定范围。。

4.初始化模式网格保留观测序号的存储数组NE

5.根据模式网格中心点E

5.1.判断待划分的模式网格中心点集

5.2.设置分割维度D，若KD树为空则设置分割维度D为经度，否则设置当前步分割维度与上一步分割维度不同，即将分割维度D交替设置为经度、纬度；

5.3.选用快速排序法对模式网格点的D维度进行排序，选取D维度上位于中位数的模式网格中心点，构造KD树的当前节点；

5.4.根据当前节点，将剩余的待划分点放入左子树点集和右子树的待划分点集中。具体规则为将D维度值小于当前节点D维度值的待划分点放入左子树中，其他的则放入右子树中。循环执行第5.1步；

6.针对观测

6.1.以KD树T为初始节点，设置当前节ND＝T；

6.2.判断当前节点ND是否为空，为空则结束搜索，否则继续搜索；

6.3.计算当前节点ND与观测位置

6.4.判断当前节点ND的分割维度D，计算当前节点ND与观测位置

6.5.当搜索到叶子节点时，搜索过程开始回溯。判断观测位置

7.根据模式网格中心点E

8.根据模式网格顶点的位置

9.循环到i+1，执行第6步，直到i＝1,…,M循环完成，执行第10步；

10.遍历模式网格保留观测序号的存储数组NE

11.结束。

实施例

1.从网站

2.将观测和模式网格的经度范围调整一致，统一为[-180°W,180°E]；

3.若观测或模式网格为循环边界，则复制经度范围在[-180°W,(-180+h)°W]之间的观测并将它们的经度值加360，复制经度范围在[(180-h)°E,180°E]之间的观测并将它们的经度值减360，其中h为指定范围。

4.初始化模式网格保留观测序号的存储数组NE

5.根据模式网格中心点E

5.1.判断待划分的模式网格中心点集

5.2.设置分割维度D，若KD树为空则设置分割维度D为经度，否则设置当前步分割维度与上一步分割维度不同，即将分割维度D交替设置为经度、纬度；

5.3.选用快速排序法对模式网格点的D维度进行排序，选取D维度上位于中位数的模式网格中心点，构造KD树的当前节点；

5.4.根据当前节点，将剩余的待划分点放入左子树点集和右子树的待划分点集中。具体规则为将D维度值小于当前节点D维度值的待划分点放入左子树中，其他的则放入右子树中。循环执行第5.1步；

6.针对观测

6.1.以KD树T为初始节点，设置当前节ND＝T；

6.2.判断当前节点ND是否为空，为空则结束搜索，否则继续搜索；

6.3.计算当前节点ND与观测位置

6.4.判断当前节点ND的分割维度D，计算当前节点ND与观测位置

6.5.当搜索到叶子节点时，搜索过程开始回溯。判断观测位置

7.根据模式网格中心点E

8.根据模式网格顶点的位置

9.循环到i+1，执行第三步，直到i＝1,…,M循环完成，执行第10步；

10.遍历模式网格保留观测序号的存储数组NE

11.结束。

对传统的穷举法和基于KD树的观测稀疏化方法进行测试。FVCOM全球海洋模式采用非结构的三角网格，其网格中心点数为688991，网格顶点数为359120。针对全球1/20°的海表温度观测资料进行稀疏化，除去陆地、湖泊、河流、海冰点后，观测数为14160103。采用湖南大学国家超级计算长沙中心的天河一号，计算节点主频为2.3GHz，内存为48G。仅使用单核，传统的穷举法处理1个观测的计算时间约为218s，基于KD树的非结构网格地球系统模式观测稀疏化方法处理14160103个观测的计算时间约为15s，加速效果显著。基于KD树的非结构网格地球系统模式观测稀疏化方法能够剔除模式网格外的观测，且能够有效去除观测冗余信息。

本发明的有益效果如下：

本发明与传统的穷举法相比，加速效果显著，能够剔除模式网格外的观测，且能够有效去除观测冗余信息。

本文所使用的词语“优选的”意指用作实例、示例或例证。本文描述为“优选的”任意方面或设计不必被解释为比其他方面或设计更有利。相反，词语“优选的”的使用旨在以具体方式提出概念。如本申请中所使用的术语“或”旨在意指包含的“或”而非排除的“或”。即，除非另外指定或从上下文中清楚，“X使用A或B”意指自然包括排列的任意一个。即，如果X使用A；X使用B；或X使用A和B二者，则“X使用A或B”在前述任一示例中得到满足。

而且，尽管已经相对于一个或实现方式示出并描述了本公开，但是本领域技术人员基于对本说明书和附图的阅读和理解将会想到等价变型和修改。本公开包括所有这样的修改和变型，并且仅由所附权利要求的范围限制。特别地关于由上述组件(例如元件等)执行的各种功能，用于描述这样的组件的术语旨在对应于执行所述组件的指定功能(例如其在功能上是等价的)的任意组件(除非另外指示)，即使在结构上与执行本文所示的本公开的示范性实现方式中的功能的公开结构不等同。此外，尽管本公开的特定特征已经相对于若干实现方式中的仅一个被公开，但是这种特征可以与如可以对给定或特定应用而言是期望和有利的其他实现方式的一个或其他特征组合。而且，就术语“包括”、“具有”、“含有”或其变形被用在具体实施方式或权利要求中而言，这样的术语旨在以与术语“包含”相似的方式包括。

本发明实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以多个或多个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。上述提到的存储介质可以是只读存储器，磁盘或光盘等。上述的各装置或系统，可以执行相应方法实施例中的存储方法。

综上所述，上述实施例为本发明的一种实施方式，但本发明的实施方式并不受所述实施例的限制，其他的任何背离本发明的精神实质与原理下所做的改变、修饰、代替、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。