人群搜索算法拟合Kriging参数的空间数据插值

胡芸,侯明勋

(上海交通大学船舶海洋与建筑工程学院,上海 200240)

0 引言

空间数据插值是利用观测数据,通过建立连续的或分段连续的表面函数模型来计算空间任意位置的现象属性数据的空间分析方法[1],在矿藏预测[2]、地质勘探[3]、气象预报[4]和农林产值预测[5]等众多领域得到广泛应用.Kriging插值法是国际上最常用的空间估值方法,其独特优势在于考虑样本的空间相关性的同时可以计算估计误差[6].决定Kriging插值精度的关键是变异函数理论模型的选取和模型参数的设置,模型参数的拟合仍存在不足[7].

Kriging变异函数模型参数的确定方法通常为自动拟合[8],自动拟合主要分为最小二乘法和加权最小二乘回归法.加权回归法相比最小二乘法考虑了数据对的权重,实现了最小二乘意义上的最优拟合[9],是Kriging插值常用的拟合方法[8,10].然而,加权回归法作为线性拟合方法,不适用于考虑变异函数变程范围外的变异值,其最优拟合仅局限于变程范围内.近年来,已有学者对Kriging理论模型参数的拟合进行相关研究.曾怀恩等提出在变程范围内使用遗传算法求解模型参数避免了传统人工拟合方法的人为影响[11],但与传统自动拟合方法相比其优势尚不明确;王英博等以预测误差作为适应度评价函数确定模型参数,减小了预测误差[7],但对于采样数据的各向异性缺乏探讨;ZHANG X等以变异函数与变异值最小差值为目标设计拟合算法提高了变异函数模型拟合精度[12],但变程椭圆方向不明确,对于空间相关性的考虑不足.以上研究从不同角度探讨了变异函数拟合的方法并提高了函数模型准确度,但针对变异函数不同区域变异值权重分配及减小误差的方法等方面的研究相对较少,为了更全面计算变异值,得到更准确的变异函数模型,本研究采用新型智能算法拟合Kriging模型参数.

人群搜索算法(seeker optimization algorithm,SOA)是一种新的基于种群的启发式随机搜索算法.该算法模拟人的随机搜索行为,完成对优化问题解的进化过程[13].SOA自提出以来,就引起国内外学者的广泛关注和应用,如李彬等采用SOA重构实际汽轮机叶片型线[14],沈俊等结合SOA实现了机载激光雷达回波信号的波形分解[15],封京梅等引进单纯形算法,提高SOA局部搜索能力,用其求解一类不可微绝对值方程[16].

本研究使用SOA算法以变异函数拟合值与观测值的加权残差平方和最小为目标建立适应度函数,将具备空间相关性的变异值分组并计算相应权重,同时考虑样本数据的各向异性,搜索迭代模型最优参数进行Kriging插值,以提高Kriging插值方法的准确度.

1 Kriging原理

Kriging法是依据协方差函数对随机场进行空间建模和预测插值的回归算法在特定的随机过程,例如固有平稳过程中,Kriging能够给出最优线性无偏估计(best linear unbiased prediction,BLUP),在地统计学中也被称为空间最优无偏估计器(spatial BLUP).

Kriging在对空间场进行预测时,将空间场视为随机过程(stochastic process)的推广,即随机场(random field).用Z(x)表示随机场的随机函数,Z(xi)(i=1,2,…,n)是其观测数据,待估块段V的真值是估计邻域内n个信息值的线性组合,即

式中,γ(x)为变异函数,μ为拉格朗日常数.

Kriging方程表明变异函数γ(x)对估计影响较大.

变异函数通常采用以变异值γ(x)及其滞后距h为对应的变异曲线表征,其理论模型为根据实验变异函数值拟合的最优理论变异函数曲线.目前各类插值软件中变异函数的拟合方法多采用加权回归法(weightedleastsquare,WLS),是一种通过对预测残差求导得到使残差最小的模型参数.

2 人群搜索算法

SOA的不确定推理行为是利用模糊系统的逼近能力,模拟人的智能搜索行为,用以建立感知和行为,即目标函数值和步长之间的联系.

在不确定推理的过程中,为了设计一个适用于大多数优化问题的模糊系统,将目标函数值按递减的顺序排序,从而把实函数值按递减的顺序排序,从而把实函数值转换成从1到S(S是种群大小)的自然数作为不确定推理的输入.设xij是j维空间的第i个种群的位置,在第t+1次迭代时,根据搜索步长αij(t)和搜索方向dij(t)对每个搜寻个体进行位置更新：

xij(t+1)=xij(t)+Δxij(t+1)

(4)

Δxij(t+1)=αij(t)dij(t)

(5)

式中：ω是惯性权值,随进化代数的增加从0.9线性递减至0.1.

搜索步长αij根据不确定推理的行为可得出：

式中：ui为搜索空间目标函数值i的隶属度,δi为高斯隶属函数参数,可由下式确定：

式中,xmin和xmax分别是同一子群中的具有最小和最大函数值的位置;t和Tmax分别为迭代次数和最大迭代次数;函数abs(·)对输入的每个元素取绝对值.

对于隶属度ui,是通过模拟人的搜索行为中的随机性,设置均匀、随机地分布在区间[ui,1]上的实数函数rand(ui,i),由不确定性推理条件得出.目标函数的模糊变量采用线性隶属函数,使隶属度直接与函数值的排列顺序成正比,即在最佳位置有最大隶属度umax=1.0,通常隶属度小于0.011 1时输出变量可忽略不计,故最小隶属度umin=0.011 1.

uij=rand(ui,1),j=1,2,…,D

(12)

ui=umax-(s-Ii)·(umax-umin)/(s-I),i=1,2,…,s

(13)

式中：ui为目标函数值i的隶属度;Ii是种群函数值按降序排列后xi(t)的序列编号;D为搜索空间维数[13].

3 SOA拟合模型参数的Kriging插值

Kriging插值的首要步骤是对已有采样点数据进行处理,即计算出各采样点对间的距离和角度后根据滞后距和容差使用点对数法进行分组,再计算出各个方向滞后距对应的变异值.

根据变异值拟合变异函数模型的过程中,传统WLS拟合会由于权重、分配方式等因素影响模拟的准确度,本文中采用的SOA求解Kriging模型参数的具体方法就是通过SOA迭代函数模型的块金值、变程和基台值,减小模拟的变异值与实际变异值的误差.

插值流程如图1所示,其中SOA拟合变异函数的具体过程如下：

1)确定搜索者种群规模,变量范围和最大迭代次数,空间维数m=3,初始化搜索者的位置.

2)计算搜索者个体的适应度函数值.Kriging变异函数模型参数确定的适应度函数表达式为

式中C为适应度;hi为滞后距;N(hi)为滞后距离分组后各组hi的点对数,w(i)为权重系数,γ*(hi)为变异函数实验值,γ#(hi)为γ*(hi)的拟合值.

3)比较个体当前适应度值与当前个体最优位置的适应度值,更新个体最优值和种群最优值.

4)搜寻策略即根据公式计算搜索个体i在每一维j的搜索方向和步长;按照公式更新搜索者位置,计算适应度函数值.

5)若达到最大迭代次数或全局最优适应度值小于设定值,停止搜索,否则回到步骤2).

4 高程估计与结果分析

4.1 理论模型参数的估计选取美国马萨诸塞州某地区527个数据点[8],滞后宽度取30 m,将最大滞后距1 650 m分为55组,考虑到变异函数存在各向异性的可能性,对变异栅格的正北-正南(N-S)、东北-西南(EN-WS)、正东-正西(E-W)和东南-西北(ES-WN)方向,容差均为 45°的4个区域进行拟合.

4.1.1 参数设置 人群搜索算法的拟合标准即适应度函数的最小值：

种群规模S=100,最大迭代次数maxgen=100,空间维数m=3,最大隶属度值Umax=0.950 0,最小隶属度值Umin=0.011 1,具体拟合流程见第3部分.

对于球状模型各参数初始设置考虑如下：块金值在理想状态下为0,本实例中变异值的高频数为25,最大滞后距的数值为1 650.在初始种群的设置上,考虑到变异函数各参数的特性以及数量级上的差异,将块金值、基台值和变程分别以0、25、1 650作为初始值,加上随机数生成大小为100的初始种群.

4.1.2SOA搜索性能探讨 为探讨SOA全局搜索性能和收敛性,以球状模型为例,现设置其他三组不同初始参数和原初始参数(见表1)进行实验,收敛曲线见图2.实验中,a、b、c组均可以快速收敛并得到全局最优解,d组收敛较快但出现早熟.实验结果表明,即使在初始参数偏差非常大的情况下,人群搜索算法也可快速收敛并得到最优解,但如果不考虑各参数特性,不设置初始值将会陷入局部搜索,根据实际变异函数模型特性设定初始值则可解决该问题.总体而言,人群搜索算法的具备收敛速度快、参数设置简单、全局搜索能力较强的特点.

表1 球状模型参数初始值

图2 SOA迭代图

4.2 拟合结果SOA对Kriging球状模型的N-S、EN-WS、E-W、ES-WN方向模型参数拟合结果如表2所示,由表2可知变异函数具有各向异性,使用SOA算法计算各向异性变异函数理论模型(图3)各向参数,确定变程椭圆长轴位于39°～219°,变异函数的理论模型块金值为3.47 m2,偏基台值为23.35 m2,42°～222°方向变程为736.82,129°～309°方向变程为275.96,各向异性比为2.67,各向异性变异函数理论模型为

表2 球状模型参数拟合结果

图3 变异函数模型

γ(h)=3.47+23.35×[1.5×(h/736.82)

-0.5×(h/736.82)3]

(16)

各向异性分离距离h的套和公式如式(17)所示,其中hc为39°～219°方向的分离距离,hd为129°～309°方向的分离距离,K为各向异性比.

使用SOA-Kriging算法进行估计的高程和估计方差分布情况如图4所示.估计误差分布较均匀,数据点过于稀疏的区域误差偏大.

图4 插值预测结果及方差分布图

同理可拟合高斯模型和指数模型的各向参数得到估计结果.

4.3 SOA与加权回归法的比较为了探讨SOA拟合变异函数的优势,比较两种方法所得模型的加权残差和(图5),并将SOA所得参数和常用的加权最小二乘法(weighted least square,WLS)所得参数进行Kriging插值,使用交叉验证计算各插值结果误差并进行比较,以平均绝对误差(mean absolute error,MAE)和均方差(mean square error,MSE)作为评价指标,并统计误差分布频率,计算结果见图6和表3.

表3 SOA和WLS拟合Kriging模型的误差比较

图5 三种模型下SOA相对于WLS各向加权残差平方和的降幅

图6 估计误差分布图

从图5可看出,三种模型下的实验结果中,在不同方向使用SOA拟合的加权残差平方和(式(14)中的C)相比WLS均有所降低,各个模型大部分方向C的降低幅度位于5%～40%之间,说明由SOA拟合变异函数相比WLS能得到与实际变异值差异更小的理论模型参数.

从图6中Kriging插值的误差频率分布情况可以看出,球状模型(sphere model,SM)基于SOA拟合的空间插值估计误差主要集中于[-9.5,7.5],基于WLS的估计误差主要集中于[-12.5,12.5],指数模型(exponential model,EM)下的误差分布主要集中于[-8.5,11.5]和[-12.5,11.5];高斯模型(Gaussian model,GM)下的误差分布主要集中于[-9.5,12.5]和[12.5,12.5].以上3种模型实验结果表明,SOA-Kriging拟合的误差分布更集中,说明SOA-Kriging算法可以减小误差离散程度和估计不确定性,给出更精准的估计结果.

由表3可见,无论是平均绝对误差还是均方误差,WLS拟合变异函数模型所得空间插值结果都比SOA要大,因此可认为,WLS拟合变异函数的不确定性更大.由表3计算得到,相比于加权最小二乘法拟合,采用SOA拟合Kriging模型参数时,在球状模型、指数模型和高斯模型下Kriging插值的绝对误差分别降低31.02%、24.02%和25.13%.

4.4 Kriging法与其他插值方法的比较为探讨SOA拟合Kriging参数的插值方法的优势,将其与反距离加权法、径向基函数插值法和最小曲率法的插值结果进行比较,见图7.从图中可以看出,反距离加权法所绘制等高线图环状地形较多与实际不符;径向基函数法的等高线较稀疏,插值结果较差;最小曲率法插值效果较差;基于SOA的Kriging法的插值结果则比较符合实际,插值效果较好.

图7 4种方法插值结果的效果图

5 结束语

本文中采用近年来出现的新型算法,人群搜索算法对Kriging变异函数模型参数进行拟合,与加权回归法拟合参数对比插值误差,并与反距离加权法、径向函数法和最小曲率的对比插值效果.研究结果表明：基于人群搜索算法的Kriging插值精度高于现广泛使用的基于加权最小二乘法拟合模型参数的Kriging方法,插值效果优于反距离加权法、径向函数法和最小曲率法,具备较好的预测准确度,是一种可行的拟合方法.