下载文档原格式
三维重心插值计算公式在计算机图形学和计算机辅助设计领域中，三维重心插值是一种常用的插值方法，用于在三维空间中对数据进行插值和补偿。
三维重心插值是基于三角形网格的，通过计算三角形的重心坐标来实现插值。
在本文中，我们将介绍三维重心插值的计算公式，并讨论其在实际应用中的重要性和优势。
三维重心插值的计算公式如下：假设有一个三角形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）领域中，三维模型的变形是一项重要的技术，可以用于模拟材料的变形和形状的调整。
三维重心插值可以通过对三角形网格进行插值，实现对三维模型的形状变形，为工程设计和仿真分析提供了有效的工具和方法。
总之，三维重心插值作为一种常用的插值方法，在计算机图形学和计算机辅助设计领域具有重要的应用价值。
㊀２０２３年８月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.。
基于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算法中起到一定的作用。
常见插值算法--拉格朗⽇插值、三次卷积插值、三次样条插值、兰克索斯插值写在前⾯本⽂简单介绍了⼏种常见的插值算法并附带了相应的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)。
本科生毕业论文题目: 拉格朗日插值和牛顿插值多项式的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引言插值法是一种古老的数学研究方法，他的产生来自与社会的生产实践活动。
