MapX遙感影像處理系統坡度和通視算法的研究與實現

摘 要：本論文描述了VC環境下開發的MapX遙感影像處理系統，基于GDAL（Geospatial Data Abstraction Library）實現坡度量算和通視分析的算法。

關鍵詞：MapX遙感影像處理系統；坡度量算；通視分析；GDAL

中圖分類號：TP751 文獻標識碼：A 文章編號：2096-4706（2018）06-0047-02

Abstract：This paper describes the MapX remote sensing image processing system developed under the environment of VC，and the algorithm based on GDAL（Geospatial Data Abstraction Library）to realize the slope calculation and the general view analysis.

Keywords：MapX remote sensing image processing system；slope measurement；visibility analysis；GDAL

0 引 言

遙感圖像處理（processing of remote sensing image data）是對遙感圖像進行輻射校正和幾何糾正、圖像整飾、投影變換、鑲嵌、特征提取、分類以及各種專題處理等一系列操作，以求達到預期目的的技術。遙感圖像處理可分為兩類：一是利用光學、照相和電子學的方法對遙感模擬圖像（照片、底片）進行處理，簡稱為光學處理；二是利用計算機對遙感數字圖像進行一系列操作，從而獲得某種預期結果的技術，稱為遙感數字圖像處理。

在VC環境下使用MapX進行遙感影像數字圖像處理的過程中，使用MapX本身自帶的功能可以實現簡單的地圖操作，但是對量算常用的坡度量算和通視分析等方法沒有提供現成的工具。本論文就是通過分析與設計實現了坡度量算和通視分析的自定義工具。

1 關鍵技術研究

1.1 高程數據

對于TIF格式的遙感影像數字圖像，我們使用GDAL庫獲取高程數據。GDAL（Geospatial Data Abstraction Library）是一個在X/MIT許可協議下的開源柵格空間數據轉換庫。GDAL提供對多種柵格數據的支持，包括Arc/Info ASCII Grid（asc），GeoTiff（tiff），Erdas Imagine Images（img），ASCII DEM（dem）等格式。它利用抽象數據模型來表達所支持的各種文件格式。它還有一系列命令行工具來進行數據轉換和處理。

1.2 坡度

坡度的計算方法有百分比法、度數法、密位法和分數法四種，其中以百分比法和度數法較為常用。（1）百分比法：表示坡度最為常用的方法，即兩點的高程差與其水平距離的百分比，其計算公式如下：坡度=（高程差/水平距離）×100%；（2）度數法：用度數來表示坡度，利用反三角函數計算而得，其公式tanα（坡度）=高程差/水平距離，所以α（坡度）=tan-1（高程差/水平距離）。這里我們選用百分比法。

1.3 點對點通視

點對點通視分析有四種常用方法：（1）Janus算法：以固定步長插視線上點的高程，高程內插采用格網單元的四個頂點進行，該方法計算點少，效率較高；（2）Dyntacs算法：通過視線在XY平面的投影與格網單元邊交點進行判斷，根據相交邊頂點高程通過反距離加權內插出交點高程，依據斜率判斷同時與否。該方法精度較高，但速度較慢；（3）ModSAF算法：該方法在Dyntacs算法上改進，不斷判斷視線與格網單元邊的交點，還要判斷視線與單元格對角線的交點，同時觀察點與目標點高程由其所在三角面上的三個頂點決定，這種方法精度更好，但計算量也隨之加大；（4）Bresenham算法：該算法在Bresenham法繪制直線的基礎上形成，高程點取視線通過格網中心點高程。這種方法速度快，但精度不高。

通過以上分析，這里我們采用Dyntacs算法。

2 系統設計

2.1 高程數據獲取

取視線點的經緯度pointX和pointY，使用GDAL庫的RasterIO（）函數獲取視線點的高程pointEleva，算法主要步驟如下：

GetDem（double pointX， double pointY， double pointEleva）

{

GDALDataset *poDataset；//GDAL數據集

GDALAllRegister（）；//注冊所有的驅動

CString tifname；//加載tif高程文件

tifname.Format（\"dem//ASTGTM2_N%dE0%d_dem.tif\"， （int）pointY， （int）pointX）；

poDataset = （GDALDataset *）GDALOpen（tifname， GA_ReadOnly）；

GDALRasterBand *poBand1；//獲取圖像波段

poBand1 = poDataset->GetRasterBand（1）；

int nImgSizeX = poDataset->GetRasterXSize（）；//獲取圖像的尺寸X

int nImgSizeY = poDataset->GetRasterYSize（）；//獲取圖像的尺寸Y

int nImgc = poDataset->GetRasterCount（）；//獲取圖像的波段個數

double trans[6]；//獲取坐標變換系數

CPLErr aaa = poDataset-> GetGeo Transform （trans）；

GDALDataType btype = poBand1->GetRaster DataType（）；

vector vecElevation；

BYTE *pafScanblock1；//開辟緩存區

pafScanblock1 = （BYTE *）CPLMalloc（sizeof（BYTE）*（1）*（1））；//建議能小則小否則會造成內存不足的情況

double XMin = trans[0]；//取X方向最小值

double YMin = trans[3]；//取Y方向最小值

double XMax = trans[0]+nImgSizeX*trans[1]+ trans[2]；//取X方向最大值

double YMax = trans[3]+trans[4]+nImgSizeY* trans[5]；//取Y方向最大值

int col = （pointX-XMin-trans[2]）/trans[1]；//獲取選取點所在行數

int row = （pointY-YMin-trans[4]）/trans[5]；//獲取選取點所在列數

poBand1->RasterIO（GF_Read， col， row， 1， 1， pafScanblock1， 1， 1， btype， 0， 0）；

pointEleva = pafScanblock1[0]；//獲取選取點高程數據

delete poDataset；

}

2.2 坡度量算

坡度量算如圖1所示，地圖上點1（p1x，p1y）和點2（p2x，p2y）的坡度算法主要步驟如下：（1）取點1高程數據GetDem（p1x，p1y，p1h）；（2）取點2高程數據GetDem（p2x，p2y，p2h）；（3）計算兩點之間的平面距離：dis=Distance（p1x，p1y，p2x，p2y）；（4）計算兩點坡度slope=（fabs（p2h-p1h）/dis）*100。

2.3 點對點通視分析

點對點通視分析如圖2所示，地圖上視點1（p1x，p1y）和目標點2（p2x，p2y）的通視算法主要步驟如下：（1）取視點1高程數據GetDem（p1x，p1y，p1h）；（2）取目標點2高程數據GetDem（p2x，p2y，p2h）（3）計算視點1與目標點2二者之間的坡度S；（4）取視點下一步（圖像的步幅為1秒）的平面坐標（pix，piy）：pix=p1x+i*cx/3600.0；piy=p1y+i*（row/col）*cy/3600.0（其他row和col表示視點和目標點之間的格網總行數值和總列數值；cx和cy表示視點和目標點之間的平面坐標的漸進關系，目標點比視點大，則為1，反之為-1；i表示步數，最大為row）；（5）取所有步點的高程值GetDem（pix，piy，pih），并計算其與目標點之間的坡度Si，若該值大于視點與目標點的坡度，則不通視。

3 結 論

通過以上研究，我們在MapX遙感影像處理系統中實現了坡度量算和通視分析自定義工具。不僅如此，在以MapX為開發平臺的其他系統也可以參考此方法進行設計與實現。

參考文獻：

[1] 李民錄，GDAL源碼剖析與開發指南 [M].北京：人民郵電出版社，2014.

[2] 張宏，溫永寧，劉愛利，等.地理信息系統算法基礎 [M].北京：科學出版社，2010.

作者簡介：李剛（1982-），男，漢族，山東棗莊人，工程師，本科。研究方向：軟件系統開發、軟件測試；王光輝（1987-），女，漢族，黑龍江依蘭人，工程師，研究生。研究方向：軟件測試、軟件質量管理。