因為多種原因，柵格(ge)數據中(zhong)經常存(cun)在一(yi)些(xie)無效(xiao)像(xiang)素(su)。在航天測繪遙感(gan)中(zhong)，遙感(gan)器獲(huo)取的(de)(de)單景(jing)影(ying)(ying)像(xiang)在經過投影(ying)(ying)坐(zuo)標(biao)變換后往(wang)(wang)往(wang)(wang)會(hui)存(cun)在黑邊（見(jian)圖(tu)(tu) 1(a)），黑邊內的(de)(de)像(xiang)素(su)為無效(xiao)像(xiang)素(su)。在航空測繪遙感(gan)中(zhong)，位于測區外的(de)(de)像(xiang)素(su)也可能(neng)被拍攝/掃描到(dao)，但是在成果交(jiao)付(fu)時往(wang)(wang)往(wang)(wang)只(zhi)需(xu)要(yao)輸出(chu)測區內的(de)(de)像(xiang)素(su)，測區外的(de)(de)像(xiang)素(su)即為無效(xiao)像(xiang)素(su)（見(jian)圖(tu)(tu) 2(a)）。柵格(ge)數據常采用二值的(de)(de)掩(yan)膜來表達像(xiang)素(su)有(you)效(xiao)或無效(xiao)。

基于像素值的掩膜計算

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

byte((b1 ne 0)*255)

設(she)值范圍在[M,N]內的為(wei)(wei)有(you)效(xiao)值，其它(ta)值為(wei)(wei)無效(xiao)值，則計算掩膜的方(fang)法為(wei)(wei)：

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

基(ji)于各個波(bo)段的(de)掩(yan)膜(mo)，可以計算(suan)(suan)出整個數據(ju)集(ji)的(de)掩(yan)膜(mo)。對于3波(bo)段真彩色影像，計算(suan)(suan)方法(fa)為：

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

在生成(cheng)的掩膜(mo)中(zhong)，255代表(biao)(biao)有效像(xiang)素，0代表(biao)(biao)無效像(xiang)素。相關示例(li)可(ke)見圖 1。

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

基于矢(shi)量邊界的掩膜計算

無論是基于測(ce)(ce)區邊(bian)(bian)界(jie)，圖幅邊(bian)(bian)界(jie)還是ROI生成掩膜(mo)，均可(ke)(ke)以歸結為基于矢量(liang)邊(bian)(bian)界(jie)的掩膜(mo)計(ji)算。如果矢量(liang)邊(bian)(bian)界(jie)對應的面狀(zhuang)幾何元素存(cun)在多(duo)(duo)個不(bu)重(zhong)疊(die)的部分，例(li)如測(ce)(ce)區由多(duo)(duo)個不(bu)聯通、不(bu)重(zhong)疊(die)的子測(ce)(ce)區構成[1]，則矢量(liang)邊(bian)(bian)界(jie)可(ke)(ke)以表達為OGC簡單特征規范定義的多(duo)(duo)多(duo)(duo)邊(bian)(bian)形（MultiPolygon），否則可(ke)(ke)以表達為多(duo)(duo)邊(bian)(bian)形（Polygon），多(duo)(duo)邊(bian)(bian)形允許帶有(you)空洞[2]。

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

文(wen)獻[3]提(ti)出一種(zhong)基于(yu)四(si)叉樹(shu)的方(fang)法(fa)，其核心思(si)想(xiang)是(shi)對(dui)(dui)柵格數(shu)據進行(xing)(xing)四(si)叉樹(shu)分(fen)(fen)(fen)(fen)(fen)割，判斷分(fen)(fen)(fen)(fen)(fen)割出的每個(ge)部分(fen)(fen)(fen)(fen)(fen)對(dui)(dui)應的矩形(xing)和矢(shi)量(liang)邊(bian)界對(dui)(dui)應的面狀幾何元(yuan)素(su)之間是(shi)否(fou)存(cun)在(zai)(zai)重(zhong)疊。如果存(cun)在(zai)(zai)重(zhong)疊，則對(dui)(dui)矩形(xing)再(zai)次進行(xing)(xing)四(si)叉樹(shu)分(fen)(fen)(fen)(fen)(fen)割，直至(zhi)分(fen)(fen)(fen)(fen)(fen)割出的每個(ge)部分(fen)(fen)(fen)(fen)(fen)所含像素(su)數(shu)目低于(yu)指定閾值。由于(yu)該算(suan)法(fa)是(shi)遞(di)歸執行(xing)(xing)的，在(zai)(zai)柵格數(shu)據行(xing)(xing)列數(shu)較(jiao)多時，對(dui)(dui)棧空間的需(xu)求(qiu)較(jiao)大(da)，容(rong)易發生棧溢(yi)出。為提(ti)高計(ji)算(suan)穩定性(xing)，可(ke)以(yi)逐(zhu)塊(kuai)(kuai)（Block）計(ji)算(suan)。對(dui)(dui)于(yu)支持按(an)塊(kuai)(kuai)存(cun)儲的格式，例(li)如Tiff，按(an)塊(kuai)(kuai)存(cun)儲柵格數(shu)據可(ke)以(yi)有效降低計(ji)算(suan)對(dui)(dui)棧空間的需(xu)求(qiu)。設每塊(kuai)(kuai)大(da)小為256*256，則四(si)叉樹(shu)最大(da)深度為8層(ceng)。對(dui)(dui)于(yu)Asc格式，雖然其不支持按(an)塊(kuai)(kuai)存(cun)儲，但(dan)該方(fang)法(fa)仍然可(ke)用，只不過一行(xing)(xing)為一個(ge)塊(kuai)(kuai)，四(si)叉樹(shu)分(fen)(fen)(fen)(fen)(fen)割變為二叉樹(shu)分(fen)(fen)(fen)(fen)(fen)割。

判(pan)斷單(dan)個像素是否在矢量邊界內(nei)，可(ke)采用如下(xia)兩種判(pan)斷方法(fa)：

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

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

圖 2是利用矢量邊界裁剪DEM數據的效果圖。

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

掩膜的表達形(xing)式

掩膜的主(zhu)要表達形式包(bao)括：

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

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

圖 3 常見GIS遙感軟件對無數據值的支持

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

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

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

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

參考文獻

[1]程曉光,魏(wei)婧,嚴(yan)明,等(deng). 應用地表覆蓋數據估算LiDAR內業工作量的方法(fa)研(yan)究(jiu)[J]. 地理信息世界(jie), 2021, 28(2): 89-94.

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

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