柵格數(shù)據(jù)有效/無效掩膜的計算與表達(dá)

因為多種原因，柵格數(shù)據(jù)中經(jīng)常存在一些無效像素。在航天測繪遙感中，遙感器獲取的單景影像在經(jīng)過投影坐標(biāo)變換后往往會存在黑邊（見圖 1(a)），黑邊內(nèi)的像素為無效像素。在航空測繪遙感中，位于測區(qū)外的像素也可能被拍攝/掃描到，但是在成果交付時往往只需要輸出測區(qū)內(nèi)的像素，測區(qū)外的像素即為無效像素（見圖 2(a)）。柵格數(shù)據(jù)常采用二值的掩膜來表達(dá)像素有效或無效。

基于像素值的掩膜計算

單個波段的掩膜的計算非常簡單。設(shè)對于某個波段，值為0代表無效值，值非0代表有效值，則利用Envi的band math功能計算掩膜的方法為：

byte((b1 ne 0)*255)

設(shè)值范圍在[M,N]內(nèi)的為有效值，其它值為無效值，則計算掩膜的方法為：

byte((b1 ge M and b1 le N)*255)

基于各個波段的掩膜，可以計算出整個數(shù)據(jù)集的掩膜。對于3波段真彩色影像，計算方法為：

byte((b1 eq 255) and (b2 eq 255) and (b3 eq 255))*byte(255)

在生成的掩膜中，255代表有效像素，0代表無效像素。相關(guān)示例可見圖 1。

圖 1 基于像素值的掩膜計算示例

基于矢量邊界的掩膜計算

無論是基于測區(qū)邊界，圖幅邊界還是ROI生成掩膜，均可以歸結(jié)為基于矢量邊界的掩膜計算。如果矢量邊界對應(yīng)的面狀幾何元素存在多個不重疊的部分，例如測區(qū)由多個不聯(lián)通、不重疊的子測區(qū)構(gòu)成[1]，則矢量邊界可以表達(dá)為OGC簡單特征規(guī)范定義的多多邊形（MultiPolygon），否則可以表達(dá)為多邊形（Polygon），多邊形允許帶有空洞[2]。

基于矢量邊界計算掩膜，可以采用矢量到柵格轉(zhuǎn)換的方法（需要確保兩種柵格數(shù)據(jù)的分辨率一致），也可以采用幾何分析的方法。下面介紹幾何分析的方法。

文獻(xiàn)[3]提出一種基于四叉樹的方法，其核心思想是對柵格數(shù)據(jù)進(jìn)行四叉樹分割，判斷分割出的每個部分對應(yīng)的矩形和矢量邊界對應(yīng)的面狀幾何元素之間是否存在重疊。如果存在重疊，則對矩形再次進(jìn)行四叉樹分割，直至分割出的每個部分所含像素數(shù)目低于指定閾值。由于該算法是遞歸執(zhí)行的，在柵格數(shù)據(jù)行列數(shù)較多時，對?？臻g的需求較大，容易發(fā)生棧溢出。為提高計算穩(wěn)定性，可以逐塊（Block）計算。對于支持按塊存儲的格式，例如Tiff，按塊存儲柵格數(shù)據(jù)可以有效降低計算對?？臻g的需求。設(shè)每塊大小為256*256，則四叉樹最大深度為8層。對于Asc格式，雖然其不支持按塊存儲，但該方法仍然可用，只不過一行為一個塊，四叉樹分割變?yōu)槎鏄浞指睢?/span>

判斷單個像素是否在矢量邊界內(nèi)，可采用如下兩種判斷方法：

(1)基于像素中心的方法。如果像素中心在矢量邊界對應(yīng)的面狀幾何元素內(nèi)，則認(rèn)為像素為有效像素，否則為無效像素；

(2)基于像素邊界的方法。計算像素邊界對應(yīng)的矩形和矢量邊界對應(yīng)的面狀幾何元素之間的重疊面積。如果重疊面積占像素面積的比例大于指定閾值p，則認(rèn)為像素為有效像素，否則為無效像素。其中，如果p取0，則意味著如果像素有任何一部分在矢量邊界內(nèi)，像素即為有效像素；如果p取0.5，且在單個像素內(nèi)的矢量邊界為直線段，則此時該判斷方法和基于像素中心的方法等價。

圖 2是利用矢量邊界裁剪DEM數(shù)據(jù)的效果圖。

圖 2 基于矢量邊界的裁剪DEM的效果，紅線為矢量邊界

掩膜的表達(dá)形式

掩膜的主要表達(dá)形式包括：

(1)掩膜文件。GDAL采用后綴為msk，但實質(zhì)為TIFF格式的柵格文件（簡稱掩膜文件）存儲掩膜數(shù)據(jù)。掩膜文件的行列數(shù)、塊大小、路徑與對應(yīng)的柵格數(shù)據(jù)相同，但文件名比柵格數(shù)據(jù)多“.msk”。掩膜文件存儲的數(shù)據(jù)類型為1字節(jié)無符號整型（下文簡稱uint_8t），0代表像素?zé)o效，非0（一般是255）代表像素有效。掩膜文件的波段數(shù)要么為1，適用于整個柵格數(shù)據(jù)，要么和柵格數(shù)據(jù)的波段數(shù)相同，第K個波段適用于柵格數(shù)據(jù)的第K個波段。

(2)無數(shù)據(jù)值（nodata value）。在GeoTiff、Asc、ERDAS IMAGINE等格式中，支持無數(shù)據(jù)值屬性。有的格式支持對單個波段設(shè)置無數(shù)據(jù)值（如GeoTiff），有的格式只支持對數(shù)據(jù)集整體設(shè)置無數(shù)據(jù)值（如Asc）。常見的GIS和遙感軟件均支持該屬性（見圖 3），在讀取波段數(shù)據(jù)時，如果某像素在某波段的值為指定的無數(shù)據(jù)值，則認(rèn)為該像素在該波段無效。在使用該方式時，需要將無效像素的值設(shè)置為無數(shù)據(jù)值。在DEM/DSM中，通常選擇一個明顯不合理的高程值作為無數(shù)據(jù)值，例如-9999 m。在數(shù)據(jù)類型為uint_8t的3波段真彩色DOM中，通常選擇0或255作為無數(shù)據(jù)值。但是該做法的風(fēng)險在于，DOM中的有效像素的波段值也可能為0或255。例如，陰影、水體等低光部分在某些波段可能值為0（見圖 4(a)），白車車頂、云等高光部分在某些波段可能值為255。所以，為確保有效像素不會出現(xiàn)無數(shù)據(jù)值，有可能需要修改像素值[3]。

圖 3 常見GIS遙感軟件對無數(shù)據(jù)值的支持

(3)alpha波段。在數(shù)據(jù)類型為uint_8t的影像中，alpha波段常被用來存儲不透明度（0代表完全透明，255代表完全不透明）。當(dāng)使用alpha波段存儲掩膜時，GDAL規(guī)定無效像素表示為0值，有效像素表示為非0值（效果見圖 4(b)）。該方法在非遙感GIS領(lǐng)域內(nèi)應(yīng)用較多。但缺點在于，一來只能用于數(shù)據(jù)類型為uint_8t的真彩色影像或灰度影像；二來會變相增加1個波段，占用較多的存儲空間。

(4)全部像素均有效。在無特殊標(biāo)記/說明的情況下，全部像素均為有效像素。

在GDAL中，按(1)、(2)、(3)、(4)的優(yōu)先順序識別掩膜。

圖 4 采用不同方式的掩膜的效果

參考文獻(xiàn)

[1]程曉光,魏婧,嚴(yán)明,等. 應(yīng)用地表覆蓋數(shù)據(jù)估算LiDAR內(nèi)業(yè)工作量的方法研究[J]. 地理信息世界, 2021, 28(2): 89-94.

[2]程曉光,陸泉源,鄭誠慧,等. 使用開源軟件的地圖分幅軟件設(shè)計與實現(xiàn)[J]. 測繪地理信息, 2021, 46(4): 122-127.

[3]佘佐明,程曉光,王艷軍. 基于GDAL開源庫的海量DOM分幅裁剪[J]. 城市勘測, 2023, No.195(1): 74-78.