改进的空间高精度曲面建模-ad算法

改进的空间高精度曲面建模-ad算法

1改进的hism-ad算法高精度曲线建模算法（ham）的研究始于20世纪80年代末。在冰斗形态建模中，首先应用于冰斗模型建模。在此基础上，我们成功应用了基于dem的多段模型理论。土壤属性空间模拟、气候因素空间模拟以及其他地理和生态领域的模拟。同时,随着HASM理论的发展与完善,各种HASM求解算法也相继提出,如预处理共轭梯度法、多重网格算法、平差算法、自适应算法及基于GPU的HASM加速算法等,但是,HASM计算效率一直没有得到明显的改善,以预处理共轭梯度算法为例,当计算规模为1000×1000时,采用当前普通计算机的时间消耗接近3小时。计算速度的限制已经成为限制HASM理论进一步应用的瓶颈。HASM计算耗时大的主要原因是该模型将曲面模拟问题转化为大规模方程组的求解,因而不管采用何种迭代法,都难以获得令人满意的计算速度。为了提高计算速度,岳天祥等(2010)根据分组平差原理构建了HASM局部算法,即HASM-AD算法。该算法将全局建立线性方程组求解问题转化为局部方程组建立求解,一定程度上提高了模型的计算速度。但是,后续的实验表明,由于该算法需要在局部窗口中判断采样点信息,当计算规模增加后重复的查找操作会使得计算时间急剧增加。为了解决上述问题,本文对现有的HASM-AD算法进行改进,在分组平差计算原理的基础上,通过在计算过程中为采样点添加索引信息,避免了重复的查找操作而大幅提高计算速度。同时,改进的HASM-AD算法在遍历独立计算单元的过程中,计算第一类、第二类基本量及克式符号等变量,降低了模型计算过程中的内存消耗。本文以全国陆地降水空间分布模拟为例,验证了改进的HASM-AD算法在空间变量模拟研究中的适用性。2改进的哈m-ad算法2.1克式符号的计算与其他HASM算法相比,HASM-AD算法可以看作是一种局部算法。该算法根据分组平差计算原理,将计算区域划分为一系列相邻的5×5分析窗口,称作独立计算单元(图1)。以中间3×3的格网点上真实值与初始值的差(称作改正数)作为未知数,将高斯方程作为约束条件建立平差方程组,并结合采样点信息(如果当前5×5分析窗口中包含采样点)建立最终的求解方程组(式1),5×5分析窗口的边界点仅参与曲面基本量及克式符号的计算。其中,vi(i=0,…,8)代表分析窗口内部解算格网点上的模拟变量的改正数,作为方程组的未知数;fi,j是格网点上初始值,由HASM模型的输入驱动场给出;是采样点值,同样由HASM模型的输入采样点信息得到;E,G,L,M,N,Γ为曲面基本量和克式符号,可由驱动场计算得到,具体的计算表达式可参见参考文献。在实际应用中,可认为一个独立计算单元中至多包含一个采样点,这样一个独立计算单元所建立的条件平差方程组是3阶(独立计算单元中不包含采样点)或4阶(独立计算单元中包含采样点)。因此,无论是采用直接法还是普通迭代法求解,计算量都很小,可以提高计算速度并减小内存消耗。当计算规模较小时,HASM-AD算法可以提高计算效率,但是实验表明,由于每个独立计算单元都需要进行采样点信息查找操作,当计算规模增大时HASM-AD算法计算耗时也急剧增加。例如,以个人计算机为例(基于VS2010平台C++程序),当计算规模为1000×1000时,预处理共轭梯度算法耗时为10774s,平差算法耗时为5741s,但是,当计算规模增加到1400×1400时,预处理共轭梯度法消耗时间为20746s,而平差算法耗时则超过预处理共轭梯度算法,达到了22174s。因此,本文对HASM-AD算法做了进一步改进,旨在降低采样点查找操作带来的时间消耗,从而整体上提高HASM模型的效率。2.2立计算单元的检索本文对HASM-AD算法作了两处改进:(1)为采样点添加索引信息,避免在每个独立计算单元求解时都进行查找操作;(2)根据HASM-AD分组算法的特点,在分组计算过程中求解建立方程组所需要的第一类、第二类基本量及克式符号。这两处改进不仅显著提高了HASM-AD算法的运算速度,也节省了计算过程中的内存空间需求,下面进行详细说明。2.2.1采样点索引创建规则图中黑色实心圆代表采样点,空心三角形则代表包含该采样点的独立计算单元中心点,其中,左上角、右下角的独立计算单元用虚线框标出。当独立计算单元中包含采样点时,为建立条件平差方程组,还需要知道采样点在3×3计算格网中的次序,以便决定条件平差方程组中第4个方程的非0系数。例如,图2中左上角的空心三角形标注的数字“8”表示以该点为中心的独立计算单元中,采样点所对应的未知变量次序号为8(序号从0开始计),其余类推。{booltag;//标识以当前格网点为中心的计算单元中是否包含采样点intcixu1;//标识采样点(如果有)在当前计算单元中的顺序(取值0-8)intcixu2;//标识采样点(如果有)在采样数组中的顺序号}Index[m×n];//m,n分别为计算区域的行号和列号在计算过程中,对采样点数组遍历一次,遍历过程中对每个格网点都建立采样点索引结构体。索引创建规则为:对于当前采样点,首先,根据其坐标信息确定该采样点在计算矩阵中的位置(例如图2中实心黑点);然后,对该采样点所在的位置极其周围的格网点(共9个)建立索引信息(图2中空心三角形)。令空心三角形位置格网点所对应的结构体数组元素的tag值为1,表示以该格网点为中心的独立计算单元内包含采样点;第二个变量cixu1按照图中所标注的进行赋值,例如,图2第一个空心三角形处的格网点的cixu1赋值为8,表示以该格网点为中心的独立计算单元中,采样点在3×3计算单元中的序号为8;第三个变量cixu2赋值为当前采样点在采样点数组中的序号,以便后面建立方程组时直接从采样点数组中获取该采样点的属性值。在后面独立计算单元遍历时,首先根据索引结构体的tag变量判断当前计算单元是否包含采样点,若包含,则分别从索引结构体的第二个、第三个变量中得到采样点在当前计算单元中的顺序,以及采样点的属性值,进而建立平差方程求解;若不包含,则直接建立不包含采样点的平差方程组进行求解。2.2.2高效阶次的高效计算本文对HASM-AD算法的另一个改进是取消计算过程中用来存储第一类、第二类基本量及克式符号的全局数组。其他HASM求解算法需要建立全局线性方程组进行求解,因此在方程组建立之前,首先需要计算每个格网点上的第一类基本量(3个)、第二类基本量(3个)及克式符号(6个),为此需要消耗的内存是计算规模的12倍。如此大的内存消耗在大规模计算中是不容忽视的。而分析HASM-AD算法的特点可以发现,该算法实质上是以5×5分析窗口遍历计算区域,遍历过程中建立条件平差方程组求解,因此,可在遍历过程中实时计算当前需要的第一类、第二类基本量,以及克式符号。此时不需要开辟数组空间,仅需定义若干变量即可。综上分析,改进后的HASM-AD算法可以按照图3的步骤进行。2.3计算速度比较HASM-AD算法包括采样点信息索引创建和独立计算单元计算两个过程,其中,创建索引的过程只需要遍历一遍采样数组,时间可以忽略不计,因此,独立计算单元计算是该算法的主要时间消耗。因为每个计算单元所求解方程组为3阶或4阶,故单个计算单元时间消耗基本一致,假设完成一个独立计算单元计算所消耗的平均时间为t,而整个计算区域包含的独立计算单元数目为M×N,M和N分别为计算区域的行数和列数,因此,最后需要的计算总时间为(M×N)×t。可以看出改进的HASM-AD算法时间消耗与计算区域的格网数成正比,即计算时间随着计算格网数目的增加而线性增加。为了验证以上分析,本文统计了不同计算规模的改进HASM-AD算法时间消耗(图4(b))。其中,实验曲面为高斯合成曲面,采样点采样率为1%并且为随机采样。作为对比,同时统计了HASM预处理共轭梯度算法在相同计算规模下的时间消耗(图4(a))。可以看出,HASM-PCG计算时间随计算规模增加呈指数形式增长,增长速度非常快,例如,当计算规模为200×200时,时间消耗为165s,当计算规模增加到1000×1000时,时间消耗达到了10744s。而改进的HASM-AD计算时间与计算规模呈近似线性关系,当计算规模为200×200时,时间消耗仅为1s,当计算规模增加到1000×1000时,时间消耗也仅为32s。由于将全局大规模线性方程组计算转化为依次计算独立计算单元中的条件平差方程组,并且改进的HASM-AD算法在独立计算单元方程组建立过程中避免了重复的采样点查找操作,而只需要进行布尔运算,因此,改进的HASM-AD算法计算效率与HASM-PCG算法相比得到了显著提高。改进的HASM-AD算法同样减小了计算过程中的内存需求。由于将第一类、第二类基本量及克氏符号改为在遍历独立计算单元过程中实时计算,因此,算法在计算过程中仅需要存储驱动场数组和存储索引信息的结构体数组(采样点信息存储空间忽略不计),并且索引结构体中的变量为1个布尔型,2个整型,所以,若假设驱动场数组的存储空间为M,则改进的HASM-AD算法的内存需求不足2M。而对于其他建立全局线性方程组的算法,以HASM-PCG算法为例,在计算过程中需要存储驱动场数据、第一类基本量、第二类基本量、克氏符号、方程组系数矩阵(用稀疏矩阵三元组格式)、方程组右端向量,以及迭代计算过程中的临时变量,内存总需求超过了20M。综上分析看出,相比其他HASM求解算法,改进的HASM-AD算法计算速度有大幅度提高,同时计算过程中的内存需求也明显减少,因此,适用于大区域、大规模数据的空间变量模拟。关于数学曲面模拟的验证可见参考文献,本文以全国降雨量模拟来验证改进后的HASM-AD算法在空间变量模拟中的精度,并说明该算法相比于其他HASM算法的效率优势。3模拟结果对比降水是陆地生物圈模型的主要驱动变量,是重要的农业气候资源。根据有限的降水观测站插值得到区域降水连续的分布是生态学、农用气候学研究的热点问题之一。目前应用于降水空间分布模拟的研究中的方法主要有数学插值方法、人工神经网络、多元回归分析等。王晨亮(2012),赵娜(2013)等将HASM引入到降水空间插值模拟中,获得了精度更好的结果。为了验证改进的HASM-AD算法的空间变量模拟精度及计算速度的优势,本文以模拟中国1951-2010年60年间的平均降水量作为验证示例。可以获取降水数据的气象站点为712个(图5),其中,选取605个站点(85%)的数据用来模拟降水分布,107个站点(15%)用作精度验证。验证站点为随机选择,共进行20次随机选择验证。最后,统计改进的HASM-AD算法及比较算法(HASM-PCG算法)的模拟结果精度及时间消耗,如表1所示。由上可见,与HASM-PCG算法相比,在模拟精度方面改进的HASM-AD算法略有提高,HASM-PCG算法模拟结果的均方根误差和平均绝对误差分别是131.6400和77.9690,而改进的HASM-AD算法模拟结果的均方根误差和平均绝对误差则分别降低到131.4991和77.5088.从模型输出结果(图6)可以看出,改进的HASM-AD算法模拟得到的中国年平均降水分布合理,符合实际情况,能够满足大尺度上水文、生态研究需求。在计算速度方面,改进的HASM-AD算法优势更为明显,对于本问题的计算规模(计算区域网格数目约为2×105),在相同的计算机环境下,HASM-PCG算法耗时1920s,而改进的HASM-AD算法仅需4s,计算速度提高了480倍。表明改进的HASM-AD算法在模拟精度不低于HASM-PCG算法的前提下能显著提高计算速度,从而降低计算时间消耗。4改进的hism-ad算法模拟结果的比较与传统的全局HASM算法不同,HASM-AD算法根据分组平差计算原理,将大规模线性方程组求解转化为求解独立计算单元中的平差方程组(3阶或4阶),因而在方程组求解阶段可以显著减少计算时间。而本文在HASM-AD算法的基础上对采样点添加了索引信息,避免了遍历独立计算单元过程中对采样点信息的重复查找,因此,大幅提高了整个模拟过程的效率。本文得出的主要结论如下:(1)改进的HASM-AD算法添加了采样点索引信息,避免了在计算过程中重复寻找计算单元中是否含有采样点,极大地提高了HASM的计算速度。不同计算规模的数值实验表明,HASM-AD算法计算时间随计算区域网格数目增加而线性增加;(2)改进的HASM-AD算法取消模型运行初始阶段的第一类、第二类基本量,以及克式符号的计算,在独立计算单元的遍历过程中实时动态计算所需要的基本量及克式符号值,显著降低了HASM模型运行的内存需求;(3)采用全国陆地降水模拟作为实例验证改进的HASM-AD算法的模拟精度及计算效率,结果表明,改进的HASM-AD算法模拟结果精度优于HASM-PCG算法的模拟结果,平均绝对误差及误差均方根均小于后者,而计算效率优势更为明显,HASM-PCG算法耗时1920s,而HASM-A