劉 旭, 趙 錚, 陳建華, 甘先霞, 張洪吉

(1.成都理工大學 地球物理學院，成都 610059;2.四川省自然資源科學研究院，成都 610015)

0 引言

隨著互聯(lián)網(wǎng)技術的快速發(fā)展，地理信息系統(tǒng)(Geographic Information System，GIS)與互聯(lián)網(wǎng)相結合形成的網(wǎng)絡地理信息系統(tǒng)(WebGIS)，為地理信息數(shù)據(jù)發(fā)布、共享和交流協(xié)作提供了良好的基礎[1]。WebGIS已經(jīng)融入到人們的日常生活和工作中，在城市交通、環(huán)境檢測[2]、地質災害等領域都有著良好地應用。但是目前WebGIS還存在著如下問題：①空間分析能力較弱，特別是涉及遙感圖像計算的空間分析；②地理信息數(shù)據(jù)發(fā)布步驟復雜，需要繁瑣的發(fā)布流程。

空間分析作為地理信息系統(tǒng)的核心功能，目前WebGIS已經(jīng)實現(xiàn)了一些基本的基于矢量的空間分析功能，包括路徑規(guī)劃、緩沖區(qū)分析、空間查詢[3]等，但是對于遙感圖像的空間分析還較為薄弱。遙感圖像具有數(shù)據(jù)量大，像素點多的特點，導致遙感圖像處理和分析通常運算復雜，占用大量運行內存，在過去服務器性能較低時，遙感圖像分析通常在桌面端完成，其運算速度取決于硬件的好壞。云計算的出現(xiàn)提高了服務器性能，其通過網(wǎng)絡集中計算資源并按需調配，達到節(jié)約經(jīng)濟地利用計算資源的目的，能夠快速地實現(xiàn)大規(guī)模計算，其效率遠高于本地進行遙感圖像分析。因此筆者針對在服務器開展遙感圖像的空間分析進行了研究。

GIS服務用于處理客戶端讀取GIS數(shù)據(jù)的請求[4]，將地理信息數(shù)據(jù)發(fā)布在服務器是WebGIS客戶端訪問數(shù)據(jù)的基礎，對于服務器的空間分析結果也應該能夠發(fā)布到服務器以便于其共享。但空間分析結果具有動態(tài)變化的特點，其數(shù)據(jù)參數(shù)隨著輸入數(shù)據(jù)不同而不同，因此需要一種能自動發(fā)布空間分析結果的方法。

GIS服務一般由GIS服務發(fā)布工具提供，GIS服務發(fā)布工具主要有ESRI的ArcGIS Server和開源的GeoServer等[5]。ArcGIS Server能夠提供數(shù)據(jù)訪問服務，以及有限的空間分析服務，但是價格昂貴，并且由于是商業(yè)軟件，功能擴展能力較弱；GeoServer能夠提供數(shù)據(jù)訪問服務，空間分析服務能力較弱，但具有軟件開源、擴展能力強[6]的優(yōu)點。

基于上述需求，筆者面向開源GIS服務器GeoServer，在其基礎上擴展空間分析能力，為客戶端提供動態(tài)空間分析服務，實現(xiàn)全自動將動態(tài)空間分析結果發(fā)布為GeoServer服務，并為了加速數(shù)據(jù)訪問速度在此基礎上對發(fā)布的數(shù)據(jù)自動切片。

1 服務自動發(fā)布與切片方法

1.1 服務器端體系結構

服務端體系結構見圖1，服務端提供了兩種服務：①處理基本業(yè)務邏輯和執(zhí)行空間分析功能的空間分析服務;②基于GeoServer的GIS服務。

圖1 體系結構設計圖

空間分析服務是由基本的Web服務器提供，負責處理普通的數(shù)據(jù)管理業(yè)務和執(zhí)行空間分析功能。空間分析服務能夠接受客戶端的HTTP請求，并將處理結果返回給客戶端；能夠調用GeoServer實現(xiàn)對分析數(shù)據(jù)的自動發(fā)布和自動切片；可以讀寫數(shù)據(jù)庫及數(shù)據(jù)文件?？蛻舳讼蚩臻g分析服務發(fā)送空間分析操作請求以及需要分析的數(shù)據(jù)，空間分析服務根據(jù)操作請求調用相應的處理函數(shù)，實現(xiàn)對數(shù)據(jù)的空間分析，并將空間分析的結果根據(jù)其類型分別存入數(shù)據(jù)庫和數(shù)據(jù)文件。在空間分析完成后，空間分析服務自動將結果文件的路徑等參數(shù)信息利用HTTP請求的方法，發(fā)送給GeoServer服務，調度GeoServer發(fā)布結果數(shù)據(jù)并將其切片。

GeoServer服務由開源GIS服務器GeoServer提供，負責各種地理信息數(shù)據(jù)發(fā)布。GeoServer服務只能接收客戶端對于獲取地理信息數(shù)據(jù)的請求；不能接收客戶端上傳的地理信息文件；不能接收客戶端發(fā)送的服務發(fā)布請求，發(fā)布服務只能通過空間分析服務調用；只具備讀取數(shù)據(jù)文件的權限，不能修改數(shù)據(jù)文件。

1.2 空間分析服務

由于GeoServer服務并不具有空間分析能力，因此需要一個單獨的提供空間分析的服務，它能夠接收客戶端請求和客戶端上傳的各種地理信息數(shù)據(jù)，將上傳的數(shù)據(jù)按照一定的組織結構有序地存放在服務器端的文件系統(tǒng)中，并根據(jù)客戶端請求對文件執(zhí)行相應的空間分析功能，得到分析后的結果。

空間分析又分為基于柵格的空間分析和基于矢量的空間分析?；跂鸥竦目臻g分析可以通過由開源地理空間數(shù)據(jù)抽象庫(Geospatial Data Abstraction Library，GDAL)編碼實現(xiàn)，其是一個在X/MIT協(xié)議下的開源柵格空間數(shù)據(jù)轉換庫，具有開源、高效、跨平臺的優(yōu)點，能夠讀取、寫入、轉換、處理大多數(shù)柵格格式文件[7]?；谑噶康目臻g分析服務，可以通過開源地理信息數(shù)據(jù)庫PostGIS或GDAL的子項目OGR(OpenGIS Simple Features Reference Implementation)庫實現(xiàn)。綜上GDAL庫相比PostGIS的空間分析功能更加強大，兼容柵格數(shù)據(jù)和矢量數(shù)據(jù)，因此我們使用Python作為基礎編程語言，結合GDAL實現(xiàn)了一系列的空間分析功能，其中包括柵格計算、線性變換、柵格分類、柵格圖像生成、空間插值等。

論文中針對于柵格數(shù)據(jù)的空間分析，具體的方法為：使用GDAL提取柵格數(shù)據(jù)波段信息，其主要代碼如下：

band= self.dataset.GetRasterBand(nband)

data = band.ReadAsArray(0, 0, self.cols, self.rows)

其中第一行是獲取波段對象；第二行是獲取波段信息，ReadAsArray將讀取的波段信息存入數(shù)組。

在獲取兩幅影像的波段信息之后，對兩幅影像的DN值求差，但值得注意的是由于圖像DN值沒有負數(shù)，因此求差結果需要去絕對值，其主要代碼如下：

data = ds.getRasterInformation(n)[0:ds.computeRows(),0:db.computeCols()-computoffset(ds,db)[0]]

result=abs(data2-data1)

其中,第一行主要是讀取影像的DN值，第二行是獲取到兩幅影像的DN值之后對其求差并取絕對值。

1.3 GeoServer服務自動發(fā)布

GeoServer是遵循OpenGIS規(guī)范實現(xiàn)的開源地圖服務器，其允許用戶對地理數(shù)據(jù)(包括矢量數(shù)據(jù)和柵格數(shù)據(jù))進行更新、刪除及插入操作，支持WMTS、WCS、WFS和WMS服務，方便用戶發(fā)布地圖數(shù)據(jù)并共享地理空間信息[8]。GeoServer是輕量級的GIS服務器，其自帶的功能滿足多數(shù)Web地圖應用開發(fā)，且與前端框架OpenLayers集成較好[9]，但是其空間分析功能相對ArcGIS Server較弱，因此需要擴展實現(xiàn)空間分析功能，以增強GeoServer的分析能力。

假定GeoServer已經(jīng)在計算機中配置完成，GeoServer沒有可視化用戶界面，其會構建一個服務，通過瀏覽器訪問IP地址對GeoServer進行管理。GeoServer手動發(fā)布服務的步驟為通過瀏覽器訪問GeoServer提供的IP地址，選擇自己需要發(fā)布的服務類型，然后填寫對應的參數(shù)，包括基本資源信息、坐標參考系統(tǒng)、覆蓋參數(shù)、波段計算細節(jié)等。參數(shù)配置過程繁瑣，且需要根據(jù)數(shù)據(jù)信息選擇參數(shù)，這些配置的選擇要求發(fā)布人員具有一定的專業(yè)知識，限制了GeoServer發(fā)布服務的靈活性。

為了提高GeoServer發(fā)布服務的靈活性、簡化相關人工工作步驟，空間分析動態(tài)生成結果的自動發(fā)布顯得十分迫切。GeoServer提供了一種REST接口，用戶可以通過這個接口檢索有關實例的信息并進行配置更改，使用REST接口接收HTTP請求，可以達成配置GeoServer服務的目的，無需再通過瀏覽器訪問Web管理界面。REST接口實質上對應一個HTTP地址，通過get、post、put、delete四種不同的請求方式達成對GeoServer實例的讀取、增加、修改、刪除的目的。REST接口為自動化發(fā)布GeoServer服務提供了基礎。

在空間分析服務接收到客戶端發(fā)送的請求數(shù)據(jù)以及發(fā)送過來的文件后，從請求數(shù)據(jù)中解析出需要執(zhí)行的空間分析功能，調用對應的封裝好的空間分析接口，對數(shù)據(jù)進行空間分析，得到空間分析結果并保存到文件系統(tǒng)中。由于發(fā)布服務需要用到發(fā)布數(shù)據(jù)的相關參數(shù)信息，包括基本資源信息、坐標參考系統(tǒng)、邊界矩形等，因此需要利用GDAL讀取待發(fā)布數(shù)據(jù)的對應基本信息。

以讀取影像數(shù)據(jù)的參數(shù)為例，獲取包括影像的波段信息、影像數(shù)據(jù)范圍以及影像仿射變換矩陣，在獲取到參數(shù)之后，按照一定組織規(guī)律生成影像的配置文件。在生成配置文件之后，根據(jù)具體的發(fā)布需求，在配置文件中加上服務配置參數(shù)，包括影像背景值、發(fā)布后影像的數(shù)據(jù)格式等，從而完成對配置文件的設置。

將配置文件綁定到HTTP請求中，向GeoServer REST接口發(fā)送請求，在發(fā)布成功之后GeoServer服務會將發(fā)送結果返回給空間分析服務，最后將結果返回給客戶端，以柵格計算為例，自動發(fā)布服務流程圖如圖2所示。

圖2 自動發(fā)布流程圖

以空間分析服務向GeoServer REST接口發(fā)送新增圖層的方法為例，利用Python發(fā)送請求的主要代碼如下：

ps=requests.post(url=url,auth=('admin','geoserver'),data=open(coveragepath+'/coverage.xml').read(), headers=headers)

該方法使用post請求向Geoserver發(fā)送數(shù)據(jù)，其中url參數(shù)指定了Geoserver服務器的IP地址，auth參數(shù)指定了Geoserver服務器的用戶名和密碼；data參數(shù)指定了配置文件的路徑以及名稱；headers參數(shù)指定了http請求的頭部信息。

1.4 GeoServer自動切片

地理信息數(shù)據(jù)往往包含大量信息，因此地理信息數(shù)據(jù)占用的空間都是比較大的，如果要一次性加載全部信息，會導致服務器壓力大、客戶端等待時間長、客戶端使用卡頓的問題。因此，可以將較大的地理信息數(shù)據(jù)采用金字塔結構(四叉樹)進行組織和存儲，切分成大量不同級別、分辨率不同的小幅數(shù)據(jù)，每提高一個瀏覽級別，就用四張更高分辨率的數(shù)據(jù)代替較低分辨率的一張數(shù)據(jù)。用戶在客戶端查看數(shù)據(jù)的時，一次僅對少量需要的數(shù)據(jù)進行請求，服務器發(fā)送響應數(shù)據(jù)即可滿足用戶需求，這樣減少了網(wǎng)絡數(shù)據(jù)量，可大大減少系統(tǒng)響應時間，改善用戶體驗。

GeoServer提供了對數(shù)據(jù)切片的功能，在客戶端加載數(shù)據(jù)時，GeoServer在實時的對數(shù)據(jù)進行切片，切片之后再根據(jù)請求范圍返回對應的數(shù)據(jù)，實時切片的方式增加了服務器的響應時間。如果能在服務發(fā)布之后自動將服務數(shù)據(jù)切片，在訪問之前預先完成切片將減小服務響應時間。

由于GeoServer REST接口并不提供執(zhí)行切片功能的接口，因此之前使用REST接口自動發(fā)布的方式在自動切片部分不再適用。對于已經(jīng)發(fā)布的數(shù)據(jù)，可以通過訪問GeoServer提供的對應IP地址在可視化頁面中交互式的操作數(shù)據(jù)(圖3)。在對數(shù)據(jù)進行切片時，需要選擇切片級數(shù)范圍、數(shù)據(jù)的坐標系、切片后的數(shù)據(jù)格式等，參數(shù)設置完成后，點擊“Submit”按鈕開始對數(shù)據(jù)進行切片。

圖3 自動切片可視化界面

分析發(fā)現(xiàn)，GeoServer管理界面的URL固定，以IP地址加問號加圖層ID構成，并且每次訪問的頁面保持一致且易于登錄，因此可以通過Web自動化測試工具(Web爬蟲)完成自動切片。Web自動化測試工具通過編寫腳本代碼，在執(zhí)行腳本時，自動生成一個模擬的瀏覽器，并訪問給定的URL，按照腳本代碼做出點擊、輸入、打開、驗證等操作。

因此，采用Python語言結合自動化測試工具Selenium，通過添加Geoserver用戶名密碼代理驗證，代理驗證通過后，通過驅動模擬瀏覽器自動訪問服務地址，將相應的參數(shù)通過代碼自動填入到對應的位置，如圖3所示的參數(shù)中Grid Set設置切片的格網(wǎng)方案、Format設置切片的格式、Zoom start和Zoom stop設置切片的級別范圍、Bounding box設置了切片的邊界范圍。在參數(shù)設置完成后，模擬點擊“Submit”按鈕實現(xiàn)數(shù)據(jù)自動切片，自動切片流程如圖4所示。

圖4 自動切片流程圖

2 實例與分析

針對河流水質動態(tài)監(jiān)測的現(xiàn)實需求，以開源GeoServer為基礎，結合上述提出的GeoServer服務自動化方法，設計、實現(xiàn)了一套面向Web的高分遙感河流水質動態(tài)監(jiān)測系統(tǒng)。

高分遙感河流水質動態(tài)監(jiān)測系統(tǒng)通過對相應河流區(qū)域遙感圖像進行反演之后，將任兩期反演結果柵格圖像上傳到服務器，在服務器進行差異性分析，從而判斷水質變化，具體流程如下。

對原始多期國產高分遙感數(shù)據(jù)作處理，包括輻射定標、大氣校正、正射校正、圖像配準、圖像融合、圖像拼接、圖像裁剪、水體提取。依據(jù)水體影像不同波段組合，反演出單波段水質指標濃度(包括葉綠素a、總磷、總氮、高錳酸鹽、懸浮物及透明度等)。由各指標反演結果計算其營養(yǎng)狀態(tài)指數(shù)，并計算綜合營養(yǎng)狀態(tài)指數(shù)。將所有反演結果、營養(yǎng)狀態(tài)指數(shù)、綜合營養(yǎng)狀態(tài)指數(shù)導出為包含濃度信息的柵格圖像數(shù)據(jù)。

客戶端將任兩期高分遙感河流水質反演結果柵格圖像上傳到服務器，并選擇自動發(fā)布、自動切片的相關參數(shù)。

利用空間分析服務對客戶端上傳的兩期結果圖像依次進行如下計算：

1)柵格計算，獲取兩期河流水質反演結果圖像的DN值，對兩者作求差計算處理。

2)線性變換，將DN值按一定的比例對數(shù)據(jù)進行線性變化。

3)柵格分類，按照一定的分類規(guī)則對生成的差異柵格圖像進行分類、著色。

4)柵格圖像生成，在分類完成后生成新的差異柵格圖像并保存至文件系統(tǒng)。

在差異柵格圖像生成基礎上，讀取其文件系統(tǒng)中柵格圖像相關參數(shù)信息(坐標系統(tǒng)、邊界矩形等一系列相關信息)，調用根據(jù)GeoServer服務自動發(fā)布原理構建的自動發(fā)布接口，完成柵格圖像發(fā)布，并根據(jù)參數(shù)決定是否自動切片，在執(zhí)行完成后服務器向客戶端返回發(fā)布結果(圖5)。

圖5 差異柵格圖像自動生成、自動發(fā)布、自動切片流程圖

基于上述原理，實現(xiàn)的面向Web的高分遙感河流水質動態(tài)監(jiān)測系統(tǒng)整體效果(圖6)，差異柵格圖像在自動切片后不同顯示級別的效果(圖7)。

圖6 系統(tǒng)整體效果

圖7 切片顯示柵格圖像效果

系統(tǒng)實現(xiàn)差異柵格圖像自動生成、自動發(fā)布、自動切片功能，取代復雜的傳統(tǒng)發(fā)布過程，具有顯著的省時、省力、自動化等特點?？蛻舳嗽L問服務，加載差異柵格圖像并疊加柵格影像數(shù)據(jù)更加真實的反應了河流水質變化(圖7)，圖7(a)～圖7(d)為四個不同展示級別的結果，四個圖像具有不同的分辨率。表明：將圖像自動發(fā)布、自動切片是提升WebGIS服務質量的有效途徑。

3 結論

這里提出了一種面向空間分析的開源GeoServer服務自動發(fā)布與自動切片方法，運用Python與GDAL構建空間分析服務，彌補了GeoServer空間分析能力不足的弱點，針對空間分析結果動態(tài)變化的特點，實現(xiàn)了結果的自動服務發(fā)布、自動切片。該方法顯著提升了基于GeoServer的WebGIS系統(tǒng)的靈活性，解決了傳統(tǒng)GIS服務需要專業(yè)人士手動發(fā)布、手動配置，無法動態(tài)發(fā)布等問題。將該方法運用到面向Web的高分遙感河流水質動態(tài)監(jiān)測系統(tǒng)，能夠有效實現(xiàn)差異柵格圖像動態(tài)生成并對其自動發(fā)布、自動切片，實現(xiàn)了水質監(jiān)測結果快速生成、持久保存、快速加載顯示的目的。