基于正交試驗(yàn)選優(yōu)的湖底地形分區(qū)插值方法

龔 懿，陳再揚(yáng)，張 維，汪 靚，王玉琳，楊小玲

(1.揚(yáng)州大學(xué) 水利科學(xué)與工程學(xué)院，江蘇 揚(yáng)州 225009； 2.江蘇省水利工程建設(shè)局，江蘇 南京 210029； 3.揚(yáng)州大學(xué) 環(huán)境科學(xué)與工程學(xué)院，江蘇 揚(yáng)州 225009； 4.重慶市農(nóng)業(yè)科學(xué)院，重慶 401329)

水文水資源、水環(huán)境數(shù)值模擬等研究需要高精度的河流或湖泊地形數(shù)字信息[1-2]。建立數(shù)字高程模型較為成熟的方法是將現(xiàn)有地形圖數(shù)字化通過內(nèi)插方法生成DEM（Digital Terrain Model），國內(nèi)常利用多種高程信息（等高線、高程點(diǎn)、骨架點(diǎn)）通過插值算法生成DEM[3]。湖底高程信息測量的主要手段有網(wǎng)絡(luò)RTK 聯(lián)合聲波探測儀[4]、水下測量機(jī)器人[5]等。然而，由于研究區(qū)域地形地貌、測量手段等多方面因素限制，湖底設(shè)置觀測點(diǎn)的數(shù)量有限，其布局也無法完全滿足研究要求[6-9]，利用觀測點(diǎn)高程數(shù)據(jù)進(jìn)行空間插值展布，成為獲得高精度湖底地形數(shù)據(jù)的主要手段之一?？臻g地形插值的方法有很多，最為常見的是克里金插值[10]和反距離加權(quán)插值[11]，其中反距離加權(quán)插值法采用插值點(diǎn)到所有空間已知高程點(diǎn)的反距離加權(quán)平均來控制插值，計(jì)算簡單，應(yīng)用最為廣泛[12-14]。很多學(xué)者將反距離加權(quán)插值法引入水下地形研究，如華祖林等[15]利用對比三角網(wǎng)格法、反距離加權(quán)插值法和克里金法，分析河道強(qiáng)彎段、干支交匯處等地形復(fù)雜區(qū)域的插值效果，并將獲得的黃浦江全段插值地形用于潮位計(jì)算，結(jié)果表明克里金法對黃浦江全段的適應(yīng)性最好、精度最高；蔣偉達(dá)等[16]研究埕島海域地形演變，利用多種插值方法進(jìn)行適應(yīng)性比較，發(fā)現(xiàn)反距離加權(quán)插值法能夠很好地生成埕島海域范圍內(nèi)的水深地形；李冕等[17]基于自適應(yīng)凸包選點(diǎn)改進(jìn)反距離加權(quán)插值法，并將其應(yīng)用于2011年的渤海無機(jī)氮評價，相比傳統(tǒng)的反距離加權(quán)插值法（Inverse Distance Weight，IDW），改進(jìn)的插值法能夠更好地反映海水污染物空間分布狀況；秦俊桃等[18]比較了反距離加權(quán)插值法等8 種空間插值方法在地下水位插值中的應(yīng)用效果，認(rèn)為反距離加權(quán)插值法較適用于西北干旱地區(qū)的地下水位插值。

然而在以往對反距離加權(quán)插值法的研究中，主要涉及利用反距離加權(quán)插值法解決實(shí)際問題，對反距離加權(quán)插值算法本身的研究相對匱乏，有關(guān)距離加權(quán)插值的長度冪值對插值影響的研究較少。Li 等[19]在反距離加權(quán)插值法公式中添加系數(shù)K以調(diào)整距離采樣點(diǎn)的權(quán)重，提高了插值精度；李翔宇等[20]研究了全國不同密度的近海岸反距離加權(quán)插值法的參數(shù)選取，近海岸插值的加權(quán)冪值選2 為最優(yōu)；張錦明等[21]研究了插值參數(shù)對插值精度的影響，認(rèn)為權(quán)指數(shù)對插值精度影響最為顯著。與普通的反距離加權(quán)插值方法相比，這些方法雖然都考慮了權(quán)指數(shù)對插值精度的影響，但對地形變化較大的湖泊，如整個湖泊給定同一冪值，存在平滑效應(yīng)影響插值效果的問題。

由于反距離加權(quán)插值法通過長度冪值來控制插值周邊采樣點(diǎn)對插值點(diǎn)的影響，地形起伏較大的湖泊冪值的選取對模擬精度會產(chǎn)生較大影響。本文提出對采樣點(diǎn)進(jìn)行方差計(jì)算，分析每個采樣點(diǎn)與其周圍采樣點(diǎn)高程的變異程度；再利用K-means 聚類，對采樣點(diǎn)以方差為指標(biāo)進(jìn)行聚類分析；對不同類別點(diǎn)的冪值采用正交試驗(yàn)表進(jìn)行選優(yōu)，以高程-庫容曲線的庫容值為衡量指標(biāo)選取最優(yōu)冪值，最終得到湖底地形形態(tài)。

1 研究方法

1.1 普通反距離加權(quán)插值法

反距離加權(quán)插值法也稱為反距離倒數(shù)乘方法，每個高程數(shù)據(jù)已知點(diǎn)對插值點(diǎn)都有影響，即為權(quán)重，反距離加權(quán)插值法主要依賴反距離的權(quán)重來控制插值。采樣點(diǎn)距離插值點(diǎn)越近將被賦予越大的權(quán)重；反之，距離越遠(yuǎn)，其權(quán)重越小[22]。將插值函數(shù)Z(x，y)定義為各已知高程點(diǎn)Z*(xi，yi)的加權(quán)平均值，設(shè)有n個已知高程數(shù)據(jù)點(diǎn)，平面坐標(biāo)為(xi，yi)，其計(jì)算式為：

式中：wi為每個已知高程點(diǎn)的長度權(quán)重；p為長度的冪值，為任意正實(shí)數(shù)；hi為第i個已知高程點(diǎn)到插值點(diǎn)間的距離（m）；hj為第j個已知高程點(diǎn)到插值點(diǎn)間的距離（m）；n為已知高程點(diǎn)的數(shù)量。

普通的IDW 插值利用采樣點(diǎn)距離插值點(diǎn)的長度確定權(quán)重，雖然每個采樣點(diǎn)都定義了冪值以考慮所有采樣點(diǎn)的影響，但僅考慮了距離因素未考慮采樣點(diǎn)的變異程度對插值的影響。因此，在地形起伏較大的區(qū)域應(yīng)用有一定的局限性。

1.2 分區(qū)域反距離加權(quán)插值方法構(gòu)建

1.2.1 地形分區(qū) 方差在統(tǒng)計(jì)學(xué)中是衡量一組數(shù)據(jù)離散程度的重要度量，選取方差作為衡量采樣點(diǎn)數(shù)據(jù)起伏程度的指標(biāo)，計(jì)算每個點(diǎn)與其周圍z個采樣點(diǎn)的高程值之間的方差。本文采用K-means 聚類算法對衡量數(shù)據(jù)起伏程度的方差進(jìn)行分析，對起伏程度相似的采樣點(diǎn)賦相同的長度冪值。選取歐式距離作為相似的測度，每個采樣點(diǎn)的方差作為聚類樣本，K-means 算法的具體步驟如下：

（1）在n個高程方差樣本中選取k個樣本作為初始聚類中心ci(i=1, 2,…,k)。

（2）計(jì)算除聚類中心外其他高程方差樣本pj到聚類中心的歐式距離：

式中：j為高程采樣點(diǎn)中除聚類中心外的其他樣本點(diǎn)，j=1, 2,…,n?1。

（3）根據(jù)聚類中心構(gòu)造k個集合，每個樣本依據(jù)與聚類中心的歐式距離，劃歸到與之最近的集合中。

（4）采用均值更新聚類中心：

式中：ni為第i個集合的樣本個數(shù)；pir為第i個聚類中第r個采樣點(diǎn)到聚類中心的距離（m）。

（5）重復(fù)步驟（2）～（4），直到滿足精度要求，聚類中心不發(fā)生改變。

1.2.2 冪值選優(yōu) 由于區(qū)域的冪值選擇數(shù)目較多，若將所有組合遍歷選優(yōu)，則工作量過大。正交試驗(yàn)方法遵循“均衡搭配”的原則，通過盡可能少的試驗(yàn)次數(shù)求解出最優(yōu)方案。因此，本文采用正交試驗(yàn)選優(yōu)的方法，通過構(gòu)造正交試驗(yàn)表選取一部分組合試驗(yàn)，然后根據(jù)試驗(yàn)值推求得到全部組合的理論最優(yōu)解。

在正交試驗(yàn)表Lp(tq)中，t為水平數(shù)，即冪值n在可行域內(nèi)的離散個數(shù)；q為正交表最多安排因素的個數(shù)，即地形劃分的區(qū)域數(shù)；p為正交表的總方案數(shù)。其中p，t，q之間的關(guān)系[23]如下：

式中：v為任意整數(shù)；Y為正交表實(shí)際因素的個數(shù)。

1.3 分區(qū)域反距離加權(quán)插值方法實(shí)現(xiàn)

與傳統(tǒng)的插值法相比，分區(qū)域反距離加權(quán)插值法據(jù)所屬類別賦予每個插值點(diǎn)特定冪值，但在實(shí)際插值中，劃分插值點(diǎn)所屬類別較為困難。因此，本文采用K 近臨算法來劃分插值點(diǎn)的類別：如果一個樣本在特征空間中m個最相似（即距離樣本距離最近的m個點(diǎn)）的樣本都屬于同一類別，則該樣本也屬于這個類別。將高程采樣點(diǎn)按照聚類分析構(gòu)建k個數(shù)據(jù)集：

式中：T1、T2、…、Tk為k個區(qū)域構(gòu)建的高程采樣點(diǎn)集合；pk為每個區(qū)域給定的冪值；xkc、ykc、zkc分別為第k個數(shù)據(jù)集中第c個采樣點(diǎn)的坐標(biāo)。

選用歐式距離作為計(jì)算插值點(diǎn)與各個數(shù)據(jù)集中采樣點(diǎn)之間的距離，按照距離遞增關(guān)系進(jìn)行排序，選取距離最小的M個點(diǎn)；計(jì)算前M個點(diǎn)在各數(shù)據(jù)集出現(xiàn)的頻率，選取出現(xiàn)頻率最高的數(shù)據(jù)集作為插值點(diǎn)的所屬區(qū)域，即選用該數(shù)據(jù)集的長度冪值作為該插值點(diǎn)的長度冪值。

2 研究區(qū)域與數(shù)據(jù)

選取太湖湖底高程進(jìn)行模擬插值。太湖地處長江三角洲的南部，是中國第三大淡水湖，橫跨江、浙兩省，北臨無錫，南靠湖州，西近宜興，東依蘇州。太湖整體呈西高東低，平均年出湖流量約75 億m3，蓄水量約44 億m3。太湖湖泊面積2 427.8 km2，水域面積為2 338.1 km2。太湖湖底高程一般為1.0 m，中東部為洼地，其地面高程為3.0～4.5 m，最低處僅2.5 m，其他區(qū)域地面高程為5.0～8.0 m。

2.1 數(shù)據(jù)基礎(chǔ)

試驗(yàn)數(shù)據(jù)來源于太湖湖底的實(shí)測數(shù)據(jù)，實(shí)測高程點(diǎn)有6 242 個，區(qū)域范圍為整個湖區(qū)，離散點(diǎn)空間密度為2.57 個/km2，測點(diǎn)的平均高程為?0.64 m。實(shí)測點(diǎn)空間分布見圖1，其數(shù)據(jù)格式為“X-Y”。由圖1 可以看出，測點(diǎn)主要分布在湖岸、人工島周圍和湖心區(qū)，其中湖的東側(cè)、東南側(cè)實(shí)測點(diǎn)分布較為集中。

圖1 實(shí)測點(diǎn)空間分布Fig.1 Spatial distribution of measured points

2.2 分區(qū)域反距離加權(quán)插值方法應(yīng)用

2.2.1 太湖采樣點(diǎn)分區(qū) 太湖水面面積較大，湖底高程采樣點(diǎn)較為分散。為更好衡量湖泊地形采樣點(diǎn)的局部變異程度，選取每個采樣點(diǎn)與其周圍最近的10 個采樣點(diǎn)為對象，進(jìn)行方差分析（見圖2）。

從圖2 可以看出，湖心區(qū)域地形趨于平緩，靠近湖岸區(qū)域有較小地形起伏，太湖的西北側(cè)、湖心島等區(qū)域地形起伏較大，這符合太湖地形的實(shí)際情況，因?yàn)樘膮^(qū)域分布多座小型島嶼。根據(jù)太湖的湖底地形起伏情況，將湖區(qū)分成3 個區(qū)域：湖心平緩區(qū)、湖心島周邊區(qū)、湖岸區(qū)，故將K-means 聚類的初始聚類中心設(shè)置為3。

圖2 太湖地形采樣點(diǎn)局部起伏分布Fig.2 Local undulation distribution map of Taihu Lake terrain sampling points

K-means 聚類將地形采樣點(diǎn)方差數(shù)據(jù)劃分為3 個區(qū)間，即[0，0.20）、[0.20，0.45）、[0.45，1.00]。

2.2.2 太湖各區(qū)域冪值選優(yōu) 試驗(yàn)因素、試驗(yàn)水平、試驗(yàn)指標(biāo)的確定與正交表的構(gòu)建如下。

（1）試驗(yàn)因素： K-means 聚類依據(jù)局部地形起伏程度將太湖湖底劃分為3 個區(qū)域，選取劃分的區(qū)域作為試驗(yàn)因素，試驗(yàn)因素?cái)?shù)量為3。

（2）試驗(yàn)水平： 選取反距離插值法的冪值作為試驗(yàn)水平。為了試驗(yàn)選優(yōu)的方便性，試驗(yàn)水平的可行區(qū)間取[1，4]，但考慮到各區(qū)域連接處的平順插值，各區(qū)域冪值不能相差太大，所以將冪值可行區(qū)間劃分為3 段：1～2、2～3、3～4，各可行區(qū)間步長取0.2，冪值水平數(shù)為5。

（3）試驗(yàn)指標(biāo)： 選取實(shí)測高程為正交試驗(yàn)的指標(biāo)。

（4）正交表選擇： 3 因素5 水平，全部組合為3×53=375，將全部方案代入分區(qū)域反距離加權(quán)插值法生成高程工作量較大。為此，采用部分試驗(yàn)選優(yōu)的方法，構(gòu)建3 個3 因素5 水平的正交試驗(yàn)表L25(53)，按正交表選擇75 種試驗(yàn)組合方案，即可獲得對應(yīng)全部試驗(yàn)組合的理論最優(yōu)區(qū)域冪值方案。試驗(yàn)組合方案詳見圖3。

本文采用留一法進(jìn)行最優(yōu)方案的確定，選取10%的實(shí)測高程數(shù)據(jù)作為驗(yàn)證點(diǎn)，其高程點(diǎn)的平均高程值為?0.635 m，然后利用規(guī)則網(wǎng)格模型內(nèi)插得出對應(yīng)點(diǎn)的高程并與之比較。采用與實(shí)測高程點(diǎn)的均方根誤差來表征空間插值結(jié)果的相對誤差，將其作為試驗(yàn)的最后評定指標(biāo)，分析結(jié)果見圖3。

圖3 L25(53)正交表區(qū)域冪值選優(yōu)方案組合Fig.3 L25(53) orthogonal table region power selection scheme combination

通過正交分析，可以獲得全部試驗(yàn)組合的理論最優(yōu)解。在冪值范圍1～2 中，各方案與實(shí)測高程點(diǎn)的平均均方根誤差為0.273 m，方案21 為最優(yōu)，與實(shí)測高程點(diǎn)的均方根誤差為0.198 m。在冪值范圍2～3 中，各方案與實(shí)測高程點(diǎn)的平均均方根誤差為0.156 m，方案21 最優(yōu)，均方根誤差為0.145 m。在冪值范圍3～4 中，各方案與實(shí)測高程點(diǎn)的平均均方根誤差為0.146 m，方案22 最優(yōu)，均方根誤差為0.144 m。通過各正交表的誤差分析，可得區(qū)域1 冪值取4.0、區(qū)域2 冪值取3.4、區(qū)域3 冪值取3.2 時，空間插值得到的高程值與實(shí)測高程值均方根誤差最小，空間插值結(jié)果最優(yōu)。

3 結(jié)果分析與討論

為考量分區(qū)域反距離加權(quán)插值法在地形插值中的效果，將本文的分區(qū)域反距離加權(quán)插值法與經(jīng)典的IDW 法、普通克里金方法、自然鄰域法對比，利用上述留一法對空間插值結(jié)果與實(shí)測高程值進(jìn)行驗(yàn)證。根據(jù)《多尺度數(shù)字高程模型生產(chǎn)技術(shù)規(guī)定》[24]，平地、山地、丘陵地和高山地的2 m 網(wǎng)格DEM 點(diǎn)位中誤差限分別為0.7、1.7、3.3 和6.7 m，表1 中均方根誤差數(shù)據(jù)表明分區(qū)反距離加權(quán)法、反距離加權(quán)法、克里金法和自然鄰域法不同地形區(qū)中誤差皆在允許范圍之內(nèi)。同時，由表1 可以看出，克里金法與反距離加權(quán)法精度相近，自然鄰域法的插值精度與克里金法、反距離加權(quán)法相比明顯提高，本文提出的分區(qū)域反距離加權(quán)法精度略高于自然鄰域法。

表1 4 種空間插值方法的誤差分析結(jié)果（高程）Tab.1 Error analysis results of four spatial interpolation methods (elevation)

以太湖湖底高程為插值對象，對比反距離加權(quán)法、克里金法、分區(qū)域反距離加權(quán)法對太湖地形空間展布結(jié)果：在湖心與岸邊交界區(qū)域、湖心島附近，反距離加權(quán)法和克里金法的空間展布都有明顯的梯度，出現(xiàn)了由湖岸向湖心高程陡降的異常突變，反距離加權(quán)法的突變更加明顯（見圖4），與實(shí)際情況不符。分區(qū)域的反距離加權(quán)法較好地反映了插值結(jié)果受采樣點(diǎn)起伏程度而表現(xiàn)出的空間變異性，能夠正確處理采樣點(diǎn)變異度較大而引起的數(shù)據(jù)異常，更加接近太湖湖底高程分布的實(shí)際情況。

圖4 太湖湖底地形空間展布Fig.4 Spatial distribution of topography at the bottom of Taihu Lake

為衡量分區(qū)域反距離加權(quán)法在湖底地形模擬中的整體效果，提取不同高程對應(yīng)的庫容值與實(shí)測庫容值[25]進(jìn)行誤差分析（表2）。4 種插值方法得到的太湖湖泊庫容與實(shí)測值都存在一定差異，自然鄰域法插值提取的庫容數(shù)據(jù)與實(shí)測數(shù)據(jù)的偏差最大，均方根誤差達(dá)到2.72 億m3，反距離加權(quán)法與克里金法的插值均方根誤差相近，均方根誤差分別為2.47 億 m3和2.48 億m3，分區(qū)域的反距離加權(quán)法誤差最小，其插值效果明顯優(yōu)于反距離加權(quán)法、克里金法和自然鄰域法。

表2 多種方法插值結(jié)果誤差分析（庫容）Tab.2 Error analysis of interpolation results by various methods (storage capacity)

4 結(jié) 語

通過K-means 聚類將太湖湖底進(jìn)行區(qū)域劃分，針對深水湖泊湖底高程起伏較大的特點(diǎn)，以插值后地形數(shù)據(jù)提取的庫容曲線與實(shí)測庫容曲線的均方根誤差為指標(biāo)，提出對湖底分區(qū)域給定冪值，并通過正交試驗(yàn)選優(yōu)的優(yōu)化方法，充分考慮了湖底高程采樣點(diǎn)的變異程度，優(yōu)化確定各區(qū)域冪值，這對提高插值精度有重要意義。

在其他條件相同的情況下，對比反距離加權(quán)法、克里金法和自然鄰域法，分區(qū)域反距離加權(quán)插值法能夠更好適應(yīng)地形起伏，插值效果良好，可以獲得比較精確的湖泊地形數(shù)據(jù)，對湖泊研究和管理具有積極作用。