幾種插值算法在洪水預(yù)報(bào)、預(yù)警發(fā)布、水資源綜合服務(wù)系統(tǒng)中的應(yīng)用

1、幾種插值算法在洪水預(yù)報(bào)、預(yù)警發(fā)布、水資源綜合信息服務(wù)中的應(yīng)用陳兵1, 韓曦1,李會(huì)會(huì)2 （1. 北京北科博研科技有限公司, 北京 100053；2.北京金水科技發(fā)展有限公司, 北京 100053）摘 要：插值算法在工程領(lǐng)域有非常重要的作用，工程測量數(shù)據(jù)不可能無窮表達(dá)對(duì)象的所有特征值，因此需要由特征點(diǎn)及連線代表趨勢(shì)特征，比如只能用小時(shí)或者分鐘點(diǎn)水位代表水位流量特征。而實(shí)際應(yīng)用過程中我們要求的數(shù)據(jù)或要求展示的效果也苛刻，因此在實(shí)際應(yīng)用中采用插值算法，來近似的表達(dá)某些特征值用于工程計(jì)算。本文重點(diǎn)描述幾種插值算法在洪水預(yù)報(bào)、預(yù)警發(fā)布、水資源綜合信息服務(wù)領(lǐng)域的應(yīng)用，并起到非常關(guān)鍵的作用。關(guān)鍵詞：直線插值

2、（線性插值）,一元三點(diǎn),一元多點(diǎn),拉格朗日,插值插值算法,洪水預(yù)報(bào),預(yù)警發(fā)布,水資源綜合信息服務(wù)1 引言插值法又稱“內(nèi)插法”，是利用函數(shù)f (x)在某區(qū)間中插入若干點(diǎn)的函數(shù)值，作出適當(dāng)?shù)奶囟ê瘮?shù)，在這些點(diǎn)上取已知值，在區(qū)間的其他點(diǎn)上用這特定函數(shù)的值作為函數(shù)f (x)的近似值，這種方法稱為插值法。如果這特定函數(shù)是多項(xiàng)式，就稱它為插值多項(xiàng)式。插值法是函數(shù)逼近的一種重要方法，是數(shù)值計(jì)算的基本課題。本文介紹的內(nèi)容有：直線插值,一元三點(diǎn)，一元多點(diǎn)。在做工程領(lǐng)域，往往得到一堆離散的數(shù)據(jù)，現(xiàn)在想用數(shù)學(xué)公式模擬這堆離散數(shù)據(jù)。怎么辦，數(shù)學(xué)家們提出了插值問題。插值問題的提法是這樣的給定一堆數(shù)據(jù)點(diǎn)(x0, y0),

3、 (x1, y1), (x2, y2).(xn, yn)，要求一個(gè)函數(shù) y = f(x) ，要求該函數(shù)經(jīng)過上面所有的數(shù)據(jù)點(diǎn)。在洪水預(yù)報(bào)、預(yù)警發(fā)布、水資源綜合信息服務(wù)領(lǐng)域的應(yīng)用中，插值被廣泛的應(yīng)用，如水位流量關(guān)系曲線、相關(guān)圖、API預(yù)報(bào)模型、動(dòng)態(tài)臨界雨量模型、趨勢(shì)分析、圖形操作等應(yīng)用，都需要插值算法的支持。相關(guān)圖實(shí)例2 插值算法2.1 直線插值直線插值是一種使用線性多項(xiàng)式進(jìn)行曲線擬合的方法，多使用在數(shù)量分析和計(jì)算機(jī)制圖方面，是內(nèi)插法的最簡單形式。兩個(gè)已知點(diǎn)之間的直線內(nèi)插法：如果兩已知點(diǎn)(x0,y0)(x1,y1),那么根據(jù) (y-y0)/(x-x0)=(y1-y0)/(x1-x0)解方程得：y=

4、y0+(x-x0)*(y1-y0)/(x1-x0)經(jīng)過擴(kuò)展，可以計(jì)算n個(gè)已知點(diǎn)的情況。在實(shí)際應(yīng)用過程中，我們通常輸入系列值，通過程序自動(dòng)找到最接近的兩點(diǎn)，然后通過直線插值查找到適合的數(shù)據(jù)，如下水文流量關(guān)系曲線中，經(jīng)常用直線插值，下面是JAVA示例代碼:/x,y是線性點(diǎn)的系列值，t是要插入的點(diǎn)（如水文或流量）public static double interL(double x, double y,double t) double iv=0;int n=x.length;int p=1;for(int i=1;i < n-1;i+) if(t <= xi) break; /查找到接

5、近X的值p+; iv=(t*yp - xp-1*yp + xp*yp-1 - yp-1*t)/(xp-xp-1);return iv;2.2 拉格朗日在數(shù)值分析中，拉格朗日插值法是以法國十八世紀(jì)數(shù)學(xué)家約瑟夫·拉格朗日命名的一種多項(xiàng)式插值方法。許多實(shí)際問題中都用函數(shù)來表示某種內(nèi)在聯(lián)系或規(guī)律，而不少函數(shù)都只能通過實(shí)驗(yàn)和觀測來了解。如對(duì)實(shí)踐中的某個(gè)物理量進(jìn)行觀測，在若干個(gè)不同的地方得到相應(yīng)的觀測值，拉格朗日插值法可以找到一個(gè)多項(xiàng)式，其恰好在各個(gè)觀測的點(diǎn)取到觀測到的值。這樣的多項(xiàng)式稱為拉格朗日（插值）多項(xiàng)式。數(shù)學(xué)上來說，拉格朗日插值法可以給出一個(gè)恰好穿過二維平面上若干個(gè)已知點(diǎn)的多項(xiàng)式函數(shù)。拉

6、格朗日中值定理又稱拉氏定理，是羅爾中值定理的推廣，同時(shí)也是柯西中值定理的特殊情形。如果函數(shù)f(x)在(a，b)上可導(dǎo)，a,b上連續(xù)，則必有一(a,b)，使得f'()*(b-a)=f(b)-f(a) 拉格朗日中值定理的幾何意義。在(a，b)上可導(dǎo)，a,b上連續(xù)是拉格朗日中值定理成立的充分條件。這個(gè)定理表達(dá)了：1.在滿足定理?xiàng)l件的前提下，函數(shù)f(x)上必有【一點(diǎn)的切線】與【f（x）在x=a，b處對(duì)應(yīng)的兩點(diǎn)（(a,f(a)和(b,f(b)點(diǎn)的連線平行）。f'()=f(b)-f(a)/（b-a），等號(hào)后為x=a，b對(duì)應(yīng)兩點(diǎn)的連線斜率，等號(hào)前為f(x)上一點(diǎn)的導(dǎo)數(shù)的值，也就是f(x)上一

7、點(diǎn)的斜率，兩斜率相等，兩線平行。這是幾何上的理解方式。拉格朗日中值定理的幾何意義2.我們將f(x)函數(shù)求導(dǎo)，得到f'(x)，眾所周知f'(x)函數(shù)記錄的其實(shí)就是【f(x)函數(shù)在每一個(gè)瞬間的變化狀態(tài)】。即，在x=x1這一瞬間f（x）進(jìn)行了程度為f'(x1)的變化，在x=x2這一瞬間f（x）進(jìn)行了程度為f'(x2)的變化。函數(shù)由f(a)變化到f(b)的過程，其實(shí)就是f'(x)函數(shù)在(a，b)區(qū)間中記錄的變化狀態(tài)的依次累加，就是對(duì)f'(x)函數(shù)在(a，b)區(qū)間的值進(jìn)行積分的過程。那么，將這一過程中所有的變化狀態(tài)的值一起取一個(gè)平均，這個(gè)平均值的數(shù)值一定在

8、f'(x)的某一點(diǎn)上出現(xiàn)過（即f'()），因?yàn)閒(x)連續(xù)，則其導(dǎo)數(shù)也連續(xù)。這個(gè)平均值乘上變化的區(qū)間（a到b）的長度 就等于這個(gè) 變化的變化量【】。即所謂的必有一，使f'()*(b-a)=f(b)-f(a)。即，【a，b區(qū)間上f（x）函數(shù)的變化量】=【a，b區(qū)間內(nèi)f（x）函數(shù)變化狀態(tài)的平均值乘以區(qū)間長度】。這是代數(shù)理解方式。/x,y是線性點(diǎn)的系列值，t是要插入的點(diǎn)（如水文或流量）public static double getValueLagrange(double x, double y, double t) int n = x.length; int i, j, k

9、, m; double z, s;/ 初值z(mì) = 0;/ 特例處理if (n < 1) return z;if (n = 1) z = y0;return z;if (n = 2) z = (y0 * (t - x1) - y1 * (t - x0) / (x0 - x1);return z;/ 開始插值i = 0;while (i < n&&xi < t ) i = i + 1;k = i - 4;if (k < 0)k = 0;m = i + 3;if (m > n - 1)m = n - 1;for (i = k; i <= m; i+

10、) s = 1;for (j = k; j <= m; j+) if (j != i) / 拉格朗日(Lagrange)插值公式s = s * (t - xj) / (xi - xj);z = z + s * yi;/ 返回結(jié)果return z;2.3 一元三點(diǎn)利用拉格朗日(Lagrange)插值公式，選取最靠近插值點(diǎn)的3個(gè)節(jié)點(diǎn)代入公式，因此對(duì)于一條曲線就有多個(gè)插值公式，采用分段插值的方法，將該曲線使用范圍內(nèi)的每一組對(duì)應(yīng)值插算出來。選取節(jié)點(diǎn)應(yīng)從小到大順序選取，曲線的兩個(gè)端點(diǎn)必須選入，其間可隨曲線曲率變化而定，端點(diǎn)附近與曲率變化較大處選密一些，反之可選稀一些，選取節(jié)點(diǎn)數(shù)量應(yīng)滿足計(jì)算精度的要

11、求且等差為佳。節(jié)點(diǎn)不得少于3個(gè)。 存在的問題：由于項(xiàng)次的有限性，造成了接點(diǎn)前后公式的非連續(xù)性，使接點(diǎn)插值不同一，節(jié)點(diǎn)間使用的插值公式不同一，給插值帶來一定的困難。/x,y是線性點(diǎn)的系列值，t是要插入的點(diǎn)（如水文或流量）static double interp3(double x, double y, double t) double v=0; int p=0; int n=x.length;for(int i=1;i < n-2;i+) if(t <= xi) break; /查找到接近X的值p+; if(p>0 && p<n-2 &&

12、(t-xp<xp+1-t) / p-; v = (t-xp+1)/(xp-xp+1)*(t-xp+2)/(xp-xp+2)*yp; v += (t-xp)/(xp+1-xp)*(t-xp+2)/(xp+1-xp+2)*yp+1; v += (t-xp)/(xp+2-xp)*(t-xp+1)/(xp+2-xp+1)*yp+2; return v; 2.4 一元多點(diǎn)在一元三點(diǎn)插值中，由于項(xiàng)次的有限性，造成了節(jié)點(diǎn)前后公式的非連續(xù)性，使節(jié)點(diǎn)插值不統(tǒng)一，節(jié)點(diǎn)間使用的插值公式不統(tǒng)一，給插值帶來一定的困難。因此通過一元多點(diǎn)（改進(jìn)型一元三點(diǎn)插值）有效解決該問題，是曲線更加連續(xù)。對(duì)于節(jié)點(diǎn)間的插值，同時(shí)采用

13、兩邊的公式計(jì)算，即使用所需要的插值的區(qū)間兩邊的節(jié)點(diǎn)，再添加更前面的一個(gè)節(jié)點(diǎn)擬合一條二次多項(xiàng)式，然后同樣方法對(duì)添加了后面的一個(gè)節(jié)點(diǎn)也擬合一條二次多項(xiàng)式，分別用這兩個(gè)公式計(jì)算插算值，然后平均，這樣區(qū)間的插值就統(tǒng)一了，不會(huì)出現(xiàn)跳躍點(diǎn)，而且可以減小單一公式帶來的誤差。同樣這樣的修正方法適用于許多水文水力計(jì)算中的一元三點(diǎn)插值公式，避免節(jié)點(diǎn)前后公式的非連續(xù)性。/x,y是線性點(diǎn)的系列值，t是要插入的點(diǎn)（如水文或流量）static double improveInterp3(double x, double y,double t) double v0=0,v=0,v1=0; int p=0,p0=0,p1=0

14、; int n=x.length;for(int i=1;i < n-2;i+) if(t <= xi) break; /查找到接近X的值p+; if(p>0 && p<n-2 && (t-xp<xp+1-t) p-; v = (t-xp+1)/(xp-xp+1)*(t-xp+2)/(xp-xp+2)*yp; v += (t-xp)/(xp+1-xp)*(t-xp+2)/(xp+1-xp+2)*yp+1; v += (t-xp)/(xp+2-xp)*(t-xp+1)/(xp+2-xp+1)*yp+2; p0=p-1;if(p0<

15、;0) p0=0;v0 = (t-xp0+1)/(xp0-xp0+1)*(t-xp0+2)/(xp0-xp0+2)*yp0; v0 += (t-xp0)/(xp0+1-xp0)*(t-xp0+2)/(xp0+1-xp0+2)*yp0+1; v0 += (t-xp0)/(xp0+2-xp0)*(t-xp0+1)/(xp0+2-xp0+1)*yp0+2;p1=p+1;if(p1>=n-3) p1=n-3;v1 = (t-xp1+1)/(xp1-xp1+1)*(t-xp1+2)/(xp1-xp1+2)*yp1; v1 += (t-xp1)/(xp1+1-xp1)*(t-xp1+2)/(xp1+

16、1-xp1+2)*yp1+1; v1 += (t-xp1)/(xp1+2-xp1)*(t-xp1+1)/(xp1+2-xp1+1)*yp1+2;return (v0+v1+v1)/3; 3 插值算法應(yīng)用示例2.1 在水位流量曲線的應(yīng)用在洪水預(yù)報(bào)中，通過水位推流量或者通過預(yù)報(bào)的流量推水位被經(jīng)常用到。實(shí)際過程中，一般通過實(shí)測的水文和流量特征值制作水文流量關(guān)系曲線，在實(shí)際應(yīng)用過程中，水位流量受到下游洪水頂托、泥沙等多種因素影響，實(shí)際上是變化的，因此采用近似的綜合線或者多條線代表水位關(guān)系曲線，如下圖所示：水文關(guān)系曲線通過設(shè)定相關(guān)的插值方式，插值方式包括（無插值、直線插值、一元三點(diǎn)、一元多點(diǎn)），拖轉(zhuǎn)曲線

17、上的相關(guān)點(diǎn)，曲線及影響到的表格數(shù)據(jù)聯(lián)動(dòng)。拖動(dòng)效果如下圖所示：利用插值點(diǎn)拖動(dòng)效果實(shí)際在預(yù)報(bào)、預(yù)警過程中，采用上面的算法進(jìn)行計(jì)算，一般采用直線插值計(jì)算水位流量關(guān)系。2.2 在API預(yù)報(bào)模型的應(yīng)用PPaR模型是三參數(shù)相關(guān)模型， API模型以流域降雨產(chǎn)流的物理機(jī)理為基礎(chǔ),以主要影響因素作參變量, 建立降雨量P與產(chǎn)流量R 之間定量的相關(guān)關(guān)系。常用的參變數(shù)有前期影響雨量Pa( 反映前期土濕) 、季節(jié)( 或用月份、周次, 反映洪水發(fā)生時(shí)間) 和降雨歷時(shí)( 或降雨強(qiáng)度) 等, 也有采用雨型、暴雨中心位置等因素。單位過程線一般簡稱單位線, 是謝爾曼( L.K.Sherman) 于1932 年提出的, 由地面徑流

18、和壤中流( 總稱直接徑流) 形成單位線。該單位線在我國洪水預(yù)報(bào)的匯流計(jì)算中應(yīng)用廣泛, 效果較好。P+PaR兩參數(shù)相關(guān)模型在PPaR三參數(shù)相關(guān)模型基礎(chǔ)上將P+Pa的和作為降雨徑流相關(guān)圖的縱坐標(biāo)，以R為橫坐標(biāo)建立其相關(guān)關(guān)系。使用時(shí)，首先計(jì)算洪水起漲時(shí)的土壤含水量Pa值，再把時(shí)段雨量序列變成累積雨量序列，用累積雨量查出累積凈雨，最后由累積凈雨轉(zhuǎn)化成時(shí)段凈雨量序列。圖附1-1 P+Pa-R關(guān)系曲線示意圖Pa由前期雨量計(jì)算，也稱前期影響雨量，是反映土壤濕度的參數(shù)。其計(jì)算公式為：當(dāng)上時(shí)段有雨量，即時(shí) (1)當(dāng)上時(shí)段無雨時(shí) (2)式中 K土壤含水量衰減系數(shù)；Pa,t、Pa,t1本時(shí)段和上時(shí)段的前期影響雨量；Pt1上時(shí)段降雨量。用Im表示土壤最大初損量，以mm計(jì)。通常Im60100mm。當(dāng)計(jì)算的Pa，t>Im時(shí)，則以Im作Pa值計(jì)算，即認(rèn)為此后的降雨不再補(bǔ)充初損量，全部形成徑流。馬良坪API方案馬良坪控制面積2.3 在相關(guān)圖預(yù)報(bào)模型的應(yīng)用根據(jù)河段上、下游斷