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

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


背景技术：

2.地球系统模式采用数值模拟的手段预报地球未来的大气、海洋、海冰、陆面等演变状态，对气候变化、防灾减灾、经济活动等具有重要意义。地球系统模式的预报效果依赖于初始场的质量，资料同化能够将模式背景场与观测资料相融合，从而改进初始场。为确保资料同化的有效性，当观测的空间分布相较于模式格点更为密集时，需要采用观测稀疏化方法去除观测冗余信息和观测误差相关性。
3.观测稀疏化方法在每个模式网格中保留至多一个观测，排除冗余观测和位于模式格点之外的观测。对于规则结构网格，模式网格采用二维数组方式存储，模式网格在东西方向和南北方向以等间距分布，能够高效实现观测稀疏化。对于非结构网格，模式网格采用一维数组方式存储，模式网格是无序排列的，目前的高效观测稀疏化方法还很少。然而，非结构网格地球系统模式能够灵活地贴合复杂不规则几何形状的岸线和陡峭的底部地形，且能够针对重点区域进行加密而无需嵌套，这些是结构网格不具备的优势。随着中小尺度过程研究的不断发展与进步，以及超级计算机计算能力不断提高带来的模式分辨率的飞速提升，工程实践对非结构网格地球系统模式观测稀疏化方法的需求也越来越迫切。
4.现有的观测稀疏化方法对于非结构网格通常采用穷举法。假设共有m个观测、n个模式网格和k个模式网格顶点，每个模式网格为p边形，即具有p个顶点。穷举法首先读取观测(i＝1,
…
,m)的经纬度位置模式网格gj中心点ej(j＝1,
…
,n)的经纬度位置模式网格顶点nk(k＝1,
…
,k)的经纬度位置(lonk,latk)，以及模式网格中心点与顶点的对应关系r
jp
(j＝1,
…
,n,p＝1,
…
,p)；然后，对于每一个观测的经纬度位置遍历模式网格gj(j＝1,
…
,n)，获取该模式网格对应的顶点r
jp
(p＝1,
…
,p)，则可获取这些顶点的位置依此即可判断观测是否位于gj内部，若位于内部则记录下来并处理下一个观测，否则循环到下一个模式网格，若不位于任一模式网格内，则直接排除最后，循环处理完所有观测(i＝1,
…
,m)后，查看每个模式网格gj(j＝1,
…
,n)中是否有多个观测，若观测数目大于1，则计算该模式网格中所有观测与模式网格中心点的距离，保留距离最小的观测，其余观测都排除。
5.可见传统的观测稀疏化方法需要对观测数m和模式网格数n进行两层循环，计算时间较长，计算效率较低。而目前典型的全球地球系统模式网格点数和观测数分别可达106～107，传统的观测稀疏化方法已经难以满足当前业务预报系统的预报时效要求。


技术实现要素：

6.如何提供高效、易行、稳定的观测稀疏化方法是本领域急需解决的技术问题。有鉴于此，本发明提出了基于kd树的非结构网格地球系统模式观测稀疏化方法。
7.本发明公开的基于kd树的非结构网格地球系统模式观测稀疏化方法，包括以下步骤：
8.s1：获取m个观测位置和n个模式网格位置数据数据，并对数据进行预处理；
9.s2：初始化模式网格保留观测序号的存储数组nej＝-1(j＝1,
…
,n)和模式网格保留观测的距离数组ndj＝-1(j＝1,
…
,n)；
10.s3：根据模式网格中心点ej(j＝1,
…
,n)的经纬度位置构建kd树t；
11.s4：针对观测的经纬度位置查找kd树t中与之距离在指定范围h内的模式网格中心点，得到l个模式网格中心点e
l
(l＝1,
…
,l)；
12.s5：根据模式网格中心点e
l
和模式网格中心点与p个顶点的对应关系r
lp
(p＝1,
…
,p)，获取模式网格顶点的位置
13.s6：根据模式网格顶点的位置判断观测是否位于模式网格中心点e
l
所对应的模式网格g
l
中，若是则计算观测位置与模式网格中心点e
l
位置的距离d：若ne
l
《0则令ne
l
＝i且nd
l
＝d，否则判断d与nd
l
的大小，若d《nd
l
则令ne
l
＝i且nd
l
＝d；否则循环到l+1，执行第s5步，直到l＝1,
…
,l循环完成，执行第s7步；
14.s7：循环到i+1执行s4步，直到i＝1,
…
,m循环完成，执行第s8步；
15.s8：遍历模式网格保留观测序号的存储数组nej(j＝1,
…
,n)，nej》0中保存的观测序号即为观测稀疏化后保留的观测序号。
16.进一步的，获取观测位置和模式网格位置数据数据包括：读取观测(i＝1,
…
,m)的经纬度位置模式网格gj中心点ej(j＝1,
…
,n)的经纬度位置模式网格顶点nk(k＝1,
…
,k)的经纬度位置(lonk,latk)，以及模式网格中心点与顶点的对应关系r
jp
(j＝1,
…
,n,p＝1,
…
,p)。
17.进一步的，所述对数据进行预处理，包括：
18.将观测和模式网格的经度范围调整一致，统一为[-180
°
w,180
°
e]；
[0019]
若观测或模式网格为循环边界，则复制经度范围在[-180
°
w,(-180+h)
°
w]之间的观测并将它们的经度值加360，复制经度范围在[(180-h)
°
e,180
°
e]之间的观测并将它们的经度值减360，其中h为指定范围。
[0020]
进一步的，步骤s3包括：
[0021]
s301：判断待划分的模式网格中心点集是否为空集，若是，则结束构建过程并执行s4步骤，否则开始构建过程；
[0022]
s302：设置分割维度d，若kd树为空则设置分割维度d为经度，否则设置当前步分割维度与上一步分割维度不同，即将分割维度d交替设置为经度、纬度；
[0023]
s303：选用快速排序法对模式网格点的d维度进行排序，选取d维度上位于中位数
的模式网格中心点，构造kd树的当前节点；
[0024]
s304：根据当前节点，将剩余的待划分点放入左子树点集和右子树的待划分点集中，循环执行步骤s301。
[0025]
进一步的，所述将剩余的待划分点放入左子树点集和右子树的待划分点集中，具体规则为将d维度值小于当前节点d维度值的待划分点放入左子树中，其他的则放入右子树中。
[0026]
进一步的，步骤s4包括：
[0027]
s401：以kd树t为初始节点，设置当前节nd＝t；
[0028]
s402：判断当前节点nd是否为空，为空则结束搜索，否则继续搜索；
[0029]
s403：计算当前节点nd与观测位置之间的距离d1，若d1小于等于指定范围h，则记录当前节点nd为符合要求的模式网格中心点e
l
；
[0030]
s404：判断当前节点nd的分割维度d，计算当前节点nd与观测位置在分割维度d上的距离d2，若d2》h意味着观测位置位于nd的左子树空间内，与左子树空间点的距离可能比与右子树空间点的距离更小，因此优先递归搜索nd的左子树，反之，优先搜索nd的右子树；搜索循环执行s402步骤；
[0031]
s405：当搜索到叶子节点时，搜索过程开始回溯；判断观测位置与当前节点nd在分割维度d上的距离是否大于指定范围h，若是，意味着当前节点nd的另一子树空间内不可能存在更近的点，则继续回溯；否则，需要执行s402步骤进一步搜索当前节点nd的另一子树空间，直到回溯到kd树根节点t结束。
[0032]
本发明的有益效果如下：
[0033]
本发明与传统的穷举法相比，加速效果显著，能够剔除模式网格外的观测，且能够有效去除观测冗余信息。
附图说明
[0034]
图1本发明的整体流程图；
[0035]
图2本发明的详细流程图；
[0036]
图3本发明的构建kd树流程图；
[0037]
图4本发明查找kd树中与观测距离在指定范围内的模式网格中心点的流程图。
具体实施方式
[0038]
下面结合附图对本发明作进一步的说明，但不以任何方式对本发明加以限制，基于本发明教导所作的任何变换或替换，均属于本发明的保护范围。
[0039]
本发明要解决的技术问题是针对现有非结构网格观测稀疏化方法存在的计算效率较低问题，提供一种高效、易行、稳定的观测稀疏化方法。当模式网格数n远大于log2n，基于kd树的观测稀疏化方法搜索效率远高于穷举法，且不受模式网格类型的影响，应用价值高、适用范围广。
[0040]
本发明的技术方案如图1所示：首先读取模式网格和观测的位置信息并处理循环
边界，其次针对非结构模式网格构建kd树，然后根据kd树查找观测所在的模式网格并保存其相对于模式网格中心点的距离，最后保留距模式网格中心点最近的观测。
[0041]
参考图2，本发明具体技术方案包括以下步骤：
[0042]
1.读取观测(i＝1,
…
,m)的经纬度位置模式网格gj中心点ej(j＝1,
…
,n)的经纬度位置模式网格顶点nk(k＝1,
…
,k)的经纬度位置(lonk,latk)，以及模式网格中心点与顶点的对应关系r
jp
(j＝1,
…
,n,p＝1,
…
,p)；
[0043]
2.将观测和模式网格的经度范围调整一致，统一为[-180
°
w,180
°
e]；
[0044]
3.若观测或模式网格为循环边界，则复制经度范围在[-180
°
w,(-180+h)
°
w]之间的观测并将它们的经度值加360，复制经度范围在[(180-h)
°
e,180
°
e]之间的观测并将它们的经度值减360，其中h为指定范围。。
[0045]
4.初始化模式网格保留观测序号的存储数组nej＝-1(j＝1,
…
,n)和模式网格保留观测的距离数组ndj＝-1(j＝1,
…
,n)，-1表明尚未找到位于模式网格中的观测；
[0046]
5.根据模式网格中心点ej(j＝1,
…
,n)的经纬度位置构建kd树t，如图3所示，具体方法是：
[0047]
5.1.判断待划分的模式网格中心点集是否为空集，若是则结束构建过程并执行第6步，否则开始构建过程；
[0048]
5.2.设置分割维度d，若kd树为空则设置分割维度d为经度，否则设置当前步分割维度与上一步分割维度不同，即将分割维度d交替设置为经度、纬度；
[0049]
5.3.选用快速排序法对模式网格点的d维度进行排序，选取d维度上位于中位数的模式网格中心点，构造kd树的当前节点；
[0050]
5.4.根据当前节点，将剩余的待划分点放入左子树点集和右子树的待划分点集中。具体规则为将d维度值小于当前节点d维度值的待划分点放入左子树中，其他的则放入右子树中。循环执行第5.1步；
[0051]
6.针对观测的经纬度位置查找kd树t中与之距离在指定范围h内的模式网格中心点，得到l个模式网格中心点e
l
(l＝1,
…
,l)，查找方法如图4所示：
[0052]
6.1.以kd树t为初始节点，设置当前节nd＝t；
[0053]
6.2.判断当前节点nd是否为空，为空则结束搜索，否则继续搜索；
[0054]
6.3.计算当前节点nd与观测位置之间的距离d1，若d1小于等于指定范围h，则记录当前节点nd为符合要求的模式网格中心点e
l
。
[0055]
6.4.判断当前节点nd的分割维度d，计算当前节点nd与观测位置在分割维度d上的距离d2，若d2》h意味着观测位置位于nd的左子树空间内，与左子树空间点的距离可能比与右子树空间点的距离更小，因此优先递归搜索nd的左子树。反之，优先搜索nd的右子树。搜索循环执行第6.2步。
[0056]
6.5.当搜索到叶子节点时，搜索过程开始回溯。判断观测位置与当前节点nd在分割维度d上的距离是否大于指定范围h。若是，意味着当前节点nd的另一子树
空间内不可能存在更近的点，则继续回溯。否则，需要执行第6.2步进一步搜索当前节点nd的另一子树空间。直到回溯到kd树根节点t结束。
[0057]
7.根据模式网格中心点e
l
和模式网格中心点与顶点的对应关系r
lp
(p＝1,
…
,p)，获取模式网格顶点的位置
[0058]
8.根据模式网格顶点的位置判断观测是否位于模式网格中心点e
l
所对应的模式网格g
l
中。若是则计算观测位置与模式网格中心点e
l
位置的距离d：若ne
l
《0则令ne
l
＝i且nd
l
＝d，否则判断d与nd
l
的大小，若d《nd
l
则令ne
l
＝i且nd
l
＝d；否则循环到l+1，执行第7步，直到l＝1,
…
,l循环完成，执行第9步；
[0059]
9.循环到i+1，执行第6步，直到i＝1,
…
,m循环完成，执行第10步；
[0060]
10.遍历模式网格保留观测序号的存储数组nej(j＝1,
…
,n)，nej》0中保存的观测序号即为观测稀疏化后保留的观测序号，其余观测即为观测稀疏化剔除的观测；
[0061]
11.结束。
[0062]
实施例
[0063]
1.从网站https://resources.marine.copernicus.eu/product-detail/sst_glo_sst_l4_nrt_observations_010_001/information下载海表温度观测数据，从观测数据中读取观测(i＝1,
…
,m)的经纬度位置从非结构网格有限体积海洋模式fvcom的全球海洋模式数据中模式网格gj中心点ej(j＝1,
…
,n)的经纬度位置模式网格顶点nk(k＝1,
…
,k)的经纬度位置(lonk,latk)，以及模式网格中心点与顶点的对应关系r
jp
(j＝1,
…
,n,p＝1,
…
,3)；
[0064]
2.将观测和模式网格的经度范围调整一致，统一为[-180
°
w,180
°
e]；
[0065]
3.若观测或模式网格为循环边界，则复制经度范围在[-180
°
w,(-180+h)
°
w]之间的观测并将它们的经度值加360，复制经度范围在[(180-h)
°
e,180
°
e]之间的观测并将它们的经度值减360，其中h为指定范围。
[0066]
4.初始化模式网格保留观测序号的存储数组nej＝-1(j＝1,
…
,n)和模式网格保留观测的距离数组ndj＝-1(j＝1,
…
,n)，-1表明尚未找到位于模式网格中的观测；
[0067]
5.根据模式网格中心点ej(j＝1,
…
,n)的经纬度位置构建kd树t，如图3所示，具体方法是：
[0068]
5.1.判断待划分的模式网格中心点集是否为空集，若是，则结束构建过程并执行第6步，否则开始构建过程；
[0069]
5.2.设置分割维度d，若kd树为空则设置分割维度d为经度，否则设置当前步分割维度与上一步分割维度不同，即将分割维度d交替设置为经度、纬度；
[0070]
5.3.选用快速排序法对模式网格点的d维度进行排序，选取d维度上位于中位数的模式网格中心点，构造kd树的当前节点；
[0071]
5.4.根据当前节点，将剩余的待划分点放入左子树点集和右子树的待划分点集中。具体规则为将d维度值小于当前节点d维度值的待划分点放入左子树中，其他的则放入
右子树中。循环执行第5.1步；
[0072]
6.针对观测的经纬度位置查找kd树t中与之距离在指定范围h内的模式网格中心点，得到l个模式网格中心点e
l
(l＝1,
…
,l)，查找方法如图4所示：
[0073]
6.1.以kd树t为初始节点，设置当前节nd＝t；
[0074]
6.2.判断当前节点nd是否为空，为空则结束搜索，否则继续搜索；
[0075]
6.3.计算当前节点nd与观测位置之间的距离d1，若d1小于等于指定范围h，则记录当前节点nd为符合要求的模式网格中心点e
l
。
[0076]
6.4.判断当前节点nd的分割维度d，计算当前节点nd与观测位置在分割维度d上的距离d2，若d2》h意味着观测位置位于nd的左子树空间内，与左子树空间点的距离可能比与右子树空间点的距离更小，因此优先递归搜索nd的左子树。反之，优先搜索nd的右子树。搜索循环执行第6.2步。
[0077]
6.5.当搜索到叶子节点时，搜索过程开始回溯。判断观测位置与当前节点nd在分割维度d上的距离是否大于指定范围h。若是，意味着当前节点nd的另一子树空间内不可能存在更近的点，则继续回溯。否则，需要执行第6.2步进一步搜索当前节点nd的另一子树空间。直到回溯到kd树根节点t结束。
[0078]
7.根据模式网格中心点e
l
和模式网格中心点与顶点的对应关系r
lp
(p＝1,
…
,3)，获取模式网格顶点的位置
[0079]
8.根据模式网格顶点的位置判断观测是否位于模式网格中心点e
l
所对应的模式网格g
l
中。若是则计算观测位置与模式网格中心点e
l
位置的距离d：若ne
l
《0则令ne
l
＝i且nd
l
＝d，否则判断d与nd
l
的大小，若d《nd
l
则令ne
l
＝i且nd
l
＝d；否则循环到l+1，执行第7步，直到l＝1,
…
,l循环完成，执行第9步；
[0080]
9.循环到i+1，执行第三步，直到i＝1,
…
,m循环完成，执行第10步；
[0081]
10.遍历模式网格保留观测序号的存储数组nej(j＝1,
…
,n)，nej》0中保存的观测序号即为观测稀疏化后保留的观测序号，其余观测即为观测稀疏化剔除的观测；
[0082]
11.结束。
[0083]
对传统的穷举法和基于kd树的观测稀疏化方法进行测试。fvcom全球海洋模式采用非结构的三角网格，其网格中心点数为688991，网格顶点数为359120。针对全球1/20
°
的海表温度观测资料进行稀疏化，除去陆地、湖泊、河流、海冰点后，观测数为14160103。采用湖南大学国家超级计算长沙中心的天河一号，计算节点主频为2.3ghz，内存为48g。仅使用单核，传统的穷举法处理1个观测的计算时间约为218s，基于kd树的非结构网格地球系统模式观测稀疏化方法处理14160103个观测的计算时间约为15s，加速效果显著。基于kd树的非结构网格地球系统模式观测稀疏化方法能够剔除模式网格外的观测，且能够有效去除观测冗余信息。
[0084]
本发明的有益效果如下：
[0085]
本发明与传统的穷举法相比，加速效果显著，能够剔除模式网格外的观测，且能够
有效去除观测冗余信息。
[0086]
本文所使用的词语“优选的”意指用作实例、示例或例证。本文描述为“优选的”任意方面或设计不必被解释为比其他方面或设计更有利。相反，词语“优选的”的使用旨在以具体方式提出概念。如本技术中所使用的术语“或”旨在意指包含的“或”而非排除的“或”。即，除非另外指定或从上下文中清楚，“x使用a或b”意指自然包括排列的任意一个。即，如果x使用a；x使用b；或x使用a和b二者，则“x使用a或b”在前述任一示例中得到满足。
[0087]
而且，尽管已经相对于一个或实现方式示出并描述了本公开，但是本领域技术人员基于对本说明书和附图的阅读和理解将会想到等价变型和修改。本公开包括所有这样的修改和变型，并且仅由所附权利要求的范围限制。特别地关于由上述组件(例如元件等)执行的各种功能，用于描述这样的组件的术语旨在对应于执行所述组件的指定功能(例如其在功能上是等价的)的任意组件(除非另外指示)，即使在结构上与执行本文所示的本公开的示范性实现方式中的功能的公开结构不等同。此外，尽管本公开的特定特征已经相对于若干实现方式中的仅一个被公开，但是这种特征可以与如可以对给定或特定应用而言是期望和有利的其他实现方式的一个或其他特征组合。而且，就术语“包括”、“具有”、“含有”或其变形被用在具体实施方式或权利要求中而言，这样的术语旨在以与术语“包含”相似的方式包括。
[0088]
本发明实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以多个或多个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。上述提到的存储介质可以是只读存储器，磁盘或光盘等。上述的各装置或系统，可以执行相应方法实施例中的存储方法。
[0089]
综上所述，上述实施例为本发明的一种实施方式，但本发明的实施方式并不受所述实施例的限制，其他的任何背离本发明的精神实质与原理下所做的改变、修饰、代替、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。