陳天偉,康傳利,陳 明,段青達(dá),韋海福

(桂林理工大學(xué) 廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室,廣西 桂林 541006)

0 引 言

格網(wǎng)DEM作為國家基礎(chǔ)地理信息基礎(chǔ), 得到了廣泛的應(yīng)用, 具有極其重要的意義。 然而其精度仍有進(jìn)一步提高的需求, 特別是在數(shù)字地形分析領(lǐng)域。 一般, 規(guī)則格網(wǎng) DEM由數(shù)學(xué)插值方法得到。 插值方法有多種, 包括多項(xiàng)式插值、 多面函數(shù)插值、 地統(tǒng)計(jì)插值、 加權(quán)平均插值等算法[1]。 由插值核函數(shù)解算插值權(quán)系數(shù), 以求得待插點(diǎn)高程值, 其高程誤差來源于原始數(shù)據(jù)誤差和模型逼近誤差兩方面, 插值模型與實(shí)際地表符合的程度決定了逼近誤差的大小。 如地學(xué)界一般認(rèn)為, 地統(tǒng)計(jì)學(xué)插值方法是在最優(yōu)線性無偏估計(jì)的準(zhǔn)則下, 對(duì)未知點(diǎn)空間屬性及插值精度進(jìn)行估計(jì)[2-3], 在平穩(wěn)假設(shè)的前提下, 理論上可以得到很高的插值精度; 然而, 地形表面的變化時(shí)常不符合平穩(wěn)假設(shè), 此時(shí)構(gòu)建的插值模型與實(shí)際不符就會(huì)顯著降低插值精度。 如由地統(tǒng)計(jì)學(xué)的楊赤中插值法求解參考點(diǎn)的權(quán)系數(shù)時(shí), 可能會(huì)出現(xiàn)負(fù)權(quán)系數(shù), 導(dǎo)致估值結(jié)果的不合理, 引起DEM建模的失真甚至錯(cuò)誤。 有國內(nèi)學(xué)者對(duì)此類負(fù)權(quán)現(xiàn)象展開研究, 認(rèn)為出現(xiàn)負(fù)值與參考點(diǎn)、 待估點(diǎn)的點(diǎn)位關(guān)系以及插值函數(shù)模型等因素相關(guān)[4]。 而筆者則認(rèn)為更主要的是由于傳統(tǒng)插值方法在構(gòu)建數(shù)學(xué)模型時(shí)存在主觀因素的干擾。 由于空間抽樣數(shù)據(jù)(如高程值)存在隨機(jī)性, 依據(jù)最大熵原理[5], 要使數(shù)學(xué)模型對(duì)樣本試驗(yàn)數(shù)據(jù)依賴程度最大, 那么熵值應(yīng)該最大, 使得插值模型更能反映空間數(shù)據(jù)的離散性, 則模型會(huì)最穩(wěn)健, 解算模型則可以避免負(fù)權(quán)系數(shù)的產(chǎn)生,筆者在文獻(xiàn)[6-7]中驗(yàn)證了上述觀點(diǎn)的正確性。 另外, 注意到地形離散點(diǎn)之間客觀存在著空間相關(guān)性, 在建模時(shí)顧及空間數(shù)據(jù)的相關(guān)特征, 應(yīng)可以使DEM模型更加符合客觀實(shí)際。

1 模型建立

1.1 最大熵模型的基本形式

最大熵方法是以空間域隨機(jī)變量的信息熵最大為準(zhǔn)則, 其數(shù)學(xué)模型包括目標(biāo)函數(shù)和約束條件。 目標(biāo)函數(shù)即熵值最大, 約束條件包括樣本點(diǎn)(參考點(diǎn))高程值的數(shù)學(xué)期望和方差等, 故建立最大熵模型[6-7]。

模型系統(tǒng)的熵值

S[p]=-(p1lnp1+p2lnp2+…+pnlnpn)/ln 2,

(1)

約束條件

(2)

其中:pi為空間域變量的產(chǎn)生概率(即權(quán)系數(shù));pi>0;i=1,2,…,n;S[p]為系統(tǒng)熵值。 約束條件第1式即0階矩約束(或稱自然矩約束); 第2式即1階矩約束: E1=E(h)是高程變量原點(diǎn)矩的期望值(隨機(jī)變量hi的原點(diǎn)矩等于相應(yīng)的樣本矩)[8], 則E(h)可利用參考點(diǎn)與估值點(diǎn)距離反比例賦權(quán), 對(duì)參與估值計(jì)算的參考點(diǎn)權(quán)系數(shù)取加權(quán)平均值算出; 第3式即2階矩約束, E2是高程變量2階中心矩的期望值, 利用高程樣本值和E1按計(jì)算方差估值的公式算出。

1.2 空間變量相關(guān)特征的表達(dá)

上述約束條件的本質(zhì)是通過期望值和方差反映研究區(qū)域高程值變量的分布中心以及離散度。根據(jù)地理學(xué)第一定律, 地貌在一定空間范圍內(nèi)具有相關(guān)性, 應(yīng)在約束條件中表達(dá)出這些特征。

在地統(tǒng)計(jì)學(xué)理論中,認(rèn)為高程屬于區(qū)域化變量,滿足2階平穩(wěn)假設(shè)情況下,在一定的范圍內(nèi),點(diǎn)和點(diǎn)之間是相關(guān)的[9],相關(guān)程度與距離成反比?？死锝鸱ú捎昧税胱儺惡瘮?shù)等指標(biāo)描述此相關(guān)特征,楊赤中法采用了綜合隨機(jī)特征函數(shù)進(jìn)行描述[10]。常用的反距離加權(quán)插值法亦是基于地理學(xué)第一定律的原理,即在一定區(qū)域,距離越近事物越相似。那么在解算DEM插值模型得到的各參考點(diǎn)權(quán)系數(shù)的大小,其實(shí)也就表達(dá)了插值點(diǎn)(待估點(diǎn))與其鄰域的參考點(diǎn)(已知點(diǎn))的相關(guān)程度，即權(quán)系數(shù)大小與參考點(diǎn)、估值點(diǎn)間距平方成反比。為了簡化最大熵模型的復(fù)雜度,在設(shè)置約束條件中僅僅增設(shè)最大最小權(quán)約束條件,即距離插值點(diǎn)最近的參考點(diǎn)權(quán)系數(shù)最大;反之,距離插值點(diǎn)最遠(yuǎn)的參考點(diǎn)權(quán)系數(shù)最小。

如圖1所示, 某待估點(diǎn)(最中心的點(diǎn))對(duì)應(yīng)的參考點(diǎn)有5個(gè),其中1#點(diǎn)距離待估點(diǎn)最遠(yuǎn),3#點(diǎn)最近,根據(jù)空間數(shù)據(jù)自相關(guān)的特性,1#點(diǎn)權(quán)系數(shù)應(yīng)該最小,3#點(diǎn)權(quán)系數(shù)應(yīng)該最大。

圖1 參考點(diǎn)、插值點(diǎn)的點(diǎn)位分布Fig.1 Distribution of reference points and interpolation point

相應(yīng)的空間相關(guān)約束條件： ①最小權(quán)約束條件p1-p2<0,p1-p3<0,p1-p4<0,p1-p5<0; ②最大權(quán)約束條件p2-p3<0,p4-p3<0,p5-p3<0。 轉(zhuǎn)換為矩陣表達(dá)式

(3)

(4)

此處,2#參考點(diǎn)在楊赤中法計(jì)算結(jié)果為負(fù)數(shù),這是因?yàn)樵擖c(diǎn)到待估點(diǎn)的方向被3#點(diǎn)屏蔽的緣故。2#、4#、5#參考點(diǎn)與待估點(diǎn)的距離介于最大、最小之間,本文不作類似的空間相關(guān)性約束,以免增加問題的復(fù)雜程度從而影響模型解算的效率。它們的權(quán)系數(shù)由以上模型自然解得?？梢?上述有關(guān)空間相關(guān)性的約束表達(dá)式相當(dāng)于給最大熵模型提供了最重要的空間相關(guān)特征信息,并進(jìn)一步縮小了可行解范圍,有利于加快模型解算收斂速度,求得最優(yōu)解。

待估點(diǎn)x0的估計(jì)量

(5)

2 模型解算

上述模型的解算本質(zhì)上屬于非線性函數(shù)約束優(yōu)化問題,即尋求全局最優(yōu)解的問題。所謂全局最優(yōu)解是指每一次解算某個(gè)未知點(diǎn)高程時(shí),先要解算搜索到的參考點(diǎn)的權(quán)系數(shù),這些權(quán)系數(shù)的解有很多組,要找到滿足約束條件且對(duì)應(yīng)熵值為最大的權(quán)系數(shù)解,進(jìn)而推算未知點(diǎn)高程。

由于目標(biāo)函數(shù)是復(fù)雜的非凸函數(shù),只能用數(shù)值計(jì)算方法進(jìn)行優(yōu)化。傳統(tǒng)的優(yōu)化算法往往求得的是局部最優(yōu)解,而不是全局最優(yōu)解,在此利用遺傳算法進(jìn)行解算,再結(jié)合最大最小權(quán)約束條件,能盡快鎖定全局最優(yōu)解。遺傳算法(GA)模仿遺傳學(xué)自然選擇進(jìn)化過程,主要是通過基因的組合交叉和變異逐步進(jìn)行問題的優(yōu)化,避免陷入局部最優(yōu)解,從而尋找到全局最優(yōu)解。該算法對(duì)于大空間、全局范圍內(nèi)求解最優(yōu)解具有較強(qiáng)的優(yōu)勢(shì)。其基本操作包括編碼、選擇初始種群、設(shè)置適應(yīng)度函數(shù)、遺傳算子操作(選擇、交叉和變異)、控制參數(shù)設(shè)定、循環(huán)終止設(shè)計(jì)[11]。GA自身通常解決無約束優(yōu)化問題,而本文建立的最大熵模型包括3個(gè)線性等式約束和n個(gè)非負(fù)約束,n即參考點(diǎn)個(gè)數(shù)。實(shí)踐中常用罰函數(shù)法,將有約束問題轉(zhuǎn)化為無約束問題。通過懲罰不可行解,使其目標(biāo)函數(shù)值增大(對(duì)求極小值問題),適應(yīng)度值減小,經(jīng)過不斷迭代,種群逐漸收斂于可行解,同時(shí)迫使無約束問題的極小點(diǎn)逐漸收斂于約束問題的極小點(diǎn)。構(gòu)造帶有懲罰項(xiàng)的適應(yīng)度函數(shù)一般有兩種: 一種是采用加法形式val(x)=f(x)+p(x), 另一種是乘法形式val(x)=f(x)·p(x)。 其中p(x)是懲罰項(xiàng)。

商業(yè)數(shù)學(xué)軟件MATLAB附帶有遺傳算法函數(shù)工具箱, MATLAB的GA主函數(shù)調(diào)用形式為[12][x, fval]=ga(ObjectiveFunction,nvars,A,b,Aeq,beq,LB,UB,ConstraintFunction,options),其輸入?yún)?shù)包括:目標(biāo)函數(shù)名、變量個(gè)數(shù)、線性不等式約束系數(shù)陣及常量、線性等式約束系數(shù)陣及常量、變量的上下界、非線性約束條件函數(shù)名、選項(xiàng)。

前述約束條件式(2)的前3項(xiàng)對(duì)應(yīng)于線性等式約束Aeq和beq,pi>0非負(fù)約束對(duì)應(yīng)于變量的上、 下界UB和LB, 而對(duì)于組成線性不等式約束系數(shù)陣A及常量b, 由空間相關(guān)約束條件式(3)、(4)得到。 一般地, 對(duì)于某個(gè)待估點(diǎn), 設(shè)參考點(diǎn)有n個(gè), 第i點(diǎn)距離待估點(diǎn)最遠(yuǎn)(權(quán)最小),第j點(diǎn)最近(權(quán)最大)。

組成線性不等式約束系數(shù)矩陣的方法是先列出最小權(quán)系數(shù)約束條件, 其不等式系數(shù)矩陣為: 在行列數(shù)均為n-1的負(fù)單位陣第i列之前插入元素為1的列向量, 即得線性不等式約束系數(shù)陣的第一部分A1, 即式(3)左邊的系數(shù)矩陣;再列出最大權(quán)系數(shù)約束條件,其不等式系數(shù)矩陣為:在行列數(shù)均為n-2的單位陣第i列之前插入元素為0的列向量(若i大于n-2,則在最后列右邊加入),再在新矩陣的第j列之前插入元素為-1的列向量(若j大于n-1,則在最后列右邊加入),即得線性不等式約束系數(shù)陣的第2部分A2,即式(4)左邊的系數(shù)矩陣;上下組合A1和A2即得A,常量b是2n-3個(gè)0元素組成的列向量,即式(3)、(4)右邊的0元素列向量。

應(yīng)注意, 若最大熵模型不存在非線性約束條件, 則在程序中指定非線性約束條件函數(shù)名為空值。

3 實(shí)例與分析

在某山地區(qū)域1∶1萬地形圖提取離散特征高程點(diǎn)60個(gè)[7],以選取某待估點(diǎn)作插值運(yùn)算為例,按極限半徑搜索該待估點(diǎn)四周的參考點(diǎn),得到5個(gè)參考點(diǎn)。分別采用楊赤中法、二次規(guī)劃法、最大熵法3種方法計(jì)算參考點(diǎn)權(quán)重系數(shù)進(jìn)行對(duì)比,結(jié)果見表1。

表1 不同插值方法參考點(diǎn)權(quán)重系數(shù)對(duì)比Table 1 Comparison of reference points weight coefficients with different interpolation methods

根據(jù)楊赤中法,不考慮負(fù)權(quán)問題即直接解算方程組(1),解得λ2=-0.033 5,結(jié)果證明負(fù)權(quán)問題是存在的。

利用二次規(guī)劃法解算權(quán)系數(shù)[13], 參照楊赤中方法算得的權(quán)系數(shù), 將負(fù)權(quán)系數(shù)改為0.001, 設(shè)權(quán)系數(shù)初始值為[0.050 0.001 0.688 0.111 0.152], 計(jì)算結(jié)果見表1中的二次規(guī)劃法。

在MATLAB平臺(tái)利用遺傳算法函數(shù)進(jìn)行編程,對(duì)最大熵模型進(jìn)行解算?；静襟E是:先利用式(2)作為目標(biāo)函數(shù), 編制適應(yīng)度評(píng)價(jià)函數(shù); 再編制主程序及確定初始值X0, 初始值即初始權(quán)系數(shù), 其值選擇恰當(dāng), 則能夠加快優(yōu)化計(jì)算的收斂速度。 根據(jù)參考點(diǎn)與估值點(diǎn)的距離反比例計(jì)算權(quán)系數(shù), 賦權(quán)公式為pi=1/di, 其中di為參考點(diǎn)i到估值點(diǎn)的距離。 對(duì)pi歸一化得到程序調(diào)用的初始權(quán)系數(shù)為0.023 9、 0.030 2、 0.403 6、 0.311 2、 0.231 1。 由此計(jì)算得矩約束條件的系數(shù)E1和E2。

在遺傳算法的主程序里,首先指定調(diào)用的目標(biāo)函數(shù),再確定約束條件的各系數(shù)矩陣,設(shè)置參數(shù)及調(diào)用遺傳算法函數(shù)的主要代碼為:

ConstraintFunction=[];

% 沒有非線性約束條件

options=gaoptimset(‘PlotFcns’,{@gaplotbestf,@gaplotmaxconstr}, ‘Display’,‘iter’);

options=gaoptimset(options,‘InitialPopulation’,X0);

[x,fval]=ga(ObjectiveFunction,nvars,A,b,Aeq,beq,LB,UB, ConstraintFunction,options)

運(yùn)行程序,可解算得各pi值即權(quán)系數(shù)得到結(jié)果(表1),可見最大熵法能消除空間權(quán)重系數(shù)負(fù)值現(xiàn)象。

3.1 權(quán)值的對(duì)比

由表1數(shù)據(jù)可見,楊赤中法存在負(fù)權(quán)現(xiàn)象(圖1)。二次規(guī)劃方法雖然能消除負(fù)權(quán)現(xiàn)象,但是得到的權(quán)值分布不合理,離待估點(diǎn)最近的參考點(diǎn)權(quán)值接近0值,遠(yuǎn)離待估點(diǎn)的參考點(diǎn)權(quán)值卻是最大的,此現(xiàn)象不符合空間數(shù)據(jù)相關(guān)的規(guī)律。負(fù)權(quán)值對(duì)應(yīng)的2#參考點(diǎn)雖然離插值點(diǎn)并不是最遠(yuǎn)的,但被靠近插值點(diǎn)的3#參考點(diǎn)所屏蔽,最大熵法求得2#參考點(diǎn)的權(quán)值為非負(fù),且大于最遠(yuǎn)的1#參考點(diǎn)的權(quán)值,這表明該法是符合空間變量相關(guān)特征的。

3.2 精度的對(duì)比

3.2.1 抽樣統(tǒng)計(jì)高程插值中誤差 對(duì)研究區(qū)域劃分12.5 m×12.5 m格網(wǎng), 利用楊赤中、 二次規(guī)劃和最大熵3種方法分別進(jìn)行格網(wǎng)點(diǎn)高程插值, 抽樣統(tǒng)計(jì)格網(wǎng)點(diǎn)的高程插值中誤差結(jié)果分別為±1.123 6、 ±2.168 6和±0.963 2 m。 可見, 最大熵法算得的估值精度最高, 滿足國家有關(guān)DEM建模精度的規(guī)范要求[14]。

圖2為上述區(qū)域采用楊赤中方法、顧及自相關(guān)特征最大熵方法生成的等高線,可見等高線回放效果有一些差異,后者反映出更多的局部地形特征細(xì)節(jié)。

圖2 不同高程插值方法生成的等高線Fig.2 Contour generated by different elevation interpolation methods

3.2.2 各種插值方法構(gòu)建DEM的高程比較 將最大熵法插值得到的網(wǎng)格點(diǎn)高程數(shù)據(jù)在ArcGIS軟件構(gòu)建DEM,并與克里金法、反距離加權(quán)法構(gòu)建的DEM對(duì)照分析如圖3所示。

圖3 不同高程插值方法的DEM柵格圖對(duì)比Fig.3 Comparison of DEM raster maps with different elevation interpolation methods

作柵格運(yùn)算,將最大熵法DEM分別減去克里金法和反距離加權(quán)法的DEM,結(jié)果如圖4所示(黑色、 灰色、 白色分別表示正最大、 最小、 負(fù)最大)。在圖4a中,高程差值絕對(duì)值最大為2 m,在起伏較大區(qū)域,最大熵方法平滑效果比較明顯, 即在山頂最大熵法高程值低于克里金法;在谷底最大熵法高程值高于克里金法;圖4b中,高程差值最大為3 m,且沒有太明顯的規(guī)律,由誤差理論分析得知該數(shù)值符合3.2.1節(jié)各種高程插值方法對(duì)中誤差的描述。

圖4 不同高程插值方法柵格相減結(jié)果Fig.4 Results of raster subtraction with different elevation interpolation methods

3.3 時(shí)間效率的對(duì)比

最大熵法、楊赤中法、二次規(guī)劃法、克里金法、反距離加權(quán)法的插值時(shí)間效率分別為3.667、2.985、3.534、3.018、2.365 s。從計(jì)算效率來看,由于附加了空間相關(guān)性約束條件、縮小了可行解范圍,相較于文獻(xiàn)[7]不帶空間相關(guān)性約束條件的最大熵法加快了模型解算收斂速度。由此可見,反距離加權(quán)法速度最快,而最大熵法運(yùn)算效率最慢,主要原因是最大熵方法解算最優(yōu)值時(shí)采用了智能算法GA,其使用隨機(jī)搜索并進(jìn)行了多次的迭代,這是其不足之處。

4 結(jié)束語

對(duì)于傳統(tǒng)的DEM格網(wǎng)插值數(shù)學(xué)模型,由于其在擬合估值核函數(shù)時(shí)存在主觀因素的干擾,模擬得出的數(shù)學(xué)模型存在不穩(wěn)健、不合理的現(xiàn)象,以至于有時(shí)出現(xiàn)負(fù)權(quán)系數(shù)。本文采用最大熵為目標(biāo)函數(shù)的建模方法,在約束中增設(shè)空間數(shù)據(jù)的相關(guān)性特征約束,最大程度地排除了主觀因素,消除了空間權(quán)重系數(shù)負(fù)值現(xiàn)象,得到的權(quán)值大小比例分布與參考點(diǎn)位分布協(xié)調(diào)一致,插值精度更高。這說明該法不僅能客觀反映出空間數(shù)據(jù)的離散性,也使得模型系統(tǒng)更加穩(wěn)健。從而DEM插值精度有了較大的提高。不足之處是,顧及空間相關(guān)特征的最大熵法模型解算總的運(yùn)算速度較慢,下一步研究將轉(zhuǎn)向減少其迭代次數(shù)或者采用解析算法以提高效率與精度。