梁志強(qiáng)，叢喜東

(1.黑龍江省自然資源權(quán)益調(diào)查監(jiān)測(cè)院，黑龍江 哈爾濱 150080；2.黑龍江省生態(tài)研究所，黑龍江 哈爾濱 150080)

1)2021年5月國(guó)家林業(yè)和草原局發(fā)布《林草濕數(shù)據(jù)與第三次全國(guó)國(guó)土調(diào)查數(shù)據(jù)對(duì)接融合技術(shù)指南》(以下簡(jiǎn)稱(chēng)“指南”)中指出，充分發(fā)揮第三次全國(guó)國(guó)土調(diào)查(以下簡(jiǎn)稱(chēng)“三調(diào)”)成果在國(guó)土空間管理中的“統(tǒng)一底版”作用，依據(jù)“三調(diào)”成果，理清林地、草地、濕地的現(xiàn)狀范圍界線，解決地類(lèi)交叉重疊問(wèn)題，融合林地、草地、濕地等資源信息，優(yōu)化國(guó)家級(jí)公益林范圍，將國(guó)家級(jí)公益林落到山頭地塊，形成與“三調(diào)”無(wú)縫銜接的林草資源圖(含濕地)，支撐森林生態(tài)效益補(bǔ)償制度，為林草生態(tài)網(wǎng)絡(luò)感知系統(tǒng)建設(shè)，開(kāi)展林草資源及生態(tài)狀況監(jiān)測(cè)，加強(qiáng)林草資源科學(xué)保護(hù)、系統(tǒng)修復(fù)和合理利用，推進(jìn)林草治理體系和治理能力現(xiàn)代化提供支撐[1]。

2)指南中明確了目的、任務(wù)、原則要求、技術(shù)流程、工作成果、質(zhì)量管理等內(nèi)容，并對(duì)落界經(jīng)營(yíng)界線和小班細(xì)化做了具體說(shuō)明。指南中指出經(jīng)營(yíng)界線存在誤差性質(zhì)偏移，無(wú)實(shí)質(zhì)改變的情況下，應(yīng)保持與三調(diào)界線重合；小班進(jìn)行分割及碎班合并，對(duì)于區(qū)劃尺度不同或區(qū)劃誤差，引起小班界線相互間里出外進(jìn)等情況，可根據(jù)影像特征和現(xiàn)地情況，對(duì)小班邊界進(jìn)行修正或合并小班，原則上應(yīng)以較高區(qū)劃精度結(jié)果為準(zhǔn)。

3)在具體的對(duì)接融合過(guò)程中由于林草濕數(shù)據(jù)與三調(diào)數(shù)據(jù)在區(qū)劃精度上存在一定的差異，常規(guī)GIS軟件，對(duì)容差范圍內(nèi)的現(xiàn)象不視為錯(cuò)誤，故存在邊線或折點(diǎn)無(wú)法重疊情況，現(xiàn)將具體問(wèn)題進(jìn)行分析，并針對(duì)問(wèn)題提出解決思路與方法，最終形成滿足對(duì)接融合要求的修正數(shù)據(jù)。

1 對(duì)接融合方法概述

1)采用GIS軟件，將林草濕圖層分別以縣為單位，進(jìn)行關(guān)鍵字賦值，用于圖斑融合后掛接屬性數(shù)據(jù)[2]，為加快圖形融合數(shù)據(jù)處理速度，可保留適量相應(yīng)圖層的影響圖斑邊界處理的重要屬性因子，其他屬性因子則在對(duì)接融合后進(jìn)行掛接，如三調(diào)圖斑保留標(biāo)識(shí)碼、地類(lèi)、權(quán)屬單位名稱(chēng)、坐落單位名稱(chēng)、圖斑面積等，林地保留單位名稱(chēng)、林班、小班、地類(lèi)等，濕地保留斑塊號(hào)、斑塊名稱(chēng)、濕地類(lèi)型等，草地保留斑塊編號(hào)、草地類(lèi)型等；

2)將統(tǒng)一空間參考的三調(diào)數(shù)據(jù)圖斑與林地、濕地、草地圖斑分別進(jìn)行聯(lián)合操作，形成相應(yīng)的對(duì)接融合圖層；對(duì)不滿足最小區(qū)劃面積條件的細(xì)碎圖斑進(jìn)行合并，合并后圖斑不可跨三調(diào)圖斑，這部分可通過(guò)技術(shù)手段自動(dòng)合并；對(duì)大于最小區(qū)劃面積的圖斑，需要采用人工方式根據(jù)實(shí)際情況進(jìn)行合并、切割、保留等處理；圖斑融合完畢后進(jìn)行空間拓?fù)洹①|(zhì)檢、修正，掛接重要屬性因子，進(jìn)行數(shù)據(jù)統(tǒng)計(jì)分析，形成空間數(shù)據(jù)、專(zhuān)題圖、統(tǒng)計(jì)分析表、分析報(bào)告等驗(yàn)收成果。

2 對(duì)接融合數(shù)據(jù)邊線及折點(diǎn)不重疊分析

在空間數(shù)據(jù)對(duì)接融合過(guò)程中，由于第三次全國(guó)國(guó)土調(diào)查數(shù)據(jù)矢量圖在GIS軟件中XY容差為0.0001m，通過(guò)分析三調(diào)數(shù)據(jù)圖斑折點(diǎn)坐標(biāo)發(fā)現(xiàn)，所有折點(diǎn)坐標(biāo)均進(jìn)行了處理，采用了4舍6入的方式，保留5位小數(shù)。而以往林草濕數(shù)據(jù)圖斑的生產(chǎn)過(guò)程中采用GIS軟件默認(rèn)XY容差為0.001m，折點(diǎn)坐標(biāo)采用默認(rèn)保留小數(shù)位，二者在容差及坐標(biāo)精度方面存在一定差異。在圖斑對(duì)接融合過(guò)程中，在容差范圍內(nèi)二者邊界雖然是一致的，但在放大比例尺到一定范圍時(shí)，二者將出現(xiàn)不完全重合的現(xiàn)象，即邊線或折點(diǎn)存在位置偏移，且該偏移值在GIS軟件的容差范圍內(nèi)，圖斑進(jìn)行編輯操作時(shí)無(wú)法實(shí)現(xiàn)國(guó)土圖斑的折點(diǎn)捕捉，從而修改操作難度大，且修改量大。如圖1對(duì)接融合圖斑與三調(diào)數(shù)據(jù)圖斑邊界及折點(diǎn)不重合現(xiàn)象示例所示。

圖1 對(duì)接融合圖斑與三調(diào)圖斑邊界及折點(diǎn)不重合示例Fig.1 Example of the Mismatch between the Boundaries and Vertices of the Butt Fusion Patch and Tritone Patch

3 解決問(wèn)題的思路和方法

針對(duì)上述對(duì)接融合數(shù)據(jù)所產(chǎn)生的情況，可以通過(guò)遍歷對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)獲取距離三調(diào)圖斑的最近折點(diǎn)或最近邊線的垂點(diǎn)的方式實(shí)現(xiàn)對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)坐標(biāo)的修改，使對(duì)接融合數(shù)據(jù)圖斑邊線及折點(diǎn)與國(guó)土數(shù)據(jù)圖斑邊線及折點(diǎn)保持完全一致。具體方法步驟如下：

1)獲取對(duì)接融合數(shù)據(jù)圖斑邊線的所有折點(diǎn)。將對(duì)接融合數(shù)據(jù)增加一個(gè)唯一值字段并賦值，用于圖斑修正后數(shù)據(jù)驗(yàn)證及屬性數(shù)據(jù)掛接。

2)獲取該對(duì)接融合數(shù)據(jù)圖斑對(duì)應(yīng)的國(guó)土圖斑邊線折點(diǎn)。對(duì)接融合數(shù)據(jù)圖斑的標(biāo)識(shí)碼(BSM)與三調(diào)數(shù)據(jù)圖斑標(biāo)識(shí)碼應(yīng)一致，可通過(guò)BSM篩選三調(diào)數(shù)據(jù)圖斑。如標(biāo)識(shí)碼不一致會(huì)導(dǎo)致修正折點(diǎn)坐標(biāo)值計(jì)算有誤。

3)獲取對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)與對(duì)應(yīng)三調(diào)數(shù)據(jù)圖斑距離最近的折點(diǎn)。逐點(diǎn)遍歷對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)，查找與該折點(diǎn)距離最近的對(duì)應(yīng)的三調(diào)數(shù)據(jù)圖斑折點(diǎn)，且該距離小于0.0005m，即視為該對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)應(yīng)與三調(diào)數(shù)據(jù)圖斑折點(diǎn)為同一點(diǎn)，對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)的坐標(biāo)用該查找到的三調(diào)數(shù)據(jù)圖斑折點(diǎn)坐標(biāo)進(jìn)行替換。

4)獲取距離對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)最近的國(guó)土數(shù)據(jù)圖斑邊線上的線段，取點(diǎn)到線段的垂點(diǎn)作為折點(diǎn)。當(dāng)未找到滿足具體誤差范圍(<0.0005m)的最近折點(diǎn)，說(shuō)明該對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)誤差限范圍內(nèi)無(wú)三調(diào)圖斑折點(diǎn)，此時(shí)需查找與該折點(diǎn)最近的三調(diào)圖斑邊線上的線段，且該點(diǎn)與線段兩個(gè)端點(diǎn)的距離與線段長(zhǎng)度的誤差也在0.0005m范圍內(nèi)，即該點(diǎn)位于線段的兩個(gè)端點(diǎn)之間。獲取該對(duì)接融合數(shù)據(jù)圖斑折點(diǎn)到三調(diào)圖斑線段的垂足點(diǎn)，且垂足點(diǎn)與折點(diǎn)距離在誤差限范圍內(nèi)，將垂點(diǎn)坐標(biāo)替換對(duì)接融合數(shù)據(jù)圖斑的折點(diǎn)坐標(biāo)，也可以獲取對(duì)接融合數(shù)據(jù)折點(diǎn)處線段與三調(diào)圖斑最近線段的交點(diǎn)作為替換坐標(biāo)值，通過(guò)分析垂點(diǎn)與交點(diǎn)誤差很小，故取垂點(diǎn)坐標(biāo)即可。

5)創(chuàng)建新的矢量圖層，將對(duì)接融合數(shù)據(jù)的每個(gè)圖斑折點(diǎn)按照重新計(jì)算的折點(diǎn)坐標(biāo)值重新生成圖斑要素。

4 軟件工具實(shí)現(xiàn)及成果驗(yàn)證

通過(guò)圖斑的折點(diǎn)修正方法實(shí)現(xiàn)對(duì)接融合數(shù)據(jù)與三調(diào)圖斑邊界的重合操作，圖斑的折點(diǎn)遍歷、點(diǎn)間距、點(diǎn)線距、垂點(diǎn)等計(jì)算采用Visual Studio開(kāi)發(fā)環(huán)境、C#語(yǔ)言、GDAL開(kāi)源技術(shù)進(jìn)行實(shí)現(xiàn)。

1)GDAL開(kāi)源技術(shù)。GDAL(Geospatial Data Abstraction Library，地理空間數(shù)據(jù)抽象庫(kù))是一個(gè)在X/MIT許可協(xié)議下的開(kāi)源空間數(shù)據(jù)(柵格數(shù)據(jù)和矢量數(shù)據(jù))轉(zhuǎn)換庫(kù)，包括讀取、寫(xiě)入、轉(zhuǎn)換、處理各種柵格和矢量數(shù)據(jù)格式[3]，可通過(guò)GDAL讀取矢量圖斑要素屬性及幾何圖形坐標(biāo)節(jié)點(diǎn)等空間及屬性信息。在Visual Studio開(kāi)發(fā)環(huán)境中，可加載GDAL的DLL類(lèi)庫(kù)文件實(shí)現(xiàn)對(duì)矢量圖層及圖斑的具體操作。

2)通過(guò)C#編碼的方式實(shí)現(xiàn)軟件功能的開(kāi)發(fā)。軟件工具的主要功能是加載對(duì)接融合和三調(diào)成果矢量數(shù)據(jù)，自動(dòng)讀取對(duì)接融合數(shù)據(jù)圖斑的邊線折點(diǎn)及對(duì)應(yīng)三調(diào)圖斑的邊線或折點(diǎn)，通過(guò)計(jì)算獲取對(duì)接融合數(shù)據(jù)圖斑的每個(gè)折點(diǎn)的修正值，并自動(dòng)創(chuàng)建新的圖斑要素，要素包含對(duì)接融合數(shù)據(jù)圖斑的關(guān)鍵字及三調(diào)圖斑標(biāo)識(shí)碼字段。對(duì)接融合數(shù)據(jù)修正后結(jié)果可通過(guò)關(guān)鍵字及BSM進(jìn)行相應(yīng)林草濕數(shù)據(jù)與三調(diào)數(shù)據(jù)的屬性連接。

軟件實(shí)現(xiàn)的部分關(guān)鍵代碼如下：

//計(jì)算點(diǎn)與點(diǎn)距離

private double GetMinDistanceP2P(double p1x, double p1y, double p2x, double p2y)

{

return Math.Pow((Math.Pow((p1x - p2x), 2) + Math.Pow((p1y - p2y), 2)), 0.5);

}

//計(jì)算點(diǎn)到線段的距離

private double GetMinDistanceP2L(double pX, double pY, double segAX, double segAY, double segBX, double segBY)

{

double product_ap_ab = (pX - segAX) * (segBX - segAX) + (pY - segAY) * (segBY - segAY);

double product_ab_ab = (segBX - segAX) * (segBX - segAX) + (segBY - segAY) * (segBY - segAY);

double r = product_ap_ab / product_ab_ab;

if (r >=1)

{

return Math.Sqrt((pX-segBX)*(pX-segBX) + (pY - segBY) * (pY - segBY));

}

else if (r <=0)

{

return Math.Sqrt((pX-segAX)*(pX-segAX) + (pY - segAY) * (pY - segAY));

}

else

{

double vector_ac_X = r * (segBX - segAX);

double vector_ac_Y = r * (segBY - segAY);

double vector_cp_X = pX - segAX - vector_ac_X;

double vector_cp_Y = pY - segAY - vector_ac_Y;

return Math.Sqrt(vector_cp_X * vector_cp_X + vector_cp_Y * vector_cp_Y);

}

}

//獲取點(diǎn)到線段的垂足坐標(biāo)，x1、y1、x2、y2為線段的端點(diǎn)坐標(biāo)

private double[] GetVerticalAixP2L (double x0, double y0, double x1, double y1, double x2, double y2)

{

double[] ret = new double[] { 0, 0 };

if (x2 - x1 == 0)

{

ret[0] = x2;

ret[1] = y0;

}

else

{

//計(jì)算斜率

double k1 = (y2 - y1) / (x2 - x1);

// 垂足x坐標(biāo)

double x = (Math.Pow(k1, 2) * x1 + (k1 * (y0 - y1)) + (x0)) / (Math.Pow(k1, 2) + 1);

ouble y = (k1) * (x - (x1)) + (y1);

ret[0] = x;

ret[1] = y;

}

return ret;

}

(3)通過(guò)修正前后數(shù)據(jù)比對(duì)分析。如圖3對(duì)接融合圖斑邊界及折點(diǎn)修正前后對(duì)比圖所示，紅色邊線為對(duì)接融合數(shù)據(jù)的修正前圖斑邊線，紫色邊線為修正后邊線，黑色邊線為三調(diào)圖斑邊線，通過(guò)對(duì)比圖可以看出修正后邊線與三調(diào)數(shù)據(jù)圖斑邊線及折點(diǎn)完全一致，并通過(guò)對(duì)某鄉(xiāng)鎮(zhèn)16000多個(gè)圖斑進(jìn)行驗(yàn)證，修正前后對(duì)接融合數(shù)據(jù)的面積、周長(zhǎng)等均與對(duì)接融合數(shù)據(jù)圖斑未修正前保持基本一致，其修正后結(jié)果滿足數(shù)據(jù)生產(chǎn)和成果質(zhì)量要求。

圖3 對(duì)接融合圖斑邊界及折點(diǎn)修正前后對(duì)比Fig.3 Contrast before and after the Patch Boundary and Vertex Correction of the Butt Fusion Pattern