基于Akima插值的地面沉降等值面服務(wù)構(gòu)建方法

彭 程,張進才,陳 勇,李 斌

(1.長安大學地質(zhì)工程與測繪學院,陜西 西安 710054; 2.河北省地質(zhì)環(huán)境監(jiān)測院,河北 石家莊 050060;3.河北省資源環(huán)境監(jiān)測與保護重點實驗室,河北 石家莊 050060)

0 引 言

隨著信息化技術(shù)的蓬勃發(fā)展，越來越多的學者將信息化技術(shù)應(yīng)用于地質(zhì)災(zāi)害防治領(lǐng)域并取得了一定的成果。在國外，先后有多名學者引入GIS技術(shù)在韓國、美國堪薩斯州、瑞士等地對滑坡、泥石流等地質(zhì)災(zāi)害進行了地質(zhì)監(jiān)測、危險性評估等研究[1]；在國內(nèi)，胡志瑞[2]依托寧夏國土資源調(diào)查監(jiān)測院建立了實時預(yù)警信息化平臺，根據(jù)實時降雨數(shù)據(jù)實現(xiàn)對全區(qū)地質(zhì)災(zāi)害的預(yù)警預(yù)報；靖常峰等人[3]探索了Web模式下降雨分布圖的自動繪制方法，在2012年北京特大暴雨中得到了應(yīng)用。這些研究充分說明了管理信息系統(tǒng)是地質(zhì)災(zāi)害防控的一把利刃[4]。

作為一種分布廣泛的緩變地質(zhì)災(zāi)害，地面沉降是地質(zhì)構(gòu)造活動和人類工程活動共同作用的結(jié)果[5]。隨著經(jīng)濟社會的高速發(fā)展，廣泛分布的地面沉降災(zāi)害給當?shù)厝嗣竦纳敭a(chǎn)安全造成了很大的損失，成為制約社會經(jīng)濟可持續(xù)發(fā)展的不利因素之一。為了滿足新形勢下區(qū)域性地面沉降防控工作的新要求，也為能在京津冀地區(qū)實施嚴格的地面沉降監(jiān)測與防控提供系統(tǒng)層面的管理信息平臺和分析技術(shù)支撐，“河北省地面沉降系統(tǒng)”項目研發(fā)了基于B/S架構(gòu)的河北地面沉降管理與分析系統(tǒng)。由于等值面能夠以直觀的方式顯示區(qū)域地面沉降漏斗的分布，在地面沉降管理與分析系統(tǒng)研發(fā)中，將其作為系統(tǒng)的核心功能進行研究。然而，繪制等值面的數(shù)據(jù)主要來自于鋪設(shè)在地表的水準點標高數(shù)據(jù)，這些數(shù)據(jù)的觀測序列不一致，即數(shù)據(jù)存在系統(tǒng)性缺失，嚴重制約著等值面模型的精度，可能會造成沉降漏斗錯位的情況。

因此，針對水準點數(shù)據(jù)系統(tǒng)性缺失的問題，本文選取具有代表性的4種插值方法從工程實踐的角度進行分析，最終應(yīng)用Akima插值[6-8]的數(shù)據(jù)插值處理方法，在繪制等值面之前完善原始數(shù)據(jù)，以提高等值面的精度。為了將桌面GIS的空間分析功能移植到Web中，本文以河北地面沉降管理與分析系統(tǒng)為基礎(chǔ)，將ArcGIS的GP服務(wù)[9-10]與Akima插值相結(jié)合，構(gòu)建符合河北數(shù)據(jù)特征的地面沉降等值面繪制服務(wù)，簡化手動分析時的復雜性，縮短分析處理時間，實現(xiàn)河北地面沉降等值面的自動繪制，對河北地面沉降防治工作有一定的參考作用。

1 全自動生成地面沉降等值面的技術(shù)路線

針對水準點原始數(shù)據(jù)系統(tǒng)性缺失問題，本文提出的全自動生成地面沉降等值面的技術(shù)流程共分為3個部分，如圖1所示：1)采用Akima插值，對存入數(shù)據(jù)庫中的標高數(shù)據(jù)進行插值，擬合出任意日期的標高，算出累計沉降量，并結(jié)合經(jīng)緯度信息構(gòu)建水準點要素圖層；2)在模型構(gòu)造器上構(gòu)建工具流模型，并依托云服務(wù)器將其發(fā)布成GP服務(wù)；3)調(diào)用GP服務(wù)，在系統(tǒng)的視圖窗口中完成等值面的渲染及圖層圖例顯示。

圖1 全自動生成等值面的GP服務(wù)構(gòu)建

1.1 Akima插值

1.1.1 Akima插值可行性論述

測量數(shù)據(jù)內(nèi)插方法有很多，如線性插值、多項式插值、樣條函數(shù)插值等，為了對各插值方法的優(yōu)劣有更清晰的認識，決定利用不同的插值方法對同一個水準點長序列做插值比較實驗，選取的水準點序列如表1所示。此外，為了保證實驗具有代表性，共選擇了高次拉格朗日插值、分段低次插值、三次樣條曲線插值及Akima插值等4種插值方法進行比較,如圖2所示。

表1 某水準點時序表(時間為2014-11-15的相對時間)

圖2(a)所示曲線為高次拉格朗日插值的實驗結(jié)果，其趨勢與另外3種插值方法具有明顯的差異，這是因為隨著階數(shù)的上升，高階多項式會發(fā)生奇異，出現(xiàn)龍格現(xiàn)象，即內(nèi)插函數(shù)在插值點與實際數(shù)據(jù)符合很好，而在插值點外出現(xiàn)較大的偏差；而圖2(b)所示曲線為分段低次插值的實驗結(jié)果，雖然曲線的總體趨勢與后2種插值相似，但曲線的光滑度很差，這是因為低階多項式的參數(shù)較少，內(nèi)插精度偏低所導致的；由于三次樣條曲線插值和Akima插值都考慮了要素導數(shù)值的效應(yīng)，因而得到的整個插值曲線都是光滑的,如圖2(c)和圖2(d)所示。但二者相比，Akima插值法所得曲線比樣條函數(shù)插值曲線更光滑、更自然。

圖2 同序列插值實驗比較結(jié)果

單從插值曲線來看，三次樣條曲線插值與Akima插值的差異不太明顯，因此，決定從圖2中選取10個點分別采用這2種方法進行內(nèi)插，剩余點作為插值結(jié)果的真實比較對象，最終得到的實驗結(jié)果如表2所示，結(jié)果表明，Akima插值的效果略優(yōu)于三次樣條曲線插值。綜上，選擇Akima插值作為原始水準點數(shù)據(jù)的插值方法。

表2 插值結(jié)果與真實值的差異結(jié)果

1.1.2 Akima插值的基本原理

Akima插值是Akima[11]于1970年提出的一種五點求導分段三次多項式插值算法，用于二維平面曲線的插值與光滑[12]。該算法利用三次多項式在每2個數(shù)據(jù)點間擬合出1條曲線，最終由這些曲線連接而成的總曲線保證一階導數(shù)連續(xù)。

給定n個型值點{(xi,yi)，i=1,…,n},按x坐標進行排列[13]，根據(jù)已知條件可見式(1)：

(1)

可以在任意2個相鄰型值點{(xk,yk)}和{xk+1,yk+1)}間唯一確定一個三次多項式，見式(2)：

a+b×(x-xk)+c×(x-xk)2+d×(x-xk)3,x∈[xk,xk+1],k=1,2,3,…,n-1

(2)

其中，a、b、c、d計算方法如式(3)：

(3)

相鄰2個型值點間線段的斜率見式(4)：

(4)

因計算切線斜率需要，補充m0=2m1-m2，m-1=2m0-m1，mn=2mn+1-mn+2，mn+1=2mn+2-mn+3。

型值點處的切線斜率見式(5)：

(5)

Akima插值算法和三次樣條函數(shù)都思考了要素導數(shù)值產(chǎn)生的效應(yīng),所以最終會得到的一條光滑的曲線。雖然用三次樣條函數(shù)插值法獲得的曲線在逼近和收斂等方面具有很大的優(yōu)勢,然而就結(jié)果而言，Akima插值法獲取的曲線比前者看起來更光順,更自然。

1.2 地理處理服務(wù)

地理處理(Geoprocessing，GP)的作用是簡化地理數(shù)據(jù)處理的流程。創(chuàng)建它的方式共有3種：1)調(diào)用ArcGIS內(nèi)置的工具箱解決單一數(shù)據(jù)處理問題；2)使用Model Builder[14-18]創(chuàng)建工具流解決復雜問題；3)使用Python腳本語言定制工具。方式雖有差別，但目的都是為了能夠自動執(zhí)行GIS空間分析任務(wù)，提高工作效率，簡化工作流程。

地理處理服務(wù)是將地理處理功能放置在Web服務(wù)器端的一種機制：客戶端將數(shù)據(jù)傳入服務(wù)器同時，服務(wù)器調(diào)用配置好的功能組件，將處理結(jié)果返回給客戶端，客戶端不需要安裝任何桌面GIS軟件，即可完成分析任務(wù)。

從ArcGIS Server 10.1開始，地理處理服務(wù)的發(fā)布摒棄了發(fā)布工具箱或發(fā)布帶有工具箱的地圖文檔等方式[19]，用戶通過運行工具流或Python腳本工具，將自己對數(shù)據(jù)的處理方式用share as發(fā)布為服務(wù)，供服務(wù)訂閱者使用。

2 利用Akima插值對原始數(shù)據(jù)進行插值處理

2.1 Akima插值代碼實現(xiàn)

在地面沉降管理過程中，經(jīng)常會遇到水準點標高數(shù)據(jù)部分缺失的情況，地面沉降每時每刻都可能在發(fā)生，由缺失數(shù)據(jù)繪制的等值面在輔助決策時很難發(fā)揮作用，為了確保地面沉降等值面繪制GP工具的精度，決定使用Akima插值對數(shù)據(jù)進行預(yù)處理。Akima插值的目的是通過不少于5組的樣本點數(shù)據(jù)，插值得到每一天對應(yīng)的水準點標高數(shù)據(jù)，利用插值后的數(shù)據(jù)最終算出任意時間段的累計沉降量。

根據(jù)Akima插值的基本原理，commons Math將算法核心封裝在AkimaSplineInterpolator類中的interpolate方法[20]中。這個方法的輸入?yún)?shù)是2組double類型的數(shù)組(每組的樣本數(shù)不得少于5)；方法返回的是由五點求導法確定的多項式樣條函數(shù)。多項式樣條函數(shù)的方法集如表3所示，其中有2種常用的方法：1)isValidPoint方法，可以用來判斷某個時間點是否在插值范圍內(nèi)；2)value方法，可以擬合出區(qū)間內(nèi)任意自變量對應(yīng)的函數(shù)值。

表3 多項式樣條函數(shù)的方法集

在使用Akima插值計算數(shù)據(jù)集的三次樣條插值時，首先需要確定輸入?yún)?shù):其中一組double數(shù)組無疑是水準點的標高數(shù)據(jù)，另一組本應(yīng)是測量日期，但是日期類型不能作為輸入?yún)?shù)，因此需要將日期數(shù)據(jù)轉(zhuǎn)換為double類型，處理方法可以參考計算機對時間的計算和顯示原則，即將水準點標高測量日期轉(zhuǎn)換為距歷元的偏移量存入第二組double數(shù)組中。

水準點標高數(shù)據(jù)預(yù)處理流程如圖3所示，一共需要遍歷2次數(shù)據(jù)庫：外層循環(huán)對數(shù)據(jù)庫中的水準點測量基本信息表(主表)進行遍歷，檢索出存入數(shù)據(jù)庫中的所有水準點的編號及經(jīng)緯度信息(繪制等值面需要經(jīng)緯度信息)；內(nèi)層循環(huán)以水準點編號為條件，對數(shù)據(jù)庫中的水準測量高程表(子表)進行遍歷，檢索出每個水準點的所有標高數(shù)據(jù)，然后以這些有限的標高及日期數(shù)據(jù)，求出每個水準點的樣條曲線函數(shù)，將用戶選擇的起始截止日期傳入函數(shù)中，就能得出每個水準點在這段時間內(nèi)的累計沉降量，最終將水準點編號、經(jīng)緯度信息及插值得到的累計沉降量作為屬性值，建立點要素集。

圖3 數(shù)據(jù)預(yù)處理流程

2.2 Akima算法可靠性檢驗

由于本文中的原始數(shù)據(jù)是離散水準點數(shù)據(jù)，沒有函數(shù)的解析式，故不能直接利用拉格朗日余項做誤差估計。Akima算法底層原理是五點求導分段三次多項式插值，即，利用三次多項式在每2個數(shù)據(jù)點間擬合出1條曲線。因此，每個數(shù)據(jù)點都包括左右2個數(shù)據(jù)區(qū)間，對同一個點的左右相鄰區(qū)間分別做Akima插值，再比較2個插值之間的誤差，就能確定代碼的可靠性。由于Akima算法的輸入樣本不能低于5組，故選擇表4中的5個點進行可靠性檢驗，最終的檢驗結(jié)果如表4所示。

表4 可靠性檢驗的實驗結(jié)果

從表4中可以直觀地看出，5個樣本點的左右區(qū)間生成的插值之差均達到了毫米級，這說明利用Java代碼實現(xiàn)的Akima算法是可靠的。

3 基于GP服務(wù)地面沉降等值面的設(shè)計與實現(xiàn)

3.1 構(gòu)建等值面繪制模型

一個模型由一個或多個過程組成。每個過程都有類似的基本結(jié)構(gòu)：輸入→函數(shù)→輸出[21]，不同模型包含的過程數(shù)可能不同，但整體的結(jié)構(gòu)保持不變。如圖4所示，模型可以是單處理模型，即只有一個輸入輸出過程；也可以通過多過程之間的組合構(gòu)建模型，即前一個過程的輸出結(jié)果作為后一個過程的輸入?yún)?shù)，按順序依次經(jīng)過多個過程。

圖4 處理過程模型分類

ArcGIS自帶的工具箱可視為最簡單的GP模型。在此基礎(chǔ)上，通過在模型構(gòu)造器(Model Builder)上設(shè)置輸入輸出參數(shù)及地理處理工具集，將簡單的GP工具串連成復雜的工作流模型供服務(wù)器統(tǒng)一調(diào)用[22]。

在構(gòu)建等值面繪制模型時，設(shè)置的輸入?yún)?shù)就是由水準點經(jīng)緯度及插值成果數(shù)據(jù)構(gòu)成的累計沉降量點要素集，而輸出參數(shù)則是以output命名的矢量圖層，兩者之間分別用下面的工具箱組合在一起，構(gòu)建的等值面繪制模型如圖5所示。

圖5 等值面繪制模型

1)反距離權(quán)重插值法(IDW)。它指的是以插值點與其鄰近樣本點間的平均距離為權(quán)重進行加權(quán)平均，離插值點越近的樣本點權(quán)重越大。與克里金插值相比，雖然克里金插值生成的柵格過渡性較好，但在某些特殊情況下(如其中存在數(shù)據(jù)為0的點)，克里金插值容易出錯。因此，選擇較為穩(wěn)定的反距離權(quán)重法將點要素類插值成柵格數(shù)據(jù)集，形成沉降量暈眩圖。

2)轉(zhuǎn)為整型工具。為了有效優(yōu)化柵格的磁盤空間，使用轉(zhuǎn)為整型工具通過截斷將柵格的每個像元值都轉(zhuǎn)換為整型[23]，輸出數(shù)據(jù)為整型柵格數(shù)據(jù)集，形成一個沉降量暈眩效果圖。

3)柵格計算器。利用Con(a,b,c)函數(shù)在柵格計算器中提前對柵格值進行歸類，將結(jié)果集輸出。由于柵格計算器沒有輸入?yún)?shù)，普通的連線是無法連接柵格和柵格計算器工具的，為了保證流程的準確性，需要在環(huán)境中設(shè)置柵格捕捉，讓輸出柵格和捕捉柵格像元大小一致。

4)柵格轉(zhuǎn)面工具。使用柵格轉(zhuǎn)面工具能夠?qū)鸥褚丶D(zhuǎn)換為面要素集。

5)裁剪工具。裁剪工具共有2個輸入?yún)?shù)：裁剪范圍和裁剪要素集。本文的裁剪范圍是河北省行政區(qū)劃圖，裁剪的要素集是上一步得到的面要素集。

3.2 發(fā)布GP服務(wù)

構(gòu)建好的地面沉降等值面模型只能在本地實現(xiàn)應(yīng)用分析，要實現(xiàn)Web環(huán)境下的服務(wù)共享，則需要將GP模型運行的結(jié)果以共享的方式進行服務(wù)發(fā)布，在發(fā)布的同時確定各類參數(shù)，比如同步或異步執(zhí)行方式：同步是所有的操作都做完，才返回給用戶結(jié)果；異步則不用等所有操作都做完，就響應(yīng)用戶請求。因此，在較為復雜的GP分析時，使用異步方式時用戶體驗較好。

本文采用ArcGIS Server作為搭載GP服務(wù)的遠程服務(wù)器。ArcGIS Server是Esri推出的支持SOA(面向服務(wù)架構(gòu))的服務(wù)端產(chǎn)品。在ArcGIS Server服務(wù)器上可以查看發(fā)布的地圖服務(wù)和地理處理服務(wù)的相關(guān)信息，如輸入輸出參數(shù)的地理坐標系統(tǒng)、字段標識符號、數(shù)據(jù)類型等，如圖6所示。其中g(shù)ridcode字段是存放各水準點累計沉降量的關(guān)鍵值，可以用于實現(xiàn)分級色彩的顯示和圖例標注的繪制。

圖6 輸出參數(shù)的相關(guān)信息

3.3 GP服務(wù)調(diào)用—ArcGIS API for JavaScript

本文采用ArcGIS API for JavaScript實現(xiàn)GP服務(wù)的調(diào)用，調(diào)用方法的選擇與發(fā)布有關(guān)，若發(fā)布時選擇同步，調(diào)用方法為excute；若發(fā)布時選擇異步，則采用submitJob提交原始數(shù)據(jù)?？紤]沉降等值面的GP服務(wù)相對復雜，因此，決定采用異步的方式，提高用戶的體驗效果。

用戶將生成沉降等值面的請求發(fā)送給Web服務(wù)器，Web服務(wù)器在接收請求之后，對存入數(shù)據(jù)庫中的插值結(jié)果數(shù)據(jù)集進行遍歷，獲取經(jīng)緯度及累計沉降量等信息，構(gòu)造JSON格式的輸入?yún)?shù)；之后調(diào)用Geoprocessor對象的submitJob提交參數(shù)，在ArcGIS Server上輸出等值面矢量數(shù)據(jù)。為了獲取ArcGIS Server上的地理處理結(jié)果，Web服務(wù)器調(diào)用dojo的循環(huán)方法對submitJob的結(jié)果集進行遍歷，獲取結(jié)果集的geometry和屬性值，根據(jù)獲取的geometry和屬性值在Web服務(wù)器上生成一個要素圖層，最后對這個要素圖層添加標注信息及符號化設(shè)置，將添加渲染之后的結(jié)果返回給客戶端，從而實現(xiàn)了等值面GP服務(wù)的調(diào)用，最終客戶端接收到的成果如圖7所示。

圖7 等值面構(gòu)建結(jié)果圖

用戶只需選擇生成等值面的起止日期，系統(tǒng)自動會通過數(shù)據(jù)庫中的原始數(shù)據(jù)生成每個水準點的樣條曲線函數(shù)，算出起止日期內(nèi)的累計沉降量，自動完成地面沉降等值面的繪制工作。用戶不需要在客戶端上安裝專業(yè)性軟件，只需通過瀏覽器向服務(wù)器發(fā)送請求，便可以一鍵式生成等值面，查看區(qū)域的沉降漏斗。

4 結(jié)束語

根據(jù)河北地面沉降項目的收官結(jié)果表明，本文提出的地面沉降等值面自動化繪圖方法在地面沉降防治工作中有著廣泛的應(yīng)用前景，所歸納出的自動化繪圖方法具有獨特的優(yōu)勢和普適性，具體歸納如下：

1)針對在地面沉降等值面開發(fā)中遇到的水準點標高樣本數(shù)據(jù)系統(tǒng)性缺失問題，本文將Akima插值方法融入到了等值面模型的創(chuàng)建過程中，以數(shù)學插值的準確性支撐模型的普適性，模型的構(gòu)建結(jié)果得到了河北省環(huán)境監(jiān)測院相關(guān)用戶的肯定。

2)本文提供了一種自動化繪制等值面的方法。該方法具有快速、高效、自動化程度高等特點，實現(xiàn)了較少的人工干預(yù)，為地面沉降管理工作提供了輔助決策，滿足部門工作者應(yīng)急需要，并能為其他行業(yè)的等值面繪制提供借鑒。

3)依托GP服務(wù)搭建的等值面繪制模型具備微服務(wù)的特征，任何引入ArcGIS API的網(wǎng)站都能以訂閱服務(wù)的方式使用模型，只需要簡單的代碼就可實現(xiàn)網(wǎng)站等值面構(gòu)建功能的開發(fā)。

4)提供的等值面自動化繪制服務(wù)具有擴展性，下一步可以借助Spark等開源集群運算框架，提高模型計算速率，增強制圖GP服務(wù)的并發(fā)承載和業(yè)務(wù)擴展能力。