基于地理處理服務(wù)的滑坡敏感性區(qū)劃研究

王衛(wèi)東，曾科，方理剛

(中南大學(xué) 土木工程學(xué)院，湖南 長沙，410075)

我國是世界上滑坡災(zāi)害最嚴(yán)重的國家之一。近10年來，我國地質(zhì)災(zāi)害預(yù)防和防治水平逐步提高，然而，對地質(zhì)災(zāi)害的預(yù)防研究較少[1]，因此，應(yīng)加強地質(zhì)災(zāi)害預(yù)警預(yù)報研究?；旅舾行詤^(qū)劃屬于滑坡空間預(yù)測的一個分支，能定量或半定量地評價區(qū)域內(nèi)滑坡敏感性等級程度，可有效應(yīng)用于防災(zāi)減災(zāi)工程[2]。地理信息系統(tǒng)(GIS)因其具有強大的地理信息處理、存儲、管理、分析和展示功能[3]，已成為分析地質(zhì)災(zāi)害的主要技術(shù)方法之一。網(wǎng)絡(luò)地理信息系統(tǒng)(WebGIS)的出現(xiàn)，使得地理信息的網(wǎng)絡(luò)共享成為可能。在 ArcGIS Server平臺上二次開發(fā)的地理處理服務(wù)(GP服務(wù))能將桌面系統(tǒng)的空間分析功能移植到Web應(yīng)用程序中，極大地增強了 Web應(yīng)用程序的地理信息處理能力[4]。在此，本文作者以主客觀綜合權(quán)重方法為原型構(gòu)建滑坡敏感性區(qū)劃 GP模型，并將 GP模型發(fā)布為網(wǎng)絡(luò)服務(wù)，實現(xiàn)基于GP服務(wù)的貴州省滑坡敏感性區(qū)劃圖自動編制。

1 GP服務(wù)的作用與實現(xiàn)

1.1 GP服務(wù)及其作用

地理信息服務(wù)(Geographic information web services, GIWS )是在網(wǎng)絡(luò)環(huán)境下，與地理信息相關(guān)的軟件功能實體，它建立了一種面向服務(wù)的地理信息共享模式，通過其接口提供服務(wù)所封裝的功能，用戶可以通過互聯(lián)網(wǎng)按需獲得和使用地理數(shù)據(jù)及計算服務(wù)，如地圖服務(wù)、空間數(shù)據(jù)格式轉(zhuǎn)換等[5]。GIWS可分為3類：地理數(shù)據(jù)服務(wù)、地理描述服務(wù)和地理處理服務(wù)[6-7]。

地理數(shù)據(jù)服務(wù)提供地理數(shù)據(jù)的管理、編輯、更新和網(wǎng)絡(luò)發(fā)布等功能，是復(fù)雜地理信息服務(wù)的基礎(chǔ)。地理描述服務(wù)包括語義描述服務(wù)和結(jié)果描述服務(wù)，語義描述服務(wù)用于描述、查詢和編目地理信息服務(wù)[8]，結(jié)果描述服務(wù)定義了結(jié)果傳輸格式和表現(xiàn)形式。地理處理服務(wù)提供GIS分析功能，強調(diào)地理數(shù)據(jù)與用戶的交互性，能充分挖掘隱藏的地理信息，其分析結(jié)果為決策提供依據(jù)[9]。GP服務(wù)所封裝的分析功能通過商業(yè)GIS軟件的二次開發(fā)實現(xiàn)，可提供研究領(lǐng)域內(nèi)的 GIS分析過程。與其他服務(wù)相比，GP服務(wù)是地理信息網(wǎng)絡(luò)服務(wù)的高級應(yīng)用和核心模塊，相當(dāng)于一般網(wǎng)絡(luò)信息系統(tǒng)中的模型庫。

本文中GP服務(wù)封裝了滑坡敏感性區(qū)劃圖編制的整個流程。專家通過互聯(lián)網(wǎng)對致災(zāi)因子的重要性進行評價后，GP服務(wù)將自動計算研究區(qū)內(nèi)各地理單元的滑坡敏感性程度(敏感性分值)，根據(jù)敏感性分值繪制滑坡敏感性區(qū)劃圖，實現(xiàn)了滑坡敏感性區(qū)劃圖編制的自動化和網(wǎng)絡(luò)化。

GP服務(wù)對地質(zhì)災(zāi)害研究的意義有：

(1) 二次開發(fā)的 GP服務(wù)封裝了地質(zhì)災(zāi)害預(yù)警預(yù)報的數(shù)學(xué)模型，一般用戶不需要知道模型的具體實現(xiàn)，只需了解它的功能和操作即可進行地質(zhì)災(zāi)害分析，降低了地質(zhì)災(zāi)害預(yù)測分析的“門檻”。如本文中，一般用戶在編制滑坡敏感性區(qū)劃圖的過程中不需要了解主客觀綜合權(quán)重方法的數(shù)學(xué)模型。

(2) 地質(zhì)災(zāi)害分析處理自動化，GP服務(wù)根據(jù)用戶輸入?yún)?shù)和預(yù)先設(shè)計的處理流程進行地理信息分析，在計算過程中不需要用戶參與[10]。本文開發(fā)的GP服務(wù)實現(xiàn)了滑坡敏感性區(qū)劃圖自動編制，簡化了手動分析時的復(fù)雜性，極大縮短了分析處理時間，提高了預(yù)警預(yù)報效率。

(3) 地理信息資源分布式部署，GP服務(wù)運行于服務(wù)器端，用戶通過互聯(lián)網(wǎng)協(xié)議獲取GP功能，本地計算機不需要安裝GIS軟件，減少了軟件購置支出。

1.2 GP服務(wù)的實現(xiàn)

GP服務(wù)的核心是GP模型。ArcGIS提供了表1所示的GP工具和可視化的GP建模環(huán)境，供專業(yè)人員進行GP模型二次開發(fā)。用戶可將GP工具組合成復(fù)雜的GIS工作流(圖1)，其中矩形方框代表GP工具，橢圓形代表輸入/輸出數(shù)據(jù)，六邊形表示中間數(shù)據(jù)，字母P標(biāo)識整個模型的輸入/輸出參數(shù)。

本文采用Python語言開發(fā)表1中的腳本工具調(diào)用GIS底層組件，用于擴展GIS功能和實現(xiàn)主客觀綜合權(quán)重方法的數(shù)學(xué)計算。GP模型在服務(wù)器端構(gòu)建之后，通過ArcGIS Server Manager發(fā)布為GP服務(wù)。

2 研究區(qū)致災(zāi)因子體系

貴州省位于北緯 103°36′～109°35′，東經(jīng) 24°37′～29°13′，占地面積176 167 km2，地處云貴高原，地形起伏大，地質(zhì)構(gòu)造復(fù)雜，歷史滑坡災(zāi)害嚴(yán)重，分布在全省絕大部分區(qū)域(圖2)。特別是每年3月至8月的雨季，強降水導(dǎo)致滑坡災(zāi)害頻發(fā)，嚴(yán)重威脅人民生命財產(chǎn)的安全。

圖1 地理處理建模示例Fig.1 Geoprocess modeling

表1 GP工具類型Table 1 Types of geoprocess tool

圖2 貴州省歷史滑坡災(zāi)害分布圖Fig.2 Landslide inventory map of Guizhou Province

控制和影響地質(zhì)災(zāi)害孕育發(fā)生的主要影響因素可以歸為3大類，即環(huán)境因素、誘發(fā)因素和歷史因素[11]。本文選取高程、坡度、巖性和地質(zhì)構(gòu)造4個環(huán)境因素，年均降雨量和鐵路建設(shè)2個誘發(fā)因素作為滑坡主要致災(zāi)因子(圖3)。同一因子取值不同時對滑坡發(fā)生的影響程度不同，因此，需要將各一級因子進一步劃分為若干二級因子。其中地質(zhì)構(gòu)造、鐵路建設(shè)分別以到地質(zhì)構(gòu)造線距離和到鐵路距離來劃分二級因子。巖性等致災(zāi)因子作為地質(zhì)變量，必須對其進行量化。由于各因子的測量單位不一致，導(dǎo)致每個因子的表現(xiàn)力不同，因此，還需要統(tǒng)一量綱。各因子分級標(biāo)準(zhǔn)、量化和量綱統(tǒng)一方法見文獻[11]。分級后的致災(zāi)因子體系如圖3所示。在因子屬性表中使用對應(yīng)的因子標(biāo)識表示各二級因子地理區(qū)域，本文重點介紹基于該致災(zāi)體系的GP模型建立過程。

3 滑坡敏感性區(qū)劃主客觀權(quán)重方法

圖3 滑坡致災(zāi)因子Fig.3 Landslide causal factors

以基于專家評價和統(tǒng)計學(xué)的主客觀權(quán)重方法[11]為原型構(gòu)建滑坡敏感性區(qū)劃GP模型。首先，依據(jù)熵值理論，由二級因子歷史滑坡災(zāi)害面積密度計算一級因子客觀權(quán)重；然后，應(yīng)用梯形模糊數(shù)方法(TFNW)，由專家打分得到二級因子主觀權(quán)重；最后，根據(jù)權(quán)重線性疊加(WLC)模型編制滑坡敏感性區(qū)劃圖。下面介紹該方法的主要計算原理，詳細(xì)計算步驟和實例見文獻[11]。

3.1 一級因子客觀權(quán)重計算

根據(jù)熵值理論[12]，某個一級因子的各個二級因子區(qū)域內(nèi)歷史滑坡災(zāi)害面積密度變化顯著，則認(rèn)為該一級因子對滑坡災(zāi)害的影響較大，即該因子權(quán)重較大。若某一級因子的所有二級因子的歷史滑坡災(zāi)害面積密度變化很小，則認(rèn)為該一級因子對滑坡災(zāi)害的影響較小，即該因子權(quán)重較小。一級因子權(quán)重計算過程如下。

由第 i個一級因子的第 j個二級因子的歷史災(zāi)點面積 ALij除以該二級因子覆蓋的區(qū)域總面積 ATij得到此二級因子的歷史滑坡災(zāi)害密度dij(見式(1))，并對其進行歸一化處理(見式(2))。

其中：n為一級因子包含的二級因子個數(shù)。根據(jù)pij求得一級因子熵值Hi(見式(3))。最后，依據(jù)式(4)計算一級因子的客觀權(quán)重ωi。

其中：m為一級因子個數(shù)，本文中取6； e = 1 /lnn，是一個常數(shù)，保證 0≤Hi≤1。并規(guī)定：若 pij=0，則ln(pij)=0。

3.2 二級因子主觀權(quán)重計算

本文采用梯形模糊數(shù)方法，通過專家打分確定二級因子權(quán)重。首先，由決策者確定第t個專家權(quán)重Rt(專家總數(shù)為k)；然后，由各專家依據(jù)事先設(shè)計好的語言評語集(見表2)，以語言變量形式對二級因子的影響進行評價。根據(jù)式(5)和式(6)計算第i個一級因子第j個二級指標(biāo)的模糊權(quán)重期望值 )(～ijEω ；最后，依據(jù)式(7)進行歸一化得到其二級因子主觀權(quán)重ijω′。

表2 語言變量和模糊數(shù)Table 2 Linguistic variables and fuzzy numbers

其中：

3.3 綜合權(quán)重及單元敏感性計算

根據(jù)權(quán)重線性組合模型(式(9))[13]，由前述iω和ijω′計算二級因子的主客觀綜合權(quán)重。

每一個地理單元的不同二級因子綜合權(quán)重之和為敏感性分值。對于某地理單元k，其敏感性分值Sk按式(10)計算，依據(jù)該值和分級標(biāo)準(zhǔn)即可得到貴州省滑坡敏感性區(qū)劃圖。

4 滑坡敏感性區(qū)劃地理處理模型構(gòu)建

4.1 模型工作空間

模型構(gòu)建前先設(shè)計模型工作空間，用于分類存儲輸入/輸出數(shù)據(jù)、臨時文件和GP工具(圖4)[14]。

其中，“臨時文件”用于存儲模型執(zhí)行過程中的中間數(shù)據(jù)；“區(qū)域GIS數(shù)據(jù)庫”包含一個Geodatabase數(shù)據(jù)庫，用于存儲和管理模型輸入/輸出數(shù)據(jù)；“滑坡敏感性區(qū)劃工具箱”即GP模型工具箱，是本文GP服務(wù)的主體部分；“貴州省滑坡敏感性區(qū)劃圖”用于存儲滑坡敏感性區(qū)劃GP服務(wù)的成果。

圖4 模型空間目錄Fig.4 Space model directory

4.2 模型工具箱

模型工具箱包括3個模型工具集，共計14個GP工具，它們分別封裝了部分地理處理功能，并組合成“滑坡敏感性區(qū)劃總模型”(圖 5)。其中，“因子權(quán)重計算工具集”提供一級因子客觀權(quán)重和二級因子主觀權(quán)重的計算功能；“敏感性區(qū)劃工具集”封裝了滑坡敏感性區(qū)劃圖編制的4個主要步驟；“輔助腳本工具集”為其他工具/模型提供輔助功能。

圖5 模型工具箱結(jié)構(gòu)Fig.5 Kit structure of model

4.3 模型輸入數(shù)據(jù)

模型輸入數(shù)據(jù)包括致因子圖層和二級因子專家評價信息，專家評價信息由客戶端提交。模型所需的致災(zāi)因子圖層如圖6所示，其中：高程、坡度、巖性、至構(gòu)造線距離、至鐵路距離圖層均來源于國家基礎(chǔ)地理信息系統(tǒng)(NFGIS-M50DB)；年平均降雨量圖層來源于貴州省氣象局；歷史滑坡災(zāi)害分布圖層來自貴州省國土資源廳多年調(diào)查和積累，并經(jīng)實地考察補測得到。所有圖層的比例尺都是1:500 000。

4.4 滑坡敏感性區(qū)劃地理處理模型

4.4.1 滑坡敏感性區(qū)劃總模型

滑坡敏感性區(qū)劃總模型由數(shù)據(jù)預(yù)處理模型、一級因子權(quán)重計算模型、二級因子權(quán)重計算模型，敏感性分值計算與區(qū)劃圖編制模型4部分組成(圖7)。數(shù)據(jù)預(yù)處理模型主要是標(biāo)準(zhǔn)化輸入數(shù)據(jù)，如大地坐標(biāo)系的統(tǒng)一。一級因子權(quán)重計算模型和二級因子權(quán)重計算模型分別封裝了3.1節(jié)和3.2節(jié)的計算過程。敏感性分值計算與區(qū)劃圖編制模型是滑坡敏感性區(qū)劃分析的中心模型，該模型綜合致災(zāi)因子主、客觀權(quán)重計算滑坡敏感性分值并編制滑坡敏感性區(qū)劃圖。

4.4.2 數(shù)據(jù)預(yù)處理模型

數(shù)據(jù)預(yù)處理模型旨在規(guī)范化因子圖層，統(tǒng)一數(shù)據(jù)格式、空間參考系和分析范圍，并將因子圖層疊加，為敏感性分析做準(zhǔn)備[15]。模型輸入數(shù)據(jù)為圖6所示的因子圖層，輸出數(shù)據(jù)為“因子疊加圖層”。因子圖層疊加過程就是屬性表融合過程，融合后屬性表中每一行對應(yīng)一個二級因子組合區(qū)域(圖8)。4.4.3 一級因子客觀權(quán)重計算模型

圖6 因子圖層Fig.6 Factors’ layers

圖7 滑坡敏感性區(qū)劃總模型Fig.7 General model for landslide-susceptibility mapping

根據(jù)3.1節(jié)構(gòu)建一級因子客觀權(quán)重計算模型。該模型輸入數(shù)據(jù)為每個致災(zāi)因子圖層(圖 6)和滑坡災(zāi)害歷史分布圖層(圖1)，輸出數(shù)據(jù)為一級因子客觀權(quán)重。

圖8 因子疊加圖層屬性表Fig.8 Attributes list of multi factors-overlay

圖9 一級因子客觀權(quán)重模型Fig.9 Weighting model for first level factor

以高程因子為列，其權(quán)重模型見圖9。熵權(quán)模型的構(gòu)建步驟為：(1) “區(qū)域因子分類”工具依據(jù)滑坡災(zāi)點歷史分布圖層，統(tǒng)計高程因子圖層中各二級因子范圍內(nèi)歷史滑坡的面積 ALij和各二級因子所覆蓋的區(qū)域面積 ATij，由式(1)計算二級因子災(zāi)害密度 dij，得到“因子災(zāi)害密度表”。(2) “密度歸一化”工具調(diào)用腳本工具，由式(2)計算pij。(3) “熵值計算”工具獲取pij并依據(jù)式(3)計算一級因子熵值 Hi。(4) “一級因子客觀權(quán)重計算”工具依據(jù)一級因子熵值和式(4)計算一級因子客觀權(quán)重。其他致災(zāi)因子熵權(quán)模型與之類似。

4.4.4 二級因子主觀權(quán)重計算模型

專家賦權(quán)界面由客戶端負(fù)責(zé)， GP模型封裝二級因子主觀權(quán)重計算過程。根據(jù)3.2節(jié)構(gòu)建二級因子主觀權(quán)重模型，見圖10。地質(zhì)災(zāi)害專家通過互聯(lián)網(wǎng)，在客戶端用語言變量評價二級因子對滑坡災(zāi)害發(fā)生的影響，若影響大，則用高、很高進行評價。該模型首先依據(jù)表 2，通過“語言變量轉(zhuǎn)換”工具將語言變量翻譯為模糊數(shù)；然后，根據(jù)式(5)和(6)計算二級因子模糊權(quán)重期望值 )(～ijEω；最后，由“權(quán)重歸一化”工具計算二級因子主觀權(quán)重。

圖10 二級因子主觀權(quán)重模型Fig.10 Weighting model for secondary level factor

4.4.5 敏感性分值計算與區(qū)劃圖編制模型

敏感性分值計算與區(qū)劃圖繪制模型可分為因子賦權(quán)、權(quán)重組合、計算敏感性分值、敏感性重分類4個步驟。前3個步驟是敏感性分值的計算，后1個步驟是敏感性區(qū)劃圖分級與繪制。

圖11 敏感性分值計算與區(qū)劃圖繪制模型Fig.11 The model of susceptibility score calculation and landslide-susceptibility mapping

圖12 貴州省滑坡災(zāi)害敏感性區(qū)劃系統(tǒng)Fig.12 Landslides susceptibility mapping system of Guizhou Province

(1) “因子賦權(quán)”工具將主觀權(quán)重和客觀權(quán)重賦給因子圖層屬性表，如圖11所示。它的輸入數(shù)據(jù)為數(shù)據(jù)預(yù)處理模型、一級因子權(quán)重計算模型和二級因子權(quán)重計算模型的輸出結(jié)果?！皺?quán)重組合”根據(jù)式(9)計算二級因子的主客觀綜合權(quán)重，得到綜合權(quán)重圖層。

(2) “計算敏感性分值”工具由式(10)累加同一地理單元內(nèi)的二級因子綜合權(quán)重(6列因子屬性值相加)，得到Sk。

(3) “敏感性重分類”工具采用等間距法，依據(jù)敏感性分值將貴州省劃分為4個區(qū)域：極高危險區(qū)、高危險區(qū)、中等危險區(qū)和低危險區(qū)，模型執(zhí)行完成后將滑坡敏感性區(qū)劃圖(圖12)返回給客戶端。

4.5 模型結(jié)果分析

模型執(zhí)行完成后，從 Geodatabase數(shù)據(jù)庫中調(diào)取滑坡敏感性區(qū)劃圖層與歷史災(zāi)害圖層進行疊加分析，計算得到各敏感級別分布面積百分比和該區(qū)內(nèi)歷史滑坡密度，如表3所示。從表3可以看出：隨著易發(fā)性的增強，歷史滑坡在該區(qū)分布的比例也增大，說明模型執(zhí)行結(jié)果具有可靠性。

表3 滑坡敏感性區(qū)劃圖結(jié)果統(tǒng)計Table 3 Results statistical for landslide-susceptibility mapping

5 滑坡敏感性區(qū)劃系統(tǒng)

在ArcGIS Server平臺下，開發(fā)滑坡敏感性區(qū)劃系統(tǒng)，將貴州省境內(nèi)滑坡的基礎(chǔ)資料、因子致災(zāi)圖層、氣象信息、敏感性區(qū)劃模型和災(zāi)害預(yù)警信息發(fā)布到互聯(lián)網(wǎng)上。系統(tǒng)界面如圖11所示。系統(tǒng)設(shè)定不同的使用權(quán)限，由防災(zāi)減災(zāi)有關(guān)人員(包括決策者、專業(yè)人士、相關(guān)人員)與GP服務(wù)進行交互，遠(yuǎn)程編制滑坡敏感性區(qū)劃圖。其他人員(如社會群眾)可在互聯(lián)網(wǎng)上獲取最新的滑坡預(yù)測信息。專家因子評價界面見圖12左下角。

6 結(jié)論

(1) 在ArcGIS的GP建模環(huán)境(ModelBuilder)下構(gòu)建滑坡敏感性區(qū)劃模型。該模型封裝了GIS空間分析功能、主客觀熵權(quán)模型和 WLC模型，實現(xiàn)了滑坡敏感性區(qū)劃圖自動編制，減少了區(qū)劃圖編制的工作量。

(2) 通過ArcGIS Server Manager將GP模型發(fā)布為網(wǎng)絡(luò)服務(wù)，實現(xiàn)了區(qū)劃圖編制的遠(yuǎn)程共享，一般地質(zhì)災(zāi)害管理人員可能通過互聯(lián)網(wǎng)由客戶端調(diào)用GP服務(wù)為本地區(qū)編制滑坡敏感性區(qū)劃圖。

[1] 殷坤龍. 滑坡災(zāi)害預(yù)測預(yù)報[M]. 北京: 中國地質(zhì)大學(xué)出版社,2004.YIN Kun-long. Landslide hazard prediction and evaluation[M].Beijing: China University of Geosciences Press, 2004.

[2] 陳永坡. 滑坡危險度區(qū)劃研究[D]. 成都: 西南交通大學(xué)土木工程學(xué)院, 2002: 4-21.CHEN Yong-bo. Study on landslide hazard zonation[D].Chengdu: Southwest Jiaotong University. School of Civil Engineering, 2002: 4-21.

[3] Lee S. Application of likelihood ratio and logistic regression models to landslide susceptibility mapping using GIS[J].Environmental Management, 2004, 34(2): 223-232.

[4] 佟晶. 基于ArcGIS Server的地理處理服務(wù)應(yīng)用研究[D]. 北京:中國地質(zhì)大學(xué)土地科學(xué)技術(shù)學(xué)院, 2010: 13-16.DONG Jing. Research on geoprocessing service application based on ArcGIS Server[D]. Beijing: China University of Geosciences. School of Land Science and Technology, 2010:13-16.

[5] Gtinther O, Miiller R. From GISystems to GIServices-Spatial computing on the interact marketplace[C]//Proceedings of the International Conference on Interoperating Geographic Information Systems. Santa Barbara. 1998:427-442.

[6] Open GIS Consortium Inc. Web service common implementation specification[EB/OL]. [2006-06-05].http://www.opengeospatial. org/specs/?page:specs.

[7] Johannes B, Theodor F, Bastian S, Bastian B. Towards a research agenda for geoprocessing services[C]//12th AGILE International Conference on Geographic Information Science. Hannover:Leibniz Universit?t, 2009: 1-12.

[8] Lutz M. Ontology-based descriptions for semantic discovery and composition of geoprocessing services[J]. Geoinformatica, 2007,11(1): 1-36.

[9] Fitzner D, Hoffmann J, Klien E. Functional description of geoprocessing services as conjunctive datalog queries[J].Geoinformatica, 2011, 15(1): 191-221.

[10] Jimenez-Peralvarez J D, Irigaray C, Hamdouni R. Building models for automatic landslide-susceptibility analysis, mapping and validation in ArcGIS[J]. Nat Hazards, 2009, 50(3): 571-590.

[11] 王衛(wèi)東. 基于 GIS的區(qū)域公路地質(zhì)災(zāi)害管理與空間決策支持系統(tǒng)研究[D]. 長沙: 中南大學(xué)商學(xué)院, 2009: 42-47.WANG Wei-dong. Regional highway geological hazard management and spatial decision support system based on GIS[D]. Changsha: Central South University. Business School,2009: 42-47.

[12] Shannon C E, Weaver W. The mathematical theory of communication[M]. Urbana: The University of Illinois Press,1947: 3-19.

[13] Baban S M J, Van Yusof K. Modelling optimum sites for locating reservoirs in tropical environments[J]. Water Resources Management, 2003, 17(1): 1-17.

[14] Enviromental Systems Research Institute Inc. ArcGIS Desktop Help[EB/OL]. [2009-04-06]. http://webhelp.esri.com/arcgisdesktop/9.3.

[15] Fitzner D, Hoffmann J, Klien E. Functional description of geoprocessing services as conjunction datalog queries[J].Geoinformatica, 2011, 15(1): 191-221.