多项式模型

多项式模型（精选5篇）

多项式模型 篇1

在过去的20年, 格子波尔兹曼方法已经成为计算流体力学中一个有效数值工具。跟传统的计算流体力学方法相比, 格子波尔兹曼方法容易实现, 可高度并行化, 并能解决复杂边界问题。更重要的是, 由于其动力学特性, 格子玻尔兹曼方法可以更直接地对复杂的流体现象建模。

格子波尔兹曼方法发展于研究纳维斯托克斯方程的格子气自动机方法[1], 作为一个无噪音可替代LGA的技术, 布尔变量被单粒子概率分布函数取代[2]。在小马赫数和小努森数时, 通过假设分布函数接近于依赖气体密度和速度本地值的局部平衡态来线性化碰撞项[3]。若干个小组通过将BGK碰撞项[4]引入格子波尔兹曼方程得到了标准的格子BGK方程[5,6]。麦克斯韦平衡态分布[5]被建议用来取代Fermi-Dirac分布[6]作为局部平衡态分布函数。在标准单松弛模型中, 可压缩效应可能会产生比较严重的数值错误。一般来说, 不可压缩可以通过使密度恒定为一个常数得到, 而在波尔兹曼方程中保持一个定值密度是比较困难的。理论上来说格子波尔兹曼方程更适合于模拟可压缩纳维斯托克斯方程, 因为在格子玻尔兹曼模拟中密度的空间变化不为零。He和Luo通过忽略可压缩效应项中马赫数的高阶部分来减弱可压缩效应[7]。

文献证明LBE可以通过对连续波尔兹曼方程使用小马赫数泰勒展开[8]或赫米特多项式展开[9]得到, 后者基于流体的动力学表达并具有完整的理论框架。Grad使用赫米特多项式提出了13阶矩来描述流体系统。受Grad工作的启示, Shan将LBGK方程在速度空间离散以便把格子波尔兹曼方法同Grad 13阶矩系统进行关联。基于这个方法, 能够得到高于纳维斯托克斯方程的波尔兹曼方程高阶近似。但是Shan的模型依然是基于BGK碰撞项的单松弛方法。

现通过泊肃叶流动和顶盖驱动方腔流两个算例比较和评价了He-Luo方法和Shan方法的性能。

1 模型和方法

1.1 格子模型

格子波尔兹曼方程使用二维九速对称格子 (D2Q9) , 在这个模型中, 空间被离散为一个方形格子, 有9个离散速度。

这里c=δx/δt, 其中δx为格子尺度, δt为格子速度。一般来说都设置为1。格子模型如图1所示。

1.2 He-Luo方法

为了弱化标准格子波尔兹曼方法中的可压缩效应, He-Luo提出了一个不可压缩格子波尔兹曼模型。

将流体的密度写为ρ=ρ0+δρ, 其中ρ0为密度中的常数部分, 而δρ则为密度中的可变部分, 将该表达式代入麦克斯韦平衡态分布函数并忽略掉含马赫数和δρ的高阶项, 就得到新的平衡态分布函数

进一步地将压力作为独立变量, 并引入一个局部压力分布函数pα=cs2fα, 这里对D2Q9模型来说是声速, 演化方程变为

其中

其中p=cs2ρ, p0=cs2ρ0。权重

压力和速度的定义为

1.3 Shan单松弛方法

与Taylor展开不同, Shan使用Hermite多项式对波尔兹曼方程展开得到LBGK方程[10]。

N阶Hermite多项式为

其权重函数为

前三个Hermite多项式为

粒子概率密度函数f在速度空间ξ使用Hermite正交多项式展开得到

其中系数a (n) 又可表达为

跟标准SRT模型类似, 宏观量可以通过概率密度函数的积分获得。

平衡态分布函数可以表达为

该表达式在温度θ=1时就变成了标准SRT模型的平衡态函数。

1.4 边界处理方法

边界处理采用Guo提出的非平衡外推方法[11]。

式中fα (eq) 是平衡态部分, O是边界节点而B是流体中临近边界的节点

这个表达式具有二阶精度。

2 数值结果和讨论

2.1 泊肃叶流动

管道中或平行平板之间的一种重要流动是泊肃叶流动。壁面速度为0, 流动区域中间的速度达到最大, 在区域0≤x≤2和0≤x≤1内一个恒定的压力梯度驱动稳态平板泊肃叶流动。速度分布的分析解为

其中h=H/2, H是两个平板之间的距离。μ是分子黏性, Δp是进口与出口之间的恒定压力降。流动的雷诺数定义为Re=uaveL/υ, 而平均速度是最大速度的2/3。在测试算例中动力黏性υ=0.1, 最大速度umax=0.01, 密度ρ=1, 长度L=200, 高度H=100, 压力降Δp通过最大速度公式给出。计算出来的雷诺数为3.33。在模拟中, 稳态收敛判据为

这里ε=1.0×10-9。该公式对整个系统求和。

图2给出了速度轮廓线, He-Luo模型和Shan模型的计算结果与分析结果都吻合的较好。

为了更细致地比较, 相对总体误差和变量相对平均值的波动率被计算出来, 变量的绝对误差 (以速度为例)

其中是分析解, 绝对值避免了在比较中出现负值。相对总体误差定义为

其中的求和是对整个系统进行的。变量相对平均值的波动率定义为 (以密度为例)

是平均密度。相对平均密度的波动率是一个代表可压缩性的量。

在表1中, Tu是模拟使用的时间, 单位是min。He-Luo模型计算出的速度u的最大绝对误差, 最小绝对误差和相对总体误差都比Shan模型的要小, 其量级为10-4。而计算出速度v的最大和最小绝对误差也都比Shan模型的要小, 其量级为10-15。HeLuo模型是恒温模型, 而Shan模型是变温模型, 统计了Shan模型中温度的绝对误差, 相对总体误差。Shan模型中温度的误差量级为10-5, 波动量级为10-4, 相对等温模型都是可以忽略的。He-Luo模型和Shan模型在计算时间上花费时间是相同的。两个模型的可压缩效应是相同的。

2.2 顶盖驱动方腔流

尽管没有分析解, 方腔 (0, 1) (0, 1) 中不可压缩黏性稳态流动多年来被用作测试新数值格式和算法的测试算例。封闭方形空间内的不可压缩流体在顶盖均匀速度平移驱动下流动。研究该流动的文献有很多, 丰富的涡现象取决于雷诺数。高度精确的基准解可以在文献中查到[12—14]。

方腔内流动由顶盖以速度U从左向右移动诱导产生。其他三个静止边界采用0速度无滑移非平衡外推处理方法。模拟中使用的顶盖均匀速度为U=0.1。顶盖驱动流的基本特性是在中心附近出现一个主涡和在下方两个角落产生两个次涡。涡中心的位置是雷诺数Re=LU/ν的函数。其中L是方腔的高度, U是顶盖平移速度, ν是运动黏性。

使用He-Luo模型和Shan模型数值计算方腔流使用的雷诺数分别为400, 1 000, 2 000和5 000, 网格节点为257 257。稳定状态当ε=1.0×10-6时得到。

图3方腔流在不同雷诺数下的流线 (A) Re=400 (B) Re=1 000 (C) Re=2 000 (D) Re=5 000, (说明:每张流线图的坐标原点都在其方形区域的左下顶点, 右侧为x轴正向, 上方是y轴正向)

图3给出了不同雷诺数下方腔流内的流线。这些图清晰地显示出不同雷诺数下的流动结构。对低雷诺数而言, 只有三个涡, 主涡靠近中心而一对次涡出现在角落里。在Re数为2 000时, 左上角出现了第三个次涡。当雷诺数达到5 000时, 右下角出现了三级涡。而且当雷诺数增加时也可以看到主涡的中心向方腔中心移动。

为了评价计算结果, 主涡和次涡中心位置被列在表2中, 其中a为Ghia, b为Vanka, c为Hou, d为He-Luo, e为Shan模型。从表格2可以看到, HeLuo模型预测的涡中心位置跟测试基准的结果很接近, Shan模型在低雷诺数时吻合还好, 高雷诺数时偏差比He-Luo模型要大。表格第三列是计算耗时, 单位是min。相比之下, He-Luo模型普遍比Shan模型快。Shan模型的计算耗时比He-Luo大, 最大时为He-Luo用时的1.76倍, 同为单松弛模型, 但是Shan模型要处理额外的温度变量, 因此消耗了更多的时间。

通过方腔中心的水平和垂直线上速度分量u, v的相对总体误差以Ghia的基准数据作为参照。从表3中可以看到, Shan模型的相对总体误差要普遍高于He-Luo模型, Shan模型在雷诺数5 000时的速度u、v和温度的相对总体误差是其在雷诺数400时的7, 2和5倍。雷诺数5 000时Shan模型的平均密度和平均温度分别是平均密度1和平均温度1/3的1.35和1.22倍。在高雷诺数顶盖驱动方腔流动中, 相比He-Luo模型Shan模型的波动比较强烈。Shan模型的平均密度波动幅度将近He-Luo模型的20到30倍, Shan模型速度u的相对总体误差随着雷诺数的增加而增加。

3 结论

本文研究了He-Luo模型和Shan模型, 测试和对比了典型的泊肃叶流动和顶盖驱动流的计算结果。仿真结果同分析解或相关实验结果也进行了对比。

对平板间由固定压力梯度驱动的稳态泊肃叶流动而言, He-Luo模型和Shan模型跟分析解吻合得都很好, 都得到了抛物线型的速度曲线。He-Lu模型计算出来的速度u的最大相对误差, 最小相对误差和相对总体误差都比Shan模型要小。这两个模型计算消耗时间是相同的。

对顶盖驱动流而言, 使用He-Luo模型和Sha模型分别计算了雷诺数为400/1 000/2 000/5 000时的情况。这些模型都能预测雷诺数2 000左右时在左上角出现的第三个二次涡和雷诺数5 000时右下角出现的三次涡。对主涡和两个二次涡的位置进行了测量并跟文献上的数据进行了对比。仿真时间也被分别记录下来。总体来说Shan模型比He-Luo消耗时间要多, 最高消耗时间达到了He-Lu的1.7倍。沿着穿过腔体中心的水平线和竖直线上速度u、v的相对总体误差也进行了统计。Shan模型的相对总体误差比He-Luo模型的大。速度u的相对总体误差随着雷诺数的增加而增加。在雷诺数5 000时, Shan模型相对平均密度和平均温度的波动比较剧烈。

总体而言, He-Luo模型的密度的波动性要小于Shan模型, 计算也要比Shan模型要快, 忽略可压缩项的优势很明显。但是He-Luo模型只能用于恒温系统, 而Shan模型同时能够计算非恒温系统, 这是其一个优势, 但是在计算恒温系统时由于温度的累计误差造成密度的波动较大。Shan模型还有待改善。

多项式模型 篇2

所谓多项式曲面拟合, 其主要目的是利用一些已知高程异常值的离散点来拟合出局部区域内的高程异常值, 通过GPS水准拟合获得所测区域的高程异常值的分布情况, 从而将GPS所测点的大地高转化为正常高。

多项式曲面拟合的数学表达式如下:

设待求点的平面坐标与待求点的高程异常值i关系式为:ξi=f (xi, yi) +ei (公式2)

上式中, 若f (x, y) =a0+a1x+a2y+a3x2+a4y2+a5xy+LL (公式3) 则称为多项式拟合。已知a0、a1、……、an是多项式曲面拟合公式的系数, 只要确定多项式拟合公式的系数, 就可以求解待定点的高程异常值。

此种方法只适用于地势起伏变化不是很大或者是地势较为平坦的测区。如果地势起伏很大或者是比较困难复杂的情况下, 运用此方法拟合误差就会比较大, 测量结果也不够精确。

当地势起伏较大时, 似大地水准面的起伏变化情况将会非常复杂更加难以拟合确定, 那么多项式曲面拟合的结果必然无法达到相应的精度标准。因此我们需要采取一些措施加以改进:

1) 从已知水准点和检验点的精度入手, 由于联测的已知水准点的精度情况直接决定了拟合后的精度, 所以要对已知水准点的精度加以检验和控制, 要使已知水准点的密度适中且尽量分布均匀, 最好能最大限度的覆盖整个测区, 这样就能够有效的减少所得高程异常值的误差, 从而提高高程异常值精度, 同时提高了正常高的精度。

2) 如果对多项式曲面拟合的方法有效合理的加入相应的改正参数或数学模型, 就可以减少复杂地形对多项式曲面拟合数学模型的影响, 就可以有效的减少复杂地形对所得高程异常值的误差影响, 就可以准确的求定在复杂地形下测区内点的高程异常值, 从而提高正常高的精度。

随着公式3中 (x, y) 幂次的不同, 多项式曲面拟合又可以分为平面拟合、二次曲面拟合、三次曲面拟合等。

1平面拟合

在地形起伏比较小的区域或地势较为平坦的区域, 可以考虑用平面逐渐逼近局部的似大地水准面的方法来拟合曲面求得高程异常值, 从而得到正常高值。在面积很小的范围内, 可以认为大地水准面近似和平面重合, 平面模型的数学表达式为:ξi=a0+a1xi+a2yi (公式4) , 平面拟合数学表达式中:x, y为平面点的坐标;a0, a1, a2为平面拟合数学模型的待定参数的系数。由已知点的正常高和GPS测得的大地高确定平面拟合模型的3个参数a0, a1, a2, 从而求得第i个拟合点的高程异常值。

2二次多项式曲面拟合

当测区内点在一定程度上布成区域面时, 可以采用二次多项式曲面拟合的数学模型来拟合局部曲面, 进而在拟合的局部曲面前提下利用多项式曲面拟合求出待测点的正常高, 也就是根据测区中的已知点的平面坐标 (x, y) 或大地坐标 (B, L) 以及相应点高程异常值ξ, 用二次多项式曲面拟合的方法来拟合局部曲面, 当多项式曲面拟合出测区的似大地水准面之后, 再通过多项式曲面拟合出该点高程异常值ξ, 最终求出待测点的正常高值。

二次曲面拟合算法的模型表达式为:ξ=f (x, y) +ei (公式5) , f (x, y) =a0+a1x+a2y+a3xy+a4x2+a5y2 (公式6) , 公式5中, e是相应的残差值, 利用一个方程就有一个对应水准联测点, 在[e2]=min条件下求出ai, 再代入公式6可求出剩余点的高程异常值。

3三次多项式曲面拟合

虽然三次多项式曲面拟合具有较好的精确性和通用性, 但是由于三次多项式曲面拟合计算过程相对困难复杂。

因此, 在精度要求不是很高的情况下, 一般不采用三次多项式曲面拟合, 三次多项式曲面拟合算法的数学模型表达式:f (x, y) =a0+a1x+a2y+a3x2+a4xy+a5y2+a6x3+a7x2y+a8xy2+a9y3 (公式7) 。

对于三次多项式曲面拟合的方法, 由于三次多项式的未知参数比较多, 而且 (x, y) 的幂次比较高, 所以在采用三次多项式曲面拟合计算高程异常值时, 计算方法和过程比较困难复杂, 因此要加倍细心。

综上所述, 相比较以上三种多项式曲面拟合方法, 如果有足够的几何水准联测拟合点和检核点, 进行高程拟合时采用二次多项式曲面拟合法能得到很好的拟合效果且更便捷。在实际工作运用中可以采用多种多项式曲面拟合的数学模型来拟合当地局部区域的似大地水准面, 通过相应的计算求出待测点的高程异常值ξ, 最终求出待测点的正常高。然后通过内符合精度、外符合精度、平均误差和拟合误差的分析, 选择最佳最优方案, 通过实践和研究证明只要使已知水准点密度适中且尽量的分布均匀, 最好能最大限度的覆盖整个测区, 这样就能够有效的减少所得高程异常值的误差, 同时也提高了正常高的精度。在满足水准高程的精度足够高的情况下, 即内符合精度、外符合精度、平均误差和拟合误差都满足的情况下, 在满足测量要求精度的前提下完全可以用GPS高程测量的成果来代替普通低等级水准测量, 且测量精度比普通低等级水准测量的精度还要高。

多项式曲面拟合是GPS高程拟合的一个热点, 有效和快速提高GPS高程转换精度, 可以提高测量精度和降低外业测量的困难程度。与传统水准测量相比, 既能满足全天候、自动化、高精度的要求, 又可以实现测量速度快、节省大量的人力和物力等优点, 从而真正的实现GPS测量的无可比拟的优越性, 使备受青睐和广泛关注的GPS有更加美好的研究价值和广阔的发展应用前景。

摘要：随着经济和科学技术的不断发展, GPS已成为测量工程领域中的一种极为重要的方法。将GPS测量得到的大地高转化为正常高通常来说是非常繁琐的, 但如果已知每个GPS点的高程异常值ξ, 那么将GPS大地高转换为正常高就会变得相对容易, 这样就可以在保证精度的前提下实现GPS高程测量代替普通水准测量。曲面拟合是实现这一过程的有效方法, 根据数学模型和原理的不同, 曲面拟合法主要分为:多项式曲面拟合法、多面函数拟合法、移动曲面拟合法等。

关键词：水准点,高程异常,多项式曲面拟合

参考文献

[1]孔祥元, 郭际明, 刘宗泉.大地测量学基础[M].武汉大学出版社, 2003.

[2]李征航, 黄劲松.GPS测量与数据处理[M].武汉大学出版社, 2005.

多项式模型 篇3

桥梁变形观测是桥梁运营期养护的重要内容,对桥梁的健康诊断和安全运营有着重要的意义。桥梁挠度是指桥面在外界荷载的作用下沿轴线的垂直位移情况［1］。桥梁挠度是评价桥梁安全性的重要指标,直接反映桥梁结构变形是否超出危险范围,过大的桥梁挠度变化不但影响行车安全,而且对桥梁的使用寿命有直接的影响,同时桥梁挠度监测是桥梁改造和加固设计的重要依据,有助于从根本上消除隐患及避免灾难性事故的发生,因此,对桥梁进行挠度监测十分必要。

对于特大跨度预应力混凝土连续钢箱梁桥,其挠度变化受到多种因素的影响［2］,如竖向预应力、温度、动荷载、混凝土收缩、徐变以及施工质量等多种因素的影响,呈动态性和非线性的特点,给桥梁挠度预测带来了很大的困难［3］。本文采用PDL模型,以南京某钢箱梁桥挠度监测数据为例,重点分析温度和车流荷载对挠度的影响,建立挠度与影响因素之间的多项式分布滞后模型,其思想是桥梁挠度同时受到自身滞后项以外其他多个因素的影响,并分布到多个时段。同时将PDL模型和灰色GM( 1,1) 模型的预测结果进行比较分析,从而更为精确地评价桥梁的健康状况。

1 模型建立

1. 1 PDL模型

PDL模型是一种动态的回归模型,不仅考虑了影响因子跨时段的影响关系,还加入因变量自身滞后项的影响,集合了时间序列模型和回归模型的优点,是一种解释能力比较强的预测模型［4］,其基本形式如下:

其中: C为常数项,yt － i是因变量的滞后项,x表示影响因子,xt － j为影响因子滞后项,αi和 βj分别为因变量滞后项和影响因子滞后项的系数,体现了各个滞后值对因变量的不同影响程度; ut为误差项; t为时间序列; m和n为滞后项数,i和j代表滞后量。影响因子可以是一个也可以是多个,若滞后期长度是一个确定的数,则称为有限滞后变量模型,若滞后期长度无限,则模型为无限滞后变量模型［5］。

1. 2 模型确定

桥梁挠度的变化受到多种因素的影响,因此影响因子的选取十分重要,如何利用现有的资料进行挠度预测值得深入研究,由于缺少日照、预应力、混凝土收缩等方面的测量资料,且各因素的复杂不确定性,本文在分析建模时仅考虑了温度和车流量荷载两个影响因子,则相应的PDL模型为:

式( 2) 中,Q表示车流量,T表示温度,βj和 γk分别为车流量和温度因子的系数,由于测量资料按月份记录,因此这里的车流量和温度都是取当月平均值。

1. 3 滞后阶数的确定

利用PDL模型进行预测时,首先要利用样本数据确定各影响因子的滞后阶数。对于有限分布滞后模型,最大滞后阶数比较难确定,而且由于样本容量有限,当滞后阶数增加时,必然使得自由度减少,而且还要考虑选取的滞后阶数之间的共线性问题［6］。由于数据的测量常常受到各种条件的限制,估计这类模型时常常会遇到数据不足的困难。滞后阶数究竟多大为好,没有明确的规定,具体需要根据样本数据以及模型来确定［7］,也可以通过一些常用的统计检验来确定,如赤池信息准则( AIC) 、施瓦茨准侧( SC) 以及调整后的修正系数。

式中是估计模型的残差平方和,k是模型中解释变量的个数,N是样本容量。在模型中逐渐增加滞后变量,扩大滞后期长度,直到模型的不再增加为止,或逐渐增加滞后变量,直到AIC和SC统计量不再降低为止。具体确定滞后阶数时,应综合考虑各方面的因素从而得到最优的结果。

1. 4 系数估计

PDL模型对系数的估计十分复杂,具体算法可参考文献[8]。本文采用阿尔蒙( Almon) 法来估计系数,即在以滞后期i为横轴、滞后系数取值为纵轴的坐标系中,如果这些滞后系数落在或近似落在一条光滑曲线上,则可以用一个关于i的次数较低的m次多项式逼近。阿尔蒙法利用有限多项式来减少待估参数的个数,消除了多重共线性及参数估计中自由度的损失,实际应用时,阿尔蒙多项式的次数通常取得较低,通常取2 或3,一般不会超过4。本文取2 次多项式来逼近,即:

2 工程应用

2. 1 工程概况

长江某大桥全长约15. 6km,其中跨江大桥长4744m,主桥为跨径648m的双塔双索面钢塔钢箱梁斜拉桥。南引桥长680m,北引桥长2780m。挠度监测点沿引桥两侧布设,上下游幅每隔50m设一个点,南引桥共布设92 个点,北引桥共布设124个点,挠度监测点布设如图1 所示。

2. 2 分析建模

大桥挠度采用水准测量的方法［9］,以南引桥S45 号点为例分别建立PDL模型和灰色GM ( 1,1)模型。首先确定滞后阶数,滞后阶数从1 开始逐渐增加,利用EVIEWS软件进行建模分析,计算各统计信息,结果如表1 所示。

表1 中R为可决系数;为调整的可决系数,系数越接近于1,说明拟合得越完美; AIC和SC值越小表明滞后阶数选取越好; F-statistic为F统计量,F越大说明模型越理想; RSS为残差平方和。当i、j、k取其他数值时,会发生自由度缺失而无法建模的情况,因此综合表1 所得信息,选取自身、车流量以及温度滞后阶数分别为2,3,3。要说明的是,这里的滞后不是说明过去的温度和车流量对现在的桥梁挠度产生影响,而是由于挠度的监测并不是每月都进行,因而仅用当月的温度和车流量来预测监测期的挠度值是不可靠的,应该考虑到挠度的变化是一个长期、累积的过程。

确定滞后阶数之后,就可利用EVIEWS进行建模,得到的预测值和GM ( 1,1) 预测结果如表2所示。

由表2 可以看出,PDL模型的预测值比GM( 1,1) 模型精度更高,与实测值非常接近,最小残差仅为0. 0027mm,相对误差最大的也只有0. 31% ,残差平方和仅有0. 0178mm2,预测效果非常理想,几乎与实际挠度曲线完全拟合。而本实验中灰色GM ( 1,1) 模型的预测效果并不理想,最大残差达到了4. 9596mm。从图2 可以很直接看出两种模型的残差比较。

( 单位: mm)

将PDL模型推广到其它挠度监测点,并与灰色模型预测结果进行比较,结果如表3 所示。

从表3 可以看出,PDL模型的预测精度均优于GM ( 1,1) 模型,说明对大桥进行挠度预测时仅考虑自身历史变化量是远远不够的,还应考虑到温度、车流量等复杂的影响因子引起的变化,而事实上这些影响因子确实对大桥挠度有一定的影响,PDL模型正是考虑到了这些因素,因而预测结果比GM ( 1,1) 模型更为可靠。

3 结论与建议

( 1) 通过对长江某大桥的挠度进行建模预测表明,用多项式分布滞后模型来预测挠度变化,所得结果与实际情况吻合较好,模型可靠性较高。

( 2) PDL模型是一种动态的参数回归模型,不仅能充分利用自身的变化规律,而且还能有效利用影响因子的历史时期数据来进行建模预测,是一种高精度的预测模型。

( 3) 在确定模型滞后阶数时,应综合考虑各种统计检验数据,以及自身样本数据的大小,避免自由度不足和多重共线性的问题,同时应该注意出现的伪回归问题。

( 4) PDL模型能综合利用各有效信息进行预测,但是所有的预测都是建立在原始观测数据基础之上的,因此要保证原始观测数据的精度,并不断补充新的观测数据,扩大样本容量,以便提高预测精度。

多项式模型 篇4

关键词：图像盲复原,运动去模糊,有理数多项式函数,自然图像梯度分布,振铃效应

在数码相机日益普及的今天, 相机抖动或者物体运动造成的图像模糊现象严重影响了使用者对数码产品的体验, 怎样从模糊图像中恢复出清晰图像成为了各国科学家研究的热点问题。模糊图像被复原成清晰图像属于图像复原, 这一复原过程被称为图像去模糊。由于模糊核函数是未知的, 故这种图像复原的过程称为模糊图像盲复原。图像去模糊技术在摄影、航空、航海、地质勘探等领域起着显著的作用。对单幅图像进行图像盲去模糊非常困难, 由于只有一个已知量 (模糊图像) , 而要求出两个未知量 (清晰图像和模糊核) , 因此该过程是一个病态问题的求解过程。求解这样的问题通常是对模糊图像或者模糊核增加一些先验约束进行正则化求解。

Fergus[1]等人利用图像的先验知识和模糊核的约束条件对模糊图像进行正则化求解。通过大量自然图像的对比试验, 提出清晰图像与模糊图像在图像的梯度分布上区别很大, 清晰图像的梯度符合长拖尾分布, 而模糊图像的梯度不符合这种特征。利用混合高斯函数作为约束函数对自然图像梯度分布进行拟合, 并采用变分贝叶斯方法估计出模糊核, 通过传统R-L算法恢复原始图像。由于图像梯度分布函数具有广泛性, 因此该方法可以得到较好的去模糊效果。但是由于混合高斯函数拟合自然图像梯度分布时, 在梯度值较大的地方往往会出现过拟合现象, 从而导致去模糊后的图像产生了严重的振铃效应。

本文对Fergus算法进一步改进, 采用有理数多项式拟合自然图像梯度分布。虽然与混合高斯函数相比, 有理数多项式函数在拟合自然图像梯度特征时在误差项平方和 (SSE) 、确定系数 (R-square) 、校正的决定系数 (Adjusted R-square) 、均方根 (RMSE) 等拟合优度标准上只具有略微的优势, 但是有理数多项式函数大大的削弱了拟合过程中的过拟合现象, 从而使模糊核的估计更准确, 减小了由此而带来的振铃效应, 对模糊图像的复原更加有效。

1 模糊图像模型

运动模糊退化图像的数学模型一般为

式中:B为输入的模糊退化图像;K为模糊核;L为清晰图像;N为成像过程中加性噪声。

图像去模糊的目标就是求解清晰图像L。一般情况下, 是在抑制噪声的情况下, 先求解模糊核K, 然后通过传统的反卷积算法求解得到清晰的图像L。设p (L) 是清晰图像的概率密度函数, p (B) 是模糊图像的概率密度函数, 而p (L|B) 表示已知模糊图像条件下, 原始清晰图像的条件概率密度。根据贝叶斯定理, 则有下式成立

由于模糊核K需要利用模糊图像B来估计, 式 (2) 可以写成

在模糊图像B已知时, 可以得到如下的关系

即可得

令E (B|L, K) =-lg (p (B|L, K) ) , E (L) =-lg (p (L) ) , E (K) =-lg (p (K) ) , 其中E (B|L, K) 是图像噪声约束项, 也是模糊图像B、清晰图像L和模糊核K之间的数据匹配约束项, E (L) 是对自然图像的特征约束项, E (K) 是对模糊核的特征约束, 设E (L, K|B) 为能量项, 则有

由于上述模型已知量B的个数小于未知量L、K的个数, 所以这是一个病态问题的求解过程, 在理论上难以获得精确的解, 因此通过L、K的大量迭代, 在实验中能够获得最优解, EK (K) 和EK (L) 为能量项, 对二者进行迭代即迭代以下两个能量

这样, 把通过贝叶斯理论建立的求最大后验概率 (MAP) 问题[2]转化成求最小能量问题。计算迭代出模糊核, 最后通过反卷积得到清晰图像。

2 有理数多项式函数拟合自然图像梯度分布特征函数

通常认为自然场景图像具有平滑性或者分块平滑性, 通过图像的梯度可以体现图像的平滑性。图像的梯度值大部分情况下比较小, 分布在梯度为零值附近, 很少一部分会出现较大的值, 这说明图像梯度是稀疏的。梯度值较大的部分, 对应的是图像中的细节。以往论文都通常假定图像梯度的分布特征服从一个高斯模型, 即

式中:∀F为图像梯度, 利用图像梯度的分布约束能够减少图像去模糊过程中的不确定性, 但是在假设中强制图像梯度分布服从高斯分布, 导致模糊核估计不准确, 恢复结果过于平滑, 因此引入较强的振铃效应[2]。

实际上, 自然场景图像梯度并不服从高斯分布, 因此利用高斯模型进行图像去模糊的效果并不是很好。本文对图像先验约束模型进一步改进, 提出用有理数多项式函数拟合自然图像梯度分布特征函数。大量图像分析统计发现, 自然场景图像的梯度分布具有一个很明显的长拖尾曲线特征[3,4]。自然场景图像具有很强的平滑性, 大部分梯度值都在零值附近, 而图像梯度分布图中很长的拖尾是由梯度值较大的点造成的, 对应的是自然场景图像中的细节。

图1是利用20幅自然图像统计的自然图像的梯度分布, 其横坐标为梯度值, 纵坐标为梯度值对应的概率的对数。

Fergus采用混合高斯函数拟合自然图像的梯度分布

式中:an, bn和cn为高斯函数参数。

图2即为混合高斯函数拟合的效果图, 其中黑点为经过计算统计得到的图像真实概率密度, untitled fit 1为混合高斯拟合函数。可以看到在拟合过程中, 混合高斯函数拟合在梯度为100~150时, 出现比较严重的过拟合现象。

为了防止过拟合现象的出现, 本文采用有理数多项式函数拟合自然图像的梯度分布

在有理数多项式函数中, 通过曲线拟合确定有理多项式中的系数p和q。调节多项式的项数可以对拟合效果进行调整, 但是考虑到项数多会产生冗余, 增加了计算复杂度和计算时间, 故本实验采用上述函数。

图3为有理数多项式函数拟合的效果图, 其中黑点为图像的真实概率密度, untitled fit 2为有理数多项式函数拟合函数。可以看到该函数克服了拟合过程中出现的过拟合现象。

为了更为详细地对比两种拟合算法, 本文计算了两种拟合函数的拟合优度

式中:Xi为原始数据;为预测值;为原始数据平均值。

如表1所示, SSE与RMSE越接近于0, 说明函数拟合越好, 数据预测越成功。决定系数R-square与校正后决定系数Adjusted R-square是通过数据的变化来表征函数拟合的好坏, 值越接近1说明拟合越好。可以看到本文的拟合算法略胜于混合高斯函数拟合算法。综合考虑本文的拟合算法必然可以带来更好的去模糊效果。

由于有理数多项式函数的分母多项式强制函数在梯度值较大值时平滑, 所以有效地避免了过拟合现象的发生, 从而使得模糊核估计得更准确, 减少振铃效应。本文通过集成学习[5,6,7], 首先利用新的拟合函数进行能量项迭代求得模糊核, 然后利用Richardson-Lucy算法[8,9]进行反卷积得到去模糊图像。

3 实验结果与分析

本文选择了人工合成图像与自然拍摄的模糊图像进行实验, 进而说明有理数多项式逼近自然图像的梯度分布作为先验约束的有效性。采用峰值信噪比PSNR (Peak Signal to Noise Ratio) 和均方误差和MSE (Mean Square Error) 这两个指标评估两种约束函数的复原结果。

图4为使用混合高斯约束去模糊的结果和有理数多项式约束去模糊的结果, 图4a为清晰图像, 图4b为人工模糊图像, 其中模糊核是均值为0, 方差为5的高斯模糊核, 图4c为使用混合高斯约束得到的去模糊图像, 图4d为使用有理数多项式约束得到的去模糊图像。在图4c中, 平坦的背景区域出现了明显的振铃效应, 这是因为在使用混合高斯函数拟合自然梯度图像时, 在大梯度值区域出现了过拟合现象, 导致模糊核估计不准确, 从而使得图像恢复失真, 出现振铃。在图4d中, 由于使用有理数多项式拟合自然图像梯度, 抑制了过拟合现象, 所以图像的细节部分保存相对完整, 且模糊核估计相对准确, 振铃效应减少。

图5为自然拍摄的模糊图像的去模糊效果图。图5a模糊图像, 图5b为使用混合高斯约束得到的去模糊图像, 图5为使用有理数多项式约束得到的去模糊图像。在图5b中, 也有明显的振铃效应。而在图5c中, 由于使用有理数多项式拟合自然图像梯度, 抑制了过拟合现象, 振铃效应得到了明显的抑制。

表2为使用有理数多项式约束的去模糊图像与混合高斯约束的去模糊图像的PSNR、MSE比较。从表中结果可以看出, 有理数多项式拟合自然图像的梯度分布作为约束在这两项客观评价上有提升, 但是在运算速度上, 由于算法本身的自然对数限制, 所以使用多项式拟合后, 迭代所需要的时间更长。

4 结论

本文提出了基于有理数多项式的图像先验模型约束, 利用有理数多项式拟合自然图像梯度分布, 通过选取合理的拟合参数, 得到拟合优度良好的拟合函数, 避免了由于混合高斯函数过拟合现象带来的模糊核估计不准确, 从而得到了较佳的去模糊效果。使用人工模糊图像与自然模糊图像的复原结果证明了本文算法的有效性和可行性。

参考文献

[1]FERGUS R, SINGH B, HERTZMANN A, et al.Removing camera shake from a single photograph[C]//Proc.ACM Trans.Graphics (TOG) .[S.l.]:ACM, 2006:787-794.

[2]SHAN Q, JIA J, AGARWALA A.High-quality motion deblurring from a single image[C]//Proc.ACM Trans.Graphics (TOG) .[S.l.]:ACM, 2008:73.

[3]KRAHMER F, LIN Y, MCADOO B, et al.Blind image deconvolution:motion blur estimation[J].IMA Preprints Series, 2006 (21) :33-35.

[4]STEFAN R, MICHAEL J.Fields of experts:freamwork for learning image priors[J].Computer Vision and Pattern Recognition, 2005 (10) :56-59.

[5]YAIR W, FREEMAN W T.What makes a good model of natural images[J].Computer Vision and Pattern Recognition, 2007 (8) :14-20.

[6]MISKIN J, MACKAY D.Ensemble learning for blind image separation and deconvolution[J].Neural Computing, 2000 (9) :123-141.

[7]OSHER S, RUDIN L I.Feature-oriented image enhancement using shock filters[J].SIAM Journal on Numerical Analysis, 1990, 27 (4) :919-940.

[8]MISKIN J W.Ensemble learning for independent component analysis[D].Cambridge:University of Cambridge, 2000.

[9]RICHARDSON W H.Bayesian-based iterative method of image restoration[J].JOSA, 1972, 62 (1) :55-59.

多项式模型 篇5

随着居民出行多样化和小汽车的使用增加, 老城区的交通问题越来越严重, 然而由于对历史文化的保护, 不可能进行大规模基础设施建设, 因此有必要研究居民出行, 特别是通勤出行。通勤者出行是指居民在家和工作地之间的空间移动, 早晚高峰出行量集中加剧了老城区交通供求矛盾。出行方式选择问题是交通规划和政策制定中的重要部分, 它能够评价城市交通方式分布的合理性, 从而制定适当的管理策略。

传统的“四阶段”法往往以交通小区作为研究对象, 难以体现个人特征的差异。随着针对个体特征模型技术的日益成熟, 国内外学者开始引入非集计建模技术[1,2,3,4]。文献[1]运用MNL模型分析通勤者出行方式选择的影响因素。文献[2]采用NL模型研究通勤者出行方式与出行链相互选择。

1 研究数据

本文研究数据来自2013年合肥巢湖市老城区的居民出行调查数据。合肥巢湖市老城区位于城市核心, 面积约5 km2。调查表主要调查通勤者日常出行时的方式选择和出行者个体特征, 如通勤者性别、职业和年龄, 以及通勤者家庭特征, 如家庭人口数、家庭收入等作为模型特性变量。调查期间一共发放调查表1 600份, 有效问卷1 232份。

根据以往建模和分析经验[5,6], 居民出行方式的选择除了受到社会家庭特征影响, 还受到出行特征的影响, 比如出行距离、出行时间等。最后, 选择影响通勤者出行方式选择的个体特征、家庭特征和出行特征作为研究变量, 每个特征所包含的变量见表1。

根据出行数据的调查统计, 用于分析的出行方式主要有自行车、电动车、公共交通和小汽车4种。4种出行方式的样本量和比例见表2。

由表2可知, 合肥巢湖市老城区电动车出行比例为38.3%, 为4种方式比例最高;其次为小汽车比例为24.3%;自行车出行和公共交通出行差别不大, 分别为20.9%和16.5%。

2 模型介绍

多项logit模型的适用范围为因变量为无序种类变量, 适合3种及3种以上因变量情况下建模, 能够得到自变量与因变量之间的关联和影响。本研究针对老城区出行方式选择, 其因变量为4种出行方式, 符合多项logit建模要求。

根据随机效用理论, 第n个个体选择第i种出行模式的效用Uin可以表示为:

式中:Vin为第n个个体选择第i种出行模式效用函数中的固定项;εin为第n个个体选择第i种活动模式效用函数中的随机项。

当Vin与其中包含的解释变量之间呈线性关系时, 可以表示为:

式中:K为解释变量个数;θ为参数矩阵;θk为第k个变量所对应的参数;xink为第n个个体选择第i种活动模式的第k个特性变量。

假设效用函数中随机项服从二重指数分布, 可以得到第n个个体选择i种活动模式的概率, 即

式中:Vjn是第n个个体选择第i种出行模式效用函数中固定项;xjnk为第n个个体选择第j种活动模式的第k个特性变量。将式 (3) 变形, 并运用极大似然估计法及牛顿-拉普松求解, 便可估计模型中的参数θ1, 。

3 模型结果与分析

用统计分析软件SPSS的multinational logistic regression模块标定通勤者出行方式选择模型。出行方式分为自行车、电动车、公共交通和小汽车。解释变量为通勤者个体特征、家庭特征和出行特征, 各个解释变量的分类和取值见表3。

SPSS在运行时通过建立0-1变量, 将J类多变量设置为J-1个二值哑变量, 并自动将各类变量中的最后一类作为参考类别。本文的建模结果见表4。

由表4可知老城区通勤者出行方式与个人特征、家庭特征和出行特征的影响关系如下:

(1) 男性在自行车、电动车和公共交通的标定值分别为-1.891、-2.126、-1.802, 均为负值, 说明男性选择这3种出行方式的概率相比女性更低。

注:以小汽车为参考类别;空白表示结果不显著。

(2) 职业变量工人在电动车的标定结果为0.24, 为正值, 说明工人倾向选择电动车作为日常的通勤方式。年龄小于40岁的通勤者在自行车的结果为正值, 表明该类出行者选择自行车的概率更大。

(3) 家庭人口较少的通勤者 (≤3人) 估计结果为-0.933, 表明人口较少的通勤者选择公共交通的概率较低。而低收入家庭 (<3万) 的模型结果为0.376, 也就是说选择电动车出行的概率较高。

(4) 家庭没有电动车和小汽车的家庭, 在自行车方式下的结果均为正值, 说明更倾向选择自行车出行;没有电动车家庭在电动车的估计值为负值-1.75, 表明该类家庭不会选择电动车出行, 与实际情况相符合;没有小汽车家庭在公共交通的标定值为2.751, 为正值, 说明没有小汽车的家庭选择公共交通的概率更大。

(5) 若通勤者不在早高峰时段出行, 这类出行者更加愿意选择电动车出行;出行时间小于20 min在电动车标定结果为正值, 而在公共交通值为负, 说明这类出行者选择电动车的意愿大于公共交通, 主要因为电动车相比公共交通在短距离有较高的优势。

(6) 根据已有的以整个城市为对象的研究结果, 个体特征、家庭特征和出行特征对出行方式的影响与本文研究结果一致, 比如文献[2]~[4]结果表明性别、职业、家庭收入等能够显著影响出行方式的选择。另外, 老城区相比城市的其他地区, 当出行距离或者出行时间较短时, 通勤者往往会更加愿意选择电动车出行, 该点与文献[3]中的结果有所差异。

4 结语

本文研究结果表明, 老城区通勤者个人特征和家庭特征对出行方式选择影响因素与以往研究结论一致。由于城区用地特征, 以及该地区居民出行特征, 居民更加倾向选择电动车。但是随着电动车的普及, 电动车往往作为出行首要考虑对象。未来随着机动化进程的加快, 除了限制小汽车出入老城区外, 还应探究针对电动车的定位, 以便制定相应的交通管理政策。

参考文献

[1]Central Transportation Planning.Transfer penalties in urban mode choice modeling (travel model improvement program) [M].Washington D.C.:Travel Model Improvement Program, 1997:30-45.

[2]栾琨, 隽志才, 宗芳.通勤者出行方式与出行链选择行为研究[J].公路交通科技, 2010, 27 (6) :107-111.

[3]宗芳, 隽志才.基于活动的出行方式选择模型与交通需求管理策略[J].吉林大学学报 (工学版) , 2007, 37 (1) :48-53.

[4]万霞, 王炜, 陈骏.居民全日出行方式选择动态模型研究[J].中国公路学报, 2012, 25 (2) :121-126.

[5]J.L.Bowman, M.E.Ben-Akiva.Activity-based disaggregate travel demand model system with activity schedules[J].Transportation Research Part A:Policy and Practice, 2001, 35:1-28.