土样条曲线的插值方法

土样条曲线的插值方法

1cu/cc值的计算机自动计算在土木工程试验中，cu和曲线系数c反映了土壤结构的特征。计算Cu/Cc的传统方法是依靠人工,根据原始数据计算出各粒径范围的含量百分数,在专用坐标纸上标出各数值点,并凭个人经验将各点连接成一条光滑的曲线,即颗粒大小分配曲线(以下简称颗分曲线),再从曲线上查取含量百分数分别为60%、30%、10%时对应的粒径值d60、d30、d10,根据公式求得Cu/Cc值。以上方法相对于电算有如下几点不足:一是计算过程繁琐,工作效率低下;二是人为影响因素多,精度较低。应用计算机取代上述人工计算和图形绘制工作,则主要是要解决如下几个问题:(1)插值函数的选择及计算边界条件的确定;(2)半对数坐标系与计算机逻辑坐标系的转换;(3)正确建立和使用三次样条插值函数。2关于颗分曲线的绘制方法在观测资料整理分析过程中,要求通过平面上已知n点(xi,yi)(i=1,2,…n)作一条光滑曲线的作图方法有多种,其中之一是拉格朗日插值公式。实践证明,拉格朗日插值对于较为规则的点分布情况,如曲线形状基本呈直线或二次抛物线,则能较好地满足要求,但对于基点较多而且分布不规则的情况,如土样的颗粒大小分配情况,做一个高次插值多项式是不理想的,因为它带有近似性,且计算也很复杂。分段线性插值和抛物线插值方法固然计算简单,但是它只能保证插值曲线的连续性,却不能保证曲线的光滑性(即各段连续点处导数的连续性不能保证)。我们很自然地想到样条插值,在施工放样中,例如克～奥溢流坝面曲线的放样,或者江河堤岸线的放样等,工程设计人员往往只知道曲线上的若干点的坐标位置,这些点称为样点,他们要用细长的样条(相当于曲线板)先把相近的几点连成曲线,接着再把另外几点连成另一曲线,并注意曲线接头处也是光滑的。这样就把所有的点都连成一条完全光滑的曲线,就称为样条曲线。颗分曲线的绘制方法也是如此,只是边界条件不同(两端的样点不是终点,而且不具有周期性)而已,需要作一些特殊的处理。样条函数有多种,其中三次样条函数最为常用,它既满足一般实际问题的要求,而建立过程相对又不太复杂,所以文中采用了三次样条函数。满足上述要求的三次样条插值函数的数学描述如下:设平面上给定n个点(xi,yi)(i=0,1,2…n)而且a=x1<x2<…<xn=b,假设通过这n个点的曲线y=f(x)在(a,b)区间上有二阶连续导数存在,在区间(a,b)上做f(x)的样条插值函数S(x),使它满足如下条件:(1)S(xi)=yi(i=0,1,2…n);(2)在区间(a,b)上存在一阶及二阶连续导数,以保证连接处曲线是光滑的;(3)在每一个子区间(xi,xi+1)(i=0,1,2…n-1)上S(x)都是三次多项式,其常用的表达式为:S(x)=[3h2i(xi+1-x)2-2h3i(xi+1-x)3]yiS(x)=[3h2i(xi+1−x)2−2h3i(xi+1−x)3]yi+[3h2i(x-xi)2-2h3i(x-xi)3]yi+1+hi[1h2i(xi+1-x)2-1h3i(xi+1-x)3]mi-hi[1h2i(x-xi)2-1h3i(x-xi)3]mi+1(1)+[3h2i(x−xi)2−2h3i(x−xi)3]yi+1+hi[1h2i(xi+1−x)2−1h3i(xi+1−x)3]mi−hi[1h2i(x−xi)2−1h3i(x−xi)3]mi+1(1)上式中,hi=xi+1-xi,mi=S′(xi)(x=xi处的一阶导数值)。3测定两端先作线性插值根据我们手工绘制颗分曲线的经验,将所有的点连接成曲线后还要作适当的延长,而且延长是顺势的,延长部分近似于直线,所以在进行插值计算前,应该在曲线的前端和末端各增加一个点。具体的方法是在原有的样点基础上,在两端先作线性插值,插值的位置是在两端延长一个数量级粒径单位,当然,超出坐标范围的点应该剔除。因此我们可以假设曲线中的x0及xn处的斜率是已知的,即S′(x0)=m0=y1-y0x1-x0(2)S′(xn)=mn=yn-yn-1xn-xn-1(3)S′(x0)=m0=y1−y0x1−x0(2)S′(xn)=mn=yn−yn−1xn−xn−1(3)(2)式和(3)式就是利用颗分试验数据建立样条插值时的边界条件。边界条件确定的恰当与否,可以从颗分曲线的效果和Cu/Cc值的计算对比中判断。4时域中的逻辑关系正确绘制颗分曲线是计算Cu/Cc的前提,而绘制曲线则首先要建立相应的坐标系。在绘制颗分曲线的坐标系中,横坐标是土样粒径大小d的半对数值,与计算机显示的逻辑坐标之间是非线性对应关系,计算机图形显示采用既定的逻辑坐标系(在本应用研究中采用MM-TWIP的图形影射模式,即在显示界面中,1/1440英寸作为一个逻辑单位),所以必须找出两者的对应关系。按照试验规范,需要表达的粒径范围为0.001～200mm,假设图形显示的横坐标范围是x0～x1,则粒径大小为d的对应的逻辑横坐标xd为:xd=x0+x1-x0log0.50.001-log0.5200(log0.5d-log0.5200)(4)xd=x0+x1−x0log0.50.001−log0.5200(log0.5d−log0.5200)(4)将各粒径值和小于其的含量百分数转换为相应的计算机图形坐标系统中的逻辑坐标后,才能进行样条函数插值等计算。5样本插值的实现5.1计算维护参数下的几何点正确建立三次样条插值函数,是正确绘制颗分曲线和计算Cu/Cc的前提。建立三次样条插值函数有多种办法,本文介绍的是参照建立拉格朗日插值公式方法的推导过程建立的,从表达式(1)中可以看出,函数建立的过程实质上是求解mi的过程,具体可分为以下三个步骤:(1)由a1=α1/2,b1=β1/2出发,按下式计算αi、βi令αi=hi-1hi-1+hiβi=3[1-αihi-1(yi-yi-1)+αihi(yi+1-yi)](5)(i=2,\:,n-1)αi=hi−1hi−1+hiβi=3[1−αihi−1(yi−yi−1)+αihi(yi+1−yi)](5)(i=2,\:,n−1)上式中αi、βi是计算的过程变量,无实际意义,以下ai、bi同;yi是第i个试验数据某粒径含量百分数转换后的纵坐标值;hi=xi+1-xi,xi是第i点粒径根据(4)式转换后的横坐标值。(2)利用下式计算ai、biai=αi2-(1-αi)αi-1bi=βi-(1-αi)bi-12-(1-αi)ai-1(i=2,3,\:,n)(6)ai=αi2−(1−αi)αi−1bi=βi−(1−αi)bi−12−(1−αi)ai−1(i=2,3,\:,n)(6)(3)按下式计算mimi=bi-aimi+1(i=n-1,n-2,\:,1)(7)mi=bi−aimi+1(i=n−1,n−2,\:,1)(7)将求得的mi代入(1)式,即可得出各个子区间的样条插值函数S(x)。5.2步长确定及相应参数确定由于三次样条函数是分段函数,即在不同的自变量区间有不同的表达式,所以在实际应用中要注意判断自变量所在的子区间。使用三次样条插值函数求解曲线上非样点的坐标值时,插值的步长通过试算法确定,即确定一个步长后,看其绘制(打印)曲线的效果,如光滑度不够,则逐步减小步长值,直到曲线光滑度达到要求为止。在本应用研究中,使用50个图形逻辑单位(相当于0.034英寸)作为样条插值的步长,经土工专家确定,其曲线精度及光滑度已完全满足要求。已知一个方向坐标值,在曲线上查取对应的另一个方向的坐标值时,如知道粒径求相应含量百分数或已知百分数求相应的粒径大小,则步长要取得足够小,以保证搜索结果值的精度。在表1、表2的电算值计算中,搜索步长为5个逻辑单位,完全保证了计算结果的精度。5.3粒径大小计算函数建立后,有三处需要使用到该函数,一是对原有的数据进行插值计算;二是针对细小颗粒土的密度计法试验时,由于原始数据中粒径大小是不确定的值,而试验报告中规定的粒径范围是固定的,所以需要利用该函数计算(相当于在颗分曲线上查取)粒径为0.075、0.05、0.005mm的颗粒含量百分数;三是计算Cu/Cc时,已知含量百分数为60%、30%、10%,计算相对应的粒径大小d60、d30、d10。在实际编程过程中,由于三次样条函数为分段函数,需要计算的值的落点区间不确定,所以需要采用遍历法。即从第一个点开始,横坐标逐次递增一个步长,然后判断坐标点所属区间,并采用相应的分段函数进行计算,直至遍历整条曲线为止。6试样数据人工或实验设计表1中是随机抽取的三组使用筛析法及密度计法试验方法所得土样试验数据。图1是利用试样1建立三次样条插值函数后,绘制的颗粒大小分布曲线和计算的土粒成分结果值的列表。表2及表3所列是表1中三组试样不均匀系数Cu和曲率系数Cc值,为方便比较,同时列出了人工计算和利用三次样条函数计算的结果。由于人工试验、计算、绘图过程中存在许多影响因素,所以不能简单地将表中任何数据判断为更接近于真实值。表中的偏差值是在假设人工计算值为理论真实值的情况下计算出来的,目的是判断电算值的可靠程度。从表2、表3的数据中可以看出,人工计算值与电算值存在偏差,但