孫欣躍

（遼寧省水資源管理集團(tuán)有限責(zé)任公司，沈陽110003）

分形插值是基于分形基礎(chǔ)理論的具體應(yīng)用研究，能很好地描述生活中廣泛存在的對象，往往具有復(fù)雜性和粗糙性。除此之外，分形插值函數(shù)靈活性強(qiáng)和穩(wěn)定性高的特點使其廣泛應(yīng)用于水文地質(zhì)、 動畫仿真、圖像處理等眾多研究領(lǐng)域［1］。

1 研究區(qū)概況

濟(jì)寧市山東省西南部的一個地級市， 濟(jì)寧市氣候四季分明，氣溫宜人。受熱帶海洋氣團(tuán)或熱帶海洋氣團(tuán)影響， 天氣炎熱多雨， 年平均降水量597～820mm；冬季多為偏北風(fēng)，受極地大陸氣團(tuán)影響，多為晴冷天氣；春秋兩季為大氣環(huán)流調(diào)整期，春季易干旱多風(fēng)，變化迅速溫暖；秋天涼爽，有時下雨。濟(jì)寧市的水資源較為充沛，年降水量大，對于防洪減災(zāi)工作難度較大。

2 計算方法和數(shù)據(jù)

2.1 計算原理

分形插值是對具有分形特征的曲線進(jìn)行的快速插值方法，其根據(jù)分形曲線的自相似性，基于迭代函數(shù)系統(tǒng)（IFS）實現(xiàn)。 首先要構(gòu)造一個迭代函數(shù)系，以此來表示采用分形的插值，令吸引子A為插值函數(shù)的圖形。對于{（xi，yi），i=0，1，···，N}∈R2，分形插值函數(shù)是 插 值 于{（xi，yi），i=0，1，···，N}的 一 個 連 續(xù) 函數(shù)f［x0，xn］→R，且其圖形就是雙曲迭代函數(shù)系統(tǒng)的吸引子A。

創(chuàng)建一個吸引子為插值函數(shù)f（x）的迭代函數(shù)系{R2：Wi，i=1，2···，N}，設(shè)函數(shù)系中的每個函數(shù)Wi是仿射變換，表示為［2］：

Wi∶R2→R2，i=1，2，…，N

2.2 降雨時間序列的提取

氣象數(shù)據(jù)由中國氣象信息中心官網(wǎng)下載， 數(shù)據(jù)類型為某氣象站點氣候資料的日值數(shù)據(jù)。 以數(shù)值形式的氣象資料為每日平均資料， 時降雨資料可基于ArcMap于時刻降雨資料雨量圖中提取， 逐小時雨量分布如圖1。

圖1 逐小時雨量分布

每幅氣象圖代表某一小時的氣象數(shù)據(jù)， 在某一具體氣象站點多幅圖的數(shù)據(jù)的提取就構(gòu)成了該點某時間區(qū)間的序列。 本文采用的數(shù)據(jù)是小時時間間隔的降雨數(shù)據(jù)。 下載2018年6月30日12：00到2018年7月5日12：00時刻雨量圖， 選取或創(chuàng)建指定氣象站點（濟(jì)寧為54326站點）的Shape 文件，提取該點（或多點）的時間序列數(shù)據(jù)。 其中點Shape數(shù)據(jù)為設(shè)定的氣象數(shù)據(jù)提取位置。通過Spitial Analyst 中的提取分析工具，提取圖像素數(shù)值，建立換算表，取圖例中的中間值，得到的換算如表1。

表1 降雨數(shù)據(jù)提取對應(yīng)

最終得到濟(jì)寧市（氣象站點58326）的120h連續(xù)降雨值，轉(zhuǎn)換為降雨強(qiáng)度單位，降雨時間序列如圖2。

圖2 120h降雨時間序列

3 研究結(jié)果

3.1 降雨數(shù)據(jù)分形特性分析

分形分析采取盒子覆蓋法計算曲線的盒維數(shù)，通過繪制單位長度和覆蓋所需的總正方形數(shù)量的雙對數(shù)曲線得到［3］。選取若干不同尺度的r，一般情況下選取離散點時間序列樣點周期的整數(shù)倍。 用邊長為r的正方形對曲線進(jìn)行覆蓋，計算不同尺度r覆蓋時間序列曲線所需單位正方形的總數(shù)N（r）。繪制ln N（r）-ln（1/r）的雙對數(shù)散點圖，采取最小二乘法估計擬合出散點斜率。 回歸直線與散點的擬合就反應(yīng)了降雨時間序列的自相似特性［4］。

根據(jù)分形維數(shù)的盒維數(shù)計算法， 采用MATLAB軟件引入這些公式，改變小正方形的邊長r，給定的單位長度r為兩個時間單位（r=2），單位長度為r=2i/N（i=1，2，…，5），統(tǒng)計正方形數(shù)目N（r）。

表2 r覆蓋統(tǒng)計

將表中的ln（1/ri）和ln（N（ri））數(shù)據(jù)值，添加趨勢擬合線，可得到趨勢擬合的回歸直線。得到的回歸方程為：y=1.012x+6.069。

圖3 時間序列曲線盒覆蓋雙對數(shù)

根據(jù)分形盒維數(shù)的定義與原理， 得知該擬合直線的斜率D即為時間序列曲線的盒維數(shù)估計數(shù)值。所以維數(shù)為1.012。 雙對數(shù)曲線的擬合相關(guān)系數(shù)R2為0.998。 從擬合直線的直觀表象來看， 對數(shù)點分布在直線的附近，擬合相關(guān)率高達(dá)99.8%，說明此次降雨過程的時間序列曲線具有很強(qiáng)的局部與整體的相似性，定量反映了該時間序列曲線具有分形特性。

3.2 降雨時間序列的分形插值計算

由于對缺少實測降雨數(shù)據(jù)， 中國氣象數(shù)據(jù)網(wǎng)站提供的地面逐小時降雨資料， 所以可依據(jù)時降雨時間序列的分形特性， 采用分形插值函數(shù)計算提取得到120個小時每兩個時刻之間的降雨量。

對時刻進(jìn)行排序，即構(gòu)造120個插值點為p0（x0，y0）=（0，0.10）、p1（x1，y1）=（1，0.10）、p2（x2，y2）=（2，0.10）···p119（x119，y119）=（119，0.10）。 上述120個點是實現(xiàn)仿射分形插值函數(shù)構(gòu)造的所需插值點，根據(jù)垂直比例因子非確定的取值方法，可初步得到各個仿射系數(shù)［5］。

以上計算過程可基于MATLAB軟件程序文件實現(xiàn)，對已知初始點進(jìn)行迭代，根據(jù)分形插值函數(shù)的計算規(guī)律，需要進(jìn)行一次迭代過程，即在原兩個數(shù)據(jù)點之間，插入118個數(shù)據(jù)點，可將原來以小時為單位的降雨時間序列細(xì)化到30s，得到插值后曲線，如圖4。

圖4 分形插值后降雨時間序列

插值得到細(xì)化時間序列將插值分形曲線圖與實測曲線圖進(jìn)行對比可看出，曲線走勢相同，變化規(guī)律基本一致。 表3是插值數(shù)據(jù)和實測數(shù)據(jù)的平均值、方差和分形維數(shù)的比較。

表3 原始數(shù)據(jù)與插值數(shù)據(jù)比較

式中 m為插值數(shù)據(jù)；x為實測數(shù)據(jù)；r為相對誤差。

3.3 分形插值與傳統(tǒng)插值方法結(jié)果對比

與傳統(tǒng)插值方法進(jìn)行對比， 可看出分形插值在插值效果上具有明顯優(yōu)越性， 為120個數(shù)據(jù)線性插值的結(jié)果。 線性插值沒有考慮已知連點間的具體變化過程，將兩點直接用直線進(jìn)行連接，而分形插值基于部分與整體的相似性， 將兩點間變化的過程具體化，反映了局部特征。拉格朗日插值也屬于一種常用的傳統(tǒng)插值方法，又稱為多項式插值。該插值方法的缺點在于插值節(jié)點個數(shù)的增加的同時， 不能確保兩個插值節(jié)點之間插值函數(shù)一定能夠很好地與被插值函數(shù)吻合。從理論角度來看，冪次數(shù)多項式需要的數(shù)據(jù)多，擬合結(jié)果比低冪次準(zhǔn)確，但從舍入誤差看，高次插值計算過程復(fù)雜，而且計算量較大，計算過程中會產(chǎn)生較大程度的誤差累積。擬合過程中出現(xiàn)了“龍格現(xiàn)象”， 使得靠近區(qū)間端點處的擬合數(shù)值很大，曲線的變化速度較高，出現(xiàn)了較大誤差。而分形插值法的誤差值較小，所得結(jié)果較為精準(zhǔn)如圖5。

圖5 線性插值后降雨時間序列

4 結(jié)語

本章從分形理論與分形的基本性質(zhì)出發(fā)， 對濟(jì)寧市2018年6月30日12：00 至2018年7月5日12：00連續(xù)120h的降雨時間序列進(jìn)行了分形性分析， 得到該降雨時間序列曲線的分形盒維數(shù)為1.012，雙對數(shù)曲線擬合度為99.8%，說明了該曲線具有很好的分形特性。進(jìn)而應(yīng)用分形插值理論，對該分形特性曲線運用迭代系統(tǒng)進(jìn)行分形插值， 根據(jù)最終所需時間序列的精確度， 確定進(jìn)行一次迭代， 實現(xiàn)每兩個小時之間118 個插值點的生成， 最終將以小時為單位的時間序列細(xì)化處理為30s為單位的新時間序列。 將源數(shù)據(jù)與插值后數(shù)據(jù)進(jìn)行誤差統(tǒng)計，平均值、均方差和分形維數(shù)的誤差分別為0.007，0.066，0.18，誤差均在可接受范圍內(nèi)， 可滿足當(dāng)前中小型城市的城市雨洪建模的降雨數(shù)據(jù)需求。