结合Markov模型和COSGSIM的插值方法
三维重心插值计算公式
三维重心插值计算公式在计算机图形学和计算机辅助设计领域中，三维重心插值是一种常用的插值方法，用于在三维空间中对数据进行插值和补偿。
三维重心插值是基于三角形网格的，通过计算三角形的重心坐标来实现插值。
在本文中，我们将介绍三维重心插值的计算公式，并讨论其在实际应用中的重要性和优势。
三维重心插值的计算公式如下：假设有一个三角形ABC，其中A、B、C分别为三个顶点，对应的数据值分别为fA、fB、fC。
现在需要在三角形内部的某一点P处进行插值，其重心坐标为(u, v, w)，则P点的插值数值可以通过以下公式计算得出：fP = u fA + v fB + w fC。
其中，u、v、w为P点的重心坐标，满足以下条件：u + v + w = 1。
u、v、w分别为P点到BC、AC、AB三条边的距离与三条边的长度之比。
三维重心插值的计算公式非常简单，但却具有广泛的应用价值。
在实际应用中，三维重心插值可以用于地形数据的插值、图像的纹理映射、三维模型的变形等方面。
下面我们将分别介绍三维重心插值在这些领域中的应用。
首先是地形数据的插值。
在地理信息系统（GIS）中，地形数据通常以离散的高程点数据的形式存在，而实际应用中需要对地形进行连续的插值和补偿。
三维重心插值可以通过对地形数据中的三角形网格进行插值，实现对地形的连续化处理，从而为地形分析和可视化提供了可靠的数据基础。
其次是图像的纹理映射。
在计算机图形学中，三维模型的表面通常需要进行纹理贴图，以增强真实感和细节。
三维重心插值可以用于在三角形网格上对纹理进行插值，从而实现对三维模型表面的纹理映射，提高了渲染效果和真实感。
最后是三维模型的变形。
在计算机辅助设计（CAD）领域中，三维模型的变形是一项重要的技术，可以用于模拟材料的变形和形状的调整。
三维重心插值可以通过对三角形网格进行插值，实现对三维模型的形状变形，为工程设计和仿真分析提供了有效的工具和方法。
总之，三维重心插值作为一种常用的插值方法，在计算机图形学和计算机辅助设计领域具有重要的应用价值。
融合VMD和XGBoost算法的GNSS高程时间序列预测方法
㊀２０２３年８月A c t aG e o d a e t i c ae tC a r t o g r a p h i c aS i n i c a A u g u s t,２０２３㊀㊀第５２卷㊀第８期测㊀绘㊀学㊀报V o l．５２,N o．８引文格式:鲁铁定,李祯,贺小星,等．融合VM D和X G B o o s t算法的G N S S高程时间序列预测方法[J]．测绘学报,２０２３,５２(８):１２３５Ｇ１２４４．D O I:１０．１１９４７/j．A G C S．２０２３．２０２２００５２．L U T i e d i n g,L IZ h e n,H E X i a o x i n g,e ta l．G N S S v e r t i c a lt i m es e r i e s p r e d i c t i o n m e t h o di n t e g r a t i n g VM D a n d X G B o o s ta l g o r i t h m s[J]．A c t aG e o d a e t i c a e tC a r t o g r a p h i c aS i n i c a,２０２３,５２(８):１２３５Ｇ１２４４．D O I:１０．１１９４７/j．A G C S．２０２３．２０２２００５２．融合V M D和X G B o o s t算法的G N S S高程时间序列预测方法鲁铁定１,２,李㊀祯１,贺小星３,周世健４１．东华理工大学测绘工程学院,江西南昌３３００１３;２．自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室,江西南昌３３００１３;３．江西理工大学土木与测绘工程学院,江西赣州３４１０００;４．南昌航空大学,江西南昌３３００６３G N S Sv e r t i c a l t i m es e r i e s p r e d i c t i o n m e t h o di n t e g r a t i n g V M D a n dX G B o o s t a l g o r i t h m sL UT i e d i n g１,２,L I Z h e n１,H EX i a o x i n g３,Z H O US h i j i a n４１．S c h o o lo fG e o d e s y a n d G e o m a t i c s,E a s tC h i n a U n i v e r s i t y o fT e c h n o l o g y,N a n c h a n g３３００１３,C h i n a;２．K e y L a b o r a t o r y o fM i n eE n v i r o n m e n t a lM o n i t o r i n g a n dI m p r o v i n g a r o u n dP o y a n g L a k e,M i n i s t r y o fN a t u r a lR e s o u r c e s, N a n c h a n g３３００１３,C h i n a;３．S c h o o l o fC i v i la n dS u r v e y i n g&M a p p i n g E n g i n e e r i n g,J i a n g x iU n i v e r s i t y o fS c i e n c e a n dT e c h n o l o g y,G a n z h o u３４１０００,C h i n a;４．N a n c h a n g H a n g k o n g U n i v e r s i t y,N a n c h a n g３３００６３,C h i n aA b s t r a c t:A i m i n g a t t h e p r o b l e m so f i m p e r f e c t f e a t u r es e l e c t i o na n d p o o r s t a b i l i t y i nt r a d i t i o n a lG N S S e l e v a t i o nt i m e s e r i e s p r e d i c t i o n m o d e l s,a c o m b i n e df o r e c a s t i n g m o d e lb a s e d o n v a r i a t i o n a l m o d e d e c o m p o s i t i o n(V M D)a n de x t r e m e g r a d i e n t b o o s t i n g(X GB o o s t)a l g o r i t h m i s p r o p o s e d．T h em o d e l o b t a i n s t h e r e c o n s t r u c t e ds i g n a l t h r o u g h m u l t i p l eV M Ds u bＧm o d e l s,a n di n p u t s i t i n t ot h eX G B o o s tm o d e la sa f e a t u r e f o r f o r e c a s t i n g o f t h eo r i g i n a l t i m es e r i e s．T ov e r i f y t h e p e r f o r m a n c eo f t h e f o r e c a s t i n g m o d e l,t h e e x p e r i m e n t s e l e c t st h ev e r t i c a l t i m es e r i e sd a t ao f４o b s e r v a t o r i e sf o r t h ef o r e c a s t i n g e x p e r i m e n t,t h e e x p e r i m e n t a l r e s u l t s s h o w t h a t t h e V M Dm o d e l c a n a c c u r a t e l y e x t r a c t t h e f e a t u r e s．C o m p a r e dw i t h t h e V M DＧC N NＧL S T Mm o d e l,t h ee x p e r i m e n t a l r e s u l t so fV M DＧX G B o o s ts h o wt h a t t h e M A Ev a l u e sa r er e d u c e db y １９．７４％~３５．９０％a n d t h e R M S E v a l u e s a r e r e d u c e d b y２２．２２％~３１．１４％．T h e f o r e c a s t i n g r e s u l t s h a v e h i g h e r s t a b i l i t y a n da r e h i g h l y c o r r e l a t e d t o t h eo r i g i n a l t i m es e r i e s,w h i c hc a nb e t t e r p r e d i c t t h eT a r g e t e dt i m e s e r i e s．T h e r e f o r e,t h e f o r e c a s t i n g m e t h o dc a nb ea p p l i e d t oG N S Sv e r t i c a l t i m e s e r i e s f o r e c a s t i n g．K e y w o r d s:V M D;X G B o o s t;G N S S;t i m e s e r i e s;f o r e c a s t i n gF o u n d a t i o ns u p p o r t:T h e N a t i o n a l N a t u r a l S c i e n c e F o u n d a t i o n o f C h i n a(N o s．４２０６１０７７;４２０６４００１;４２１０４０２３);T h e N a t i o n a l N a t u r a l S c i e n c e F o u n d a t i o n o f J i a n g x i,C h i n a(N o s．２０２０２B A B L２１３０３３;２０２０２B A B２１２０１０);T h eJ i a n g x iU n i v e r s i t y o fS c i e n c ea n dT e c h n o l o g y H i g hＧl e v e lT a l e n tR e s e a r c hS t a r t u p P r o j e c t(N o．２０５２００１００５６４);Y o u t h T a l e n tP l a n o f S c i e n c ea n dT e c h n o l o g y T h i n k T a n ko f C h i n aA s s o c i a t i o n f o r S c i e n c ea n dT e c h n o l o g y i n２０２２摘㊀要:针对传统G N S S高程时间序列预测模式存在特征选取不完善㊁稳定性差等问题,本文提出了一种融合V M D和X G B o o s t算法的预测模型.该模型通过多个V M D子模型得到重构信号,再将其作为特征输入X G B o o s t模型中进行原始时间序列的预测.为了验证预测模型的性能,试验选取４个观测站高程时间序列数据进行预测试验,试验结果表明,V M D模型能够准确地提取特征信息.与V M DＧC N NＧL S T M模型相比,V M DＧX G B o o s t模型预测结果的M A E值降低了１９．７４％~３５．９０％,R M S E值降低了２２．２２％~３１．１４％,预测结果具有更高的稳定性且与原始时间序列呈较强相关性,可以较好地预测出目标时间序列.因此,该预测模型可应用于G N S S高程时间序列预测.关键词:V M D;X G B o o s t;G N S S;时间序列;预测Copyright©博看网. All Rights Reserved.A u g u s t２０２３V o l．５２N o．８A G C S h t t p:ʊx b．c h i n a s m p．c o m 中图分类号:P２２８㊀㊀㊀㊀文献标识码:A㊀㊀㊀㊀文章编号:１００１Ｇ１５９５(２０２３)０８Ｇ１２３５Ｇ１０基金项目:国家自然科学基金(４２０６１０７７;４２０６４００１;４２１０４０２３);江西省自然科学基金(２０２０２B A B L２１３０３３;２０２０２B A B２１２０１０);江西理工大学高层次人才科研启动项目(２０５２００１００５６４);２０２２年度中国科协科技智库青年人才计划㊀㊀近３０年来,G N S S基准站不断积累的时间序列数据为大地测量和地球动力学研究提供了宝贵的数据基础[１Ｇ５].这些数据可以有效地反映由地球物理效应引起的长期变化趋势和非线性变化[６].因此,对G N S S坐标时间序列进行分析有助于监测地壳板块运动[７Ｇ９]㊁大坝或桥梁变形监测[１０Ｇ１３]㊁全球或区域坐标系维护[１４Ｇ１６]等领域的发展.随着G N S S相关技术的发展,G N S S数据已被可靠地应用于区域陆地运动的研究中,其中G N S S高程时间序列数据可以为研究人员进行区域垂直陆地运动分析提供有效的数据参考[１７].因此通过分析G N S S高程时间序列,可以预测连续时间点的高程[１８Ｇ１９],为判断运动趋势提供重要依据.已有研究发现,G N S S坐标时间序列中垂直方向上的噪声分量通常大于水平方向,且噪声组合模型丰富[２０Ｇ２１].在目前G N S S坐标时间序列建模中,研究人员通常以趋势项㊁周期项和噪声解释G N S S坐标时间序列的构造成分.但是噪声不同于趋势项和周期项,并不具有时间特征.因此,目前构建顾及噪声影响的高精度G N S S高程时间序列预测模型困难较大[２２Ｇ２３].随着G N S S坐标时间序列预测研究的深入,出现了基于信号分解模型对G N S S时间序列数据进行分解,然后对各分量逐一预测,最后等权相加得到预测结果的预测模式[２４].该预测模式下的模型虽然可以有效地进行G N S S坐标时间序列预测研究,但依旧存在着需要解决的问题:①该预测模式存在两个误差来源,即分解时间序列造成的误差和预测存在的误差.②有限次分解后得到的子时间序列相加结果普遍小于原始时间序列,预测模型对于不同子时间序列的预测结果误差可能存在两个方向误差,导致预测模型稳定性较差.③分解后存在数据尺度较小的子时间序列,预测模型对于数据尺度较小的时间序列预测效果相对欠佳,导致预测模型精度下降.随着人工智能技术的发展,机器学习算法愈加受到研究人员的青睐,越来越多强大的算法被应用于不同的领域.机器学习是一种可以通过手动输入特征进行时间序列预测的有效方法,其预测结果具有较强的解释性,并且机器学习算法对于预测结构化数据表现突出,其中基于决策树的算法在时间序列数据预测领域表现更加优越[２５Ｇ２６].文献[２７]提出的极端梯度提升算法在多个领域的目标检测和预测研究中得到了应用并取得了显著的效果[２８Ｇ３０].虽然X G B o o s t算法展现出强大的时间序列数据预测能力,但是X G B o o s t算法在提取非平稳时间序列特征方面能力欠佳.基于上述问题,本文提出一种融合变分模态分解(v a r i a t i o n a lm o d ed e c o m p o s i t i o n,VM D)和X G B o o s t算法的G N S S高程时间序列预测方法(VM DＧX G B o o s t模型).首先,VM DＧX G B o o s t 模型叠加由VM D模型分解得到的分量以得到重构信号,摒弃传统预测模型分步预测的模式,以消除分解时间序列带来的误差.然后,VM DＧX G B o o s t模型以重构信号作为特征,取代传统预测模型以时间作为特征的方法,达到削弱噪声对于预测精度影响的目的.最后,以重构信号作为特征的预测模式在为X G B o o s t模型提供有效特征的前提下,弥补X G B o o s t算法提取非平稳时间序列特征能力欠佳的问题,提升X G B o o s t模型预测能力.１㊀V M DＧX G B o o s t模型原理１．１㊀V M D算法VM D算法是一种时频分析算法,能够将信号一次性分解成多个单分量调幅调频信号,避免了迭代过程中遇到的端点效应和虚假分量问题[３１].该算法可以有效处理非线性㊁非平稳的G N S S高程时间序列,但由于它对噪声敏感,处理存在噪声的G N S S高程时间序列时,可能出现模态混叠现象.VM D的分解过程即变分问题的求解过程,在该算法中,本征模态函数(i n t r i n s i c m o d e f u n c t i o n,I M F)被定义为一个有带宽限制的调幅Ｇ调频函数,VM D算法的功能便是通过构造并求解约束变分问题,将原始信号分解为指定个数的I M F分量.６３２１Copyright©博看网. All Rights Reserved.第８期鲁铁定,等:融合V M D 和X G B o o s t 算法的G N S S 高程时间序列预测方法假设将一个信号分解为K 个I M F 分量,VM D 算法分解的具体流程如下[３２]:(１)通过H i l b e r t 变换,得到每个模态分量μK (t )的解析信号,进而得到其单边频谱为δ(t )＋j πt éëêêùûúú μK (t )(１)(２)对各模态解析信号预估一个中心频率e －jωK t,将每个模态的频谱调制到相应的基频带δ(t )＋j πt æèçöø÷ μK (t )éëêêùûúúe －j ωK t(２)(３)计算上述解调信号梯度平方L 的范数,估计出各模态信号带宽,受约束的变分问题为m i n{μK },{ωK }ðKd t δ(t )＋j πt æèçöø÷ μK (t )éëêêùûúúe －j ωK t ２２{}(３)s ．t ．㊀ðKμK ＝f(４)式中,μK 代表分解得到的K 个I M F 分量;ωK 表示各模态对应的中心频率.为了求解该约束性变分问题,引入二次惩罚因子α和拉格朗日乘法算子λ(t ),将约束性变分问题变为非约束性变分问题.扩展的拉格朗日表达式为L ({μK },{ωK },λ)＝αðK∂t δ(t )＋j πt æèçöø÷ μK (t )éëêêùûúúe －j ωK t ２２＋f (t )－ðKμK (t )２２＋‹λ(t ),f (t )－ðKμK (t )›(５)式中,α为二次惩罚因子;λ(t )为拉格朗日乘法算子.其中,α可在高斯噪声存在的情况下,保证信号的重构精度,通常取拉格朗日算子使得约束条件保持严格性.利用乘法算子交替方向法解决以上无约束变分问题,通过交替更新μn ＋１K㊁ωn ＋１K和λn ＋１寻求扩展拉格朗日表达式的 鞍点.１．２㊀X G B o o s t 算法X G B o o s t 算法是一个开源机器学习项目.该算法在水文时间序列预测㊁环境质量指标预测和网络异常入侵检测[３３]等研究有着较好的应用,但在G N S S 坐标时间序列相关研究未见报道.X G B o o s t 算法属于集成学习中的B o o s t i n g分支,其每一次的计算都是为了减少上一次的残差,进而在负梯度方向上建立一个新的树模型,即前面决策树的训练和预测效果会影响建立下一棵树模型时的样本输入.与同属于树模型算法的G B D T 算法相比,X G B o o s t 算法更为高效,并且在算法上也进行了改进.X G B o o s t 算法使用了二阶的泰勒展开式逼近目标函数的泛化误差部分,简化了目标函数的计算;X G B o o s t 算法通过在目标函数中加入正则项降低模型预测的波动性及改善模型过拟合现象.X G B o o s t 算法具体流程如下:假设有一个数据集A ＝{(x i ,y i ):i ＝１,２, ,n ,x i ɪm,y i ɪ},其中含有n 个观测值,每个观测值有m 个特征和一个对应的变量y .定义一个值^y i ,并通过广义模型表示为^y i ＝ϕ(x i )＝ðKK ＝１f k (x i )(６)式中,f k 表示一棵决策树,f k (x i )表示第k 棵树赋予第i 个观测值的分数.同时,使用函数f k 时,下述的正则化目标函数应该被最小化为L (ϕ)＝ðil (y i ,^y i )＋ðkΩ(f k )(７)式中,l 为损失函数.为了防止模型过于复杂,模型将惩罚项Ω设置为Ω(f k )＝γT ＋１２λw ２(８)式中,γ为控制惩罚项中枝叶数量T 的参数;λ为控制惩罚项中枝叶重量w 的参数.设置Ω(f k )项不仅可以简化算法生成的模型,且可以防止模型过于拟合.X G B o o s t 算法采用迭代法最小化目标函数.模型通过增加f j 项在第j 次迭代得到减小的目标函数为L j＝ðni ＝１l (y i ,^y (j －１)i ＋f j (x i ))＋Ω(f j )(９)式(９)可以通过泰勒展开进行简化,并且可以推导出树从给定节点分裂后的损失减少公式为L s p l i t ＝１２(ði ɪI Lg i )２ði ɪI L h i ＋λ＋(ði ɪI Rg i )２ði ɪI R h i ＋λ－(ði ɪI g i )２ði ɪIh i ＋λéëêêêùûúúú－γ(１０)式中,I 是当前节点中可用观测数据的一个子集;I L 和I R 分别为分裂后左右节点中可用观测数据的一个子集.函数g i 和h i 定义分别为g i ＝∂^y (j －１)i l (y i ,^y (j －１)i )(１１)h i ＝∂２^y (j －１)il (yi ,^y (j －１)i )(１２)从推导L s pl i t 的公式中找到任意给定节点的最佳分裂,这个函数只依赖于损失函数和正则化７３２１Copyright ©博看网. All Rights Reserved.A u gu s t ２０２３V o l ．５２N o ．８A G C S h t t p :ʊx b ．c h i n a s m p ．c o m 参数γ.同时,基于式(１０),X G B o o s t 算法可以优化任何损失函数,并且提供一阶和二阶梯度.１．３㊀V M D ＧX G B o o s t 模型VM D ＧX G B o o s t 模型通过使用VM D 算法构造特征取代X G B o o s t 算法提取特征的模块,从而将VM D 算法和X G B o o s t 算法进行融合.G N S S 高程时间序列数据通常为一维时间序列数据,具有统一的时间间隔.将G N S S 高程数据按照时间顺序一维排列为X １,X ２, ,X n (１３)通过VM D 算法将时间序列分解为K 个I M F 分量,分解结果可表达为X １I M F １X ２I M F １ X n I M F １X １I M F ２X ２I M F ２ X n I M F ２⋮⋮⋮X １I M F KX ２I M F KX n I M F Kéëêêêêêùûúúúúú(１４)式中,X n I M F K 表示X n 在第K 个I M F 分量中的值.已有的分部预测模式中,研究人员将通过预测模型对各个分量时间序列进行预测然后叠加得到预测结果,但由于时间序列没有被完全分解,导致预测结果存在二次误差.为了避免该情况发生,VM D ＧX G B o o s t 模型通过将各I M F 分量叠加生成与原始时间序列高度相关的新时间序列,并通过设置m 个各不相同的K 值得到m 个新时间序列,m 个新时间序列可整合为X １VM D １X １VM D ２X １V MD m X ２VM D １X ２VM D ２ X ２VM D m ⋮⋮⋮X n VM D １X n VM D ２X n VM D m éëêêêêêùûúúúúú(１５)式中,X n VM D m 表示当K ＝m 时,X n 由VM D 算法分解后叠加所得到的值.然后将整合的m 维时间序列加入原始时间序列生成一个(m ＋１)维时间序列X １VM D１X １VM D ２X １VM D m X １X ２VM D １X ２VM D ２X ２VM D mX ２⋮⋮⋮⋮X n VM D １X n VM D ２X n VM D m X n éëêêêêêùûúúúúú(１６)在X G B o o s t 模型中,将生成的(m ＋１)维时间序列中的前m 列时间序列数据作为特征取代X G B o o s t 模型提取特征步骤;将最后一列时间序列数据即原始时间序列数据作设置为目标序列进行预测,从而得到预测结果.１．４㊀精度评价指标本文使用平均绝对误差(m e a n a b s o l u t ee r r o r ,MA E )和均方根误差(r o o t m e a ns qu a r e e r r o r ,R M S E )作为模型预测精度的评价指标[２０].MA E 与R M S E 分别为MA E ＝１n ðni ＝１(y i －^yi )(１７)R M S E ＝１n ðni ＝１(y i －^y i )２(１８)式中,y i 为原始值;^y i 为预测值.MAE 和R M S E 的值越小,代表模型的预测精度越高,更适用于该时间序列.反之,则代表模型的预测精度越低,在该时间序列中适用性较差.通过引入皮尔森相关系数判断预测时间序列和原始时间序列之间的相关性.皮尔森相关系数将时间序列视为变量,从而计算两个时间序列的相关性,皮尔森相关系数可表达为R Y ^Y＝E (Y ^Y )－E (Y )E (^Y )σY σ^Y(１９)式中,Y 和^Y 分别代表进行皮尔森相关系数计算时将原始时间序列和预测时间序列视为的变量.式(１９)值域为[－１,１],R Y ^Y ＞０,则两时间序列正相关;R Y ^Y ＜０,则两时间序列负相关;R Y ^Y ＝０,则两时间序列不相关.同时,R Y ^Y 绝对值越接近１,两时间序列相关性越强.２㊀数据与试验２．１㊀数据选取本文试验数据均来自中国地震局G N S S 数据产品服务平台,试验数据集包含了中国大陆构造环境监测网络中B J G B 站(１１９．１ʎE ,４０．６ʎN )㊁H E Z J 站(１１４．８ʎE ,４０．８ʎN )㊁X J S S 站(９０．２ʎE ,４２．８ʎN )和Y N R L 站(９７．８ʎE ,２４．０ʎN )共４个观测站各１０９５个历元的单日解高程数据.试验按照２ʒ１的比例划分训练集和测试集,即训练集包含各测站７３０个历元数据,测试集包含各测站３６５个历元数据.２．２㊀构建V M D ＧX G B o o s t 模型VM D ＧX G B o o s t 模型预测流程如图１所示.基于VM D 模型改进的G N S S 高程时间序列预测模型(VM D ＧX G B o o s t 模型)具体预测步骤如下:(１)数据准备.G N S S 坐标时间序列是通过实际观测或求解得到的,它应在周㊁天㊁小时㊁秒等维度具有一致性.本文选取４个G N S S 观测站单日解高程时间序列数据作为试验数据.８３２１Copyright ©博看网. All Rights Reserved.第８期鲁铁定,等:融合V M D和X G B o o s t算法的G N S S高程时间序列预测方法图１㊀VM DＧX G B o o s t模型流程F i g．１㊀F l o w c h a r t o fVM DＧXG B o o s tm o d e l㊀㊀(２)通过V M D模型构造子序列.首先通过V M D模型进行时间序列分解,然后将分解得到的时间序列叠加得到子时间序列.试验中该步骤通过设置不同的K值得到５个子时间序列.为了保证试验的科学性,训练集和测试集数据需要分别处理,试验按照２ʒ１的比例划分训练集和测试集.(３)构造数据集.将子时间序列和原始时间序列放入同一数据集,将子时间序列作为原始时间序列的特征.(４)X G B o o s t模型预测.首先将试验数据集中的训练集输入X G B o o s t模型中,同时为了更好地体现模型的泛化能力,在训练集中随机抽取２０％个历元进行预测并通过五折交叉验证得到最终预测模型及输入特征的特征评价结果.然后将试验数据集输入训练好的X G B o o s t模型中.(５)统计V M DＧX G B o o s t模型预测结果.通过设置不同的预测时间跨度统计该时间跨度内预测精度指标值,以评判模型预测的稳定度.试验共设置了７㊁３０㊁９０㊁１８０㊁３６５d共５个时间跨度进行统计.３㊀试验结果与分析３．１㊀试验结果３．１．１㊀信号重构结果试验中K值的选取会直接影响特征集的构成,从而影响预测结果.K值较大时,相邻模态分量的中心频率则会相距较近,导致模态重复或产生额外的噪声,因此试验中K值取值不应过大.若追求预测效率,研究人员可以缩小K值的取值范围,减轻模型载荷;若追求更高的预测精度,研究人员可以适当扩大K值的取值范围,但并不是K值越大就必定会取得更高精度的预测结果.因此,顾及模型载荷会随着K值的增大而呈梯度上升的问题,试验将VM D模型的参数K 分别设置为４㊁５㊁６㊁７㊁８,从而得到不同的５个子时间序列.图２为Y N R L站信号重构结果.图２㊀Y N R L站信号重构结果F i g．２㊀S i g n a l r e c o n s t r u c t i o n r e s u l t s o fY N R Ls t a t i o n图２包含了５个VM D子模型信号重构结果和Y N R L站原始时间序列数据,每条曲线包含１０９５个历元.蓝色曲线为Y N R L站原始时间序列;橘黄色曲线为K＝４时,VM D模型信号重构结果;绿色曲线为K＝５时,VM D模型信号重构结果;红色曲线为K＝６时,VM D模型信号重构结果;紫色曲线为K＝７时,VM D模型信号重构结果;棕色曲线为K＝８时,VM D模型信号重构结果;红色直线为训练集和测试集分界线,左侧为训练集,右侧为测试集.３．１．２㊀预测结果根据试验设置,试验通过训练７３０个历元数据预测３６５个历元数据.图３为Y N R L站VM DＧX G B o o s t模型预测结果.图３㊀Y N R L站VM DＧX G B o o s t模型预测结果F i g．３㊀VM DＧXG B o o s t m o d e lf o r e c a s t i n g r e s u l t sa tY N R Ls t a t i o n９３２１Copyright©博看网. All Rights Reserved.A u gu s t ２０２３V o l ．５２N o ．８A G C S h t t p :ʊx b ．c h i n a s m p ．c o m 图３中蓝色曲线为Y N R L 站原始时间序列数据,共１０９５个历元;橘黄色曲线为训练集中VM D ＧX G B o o s t 模型拟合结果,共７３０个历元;绿色曲线为测试集中VM D ＧX G B o o s t 模型预测结果,共３６５个历元;红色直线为训练集和测试集分界线.３．２㊀试验分析３．２．１㊀重构信号质量分析试验以MA E 值和R M S E 值作为指标评价VM D 子模型重构信号的质量,通过皮尔森相关系数R 评价重构信号和原始时间序列间相关性.同时试验通过在训练集中对２０％个随机历元进行预测,得到不同子时间序列在预测模型中的得分情况,从而进行重构信号质量分析.表１为４个观测站重构信号精度.表１㊀各观测站重构信号精度T a b ．１㊀R e c o n s t r u c t i o n s i g n a l a c c u r a c y of e a c ho b s e r v a t i o n s t a t i o n测站KMA E /mm R M S E /mmRB J G B４２．６７３．４９０．９０５２．２４３．０００．９３６２．０１２．６７０．９４７１．６６２．１７０．９７８１．４３１．８７０．９８H E Z J４２．５９３．４００．９０５２．１８２．８９０．９３６１．８５２．４００．９６７１．６８２．１７０．９７８１．３４１．７５０．９８X J S S４２．７０３．６６０．９０５２．２９３．０６０．９３６１．９９２．６７０．９５７１．７３２．２８０．９７８１．５３２．０１０．９７Y N R L４３．１８４．１１０．９３５２．７７３．５９０．９５６２．４０３．０７０．９６７２．１０２．７００．９７８１．８０２．３００．９８由表１可知,随着K 值的增大,重构信号越接近原始时间序列,且重构信号均与原始时间序列呈极强相关性.虽然K 越大,重构信号误差越小,但当重构信号作为特征进行预测时并非与原始时间序列相关性越高,对于预测的贡献越大.通过X G B o o s t 算法内置的特征评价模块统计后,各观测站重构信号在训练集中得分情况如图４ 图７所示.图４㊀B J G B 站重构信号得分情况F i g ．４㊀S c o r e o f r e c o n s t r u c t e d s i gn a l o fB J G Bs t a t i on 图５㊀H E Z J 站重构信号得分情况F i g ．５㊀S c o r e o f r e c o n s t r u c t e d s i gn a l o fH E Z J s t a t i on 图６㊀X J S S 站重构信号得分情况F i g ．６㊀S c o r e o f r e c o n s t r u c t e d s i gn a l o fX J S Ss t a t i o n ㊀㊀由图４ 图７可知,当K ＝４或８时,重构信号得分最高,即对于模型预测提供了更为重要的帮助,这也证明了试验将K 的取值范围设置为４~８的合理性.３．２．２㊀预测结果分析C N N ＧL S TM 模型在多个领域的时间序列预测研究中得到良好的应用[３４Ｇ３６],因此,试验以VM D ＧX G B o o s t 模型的构建思路构建VM D ＧC N N ＧL S T M 模型作为参照模型.０４２１Copyright ©博看网. All Rights Reserved.第８期鲁铁定,等:融合V M D和X G B o o s t算法的G N S S高程时间序列预测方法图７㊀Y N R L站重构信号得分情况F i g．７㊀S c o r e o f r e c o n s t r u c t e d s i g n a l o fY N R Ls t a t i o n㊀㊀试验以MA E值和R M S E值作为指标评价VM DＧX G B o o s t模型的预测精度,通过皮尔森相关系数R评价预测时间序列和原始时间序列间相关性.在进行４个观测站高程时间序列预测试验时,通过训练７３０个历元数据,预测３６５个历元的数据.表２为各观测站预测结果精度.表２㊀各观测站预测结果精度T a b．２㊀F o r e c a s t i n g a c c u r a c y o f e a c ho b s e r v a t i o n s t a t i o n测站模型预测天数/dMA E/mmR M S E/mm RB J G B VM DＧX G B o o s tVM DＧC N NＧL S T M ㊀７０．９７１．４３０．８７３００．９９１．２８０．９８９０１．３４１．６９０．９６１８０１．９１２．３９０．９４３６５１．６３２．１００．９７３６５２．１０２．７００．９２H E Z J VM DＧX G B o o s tVM DＧC N NＧL S T M ㊀７０．９２１．１２０．９４３０１．０２１．２３０．９７９０１．３４１．６７０．９６１８０２．４１３．０９０．９５３６５１．７５２．４１０．９７３６５２．７３３．５００．９２X J S S VM DＧX G B o o s tVM DＧC N NＧL S T M ㊀７２．７６２．８４０．７７３０２．０６２．３２０．９５９０１．９２２．２１０．９４１８０２．２９２．６９０．９４３６５１．８７２．３３０．９７３６５２．３３３．１３０．９２Y N R L VM DＧX G B o o s tVM DＧC N NＧL S T M ㊀７１．６２１．９２０．９４３０１．２２１．５９０．９５９０１．３９１．７１０．９６１８０２．６３３．６８０．９４３６５２．１１２．９８０．９８３６５３．０２４．０６０．９５通过对比VM DＧX G B o o s t模型和VM DＧC N NＧL S T M模型的预测结果可以得出,在进行长时间跨度的预测时,V M DＧX G B o o s t模型精度更高.相较于V M DＧC N NＧL S T M模型,在B J G B㊁H E Z J㊁X J S S㊁Y N R L４个站,V M DＧX G B o o s t模型预测结果的MA E值分别降低了２２．３８％㊁３５．９０％㊁１９．７４％㊁３０．１３％,R M S E值分别降低了２２．２２％㊁３１．１４％㊁２５．５６％㊁２６．６０％.因此,在进行G N S S 高程时间序列预测时,X G B o o s t模型更为稳定.由表２可知,当以VM D模型的重构信号作为特征时,X G B o o s t算法在不同时间跨度的时间序列预测中表现出良好的预测能力.通过对比VM DＧX G B o o s t模型在４个观测站不同时间跨度预测结果的精度可以得出,VM DＧX G B o o s t模型具有较强的稳定性,虽然预测精度会随着时间跨度的增长出现波动,但仍然保持着较高的预测精度,且预测结果与原始时间序列间相关性略有提升.因此,VM DＧX G B o o s t模型对于观测站高程方向上的运动趋势和变化规律拥有良好的预测效果.由表２可知,随着预测天数的增加,皮尔森相关系数R呈现先增大后减小再增大的现象.这一现象与预测时间跨度和预测误差存在关联性.首先,R第一次增大出现在预测时间跨度由７变为３０时,这是由于当预测时间跨度为７时,预测点数较少,每个点对于相关系数影响较大;在预测时间跨度由３０~３６５变化时,R出现先减小后增大的现象,通过对比预测误差的变化可以得出,模型在预测第９０~１８０个历元的数据时误差较大,从而导致相关系数降低,但由于预测时间跨度的增长,大部分历元的预测精度较高,所以虽然出现减少的想象,但不会低于初始值,然后随着预测时间跨度的增长,模型依旧保持着良好的预测精度,削弱了预测较差部分对于R的影响,R呈增大趋势.此外,通过对比X G B o o s t模型在４个测站的预测结果可以得出,模型在第９０~１８０个历元区间出现较大预测误差的因素并非是时间跨度,而是时间序列的峰值在该区间占比较高.因此, VM DＧX G B o o s t模型对于峰值的预测精度有待提升.表３为H E Z J站预测精度最差的１０个历元.由表３可知,H E Z J站预测精度最差的１０个历元分布在第１０６~１６０个历元,验证了在H E Z J 站VM DＧX G B o o s t模型预测精度出现下滑的区间.图８为VM DＧX G B o o s t模型在H E Z J站的预测结果.１４２１Copyright©博看网. All Rights Reserved.A u gu s t ２０２３V o l ．５２N o ．８A G C S h t t p :ʊx b ．c h i n a s m p ．c o m 表３㊀H E Z J 站预测精度最差的１０个历元T a b ．３㊀T h e t e ne p o c h sw i t h t h ew o r s t f o r e c a s t i n g a c c u r a c yo fH E Z J s t a t i o nm测站历元真实值预测值绝对误差H E Z J１３２０．０１２５４０．０１８８２０．００６２８１５１０．０１５９００．０２２１９０．００６２９１６００．００３４３０．００９８８０．００６４５１３１０．０１２０７０．０１８８２０．００６７５１５８－０．００１５００．００５３４０．００６８４１３０－０．００１１８０．００５９９０．００７１７１０６０．０１４９５０．０２２１９０．００７２４１４４０．０１４１４０．０２２１９０．００８０５１４９０．０１０４５０．０１８８７０．００８４２１３６０．０１３４９０．０２２１９０．００８７０图８㊀H E Z J 站预测结果F i g ．８㊀F o r e c a s t i n g re s u l t s o fH E Z J s t a t i o n 图８中,蓝色曲线为H E Z J 站原始原始时间序列;绿色曲线为VM D ＧX G B o o s t 模型预测结果;１０个红色星号标记为VM D ＧX G B o o s t 模型在H E Z J 站预测精度最差的１０个点.由图８可以看出,预测结果与原始时间序列保持较好的一致性,对于原始时间序列的上升或下降趋势有着良好的预见性.根据红色星号标记分布规律可以得到,预测精度最差的１０个点虽然不是原始时间序列的峰值点,但都分布在时间序列峰值附近,因此,VM D ＧX G B o o s t 模型预测误差主要来源于峰值区域.上述试验说明了模型在高相关性特征的作用下可以得到高精度的结果,但上述试验使用目标时间序列构造特征.因此为了完善试验设置,试验通过VM D 算法重构与目标时间序列邻近的观测站数据,将重构得到的低噪声时间序列作为特征进行目标时间序列的预测.表４为４个观测站的预测结果精度.表４中,X J S S 和Y N R L 站预测结果的MA E值和R M S E 值较大的原因是X J S S 和Y N R L 站真实观测值的绝对值较大.对比表２和表４可以看出,VM D ＧX G B o o s t 模型在４个测站的预测精度均有降低,这是因为试验初期为了验证模型的性能,试验通过目标序列构造特征,使得特征中含有目标时间序列信息.后续试验中,为了完善预测试验的严谨性,试验通过邻近观测站获得特征信息,所以预测精度出现下滑,这说明了特征与目标时间序列的相关性是影响预测精度的因素之一,构造与目标时间序列相近的特征是VM D ＧX G B o o s t 模型的关键.表４㊀４个观测站的预测结果T a b ．４㊀F o r e c a s t i n g re s u l t s of f o u rG N S S s t a t i o n s mm测站特征观测站MA E R M S E RB J G B T J B D３．２７４．３２０．８３H E Z JH E C C ３．１６４．０８０．８６X J S S X J B E ４．４４５．５４０．７５Y N R LY N T H ６．４６８．３２０．８３４㊀结㊀论本文提出了一种融合VM D 和X G B o o s t 算法的G N S S 高程时间序列预测方法,即VM D ＧX G B o o s t 模型,给出了预测模型的数据处理策略,深入研究了VM D ＧX G B o o s t 模型的特点和优势.本文首先分析了传统时间序列预测模式存在的问题,在此基础上,构建了VM D ＧX G B o o s t 模型;然后评定了通过VM D 算法取代X G B o o s t 算法特征提取模块的有效性,并根据特征得分情况对K 值的选取进行分析;最后基于VM D ＧX G B o o s t 模型进行G N S S 高程时间序列预测试验.结果表明,在不同的时间跨度预测中,VM D ＧX G B o o s t 模型保持着良好的稳定性,预测结果与原始时间序列相关性较强;VM D ＧX G B o o s t 模型具有较高的精度,以VM D ＧC N N ＧL S T M 模型作为对比,精度提升可以达到１９．７４％~３５．９０％.特征相关性是影响VM D ＧX G B o o s t 模型预测精度的因素之一,多源的地学数据同样可以为模型预测提供有效特征[３７].因此,在后续研究中需要针对如何科学有效地获取高相关性特征展开进一步研究,以提高预测框架的适用性.参考文献:[１]㊀D E N GL i a n s h e n g ,J I A N G W e i p i n g,L IZ h a o ,e t a l ．A s s e s s Ｇ２４２１Copyright ©博看网. All Rights Reserved.第８期鲁铁定,等:融合V M D和X G B o o s t算法的G N S S高程时间序列预测方法m e n t o f s e c o n dＧa n dt h i r dＧo r d e r i o n o s p h e r i ce f f e c t so nr eＧg i o n a l n e t w o r k s:c a s e s t u d y i n C h i n a w i t h l o n g e rC MO N O CG P S c o o r d i n a t e t i m e s e r i e s[J]．J o u r n a l o f G e o dＧe s y,２０１７,９１(２):２０７Ｇ２２７．[２]㊀WU W e i w e i,WUJ i c a n g,M E N GG u o j i e．As t u d y o f r a n kd e f e c t a n d n e t w o r k e f f e c t i n p r o c e s s i n g t h eC MO N O Cn e tＧw o r ko nB e r n e s e[J]．R e m o t eS e n s i n g,２０１８,１０(３):３５７．[３]㊀姜卫平,王锴华,李昭,等．G N S S坐标时间序列分析理论与方法及展望[J]．武汉大学学报(信息科学版),２０１８,４３(１２):２１１２Ｇ２１２３．J I A N G W e i p i n g,W A N G K a i h u a,L IZ h a o,e ta l．P r o s p e c ta n d t h e o r y o fG N S Sc o o r d i n a t et i m es e r i e sa n a l y s i s[J]．G e o m a t i c s a n d I n f o r m a t i o nS c i e n c eo fW u h a nU n i v e r s i t y,２０１８,４３(１２):２１１２Ｇ２１２３．[４]㊀Y A OY i b i n,Y A N GY u a n x i,S U N H e p i n g,e t a l．G e o d e s yd i s c i p l i n e:p r o g re s s a n d p e r s p e c t i v e[J]．J o u r n a l o fG e o d e s y a n dG e o i n f o r m a t i o nS c i e n c e,２０２１,４(４):１Ｇ１０．[５]㊀R E N Y i n g y i n g,L I A N L i z h e n,WA N GJ i e x i a n．A n a l y s i s o fs e i s m i cd e f o r m a t i o nf r o m g l o b a lt h r e eＧd e c a d e G N S Sd i s p l a ce m e n t s:i m p l i c a t i o n sf o rat h r e eＧd i m e n s i o n a l e a r t hG N S S v e l o c i t y f i e l d[J]．R e m o t e S e n s i n g,２０２１,１３(１７):３３６９．[６]㊀W A N G J i a n,J I A N G W e i p i n g,L I Z h a o,e t a l．An e w m u l t iＧs c a l e s l i d i n g w i n d o w L S T Mf r a m e w o r k(M S S WＧL S T M):a c a s e s t u d y f o rG N S St i m eＧs e r i e s p r e d i c t i o n[J]．R e m o t eS e n s i n g,２０２１,１３(１６):３３２８．[７]㊀S T A L L E R A,ÁL V A R E ZＧGÓM E Z JA,L U N A MP,e t a l．C r u s t a lm o t i o na n dd e f o r m a t i o ni nE c u a d o r f r o m c G N S St i m e s e r i e s[J]．J o u r n a l o f S o u t hA m e r i c a nE a r t hS c i e n c e s,２０１８,８６:９４Ｇ１０９．[８]㊀H O B B SB,O R DA．N o n l i n e a r d y n a m i c a l a n a l y s i s o fG N S Sd a t a:q u a n t i f i c a t i o n,p re c u r s o r sa n ds y n c h r o n i s a t i o n[J]．P r o g r e s s i n E a r t ha n d P l a n e t a r y S c i e n c e,２０１８,５(１):１Ｇ３５．[９]㊀X U K,H E R,L IK,e ta l．S e c u l a rc r u s t a ld e f o r m a t i o nc h a r a c t e r i s t i c s p r i o rt ot h e２０１１T o h o k uＧO k ie a r t h q u a k ed e t e c t e d f r o m G N S Sa r r a y,２００３ ２０１１[J]．A d v a n c e s i nS p a c eR e s e a r c h,２０２１．[１０]㊀X IR,J I A N G W,M E N GX,e t a l．R a p i d i n i t i a l i z a t i o nm e t h o di n r e a lＧt i m e d e f o r m a t i o nm o n i t o r i n g o f b r i d g e sw i t h t r i p l eＧf r e q u e n c y B D Sa n d G P S m e a s u r e m e n t s[J]．A d v a n c e s i nS p a c eR e s e a r c h,２０１８,６２(５):９７６Ｇ９８９．[１１]㊀C H E NQ u s e n,J I A N G W e i p i n g,M E N GX i a o l i n,e t a l．V e r t i c a ld e f o r m a t i o n m o n i t o r i n g o ft h es u s p e n s i o n b r i d g et o w e ru s i n g G N S S:a c a s es t u d y o f t h e f o r t hr o a db r i d g e i nt h eU K[J]．R e m o t eS e n s i n g,２０１８,１０(３):３６４．[１２]㊀X I NJ i n g z h o u,Z H O UJ i a n t i n g,Y A N GS,e t a l．B r i d g e s t r u cＧt u r e d e f o r m a t i o n p r e d i c t i o n b a s e d o nG N S S d a t a u s i n g K a lＧm a nＧA R I MAＧG A R C H m o d e l[J]．S e n s o r s,２０１８,１８(１):２９８．[１３]㊀Z HA N G R u i c h e n g,G A O C h e n g f a,P A N S h u g u o,e ta l．F u s i o no fG N S Sa n ds p e e d o m e t e rb a s e do n VM Da n d i t sa p p l i c a t i o n i nb r i d g e d e f o r m a t i o nm o n i t o r i n g[J]．S e n s o r s,２０２０,２０(３):６９４．[１４]㊀A L T AM I M I Z,R E B I S C HU N GP,MÉT I V I E RL,e t a l．I T R F２０１４:an e wr e l e a s eo f t h e I n t e r n a t i o n a lT e r r e s t r i a lR e f e r e n c eF r a m e m o d e l i n g n o n l i n e a rs t a t i o n m o t i o n s[J]．J o u r n a l o fG e o p h y s i c a lR e s e a r c h:S o l i dE a r t h,２０１６,１２１(８):６１０９Ｇ６１３１．[１５]㊀L A H T I N E NS,J I V A L LL,HÄK L I P,e t a l．D e n s i f i c a t i o n o f t h e I T R F２０１４p o s i t i o na n dv e l o c i t y s o l u t i o n i n t h eN o r d i ca n dB a l t i cc o u n t r i e s[J]．G P S S o l u t i o n s,２０１９,２３(４):１Ｇ１３．[１６]㊀L I Z,C H E N W,D A M T,e t a l．C o m p a r a t i v e a n a l y s i s o f d i fＧf e r e n t a t m o s p h e r i c s u r f a c e p r e s s u r em o d e l sa n dt h e i r i mＧp a c t s o nd a i l y I T R F２０１４G N S Sr e s i d u a l t i m es e r i e s[J]．J o u r n a l o fG e o d e s y,２０２０,９４(４):１Ｇ２０．[１７]㊀K OWA L C Z Y KK,P A J A KK,W I E C Z O R E KB,e t a l．A na n a l y s i s o f v e r t i c a l c r u s t a lm o v e m e n t s a l o n g t h eE u r o p e a nc o a s t f r o ms a t e l l i t e a l t i m e t r y,t ide g a u g e,G N S S a n d r a d a ri n t e r f e r o m e t r y[J]．R e m o t e S e n s i n g,２０２１,１３(１１):２１７３．[１８]㊀鲁铁定,李祯．基于P r o p h e tＧX G B o o s t模型的G N S S高程时间序列预测[J]．大地测量与地球动力学,２０２２,４２(９):８９８Ｇ９０３．L U T i e d i n g,L I Z h e n．P r e d i c t i o no fG N S Sv e r t i c a l c o o r d iＧn a t et i m es e r i e sb a s e d o n p r o p h e tＧX G B o o s t m o d e l[J]．J o u r n a lo f G e o d e s y a n d G e o d y n a m i c s,２０２２,４２(９):８９８Ｇ９０３．[１９]㊀L I Z h e n,L UT i e d i n g．P r e d i c t i o n o fm u l t i s t a t i o nG N S S v e r t i c a lc o o rd i n a te t i m e s e r i e s b a s e do nX G B o o s t a l g o r i t h m[M]ʊL e c t u r e n o t e s i n e l e c t r i c a l e n g i n e e r i n g．S i n g a p o r e:S p r i n g e rN a t u r eS i n g a p o r e,２０２２:２７５Ｇ２８６．[２０]㊀L IW e n h a o,L IF e i,Z H A N GS h e n g k a i,e t a l．S p a t i o t e mＧp o r a lf i l t e r i n g a n d n o i s e a n a l y s i s f o r r e g i o n a l G N S Sn e t w o r k i n A n t a r c t i c a u s i n g i n d e p e n d e n t c o m p o n e n ta n a l y s i s[J]．R e m o t eS e n s i n g,２０１９,１１(４):３８６．[２１]㊀N I S T O RS,S U B A NS,MA C I U K K,e t a l．A n a l y s i so f n o i s e a n dv e l o c i t y i n G N S SE P NＧr e p r o２t i m es e r i e s[J]．R e m o t eS e n s i n g,２０２１,１３(１４):２７８３．[２２]㊀贺小星,花向红,鲁铁定,等．时间跨度对G P S坐标序列噪声模型及速度估计影响分析[J]．国防科技大学学报,２０１７,３９(６):１２Ｇ１８．H EX i a o x i n g,HU AX i a n g h o n g,L UT i e d i n g,e t a l．E f f e c to f t i m e s p a no nG P St i m e s e r i e sn o i s em o d e l a n dv e l o c i t ye s t i m a t i o n[J]．J o u r n a l o fN a t i o n a lU n i v e r s i t y o fD ef e n s eT e c h n o l o g y,２０１７,３９(６):１２Ｇ１８．[２３]㊀李威,鲁铁定,贺小星,等．基于P r o p h e tＧR F模型的G N S S 高程坐标时间序列预测分析[J]．大地测量与地球动力学,２０２１,４１(２)１１６Ｇ１２１．L IW e i,L U T i e d i n g,H E X i a o x i n g,e t a l．P r e d i c t i o na n da n a l y s i s o fG N S Sv e r t i c a l c o o r d i n a t e t i m es e r i e sb a s e do np r o p h e tＧR F m o d e l[J]．J o u r n a l o fG e o d e s y a n dG e o d y n a mＧi c s,２０２１,４１(２)１１６Ｇ１２１．[２４]㊀T A O R u i,L U T i e d i n g,C H E N G Y u a n m i n g,e t a l．A n i mＧ３４２１Copyright©博看网. All Rights Reserved.。
基于MFC和OpenGL的快速填充等值线实现
基于M FC和OpenGL的快速填充等值线实现孙晓燕1,伍增贵2(1.石油大学计算机与通信工程学院,山东东营257061;2.石油大学石油工程学院,山东东营257061)摘 要:通过等值线的绘制与填充过程分离的编程思想,先用OpenG L硬件加速的光栅化技术实现区域快速填充,然后利用MFC下的G D I绘图功能在相同的区域上绘制等值线。
该方法不涉及到复杂的算法,用很简单的代码就能实现与商业软件视觉效果相媲美的等值线填充效果,且适用于离散区域为任意形状的多边形网格系统。
关键词:等值线填充;OpenG L;可视化;MFC;调色板中图法分类号:TP391141 文献标识码:A 文章编号:100123695(2005)0320169202 Realizati on of Fast Filling Cont ours Based on MFC and OpenG LS UN Xiao2yan1,WU Zeng2gui2(1.College of Co m puter&Co mm unication Engineering,U niversity of Petroleum(East China),D ongying Shandong257061,China;2.College of Petroleum Engineering,U niversity of Petroleum(East China),D ongying Shandong257061,China)Abstract:This paper p resents a p r ogra mm ing idea of dra wing and filling cont ours res pectively.First the zone is filled by ta2 king advantage of rap id raster dis p lay of OpenG L,then on the sa me zone cont ourlines are generated by using G D I functi ons of MFC.W ithout any comp licated algorith m,the p r ogra m is easily realized,but p r ovides excellent visual effects comparable with commercial s oft w are.I n additi on,the method can be app lied t o vari ous grid syste m s.Key words:Cont ours Fill;OpenG L;V isualizati on;M FC;Palette 在工程计算分析中常常涉及到大量的参数场,如温度场、压力场、浓度场、速度场、应力场等数据场的分析和处理。
如何在马尔可夫链蒙特卡洛中处理缺失数据问题(Ⅲ)
马尔可夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）是一种用于随机模拟的方法，它在许多领域中都有广泛的应用，包括机器学习、统计学和物理学等。
在实际应用中，我们经常会遇到数据缺失的情况，这就需要针对缺失数据问题对MCMC进行适当的处理。
本文将讨论如何在MCMC中处理缺失数据问题。
首先，我们需要了解MCMC的基本原理。
MCMC是一种基于蒙特卡洛方法的统计推断技术，它通过构建一个马尔可夫链，从而可以对复杂的概率分布进行抽样。
在MCMC中，我们通常会使用马尔可夫链的转移核函数来生成样本，并利用这些样本来近似计算我们感兴趣的分布的期望值和方差等统计量。
然而，当我们的数据中存在缺失值时，MCMC的应用就会变得更加复杂。
因为在缺失数据的情况下，标准的MCMC算法可能会导致样本的偏误，从而影响我们对目标分布的估计。
因此，我们需要对MCMC进行适当的修改和调整，以解决缺失数据带来的问题。
一种处理缺失数据的方法是使用Gibbs采样。
Gibbs采样是MCMC的一种特殊形式，它可以有效地处理缺失数据，并且在实际应用中得到了广泛的应用。
在Gibbs采样中，我们将缺失的数据视为未知参数，并通过条件分布来进行采样。
通过交替地对每个缺失变量进行采样，我们可以逐步地减小参数空间，从而得到对未知参数的估计。
此外，我们还可以利用MCMC算法中的元算法来处理缺失数据。
元算法是一种用于加速MCMC收敛的技术，它可以有效地探索参数空间，并减少样本的自相关性。
在处理缺失数据时，我们可以利用元算法来优化参数的转移核函数，从而提高MCMC算法的采样效率。
通过优化参数的转移核函数，我们可以更好地利用数据中已有的信息，从而得到更准确的估计结果。
除了上述方法外，我们还可以考虑在MCMC中引入辅助变量来处理缺失数据。
辅助变量是一种在统计模型中引入的人工变量，它可以帮助我们对缺失数据进行建模，并且在MCMC算法中起到一定的作用。
克里格空间插值法ppt课件
4．高斯模型（Gaussian model） 变程为 。
1.9 理论变异函数模型
图是球状模型、指数模型和高斯模型的比较，可以看出，球状模型的变程最小，指数的模型变程最大，高斯模型的变程介于二者之间。球状模型和指数模型过原点存在切线，高斯模型则没有。
1.9 理论变异函数模型
3．指数模型(Exponential model) 其中，d是控制方程空间范围的距离参数。这里，仅在无穷远处相关性完全消失。变程为3d。指数模型在统计理论中地位重要，它表示了空间随机性的要素，是一阶自回归和马尔可夫过程的半方差函数。作为自相关函数，它们是采样设计有效性的理论基础。
1.4邻域函数的统计函数及其意义
摄影测量得到的正射航片或卫星影象； 卫星或航天飞机的扫描影象； 野外测量采样数据，采样点随机分布或有规律的线性分布（沿剖面线或沿等高线； 数字化的多边形图、等值线图；
1.5 空间插值的数据源
图1 各种不同的采样布置方式
1.6 采样布置方式
1.8 方差变异函数
2）曲线从较低的方差值升高，到一定的间隔值时到达基台值，这一间隔称为变程（range）。在理论函数模型中，变程用a表示。 变程是半方差函数中最重要的参数，它描述了该间隔内样点的空间相关特征。在变程内，样点越接近，两点之间相似性、即空间上的相关性越强。很明显，如果某点与已知点距离大于变程，那么该点数据不能用于数据内插（或外推），因为空间上的自相关性不复存在。 变程的高低取决于观测的尺度，说明了相互作用所影响的范围。不同的属性，其变程值可以变化很大。
1.2.2局部插值方法 分类
1.4邻域函数的统计函数及其意义
众数（majority)：邻域中出现频率最高的数值 最大值(max)：邻域中最大的数值 最小值(min)：邻域中最小的数值 中位数(median)：邻域中数值从小到大排列后位于中间的数 平均值(mean)：邻域中数值的算术平均 频率最小数(minority)：邻域中出现频率最小的数值 范围(range)：邻域中数值的范围，最大值与最小值之差 标准差(std)：邻域中数值的标准差 和(sum)：邻域中数值的和 变异度(varity)：邻域中不同数值的个数
基于CS算法的Markov模型及收敛性分析
（ｃｏ ｌ ｆ ｃ ｎ ｅＸｉ ｎＰ ｌｔｃｎｃＵ ｉｅｓｙ Ｘｉ ｎ７０ ４ ， ｍａ Ｓｈ ｏ ｉ ｃ， ’ ｏｙｅｈ ｉ ｎｖｒｉ ， ’ １０ ８Ｃｈ ） ｏＳ ｅ ａ ｔ ａ
［ ｓ ａ ｔ ｒｅ ｅｆｃ ｅ ｏｖ ｒｅｃ ｅｒ ｆ ｃ ｏ ｅｒｈＣ ） ｌｏｉ ｍ，ｈ ｒｏ ｈ ｉ ｍｏ ｅ ｏ ｔｅ ｌｏｉ ｍ ｓａｌｈ ｄ Ａｂ ｔ ｃｌＩ ｏｄ ｒｏｐｒｅｔｈ ｎ ｅｇｎ ｅ ｏｙｏ Ｃｕｋ ｏＳａｃ（ Ｓ ａｇ ｒ ｒ ｎ ｔ ｔ ｃ ｈ ｔ ｈ ｔ ｔｅ ｋ ｖｃ ａｎ ｄ ｌ ｆ ａｇ ｒ Ｍａ ｈ ＣＳ ｈ ｓ ｔ ｉ ｅｔｂｉ ｅ ｓ
ａ ｄ ｔ ｅ ｐ ｏ ｅ ｔ ｅ ｌ ｉ ｄ ａ ｄ ｈ ｍ ｏ ｅ ｅ ｕ ｆＭ ａ ｋ ｖ ｃ ａ ｎ ｉ ｎ ｌ ｚ ｄ ｎ ｒ ｐ ｒｙ ｏｆｔ ｉ ｔ ｎ ｏ ｈ ｈ ｍ ｅ ｇ ｎ ｏ ｓ ｏ ｒ ｏ ｈ ｉ ｓ ａ ａ ｙ ｅ ．Ｏｎ ｔ ｅ ｂ ｓ ｓ ｏ ｉ，ｔ ｒ ｕ ｈ ｔ ｅ ａ ａ ｙ ｉ ｆｔ ｅ ｓ ａｅ ｔａ ｓｔ ｎ ａ ｉ ｆｔ ｓ ｈ ｏ ｇ ｎ ｌ ｓ ｓ ｏ ｈ ｔ ｔ ｒ ｉ ｏ ｈ ｈ ｈ ｎ ｉ
析鸟 窝位置 的群体 状态 转移过程 ，指出随机序列将进入最优状态集 ，同时证明 ＣＳ算法满足随机搜索算法全局收敛 的 ２个条件 。通过仿真 实验验证 Ｃ Ｓ算法 可收敛于全局最优 ，从而确保 ｃ ｓ算法的全局收敛性。
关健词 ：启发式算法 ；布谷鸟搜索 ；Ｍａｋｖ链 ；状态转移 ；全局收敛性 ｒｏ
［ ｙｗｏｄ ｌｈｕｉ ｉ ａ ｏｉ ｍ； ｕ ｋ ｏＳ ａｃ （ Ｓ； ｒｋｖｃ ａｎ ｓ ｔ ｔ ｎｉｏ ； ｌｂ ｌ ｏ ｖｒｅｃ Ｋｅ ｒ ｓ ｅ ｓ ｃ ｌ ｒ ｈ Ｃ ｃｏ ｅｒｈＣ ）Ｍａ ｏ ｈ ｉ；ｔ ｅｒ ｓ ｉｎ ｇｏ ａ ｎ ｅｇｎ ｅ ｒ ｔ ｇ ｔ ａ ａ ｔ ｃ ＤＯＩ １．９ ９ ０ —４ ８２ １．１ ５ ｊｓ ０
常见插值算法--拉格朗日插值、三次卷积插值、三次样条插值、兰克索斯插值
常见插值算法--拉格朗⽇插值、三次卷积插值、三次样条插值、兰克索斯插值写在前⾯本⽂简单介绍了⼏种常见的插值算法并附带了相应的python代码，本⽂公式使⽤latex编写，如有错误欢迎评论指出，如果谁知道如何修改latex字号也欢迎留⾔关于⼀维、⼆维和多维插值三次卷积插值、拉格朗⽇两点插值（线性插值）、兰克索斯插值在⼆维插值时改变x和y⽅向的计算顺序不影响最终结果，这三个也是图像缩放插值时常⽤的插值算法，⽽其他插值在改变计算顺序时会产⽣明显差异，多维的情况笔者没有尝试，读者可以⾃⾏尝试或推导最近邻插值法（Nearest Neighbour Interpolation）在待求像素的四邻像素中，将距离待求像素最近的像素值赋给待求像素p_{11}p_{12}pp_{21}p_{22}python代码1def NN_interpolation(srcImg, dstH, dstW):2 scrH, scrW, _ = srcImg.shape3 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)4for i in range(dstH - 1):5for j in range(dstW - 1):6 scrX = round(i * (scrH / dstH))7 scrY = round(j * (scrW / dstW))8 dstImg[i, j] = srcImg[scrX, scrY]9return dstImg拉格朗⽇插值（Lagrange Interpolation）拉格朗⽇插值法需要找到k个p_i(x)函数，使得每个函数分别在在x_i处取值为1，其余点取值为0，则y_ip_i(x)可以保证在x_i处取值为y_i，在其余点取值为0，因此L_k(x)能恰好经过所有点，这样的多项式被称为拉格朗⽇插值多项式，记为L_k(x)=\sum_{i=1}^ky_ip_i(x)p_i(x)=\prod_{j \neq i}^{1 \leq j \leq k}\frac{x-x_j}{x_i-x_j}以四点即三次图像插值为例，因为横坐标间隔为1，则设四个点横坐标为-1、0、1和2，可得p_1(x)、p_2(x)、p_3(x)和p_4(x)假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得拉格朗⽇函数如下图所⽰，待插值点横坐标范围为[0,1]在K=2时在k=2时，也被称为线性插值通⽤公式p_1=\frac{x-x_2}{x_1-x_2}p_2=\frac{x-x_1}{x_2-x_1}\begin{align} L_2x &= p_1y_1+p_2y_2 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}pp_{21}p_{22}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_2x &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \\ &= (x_2-x)y_1+(x-x_1)y_2 \nonumber \\ &= (1-dx)y_1+dxy_2 \nonumber \\ &= (y_2-y_1)dx+y_1 \nonumber \end{align}L_2'x=y_2-y_1在K=3时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\begin{align} L_3x &= p_1y_1+p_2y_2+p_3y_3 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1+\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2+\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}p_{13}p_{21}p_{22}p_{23}pp_{31}p_{32}p_{33}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_3x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \\ &= \frac{-dx(1-dx)}{(-1)\cdot(-2)}y_1 + \frac{-(1+dx)(1-dx)}{1\cdot(-1)}y_2 + \frac{(1+dx)dx}{2\cdot 1}y_3 \nonumber \\ &= (\frac{1}{2}d^2x-\frac{1}{2}dx)y_1 - (d^2x-1)y_2 + (\frac{1}{2}d^2x+\frac{1}{2}dx)y_3 \nonumber \\ &= d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{2}y_1+\frac{1}{2}y_3)+y_2 \nonumber \end{align}L_3'x=dx(y_1-2y_2+y_3)+(\frac{1}{2}y_3-\frac{1}{2}y_1)在K=4时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}p_4=\frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}\begin{align} L_4x &= p_1y_1+p_2y_2+p_3y_3+p_4y_4 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3} {x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4\nonumber \end{align}图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_4x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4 \nonumber \\ &= \frac{dx[-(1-dx)][-(2-dx)]}{(-1)\cdot(-2)\cdot(-3)}y_1 + \frac{(1+dx)[-(1-dx)][-(2-dx)]}{1\cdot(-1)\cdot(-2)}y_2 + \frac{(1+dx)dx[-(2-dx)]}{2\cdot 1\cdot(-1)}y_3 + \frac{(1+dx)dx[-(1-dx)]}{3\cdot 2\cdot 1}y_4 \nonumber \\ &= \frac{d^3x-3d^2x+2dx}{-6}y1 + \frac{d^3x-2d^2x-dx+2}{2}y_2 + \frac{d^3x-d^2x-2dx}{-2}y_3 + \frac{d^3x-dx}{6}y_4 \nonumber \\ &= d^3x(-\frac{1}{6}y_1+\frac{1}{2}y_2-\frac{1} {2}y_3+\frac{1}{6}y_4)+d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4)+y_2 \nonumber \end{align}\begin{align} L_4'x &= d^2x(-\frac{1}{2}y_1+\frac{3}{2}y_2-\frac{3}{2}y_3+\frac{1}{2}y_4)+dx(y_1-2y_2+y_3)+(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4) \nonumber \\ &= -[\frac{1}{2}d^2x(y_1-3y_2+3y_3-y_4)-dx(y_1-2y_2+y_3)+\frac{1}{6}(2y_1+3y_2-6y_3+y_4)] \nonumber \end{align}python代码插值核计算的时候乘法和加减法计算的顺序不同可能会导致结果存在细微的差异，读者可以⾃⾏研究⼀下1class BiLagrangeInterpolation:2 @staticmethod3def LagrangeInterpolation2(x, y1, y2):4 f1 = 1 - x5 f2 = x6 result = y1 * f1 + y2 * f27return result89 @staticmethod10def LagrangeInterpolation3(x, y1, y2, y3):11 f1 = (x ** 2 - x) / 2.012 f2 = 1 - x ** 213 f3 = (x ** 2 + x) / 2.014 result = y1 * f1 + y2 * f2 + y3 * f315return result1617 @staticmethod18def LagrangeInterpolation4(x, y1, y2, y3, y4):19 f1 = - (x ** 3 - 3 * x ** 2 + 2 * x) / 6.020 f2 = (x ** 3 - 2 * x ** 2 - x + 2) / 2.021 f3 = - (x ** 3 - x ** 2 - 2 * x) / 2.022 f4 = (x ** 3 - x) / 6.023 result = y1 * f1 + y2 * f2 + y3 * f3 + y4 * f424return result2526def biLag2_2(self, srcImg, dstH, dstW):27 dstH, dstW = int(dstH), int(dstW)28 srcH, srcW, _ = srcImg.shape29 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')30 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)31for dstY in range(dstH):32for dstX in range(dstW):33for channel in [0, 1, 2]:34# p11 p1235# p36# p21 p2237# 储存为 p(y, x)38 p = [dstY * srcH / dstH, dstX * srcW / dstW]39 p11 = [math.floor(p[0]), math.floor(p[1])]40 p12 = [p11[0], p11[1] + 1]4142 p21 = [p11[0] + 1, p11[1]]43 p22 = [p21[0], p12[1]]4445 diff_y, diff_x = p[0] - p11[0], p[1] - p11[1]46 r1 = grangeInterpolation2(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel])47 r2 = grangeInterpolation2(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel])4849 c = grangeInterpolation2(diff_y, r1, r2)5051 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)52return dstImg5354def biLag3_3(self, srcImg, dstH, dstW):55 dstH, dstW = int(dstH), int(dstW)56 srcH, srcW, _ = srcImg.shape57 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')58 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)59for dstY in range(dstH):60for dstX in range(dstW):61for channel in [0, 1, 2]:62# p11 p12 p1363#64# p21 p22 p2365# p66# p31 p32 p3367# 储存为 p(y, x)68 p = [dstY * srcH / dstH, dstX * srcW / dstW]69 p22 = [math.floor(p[0]), math.floor(p[1])]70 p21 = [p22[0], p22[1] - 1]71 p23 = [p22[0], p22[1] + 1]7273 p11 = [p21[0] - 1, p21[1]]74 p12 = [p11[0], p22[1]]75 p13 = [p11[0], p23[1]]7677 p31 = [p21[0] + 1, p21[1]]78 p32 = [p31[0], p22[1]]79 p33 = [p31[0], p23[1]]8081 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]82 r1 = grangeInterpolation3(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel])83 r2 = grangeInterpolation3(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel])84 r3 = grangeInterpolation3(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel]) 8586 c = grangeInterpolation3(diff_y, r1, r2, r3)8788 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)89return dstImg9091def biLag4_4(self, srcImg, dstH, dstW):92 dstH, dstW = int(dstH), int(dstW)93 srcH, srcW, _ = srcImg.shape94 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')95 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)96for dstY in range(dstH):97for dstX in range(dstW):98for channel in [0, 1, 2]:99# p11 p12 p13 p14100#101# p21 p22 p23 p24102# p103# p31 p32 p33 p34104#105# p41 p42 p43 p44106# 储存为 p(y, x)107 p = [dstY * srcH / dstH, dstX * srcW / dstW]108 p22 = [math.floor(p[0]), math.floor(p[1])]109 p21 = [p22[0], p22[1] - 1]110 p23 = [p22[0], p22[1] + 1]111 p24 = [p22[0], p22[1] + 2]112113 p11 = [p21[0] - 1, p21[1]]114 p12 = [p11[0], p22[1]]115 p13 = [p11[0], p23[1]]116 p14 = [p11[0], p24[1]]117118 p31 = [p21[0] + 1, p21[1]]119 p32 = [p31[0], p22[1]]120 p33 = [p31[0], p23[1]]121 p34 = [p31[0], p24[1]]122123 p41 = [p21[0] + 2, p21[1]]124 p42 = [p41[0], p22[1]]125 p43 = [p41[0], p23[1]]126 p44 = [p41[0], p24[1]]127128 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]129 r1 = grangeInterpolation4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 130 r2 = grangeInterpolation4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 131 r3 = grangeInterpolation4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 132 r4 = grangeInterpolation4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 133134 c = grangeInterpolation4(diff_y, r1, r2, r3, r4)135136 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)137return dstImg三次卷积插值法（Cubic Convolution Interpolation）使⽤上图中的卷积核进⾏加权平均计算，卷积核为u(s)，四个等距（距离为1）的采样点记为x_0、x_1、x_2和x_3，采样数值记为y_0、y_1、y_2和y_3，且保证四个点均在[-2,2]区间上，计算得到g(x)，假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得三次卷积插值函数如下图所⽰，待插值点横坐标范围为[0,1]公式推导设u(s)=\begin{cases} A_1|s|^3+B_1|s|^2+C_1|s|+D_1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\because函数在s=0,1,2处连续\therefore\begin{cases} 1=u(0^+)=D_1 \\ 0=u(1^-)=A_1+B_1+C_1+D_1 \\ 0=u(1^+)=A_2+B_2+C_2+D_2 \\ 0=u(2^-)=8A_2+4B_2+2C_2+D_2 \end{cases} (1)\because函数在s=0,1,2处导函数连续\therefore\begin{cases} u'(0^-)=u'(0+) \\ u'(1^-)=u'(1+) \\ u'(2^-)=u'(2+)\end{cases} \Rightarrow \begin{cases} -C_1=C_1 \\ 3A_1+2B_1+C_1=3A_2+2B_2+C_2\\ 12A_2+4B_2+C+2=0 \end{cases} ~~~~ (2)联⽴⽅程组(1)(2)，设A_2=a，解得\begin{cases} A_1=a+2 \\ B_1=-(a+3) \\ C_1=0 \\ D_1=1 \\ A_2=a \\ B_2=-5a \\ C_2=8a \\ D_2=-4a \end{cases}\Rightarrow u(s)=\begin{cases} (a+2)|s|^3-(a+3)|s|^2+1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2\\ 1, &s=0 \\ 0, &otherwise \end{cases}\because g(x)=\sum_kC_ku(s+j-k), ~~~~k=j-1,j, j+1,j+2且0<s<1⼜\because \begin{cases}\begin{align} u(s+1)&=as^3-2as^2+as \nonumber \\ u(s)&=(a+2)s^3-(a+3)s^2+1 \nonumber \\ u(s-1)&=-(a+2)s^3+(2a+3)s^2-as \nonumber \\ u(s-2)&=-as^3+as^2 \nonumber \end{align}\end{cases}\begin{align} \therefore g(x) &= C_{j-1}u(s+1)+C_{j}u(s)+C_{j+1}u(s-1)+C_{j+2}u(s-2) \nonumber \\ &= C_{j-1}(as^3-2as^2+as)+C_j[(a+2)s^3-(a+3)s^2+1]+C_{j+1}[-(a+2)s^3+ (2a+3)s^2-as]+C_{j+2}[-a^3+as^2] \nonumber \\ &= s^3[aC_{j-1}+(a+2)C_j-(a+2)C_{j+1}-aC_{j+2}]+s^2[-2aC_{j-1}-(a+3)C_j+(2a+3)C_{j+1}+aC_{j+2}]+s[aC_{j-1}-aC_{j+1}]+C_j \nonumber \end{align} ~~(3)f在x_j处泰勒展开得到f(x)=f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+\cdots\therefore \begin{cases} f(x_{j+1})=f(x_j)+f'(x_j)(x_{j+1}-x_j)+\frac{1}{2}f''(x_j)(x_{j+1}-x_j)^2+\cdots \\ f(x_{j+2})=f(x_j)+f'(x_j)(x_{j+2}-x_j)+\frac{1}{2}f''(x_j)(x_{j+2}-x_j)^2+\cdots \\ f(x_{j-1})=f(x_j)+f'(x_j)(x_{j-1}-x_j)+\frac{1}{2}f''(x_j)(x_{j-1}-x_j)^2+\cdots \end{cases}令x_{j+1}-x_j=h\because x_{i+1}=x_i+1\therefore x_{j+2}-x_j=2h，x_{j-1}-x_j=-h\therefore \begin{cases} f(x_{j+2})=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+\cdots \\ f(x_{j+1})=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \\ f(x_{j-1})=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \end{cases}\therefore \begin{cases} c_{j-1}=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3) \\ c_j=f(x_j) \\ c_{j+1}=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3)\\ c_{j+2}=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+o(h^3) \end{cases} ~~ (4)将(4)代⼊(3)，得g(x)=-(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3+[(6a+3)hf'(x_j)+\frac{4a+3}{2}h^2f''(x_j)]s^2-2ahf'(x_j)s+f(x_j)+o(h^3)\because h=1,s=x-x_J\therefore sh=x-x_j\begin{align}\therefore f(x)&= f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+o(h^3) \nonumber \\ &= f(x_j)+f'(x_j)sh+\frac{1}{2}f''(x_j)s^2h^2+o(h^3) \nonumber \end{align}\therefore f(x)-g(x)=(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3-(2a+1)[3hf'(x_j)+h^2f''(x_j)]s^2+[(2a+1)hf'(x_j)]s+o(h^3)\because 期望f(x)-g(x)趋于0\therefore 2a+1=0 \Rightarrow a=-\frac{1}{2}\therefore u(s)=\begin{cases} \frac{3}{2}|s|^3-\frac{5}{2}|s|^2+1, &0<|s|<1 \\ -\frac{1}{2}|s|^3+\frac{5}{2}|s|^2-4|s|+2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\therefore g(s)=s^3[-\frac{1}{2}c_{j-1}+\frac{3}{2}c_j-\frac{3}{2}c_{j+1}+\frac{1}{2}c_{j+2}]+s^2[c_{j-1}-\frac{5}{2}c_j+2c_{j+1}-\frac{1}{2}c_{j+2}]+s[-\frac{1}{2}c_{j-1}+\frac{1} {2}c_{j+1}]+c_j图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}python代码1class BiCubicConvInterpolation:2 @staticmethod3def CubicConvInterpolation1(p0, p1, p2, p3, s):4# ⽤g（s）公式计算，已经将四个u（s）计算完毕并整理5# as^3 + bs^2 + cs + d6 a = 0.5 * (-p0 + 3.0 * p1 - 3.0 * p2 + p3)7 b = 0.5 * (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3)8 c = 0.5 * (-p0 + p2)9 d = p110return d + s * (c + s * (b + s * a))1112 @staticmethod13def CubicConvInterpolation2(s):14# ⽤u（s）公式计算15 s = abs(s)16if s <= 1:17return 1.5 * s ** 3 - 2.5 * s ** 2 + 118elif s <= 2:19return -0.5 * s ** 3 + 2.5 * s ** 2 - 4 * s + 220else:21return 02223def biCubic1(self, srcImg, dstH, dstW):24# p11 p12 p13 p1425#26# p21 p22 p23 p2427# p28# p31 p32 p33 p3429#30# p41 p42 p43 p4431 dstH, dstW = int(dstH), int(dstW)32 scrH, scrW, _ = srcImg.shape33 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')34 dstImg = np.zeros((dstH, dstW, 1), dtype=np.uint8)35for dstY in range(dstH):36for dstX in range(dstW):37for channel in [0]:38 y = dstY * scrH / dstH39 x = dstX * scrW / dstW40 y1 = math.floor(y)41 x1 = math.floor(x)4243 array = []44for i in [-1, 0, 1, 2]:45 temp = self.CubicConvInterpolation1(srcImg[y1 + i, x1 - 1, channel],46 srcImg[y1 + i, x1, channel],47 srcImg[y1 + i, x1 + 1, channel],48 srcImg[y1 + i, x1 + 2, channel],49 x - x1)50 array.append(temp)5152 temp = self.CubicConvInterpolation1(array[0], array[1], array[2], array[3], y - y1)53 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)5455return dstImg5657def biCubic2(self, srcImg, dstH, dstW):58# p11 p12 p13 p1459#60# p21 p22 p23 p2461# p62# p31 p32 p33 p3463#64# p41 p42 p43 p4465 dstH, dstW = int(dstH), int(dstW)66 scrH, scrW, _ = srcImg.shape67 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')68 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)69for dstY in range(dstH):70for dstX in range(dstW):71for channel in [0, 1, 2]:72 y = dstY * scrH / dstH73 x = dstX * scrW / dstW74 y1 = math.floor(y)75 x1 = math.floor(x)7677 array = []78for i in [-1, 0, 1, 2]:79 temp = 080for j in [-1, 0, 1, 2]:81 temp += srcImg[y1 + i, x1 + j, channel] * self.CubicConvInterpolation2(x - (x1 + j))82 array.append(temp)8384 temp = 085for i in [-1, 0, 1, 2]:86 temp += array[i + 1] * self.CubicConvInterpolation2(y - (y1 + i))87 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)8889return dstImg三次样条插值在n-1个区间上寻找n-1个三次曲线，使其满⾜相邻曲线在端点处值相等、⼀阶导数相等，⼆阶导数相等，在加以边界条件后可得每个曲线的⽅程，然后沿x轴依次偏移对应的距离即可得到插值结果，如仅需要特定范围内的结果，则可以⼤幅减少计算量公式推导设S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3, ~~~~i=0,1,...,n-1则 \begin{cases} S_i'(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i)^2\\ S_i''(x)=2c_i+6d_i(x-x_i)\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1设h_i(x)=x_{i+1}-x_i，可得\begin{cases} S_i(x)=a_i+b_ih_i+c_ih_i^2+d_ih_i^3\\ S_i'(x)=b_i+2c_ih_i+3d_ih_i^2\\ S_i''(x)=2c_i+6d_ih_i\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1\because S_i(x)过点(x_i,y_i)\therefore S_i(x)=a_i=y+i, ~~~~i=0,1,...,n-1 ~~~~~~(1)\because S_i(x)与S_{i+1}(x)在X_{i+1}处相等\therefore S_i(x_{i+1})=S_{i+1}(x_{i+1})\Rightarrow a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1}, ~~~~i=0,1,...,n-2~~~~~~(2)\because S_i'(x)与S_{i+1}'(x)在X_{i+1}处相等\therefore S_i'(x)-S_{i+1}'(x)=0\Rightarrow b_i+2c_ih_i+3d_ih_i^2-b_{i+1}=0~~~~~~(3)\because S_i''(x)与S_{i+1}''(x)在X_{i+1}处相等\therefore S_i''(x)-S_{i+1}''(x)=0\Rightarrow 2c_i+6d_ih_i-2c_{i+1}=0, ~~~~i=0,1,...,n-2~~~~~~(4)设m_i=S_i(x_i)=2c_i，即c_i=\frac{1}{2}m_i, ~~~~i=0,1,...,n-1~~~~~~(5)将(5)代⼊(4)，得2c_i+6d_ih_i-2c_{i+1}=0\Rightarrow m_i+6h_id_i-m_{i+1}=0\Rightarrow d_i=\frac{m_{i+1}-m_i}{6h_i}, ~~~~i=0,1,...,n-2~~~~~~(6)将(1)(5)(6)代⼊(2)，得\begin{align} &a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1} \nonumber \\ \Rightarrow&y_i+b_ih_i+\frac{1}{2}m_ih_i^2+\frac{m_{i+1}-m_i}{6h_i}h_i^3=y_{i+1} \nonumber \\\Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{2}m_ih_i-\frac{1}{6}(m_{i+1}-m_i)h_i \nonumber \\ \Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i, ~~~~i=0,1,...,n-2~~~~~~(7) \nonumber \end{align}将(5)(6)(7)代⼊(3)，得\begin{align} &\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+2\cdot\frac{1}{2}m_ih_i+3\frac{m_{i+1}-m_i}{6h_i}h_i^2-(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{1}{3}m_{i+1}h_{i+1}-\frac{1}{6}m_{i+2}h_{i+1})=0 \nonumber \\ \Rightarrow&\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+m_ih_i+\frac{1}{2}(m_{i+1}-m_i)h_i-\frac{y_{i+2}-y_{i+1}}{h_{i+1}}+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=0 \nonumber \\ \Rightarrow&m_ih_i(-\frac{1}{3}+1-\frac{1}{2})+m_{i+1}h_i(-\frac{1}{6}+\frac{1} {2})+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&\frac{1}{6}(m_ih_i+2m_{i+1}h_i+2m_{i+1}h_{i+1}+m_{i+2}h_{i+1})=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&m_ih_i+2m_{i+1}(h_i+h_{i+1})+m_{i+2}h_{i+1}=6(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}}), ~~~~i=0,1,...,n-2~~~~~~(8) \nonumber \end{align}由(8)可知i=0,1,...,n-2，则有m_0,m_1,...,m_n，需要两个额外条件⽅程组才有解⾃然边界(Natural)m_0=0，m_n=0\begin{bmatrix} \tiny 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}固定边界(Clamped)\begin{align} &\begin{cases} S_0'(x_0)=A\\ S_{n-1}'(x_n)=B \end{cases} \nonumber \\ \Rightarrow&\begin{cases} b_0=A\\ b_{n-1}+2c_{n-1}h_{n-1}+3d_{n-1}h_{n-1}^2=B\end{cases} \nonumber \\ \Rightarrow&\begin{cases} A=\frac{y_1-y_0}{h_0}-\frac{h_0}{2}m_0-\frac{h_0}{6}(m_1-m_0)\\ B=\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{1}{3}m_{n-1}h_{n-1}+m_{n-1}h_{n-1}+\frac{1}{2}m_nh_{n-1}-\frac{1}{2}m_{n-1}h_{n-1} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 2h_0m_0+h_0m_1=6(\frac{y_1-y_0}{h_0}-A)\\ h_{n-1}m_{n-1}+2h_{n-1}m_{n}=6(B-\frac{y_n-y_{n-1}}{h_{n-1}}) \end{cases} \nonumber \\ \end{align}\begin{bmatrix} 2 & 1 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 1 & 2 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ B-\frac{y_n-y_{n-1}}{h_{n-1}} \end{bmatrix}⾮节点边界(Not-A-Knot)\begin{align} &\begin{cases} S_0'''(x_1)=S_1'''(x_1)\\ S_{n-2}'''(x_{n-1})=S_{n-1}'''(x_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 6\cdot\frac{m_1-m_0}{6h_0}=6\cdot\frac{m_2-m_1}{6h_1}\\ 6\cdot\frac{m_{n-1}-m_{n-2}}{6h_{n-2}}=6\cdot\frac{m_n-m_{n-1}}{6h_{n-1}} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} h_1(m_1-m_0)=h_0(m_2-m_1)\\ h_{n-1}(m_{n-1}-m_{n-2})=h_{n-2}(m_n-m_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} -h_1m_0+(h_1+h_0)m_1-h_0m_2=0\\ -h_{n-1}m_{n-2}+(h_{n-1}+h_{n-2})m_{n-1}-h_{n-2}m_n=0 \end{cases} \nonumber \\ \end{align}\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 &2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & -h_{n-1} & h_{n-1}+h_{n-2} & -h_{n-2} \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}在n=4时通⽤公式⾃然边界\begin{bmatrix} 1 & 0 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}固定边界\begin{bmatrix} 2 & 1 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ B-\frac{y_3-y_2}{h_2} \end{bmatrix}⾮节点边界\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & -h_2 & h_1+h_2 & -h_1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}图像插值x_{i+1}-x_i=1 \Rightarrow h_i(x)=1在n=4时，即四个点时如下所⽰p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}⾃然边界（可⽤TDMA或化简计算）\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}固定边界（只能⽤TDMA计算）\begin{bmatrix} 2 & 1 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} y_1-y_0-A\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ y_2-y_3+B \end{bmatrix}⾮节点边界（只能化简计算）\begin{bmatrix} -1 & 2 & -1 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & -1 & 2 & -1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}python代码1class BiSplineInterpolation:2 @staticmethod3 def TDMA(a, b, c, d):4 n = len(d)56 c[0] = c[0] / b[0]7 d[0] = d[0] / b[0]89for i in range(1, n):10 coef = 1.0 / (b[i] - a[i] * c[i - 1])11 c[i] = coef * c[i]12 d[i] = coef * (d[i] - a[i] * d[i - 1])1314for i in range(n - 2, -1, -1):15 d[i] = d[i] - c[i] * d[i + 1]1617return d1819 @staticmethod20 def Simplified_Natural4(y1, y2, y3, y4):21 # 四点⾃然边界化简公式22 d1 = y1 + y3 - 2 * y223 d2 = y2 + y4 - 2 * y32425 k0 = 026 k1 = (4 * d1 - d2) * 0.427 k2 = (4 * d2 - d1) * 0.428 k3 = 02930return [k0, k1, k2, k3]3132 @staticmethod33 def Simplified_Not_A_Knot4(y1, y2, y3, y4):34 # 四点⾮节点边界化简公式35 d1 = y1 + y3 - 2 * y236 d2 = y2 + y4 - 2 * y33738 k0 = 2 * d1 - d239 k1 = d140 k2 = d241 k3 = 2 * d2 - d14243return [k0, k1, k2, k3]4445 # TDMA矩阵说明46 # a0 和 c3 没有实际意义，占位⽤47 # a0 [b0 c0 00 ] [x0] [d0]48 # [a1 b1 c1 0 ] [x1] = [d1]49 # [0 a2 b2 c2] [x2] [d2]50 # [00 a3 b3] c3 [x3] [d3]5152 def SplineInterpolationNatural4(self, x, y1, y2, y3, y4):53 # ⽤TDMA计算54 # matrix_a = [0, 1, 1, 0]55 # matrix_b = [1, 4, 4, 1]56 # matrix_c = [0, 1, 1, 0]57 # matrix_d = [0, 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 0]58 # matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)5960 # 化简计算61 matrix_x = self.Simplified_Natural4(y1, y2, y3, y4)6263 a = y264 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.065 c = matrix_x[1] / 2.066 d = (matrix_x[2] - matrix_x[1]) / 6.06768 s = a + b * x + c * x * x + d * x * x * x69return s7071 def SplineInterpolationClamped4(self, x, y1, y2, y3, y4):72 # 仅有TDMA计算，⽆法化简73 A, B = 1, 17475 matrix_a = [0, 1, 1, 1]76 matrix_b = [2, 4, 4, 2]77 matrix_c = [1, 1, 1, 0]78 matrix_d = [6 * (y2 - y1 - A), 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 6 * (B - y4 + y3)]79 matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)8081 a = y282 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.083 c = matrix_x[1] / 2.084 d = (matrix_x[2] - matrix_x[1]) / 6.08586 s = a + b * x + c * x * x + d * x * x * x87return s8889 def SplineInterpolationNotAKnot4(self, x, y1, y2, y3, y4):90 # ⽆法使⽤TDMA计算91 matrix_x = self.Simplified_Not_A_Knot4(y1, y2, y3, y4)9293 a = y294 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.095 c = matrix_x[1] / 2.096 d = (matrix_x[2] - matrix_x[1]) / 6.09798 s = a + b * x + c * x * x + d * x * x * x99return s100101 def biSpline4(self, srcImg, dstH, dstW):102 dstH, dstW = int(dstH), int(dstW)103 srcH, srcW, _ = srcImg.shape104 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')105 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)106for dstY in range(dstH):107for dstX in range(dstW):108for channel in [0, 1, 2]:109 # p11 p12 p13 p14110 #111 # p21 p22 p23 p24112 # p113 # p31 p32 p33 p34114 #115 # p41 p42 p43 p44116 # 储存为 p(y, x)117 p = [dstY * srcH / dstH, dstX * srcW / dstW]118 p22 = [math.floor(p[0]), math.floor(p[1])]119 p21 = [p22[0], p22[1] - 1]120 p23 = [p22[0], p22[1] + 1]121 p24 = [p22[0], p22[1] + 2]122123 p11 = [p21[0] - 1, p21[1]]124 p12 = [p11[0], p22[1]]125 p13 = [p11[0], p23[1]]126 p14 = [p11[0], p24[1]]127128 p31 = [p21[0] + 1, p21[1]]129 p32 = [p31[0], p22[1]]130 p33 = [p31[0], p23[1]]131 p34 = [p31[0], p24[1]]132133 p41 = [p21[0] + 2, p21[1]]134 p42 = [p41[0], p22[1]]135 p43 = [p41[0], p23[1]]136 p44 = [p41[0], p24[1]]137138 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]139 r1 = self.SplineInterpolationNatural4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 140 r2 = self.SplineInterpolationNatural4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 141 r3 = self.SplineInterpolationNatural4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 142 r4 = self.SplineInterpolationNatural4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 143144 c = self.SplineInterpolationNatural4(diff_y, r1, r2, r3, r4)145146 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)。
毕业设计--一氧化碳报警器
一氧化碳报警器摘要随着工业的发展和汽车尾气的大量排放，空气遭到了严重的污染，同时在家庭生活中，煤气的不完全燃烧也将产生大量的一氧化碳气体，而一氧化碳气体的化学性质比较稳定，在室内通风条件不大好的情况下，可能引起人体中毒，产生致命后果。
因此需要对大气中的一氧化碳气体进行监测。
本文主要对一氧化碳报警器的原理进行了分析，其中重点分析了一氧化碳传感器探测信号的模式分类及识别，其原理是采用温度调制模式，基于统计学原理，采用贝叶斯公式对各模式进行识别。
对于一氧化碳传感器的信号处理方法，主要介绍了快速傅立叶变换方法(FFT)和离散小波变换方法(DWT)。
同时还对目前一氧化碳监测中存在的主要问题进行了初步的分析。
最后，根据一氧化碳传感器的原理，针对家庭中一氧化碳气体的监测，初步确定一氧化碳传报警器的设计要求，根据当前一氧化碳传感器的发展状况，选用MOTOROLA公司的MGS1100芯片进行一氧化碳报警器的设计。
同时对MGS1100传感器的结构及工作模式进行了介绍。
根据MGS1100芯片的结构及工作模式，设计了报警器电路，并对报警器的控制软件进行了初步的设计，最后完成了对一氧化碳报警器的调试，调试结果表明该报警器是可以在实际中使用。
关键词一氧化碳报警器；温度调制模式；传感器Co sannuncitorAbstractWith the developing of industry and exhaust gas expelling of cars, the atmosphere is polluting. At the same time, the incomplete inflammation of gas will produce much CO in family, and CO is a kind of gas which is stabilization in chemic character, will bring poisoning to person when the house is stuffiness. So it is need to monitor the CO in house.This paper analyzes the principle of CO sannunciator, focus on the analysis of the model classification and recognition of CO sensor, the operational principle of CO sensor is based on statistics and temperature moderate model. The signal process method of CO sensor include FFT (Fast Fourier Transform) and DWT (Discreet Wave Transform). At the end, introduce the main problem of CO monitor.At the end, based on the operational principle of CO sensor, according to the requirement of CO monitor in house, decide the parameters of CO sannunciator. According to the development of CO sensor, select MGS1100 of MOTOLORA to design the CO sannunciator. At the same time, introduce the structure and operational principle of MGS1100, and design the alarm circuit and control soft. At the end, finish the debugging of CO sannunciator, the result show the CO sannunciator can be used in realism.Keywords CO sannunciator;temperature moderate model; senso目录摘要 (I)Abstract (II)第1章绪论 (1)1.1 课题背景 (1)1.2研究的目的和意义 (1)1.3一氧化碳报警器的发展状况 (2)第2章一氧化碳传感器信号处理技术 (5)2.1气体传感器阵列信号处理技术 (5)2.2信号预处理 (5)2.3模式分类，识别和量化 (7)2.4气体传感器温度调制及信号处理技术 (10)2.5温度调制模式 (10)2.6信号处理方法 (12)2.7 存在的问题 (12)2.8 本章小结 (13)第3章一氧化碳报警器的设计 (14)3.1器件的结构原理 (14)3.2器件的工作模式 (15)3.3一氧化碳报警器的设计 (16)3.4电路原理 (17)3.5控制过程 (18)3.5.1对于温度变化的特性处理(假设CO的浓度为60ppm)。
等距曲线也称为平行或位差曲线,它是基曲线沿法向距离为d的点的轨迹,
4.2.1 曲线上动点 Qj (u,α ) 到型值点连线的距离计算...........................32 4.2.2 曲线上动点 Qj (u,α ) 到型值点连线的距离控制...........................35
In chapter one, we firstly review the phylogeny of the offset curves, then summarize PH curves and OR curves concerning the inherited geometric structure of offset curves.
关键词: 等距曲线 逼近算法 奇异混合 离散公式 误差控制
I
Abstract
ABSTRACT
Offset curves, also called parallel or potential difference curves, are defined as locus of the points which are at constant distance d along the normal from the base curves. It has been a hot topic on computer aided geometric design in recent ten years.
Keyword: offset curves; approximation algorithm; singular blending; subdivision formula; error control
拉格朗日插值和牛顿插值多项式的C程序算法毕业论文
本科生毕业论文题目: 拉格朗日插值和牛顿插值多项式的C程序算法原创性声明本人郑重声明: 所提交的学位论文是本人在导师指导下, 独立进行研究取得的成果. 除文中已经注明引用的内容外, 论文中不含其他人已经发表或撰写过的研究成果, 也不包含为获得**大学或其他教育机构的学位证书而使用过的材料. 对本文的研究做出重要贡献的个人和集体, 均已在文中以明确方式标明. 本人承担本声明的相应责任.学位论文作者签名: 日期指导教师签名: 日期目录拉格朗日插值多项式的C程序算法 (1)1引言 (1)1.1插值问题的提出 (1)1.2插值法 (2)1.3插值法思想 (2)2拉格朗日插值法 (3)2.1拉格朗日插值法的由来 (3)2.2n次插值基函数 (4)2.3拉格朗日插值多项式 (4)3牛顿插值法 (5)3.1均差： (5)3.2牛顿插值多项式： (6)4C程序设计 (7)4.1算法设计： (7)4.2程序源码编写 (8)5程序检测 (12)5.1对拉格朗日插值的检测 (12)5.2对牛顿插值的检测 (13)总结 (15)参考文献 (16)致谢 (17)摘要本论文着重研究了用C语言编写程序计算拉格朗日插值和牛顿插值的方法。
在前人已有的研究成果的基础上，首先介绍了拉格朗日插值和牛顿插值的思想和方法，通过添加可以循环计算功能和输入非法数值时的纠错功能，改进了已有文献的方法，对其进行了推广，使之更加的合理和完美，并且通过实际的例子进行了具体的验证。
最后，总结了一下本论文的主要研究成果和应用前景。
关键词：拉格朗日插值，牛顿插值，C算法，精确解AbstractThis article discuss the method to calculate Lagrange interpolation and Newton interpolation with C program. Base on the results of predecessors' research, firstly, this article introduces the thoughts and methods of Lagrange interpolation and Newton interpolation. Improving the old method by adding functions which can repeatedly computing interpolation and correct illegal data. Then spreading it and making it more reasonable and perfect, checking it with some examples. Finally, summing up the main results of this article and application prospect.Key words：Lagrange interpolation; Newton interpolation ; C program;拉格朗日插值多项式的C 程序算法1引言插值法是一种古老的数学研究方法，他的产生来自与社会的生产实践活动。
空间插值分析算法综述
2 专论•综述 Special Issue
森多边形的每个顶点都是每个三角形的外接圆心[1]. 泰 森多边形也被称为 Voronoi 图或 dirichlet 图[3].
A
B
D
C
图 1 泰森多边形示意图
泰森多边形利用离散观测点的值对该点所在的区
域进行赋值, 得到的结果往往是数值的变化只发生在
多边形的边界上, 而多边形内部的数值则是均匀、同 质的[9]. 其数学表达式为:
Ve = Vi
(1)
其中, Ve 表示待插值点的距离, Vi 表示 i 示点的离散观 测值. i 点必须满足如下条件:
dei = min(de1, de2, · · · , den)
(2)
其中, dij 表示点 i(xi,yi) 与点 j(xj,yj) 间的欧几里德距离[8]. 泰森多边形的关键是将离散观测点合理地连接到
引用格式: 李海涛,邵泽东.空间插值分析算法综述.计算机系统应用,2019,28(7):1–8. /1003-3254/6988.html
一种基于Markov模型和协同克里格的插值方法
一种基于Markov模型和协同克里格的插值方法Markov模型是一种采用概率论的建模方法，通过对数据中出现的模式进行分析，从而创建一个可以预测未来事件的模型。
协同克里格是一种空间统计方法，用于空间数据的建模和模拟。
这两种方法的结合，可以创建一种有效的插值方法。
基于Markov模型和协同克里格的插值方法，可以通过对空间数据进行分析和建模，找到模式和趋势，从而预测未知位置的值。
这个方法的基本思路是在空间数据上进行Markov模型的建模，然后使用协同克里格方法来实现插值。
首先，我们需要将原始数据进行预处理，包括去掉异常值、插值缺失数据等。
然后，使用Markov模型对数据进行建模，提取数据中的空间相关性和趋势信息。
Markov模型是一个基于历史状态的随机模型，它使用当前状态来预测未来状态。
通过对数据进行分析，我们可以得到当前状态和下一个状态之间的转移概率，利用这些概率，我们可以预测未来的状态。
在这种插值方法中，我们利用Markov模型预测未知位置的值。
然后，使用协同克里格方法来优化插值。
协同克里格方法是一种多变量插值方法，它可以在不同空间尺度上进行空间插值，通过对空间变异性进行建模，将多个空间变量进行插值和合并。
使用协同克里格方法可以提高插值精度和可靠性。
最后，我们将两种方法结合，基于Markov模型进行预测，然后使用协同克里格方法进行插值。
这个方法可以应用于许多地理学领域，例如环境科学、社会经济研究、水资源管理等。
总之，基于Markov模型和协同克里格的插值方法可以准确地预测未知位置的值，同时还可以提高插值精度和可靠性。
这种方法在各种空间数据分析领域中都有广泛的应用。
mike模型插值方法
mike模型插值方法
Mike模型插值方法是一种地质建模方法，通常用于地质数据的插值和建模。
这种方法是基于地质学家Michael Hohn的研究成果而得名的。
Mike模型插值方法主要用于处理地质数据，例如岩性、孔隙度、渗透率等地质参数的空间变化。
这种方法的核心思想是利用已知的地质数据点来推断未知区域的地质特征，从而实现对整个区域的地质特征进行建模和预测。
Mike模型插值方法通常包括以下几个步骤：
1. 数据采集和准备，首先需要收集并准备地质数据，包括采样点的坐标、地质参数数值等信息。
2. 插值方法选择，根据数据的特点和空间分布情况，选择合适的插值方法，常见的插值方法包括克里金插值、反距离加权插值、样条插值等。
3. 插值参数设定，根据实际情况，设置插值方法所需的参数，如克里金插值中的变程、插值阶数等。
4. 插值计算，利用选定的插值方法和参数对地质数据进行插值计算，推断未知区域的地质特征。
5. 结果分析和验证，对插值结果进行分析和验证，评估插值模型的准确性和可靠性，可以通过交叉验证等方法进行验证。
Mike模型插值方法在地质建模领域得到了广泛的应用，可以帮助地质工作者更好地理解地下地质结构和特征，为资源勘探和开发提供重要的支持。
该方法的优点在于能够利用有限的地质数据对整个区域进行建模，但也需要注意插值结果的可靠性和局限性，尤其是在数据稀疏或不均匀分布的情况下，需要谨慎处理和分析插值结果。
总的来说，Mike模型插值方法是一种重要的地质建模方法，对于理解地质特征和预测地下地质结构具有重要意义。
狄利克雷过程混合模型的马尔可夫链蒙特卡洛方法 -回复
狄利克雷过程混合模型的马尔可夫链蒙特卡洛方法-回复狄利克雷过程混合模型(Dirichlet Process Mixture Model, DPMM)是概率机器学习领域中常用的非参数贝叶斯模型，用于聚类、密度估计等任务。
而马尔可夫链蒙特卡洛方法(Markov Chain Monte Carlo, MCMC)是一类重要的随机模拟方法，广泛应用于贝叶斯统计推断中。
本文将详细介绍狄利克雷过程混合模型的马尔可夫链蒙特卡洛方法。
第一部分：狄利克雷过程混合模型简介首先，我们来回顾一下狄利克雷过程(Dirichlet Process, DP)的定义。
对于一个随机概率分布G，如果对于任意的有限划分(A1, ..., An)且n趋于无穷大时，G(A1), ..., G(An)的联合分布满足以下特性：(G(A1), ..., G(An))~Dir(αG(A1), ..., αG(An))其中，α>0是一个超参数，Dir(·)表示狄利克雷分布。
我们称G是一个狄利克雷过程。
对于狄利克雷过程混合模型(DPMM)，假设我们有一组观测数据X={x1, ..., xn}，我们的目标是通过数据来估计未知的总体分布。
DPMM的核心思想是，将总体分布看作是一系列由狄利克雷过程生成的混合分布的随机样本。
具体来说，对于每个样本xi，我们选择它属于哪一个混合分布，由具有无穷多个成分的狄利克雷过程决定。
每个混合分布都有一个成分参数θ，而狄利克雷过程决定了成分参数的概率分布。
第二部分：马尔可夫链蒙特卡洛方法简介接下来，我们来简要介绍一下马尔可夫链蒙特卡洛方法(MCMC)。
MCMC 是一类随机模拟方法，用于从复杂的分布中采样。
对于贝叶斯统计推断，MCMC方法可以用于近似计算后验分布的积分，从而实现参数估计和模型选择。
MCMC的核心思想是通过构建一个马尔可夫链，使得该链的平稳分布正好是目标分布。
在MCMC算法中，马尔可夫链的状态转移是按照特定的转移概率进行的，从而保证经过足够的迭代后，我们可以在马尔可夫链的平稳状态下得到目标分布的采样。
clips插值法实现代码
clips插值法实现代码【原创版】目录1.插值法概述2.CLIPS 插值法介绍3.CLIPS 插值法的实现代码4.CLIPS 插值法在实际应用中的优势5.总结正文1.插值法概述插值法是一种通过已知的数据点来预测未知数据点的方法。
在数学和计算机科学中，插值法被广泛应用于函数逼近、数据分析以及图像处理等领域。
插值法的基本原理是根据已知的数据点，构建一个函数，使得这个函数在未知数据点处能够尽可能地接近实际值。
2.CLIPS 插值法介绍CLIPS 插值法，全称为 Cubic Locally Interpolating Polynomial Smoothing，即三次局部多项式插值平滑法，是一种基于局部多项式插值的三维空间数据插值方法。
CLIPS 插值法通过对空间数据进行局部多项式插值，能够在保持数据原貌的同时，实现数据的平滑处理。
这种方法在计算机图形学、地形建模和虚拟现实等领域有着广泛的应用。
3.CLIPS 插值法的实现代码CLIPS 插值法的实现主要分为以下几个步骤：（1）数据预处理：将输入的数据点进行预处理，包括去除异常值、计算邻域平均值等。
（2）构建局部多项式：以每个数据点为基点，构建一个局部多项式。
局部多项式的阶数通常为三次，即三次多项式插值。
（3）插值计算：利用构建的局部多项式，对输入的数据点进行插值计算，得到预测的数据点。
（4）平滑处理：对预测的数据点进行平滑处理，以实现数据的平滑过渡。
下面是一个简单的 Python 代码示例：```pythonimport numpy as npdef clip_interpolation(x, y, z, order=3):# 数据预处理x_mean = np.mean(x)y_mean = np.mean(y)z_mean = np.mean(z)# 构建局部多项式x_poly = np.poly1d(np.polyfit(x, z, order))y_poly = np.poly1d(np.polyfit(y, z, order))# 插值计算x_interp = np.interp(x, x_mean, x_poly(x_mean))y_interp = np.interp(y, y_mean, y_poly(y_mean))# 平滑处理z_interp = np.interp(x, y, x_interp * y_interp)return z_interp# 示例数据x = np.array([0, 1, 2, 3, 4])y = np.array([0, 1, 2, 3, 4])z = np.array([5, 7, 9, 11, 13])# 插值结果z_interp = clip_interpolation(x, y, z)```4.CLIPS 插值法在实际应用中的优势CLIPS 插值法在实际应用中具有以下优势：（1）保持数据原貌：CLIPS 插值法能够较好地保持原始数据的局部结构和特征。
基于Markov模型与轨迹相似度的移动对象位置预测算法
基于Markov模型与轨迹相似度的移动对象位置预测算法摘要：针对低阶Markov模型预测精度较差，以及多阶Markov 模型预测稀疏率高的问题，提出一种基于Markov模型与轨迹相似度（MMTS）的移动对象位置预测算法。
该方法借鉴了Markov模型思想对移动对象的历史轨迹进行建模，并将轨迹相似度作为位置预测的重要因素，以Markov预测模型的预测结果集作为预测候选集，结合相似度因素得出最终预测结果。
实验结果表明，与k阶Markov模型相比，该方法的预测性能不会随着训练样本大小及阶数k的变化受到很大的影响，并且在大幅降低k阶Markov模型预测稀疏率的同时将预测精度平均提高了8%以上。
所提方法不仅解决了k 阶Markov模型的预测稀疏率高及预测精度不足的问题；同时提高了预测的稳定性。
关键词：轨迹相似度；位置预测；移动对象；马尔可夫模型；稀疏性0 引言随着移动便携设备的普及和定位技术的发展，基于位置的服务得到了广泛的应用。
当前基于位置的服务主要集中于定位查询及位置共享，这些应用主要关注于给用户提供有关当前位置的服务。
为了使服务更具前瞻性，移动对象位置预测技术越来越受到国内外研究学者的关注[1]。
移动对象位置预测具有高度的研究价值和广阔的应用领域，如果能够实时动态地预测出用户将来的位置，可以提供更加灵活的预报服务，例如：推荐服务[2]、提醒服务、广告推送、智能化交通服务等。
一般来说，当前存在的移动对象位置预测方法可以分为基于欧氏空间和基于受限路网两类位置预测方法。
由于在实际情形中，基本所有的移动对象的运动都受限于实际路网条件，所以对于具有空间运动随意性的基于欧氏空间的位置预测方法在实际应用中具有很大的局限性。
而在基于受限路网的位置预测研究中，应用较广泛的预测方法又可分为基于序列分析模式的预测方法[3-5]和基于马尔可夫模型的统计模型预测方法[6-10]。
文献[6]中使用混合Markov模型用于移动路径的预测，减少了2阶Markov预测模型的空间代价，但其并没有有效地提高预测精度，而且对历史信息考虑不足；文献[7-8]使用简单的Markov模型来进行位置的预测，该模型仅考虑当前时刻位置对将来位置的影响，导致预测精度过低；文献[9]中提出了一种可变阶的局部匹配的Markov预测模型，减小了预测的稀疏率，但在其预测过程会存在大量退化为简单的Markov模型的情况；此外，SPM（Sampled PatternMatching）算法[10]、PPM（Prediction by Partial Matching）算法[11]等，也都是在Markov模型的基础上进行轨迹预测的研究。
结合Markov模型和COSGSIM的插值方法
结合Markov模型和COSGSIM的插值方法
汪彦龙;刘金华;张挺
【期刊名称】《计算机辅助设计与图形学学报》
【年(卷),期】2010(022)010
【摘要】针对协同序贯高斯模拟在确定插值公式中的权值时出现的交叉矩阵不稳定问题,提出一种基于Markov模型和协同序贯高斯模拟的插值方法.利用Markov 模型的屏蔽效应,假设待模拟硬数据可以屏蔽其位置以外的硬数据对待模拟位置处软数据的影响,对协同序贯高斯模拟进行逼近时只保留与待模拟硬数据同位置的软数据,而不考虑其他位置的软数据,从而解决了交叉矩阵不稳定的问题.实验结果证明,该方法插值误差小,模拟效果好.
【总页数】6页(P1715-1720)
【作者】汪彦龙;刘金华;张挺
【作者单位】中国科学技术大学电子工程与信息科学系,合肥,230027;浙江传媒学院电子信息学院,杭州,310018;浙江传媒学院电子信息学院,杭州,310018;C4ISR国防科技重点实验室,南京,210007
【正文语种】中文
【中图分类】TP391.41
【相关文献】
1.一种基于改进Markov模型的插值方法 [J], 杜奕;张挺;卢德唐;李道伦
2.基于CLUE-S改进模型和Markov模型的区域土地利用变化多情景模拟 [J], 刘
玉湖;张明;陈大凯
3.基于CLUE-S模型和Markov模型的土地利用变化模拟预测\r——以德阳市为例[J], 邓元杰;潘洪义;陈丹;蒋贵国;曲比伟石;孙嘉璐
4.指数平滑结合条件Markov预测模型研究\r——基于新疆货运量预测 [J], 吕冬梅
5.Markov随机游走和高斯混合模型相结合的运动目标检测算法 [J], 李菊;李克清;苏勇刚
因版权原因，仅展示原文概要，查看原文内容请购买。
cokriging插值法
cokriging插值法是一种基于kriging插值方法的一种扩展，可以用于处理多个相关的目标变量。
下面是一个简单的介绍，包括多个O 的目标、O下边具体的KR1、KR2，以及子指标。
1. 多个O的目标cokriging插值法的目标是估计一个空间域内的多个目标变量的值。
这些目标变量可以是任何需要插值的变量，如温度、湿度、海拔等等。
在cokriging插值中，这些目标变量被统称为“目标”（objective）。
2. KR1和KR2在cokriging插值法中，需要用到两个或更多的kriging模型。
这些模型被称为“kriging权重函数”（kriging weights）或简称“KR”。
KR1和KR2是两个常用的kriging权重函数，它们分别表示第一个和第二个目标变量的kriging权重。
KR1和KR2都是基于空间域内的样本数据和它们之间的空间自相关结构来计算的。
KR1和KR2的取值范围是0到1之间，表示每个样本点在估计目标变量时的权重大小。
3. 子指标除了KR1和KR2之外，cokriging插值法还需要用到一些子指标（sub-objective）。
这些子指标是除了KR1和KR2之外的其他需要考虑的因素，如地形、气候、土壤类型等等。
子指标可以用来调节kriging 权重函数的影响程度，以更好地拟合实际情况。
在cokriging插值法中，每个子指标都需要根据实际情况进行定义和计算。
不同的子指标可能会有不同的计算方法和影响程度，需要根据具体情况进行调整和优化。
总之，cokriging插值法是一种基于kriging插值方法的扩展，可以用于处理多个相关的目标变量。
通过定义多个O的目标、具体的KR1、KR2和子指标，可以更好地拟合实际情况，提高插值精度。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
ｐ ｏ ｉ ｅ ａ ｋｏ ｒ ｖ ｄ ｄ ｂｙ ａ Ｍ ｒ ｖ ｍｏｄ ｌ ｈｅ ｈａ ｄ ｄ ｔ ｅ ，ｔ ｒ ａ ｕｍ ｏ ｂ ｉ ｔ ｅ ｓｍｕｌ ｔ ｄ ｓ ｒ ｅ ｈｅ ｉ ｆｕｅ ｅ ｏ ｎｙ ｏｔ ｒ ｈ ｒ ａ ｅ ｃ ｅ ｎｓ ｔ ｎ ｌ ｎｃ ｆ ａ ｈｅ ａ ｄ ｄ ｔ ｎ ｔ ｏ ｔｄ ｔ ｍ ｔｔ ａ ｏ ａ ｉ ｎ．Ｔｈｅ ｅ ｏｒ ａ ａ ｏ ｈｅ ｓ ｆ ａ ｕ ａ ｈｅ ｓ ｍｅ １ ｃ ｔｏ ｒ ｆ ｅ，ｔ ｎｓ ａ ｉｉｙ ｏ ｒ ｓ ｔ ｉ ａ ｅ ｄ ｓ ｏ ｖ ｄ ｈｅｉ ｔ ｂ ｌｔ ｆｃ ｏ ｓｍａ ｒｘ ｃ ｎ ｂ ｉ ｓ ｌ ｅ ｂ ｒ ｓ ｒ ｉ ｎｌ ｈｅ ｓ ｔ ｄａ ｕ ｔ ｔ ｈａ ｈｅ ｓ ｍ ｅ ｏ ａ ｉ ｎ ｗｉｈ ｔ ｈ ｒ ａ ｕ ｔ ｂ ｓｍｕｌ ｔ ｄ， ｙ ｐ ｅ ｅ ｖ ｎｇ ｏ ｙ ｔ ｏｆ ｔ ｍ ｈａ ｓ ｔ ａ １ ｃ ｔｏ ｔ ｈｅ ａ ｄ ｄ ｔ ｍ ｏ ｅ ｉ ａｅ ｗｈｉｅ ｉ ｎｏ ｉ ｇ ｔ ｍ ｐ ｃ ｙ ｏ ｈｅ ｏ ｔｄａ ａ ｗｈｅ ｐｐ ｏｘ ｍａ ｉ ｇ ｔ ｌ ｇ ｒｎ ｈｅ ｉ ａ ｔｂ ｔ ｒ ｓ ｆ ｔ ｎ ａ ｒ ｉ ｔｎ ｈｅＣＯＳＧＳ Ｍ ．Ｅｘｐ ｒｍｅ ｔ ｌｒ ｓ ｌｓ Ｉ ｅ ｉ ｎ ａ ｅ ｕ ｔ ｓ ｏ ｔ ａ ｈｅ ｄ ｖ ａ ｉ Ｓ ｓ ｌ ａ ｈ ｉ ｈ ｗ ｈ ｔｔ ｅ ｉ ｔｏｎ ｉ ｍａ ｌ ｎｄ ｔ ｅ ｓｍｕｌ ｔ ｄ ｒ ｓ ｔｉ ｏ ａ ｅ ｅ ｕｌ Ｓ ｇｏ ｄ．
合肥
２０２） ３ ０ ７
杭 州 ３ ０ １ ） １ ０ ８ 南 京 ２ ０ ０ ） １０ ７
摘
要 ： 对 协 同序 贯 高 斯 模 拟 在 确 定 插 值 公 式 中 的权 值 时 出现 的 交 叉 矩 阵 不 稳 定 问题 ， 出一 种 基 于 Ｍａｋ ｖ模 针 提 ｒｏ
型 和 协 同 序 贯 高 斯模 拟 的 插 值 方 法 ． 用 Ｍａｋ ｖ 型 的 屏 蔽 效 应 ， 设 待 模 拟 硬 数 据 可 以 屏 蔽 其 位 置 以 外 的硬 数 利 ｒｏ 模 假 据 对待 模 拟位 置处 软 数据 的影 响 ， 协 同序 贯 高 斯 模 拟 进 行 逼 近 时 只保 留 与 待 模 拟 硬 数 据 同位 置 的 软 数 据 ， 不 考 对 而 虑 其他 位 置 的 软数 据 ， 而 解 决 了交 叉 矩 阵 不 稳 定 的问 题 ． 从 实验 结果 证 明 ， 方 法 插 值 误 差 小 ， 拟效 果好 ． 该 模
结 合 Ｍａ ｋ ｖ模型 和 Ｃ ＧＳＭ 的 插 值 方 法 ｒｏ ＯＳ Ｉ
汪 ， 彦龙 刘金华∞ 张 ， 挺
”（ 国科 学 技 术 大 学 电 子 工程 与信 息 科 学 系 中 （ 江 传 媒 学 院 电 子 信息 学 院 浙 ”（ Ｓ 国防 科 技 重 点 实验 室 ｃ ＩＲ
（ ｎｎ ｍａｌｕ ｔ． ｄ ． ｎ １ ｉ＠ ｅ ｉ ｓｃ ｅ ｕ ｃ ） ．
ｗ ｅｇ ｓｏｆ ｔ ｅ ｉ ｅ ｐｏ ａ ｉ ｑ ａ ｉ ｎ ｉ ｉ ｈｔ ｈ ｎｔ ｒ ｌ ｔｏｎ ｅ ｕ ｔｏ ｎ ＣＯ Ｓ ＳＩ ． Ａ ｃ ｏ ｄ ｎ ｔ ｔ ｈ ｏｔ ｓｓ ｏ ｃ ｅ ｎ ｎｇ ｅ ｆ ｃ Ｇ Ｍ ｃ ｒ ｉ ｇ ｏ ｈｅ ｙｐ ｈｅ ｉ ｆ ｓ ｒ ｅ ｉ ｆｅ ｔ
。（ ｔｏ ａ ｙ Ｌ ｂ ｒ ｔｒ ｆ Ｓ ｉｎ ｅａ ｄ Ｔｅｈ ｏ ｏ ｙ ｏ Ｓ ’ Ｎａ ｉｎ ｌ Ｋｅ ａ ｏａ ｏ ｙ ｏ ｃｅ ｃ ｎ ｃ ｎ ｌｇ ｎ Ｃ Ｉ Ｒ。Ｎａ ｊｎ ２ ０ ０ ） ｎ ｉｇ １ ０ ７
Ａｂ ｔ ａ ｔ ｓｒ ｃ ：Ａｎ ｎｔ ｒ ａ ｉ ｎ ｉ ｅ ｐｏｌｔｏ ｍｅ ｈ ｅ ｌ ｉｉ ｂｏ ｈ ｈｅ ｔ ｏｄ ｘｐ ｏ ｔｎｇ ｔ ｔ Ｍ ａ ｋｏ ｍｏ ｌ ｎｄ ｅ ｕｅ ｔａ Ｇａ ｓ ｉ ｎ Ｏ ｒ ｖ ｄｅ ａ ｓ ｑ ｎ ｉ ｌ ｕ ｓａ Ｃ — ｓｍｕｌ ｔｏ （ ｉ ａ ｉ ｎ ＣＯＳ ＧＳＩ ） ｉ ｐｒ ｐ ｅ ｏ ｏ ｖｅ ｈｅ ｎｓ ａ ｉｉｙ ｆ ｃ ｏ ｓ ａ ｒｘ Ｍ ｓ ｏ ｏｓ ｄ ｔ ｓ ｌ ｔ ｉ ｔ ｂ ｌｔ ｏ ｒ ｓ ｍ ｔ ｉ ｗｈｅ ｄｅ ｅ ｍｉ ｉ ｔ ｎ ｔ ｒ ｎ ｎｇ ｈｅ
Ｗ ａ ｇＹａ ｌｎ ，Ｌ ｕ Ｊｎ ｕ ，ａ ｄ Ｚ ａ ｇ Ｔｉｇ ｎ ｎｏ ｇ’ ｉ ｉｈ ａ ｎ ｈ ｎ ｎ 。
（Ｄ ｐ ｒｍｅ ｔ ｆ Ｅ ｅｔｏ ｉ Ｅｎ ｉ ｅｒｎ ｎ ｎ ｏ ｍａ ｉｎＳ ｉ ｃ ，Ｕ ｉｅ ｓｔ ｆ Ｓ ｉｎ ｅａ ｄ Ｔ ｃ ｎ ｌｇ ｈ ｎ ，Ｈ８ ｅ ２ ０ ２ ） ｅ ａ ｔ ｎ ｌｃｒ ｎｃ ｇ ｎ ｅ ｉｇ ａ ｄ Ｉ ｆ ｒ ｔ ｃｅ ｅ ｎ ｖ ｒｉ ｏ ｃｅ ｃ ｎ ｅｈ ｏｏ ｙ ｏ Ｃ ｉ ａ ｏ ｏ ｎ ｙ ｆ ｆ ｉ ３ ０ ７ （ ｐ ｒｍｅ ｔ ｅ ｔｏ ｉ ａ ｄ Ｉ ｆ ｒ ｔｏ ｎ ｉｅ ｒｎ Ｄｅ ａ ｔ ｎ Ｅｌ ｒ ｎ ｃ ｎ ｎ ｏ ｍａ ｉｎ Ｅ ｇ ｎ ｅ ｉ ｇ，Ｚ ｅｉ ｎ ｎ ｔ ｕ ｅ ｏ ｆ ｏ ｃ ｈ ｊ ａ ｇ Ｉ ｓｉ ｔ Ｃ ｍｍｕ ｉａ ｉｎａ ｄ Ｍ ｅ ｉ ｔ ｏ ｆ ｎ ｃ ｔ ｎ ｄ ａ，Ｈａ ｇ ｈ ｕ ３ ０ １ ） ｏ ｎ ｚ ｏ １ ０ ８
关 键 词 ：插 值 ； 同序 贯 高斯 模 拟 ； ｒｏ 协 Ｍａｋ ｖ模 型 ； 简单 克里 格 ； 局 协 同克 里 格 全
中 图法 分 类 号 ：Ｔ ３ １ ４ Ｐ ９．１
Ａｎ Ｉ ｅ ｐ ｌ ｔｏ ｅ ｈｏ ｉ ｇ ｔ ｅ Ｍ ａ ｋ ｖ Ｍ ｏ ｌａ ｄ ＣＯＳ ｎｔ ｒ ｏ ａ ｉ ｎ Ｍ ｔ ｄ Ｕｓ ｎ ｈ ｒｏ ｄｅ ｎ ＧＳ Ｍ Ｉ
第２ ２卷 第 １ 期 Ｏ ２１ Ｏ Ｏ年 １ Ｏ月
计算 机辅 助设计 与 图形 学学报
ｍｐ ｔｒＧｒｐ ｉｓ Ｊ ｕ ｎ ｌ ｆ ｍｐ ｔｒＡｉｅ ｓｇ ＆ Ｃｏ ｕ ｅ ａ ｈｃ ｏ ｒ ａ ｏ Ｃｏ ｕ ｅ — ｄ ｄＤｅｉｎ
Ｖｏ ． ２ ＮＯ １ １２ ． ０
０ｃ ． ２ ０ ｔ Ｏ１