李 想，高 磊

(1.江蘇星月測(cè)繪科技股份有限公司南京雨花臺(tái)分公司，江蘇 南京 210000；2.江蘇省測(cè)繪工程院，江蘇 南京 210013)

0 引 言

近十年來(lái)，國(guó)產(chǎn)遙感衛(wèi)星發(fā)射數(shù)量不斷增加，空間、時(shí)間和光譜分辨率不斷提升，用戶(hù)對(duì)高分辨率遙感數(shù)據(jù)的需求日益提升。目前我國(guó)已有包含高分、資源、環(huán)境減災(zāi)等多個(gè)系列在內(nèi)的陸地觀(guān)測(cè)衛(wèi)星近30顆，形成了多尺度、多觀(guān)測(cè)要素、多級(jí)別產(chǎn)品的業(yè)務(wù)化對(duì)地觀(guān)測(cè)體系，衛(wèi)星遙感影像產(chǎn)品分發(fā)服務(wù)機(jī)制也日趨完善[1-3]。

終端用戶(hù)訂購(gòu)獲取的衛(wèi)星遙感影像產(chǎn)品除影像數(shù)據(jù)外，通常還包括元數(shù)據(jù)、成像幾何模型參數(shù)文件、快視圖等[4]。其中快視圖是柵格影像經(jīng)拉伸增強(qiáng)后的真彩色或假彩色圖片，用于直接顯示影像數(shù)據(jù)呈現(xiàn)的地表環(huán)境狀態(tài)，如地理范圍、云雪覆蓋、水陸分布等，方便用戶(hù)快速了解影像質(zhì)量情況進(jìn)行選擇。據(jù)統(tǒng)計(jì)，早期的國(guó)產(chǎn)陸地衛(wèi)星遙感影像產(chǎn)品中的快視圖都沒(méi)有提供定位信息，近年來(lái)高分系列的PMS、資源一號(hào)02d、北京二號(hào)等少數(shù)產(chǎn)品的快視圖逐步增加了定位文件，其它如高分一號(hào)和高分六號(hào)的WFV、資源三號(hào)(01、02、03星)、高分七號(hào)、高景一號(hào)等遙感影像產(chǎn)品中快視圖仍然只是一張普通的JPG格式圖片，缺少定位信息，無(wú)法結(jié)合行政區(qū)劃或研究區(qū)矢量邊界范圍進(jìn)行針對(duì)目標(biāo)區(qū)域的快速瀏覽和數(shù)據(jù)篩選，在實(shí)際生產(chǎn)和應(yīng)用項(xiàng)目中較為不便。遙感產(chǎn)品的元數(shù)據(jù)文件記錄了影像產(chǎn)品的成像時(shí)間、太陽(yáng)角度、觀(guān)測(cè)角度、四至坐標(biāo)、像元大小或影像尺寸等基本信息[5]。本文在分析圖像坐標(biāo)定位文件結(jié)構(gòu)的基礎(chǔ)上，探討了利用元數(shù)據(jù)提供的影像數(shù)據(jù)坐標(biāo)信息計(jì)算定位參數(shù)的方法，并通過(guò)Python編程語(yǔ)言結(jié)合空間數(shù)據(jù)轉(zhuǎn)換庫(kù)(Geospatial Data Abstraction Library，GDAL)開(kāi)發(fā)了快視圖定位文件生成工具，以期為相關(guān)研究和生產(chǎn)作業(yè)人員提供參考。

1 研究方法

本文的研究數(shù)據(jù)源是幾種主要的國(guó)產(chǎn)中高分辨率光學(xué)遙感影像產(chǎn)品，其中高分、資源和環(huán)境等系列影像來(lái)源于中國(guó)資源衛(wèi)星應(yīng)用中心，高景影像由北京航天世景信息技術(shù)有限公司提供。

1.1 定位文件結(jié)構(gòu)分析

圖像以柵格數(shù)據(jù)形式存儲(chǔ)，柵格圖像中的每個(gè)像元都具有一個(gè)對(duì)應(yīng)于特定行號(hào)和列號(hào)的像素坐標(biāo)。要在正確的地理位置顯示圖像，需要建立特定的變換關(guān)系，以將像元坐標(biāo)轉(zhuǎn)換為真實(shí)世界坐標(biāo)，通常使用6參數(shù)構(gòu)建的仿射變換關(guān)系進(jìn)行轉(zhuǎn)換[6]：

(1)

寫(xiě)為方程組的形式：

x′=Ax+By+C

(2)

y′=Dx+Ey+F

(3)

式中，x、y分別為像素坐標(biāo)中的列號(hào)和行號(hào)，x′、y′為該像元對(duì)應(yīng)的地理坐標(biāo)；A、E分別為像素在地圖坐標(biāo)系X和Y方向上的尺寸，采用地圖單位；B和D為旋轉(zhuǎn)項(xiàng)；C和F分別左上角像素中心點(diǎn)的地圖坐標(biāo)。由于圖像坐標(biāo)系和地圖坐標(biāo)系的起點(diǎn)分別位于左上角和左下角，Y正方向相反，因此E為負(fù)值。

仿射變換信息可以存儲(chǔ)在圖像文件的文件頭內(nèi)部，也可以存儲(chǔ)在單獨(dú)的普通文本文件中，這個(gè)文本文件就是圖像坐標(biāo)定位文件，也稱(chēng)世界文件(world file)[7]，仿射變換的6個(gè)參數(shù)以十進(jìn)制數(shù)字的形式，按A、D、B、E、C、F的順序，依次分行存儲(chǔ)在文件中。坐標(biāo)定位文件的主文件名與圖像文件相同，后綴名分兩種：在圖像文件的擴(kuò)展名最后添加字母w，或使用圖像文件擴(kuò)展名中的第一個(gè)和第三個(gè)字符，最后再添加字母w，如jpg格式圖像的坐標(biāo)定位文件后綴名可以是jpgw或者jgw，通常后者更為常用。圖1為2021年1月11日某景高分一號(hào)PMS影像的坐標(biāo)定位文件內(nèi)容，產(chǎn)品號(hào)為L(zhǎng)1A0005387434。

圖1 圖像坐標(biāo)定位文件示例

1.2 定位參數(shù)解算原理

圖像定位校正與仿射變換參數(shù)的關(guān)系如圖2所示，X-O-Y為地圖坐標(biāo)系，影像產(chǎn)品的元數(shù)據(jù)中記錄了影像四角的坐標(biāo)信息，以此可以計(jì)算影像在X方向的距離w和Y方向的距離h，均為地圖單位：

圖2 圖像定位校正與仿射變換參數(shù)關(guān)系示意圖

ΔX1=xur-xul

(4)

ΔY1=yur-yul

(5)

(6)

同理，

ΔX2=xur-xlr

(7)

ΔY2=yur-ylr

(8)

(9)

并計(jì)算得到：

(10)

Rx=w/col

(11)

Ry=h/row

(12)

式中，xur、yur、xul、yur分別為圖像右上角和左上角的地理坐標(biāo)，row、col為圖像行數(shù)和列數(shù)，Rx、Ry為像元在X和Y方向上的尺寸，t為圖像旋轉(zhuǎn)角度，以X軸逆時(shí)針?lè)较驗(yàn)檎?。由此可以得到?/p>

A=Rx×cost

(13)

B=Ry×sint

(14)

D=Rx×sint

(15)

E=-Ry×cost

(16)

(17)

(18)

2 編程實(shí)現(xiàn)

遙感影像產(chǎn)品按處理的程度進(jìn)行分級(jí)，其中L0級(jí)和L1級(jí)產(chǎn)品未經(jīng)幾何校正，地理坐標(biāo)信息存儲(chǔ)在成像幾何模型參數(shù)文件中，影像文件本身不帶地理信息，大部分國(guó)產(chǎn)衛(wèi)星遙感影像如高分、資源系列等都按L1級(jí)產(chǎn)品分發(fā)；L2級(jí)產(chǎn)品是指經(jīng)過(guò)系統(tǒng)幾何校正的影像產(chǎn)品，影像文件帶有地理坐標(biāo)，如環(huán)境、中巴系列等[8]。為適應(yīng)不同級(jí)別遙感產(chǎn)品快視圖定位的需要，本文利用Python編程語(yǔ)言和GDAL轉(zhuǎn)換庫(kù)，從元數(shù)據(jù)解析出發(fā)，結(jié)合影像四角地理坐標(biāo)和快視圖圖像坐標(biāo)構(gòu)建控制點(diǎn)，快速計(jì)算仿射變換參數(shù)并生成坐標(biāo)定位文件，主要流程如圖3所示。

圖3 坐標(biāo)定位文件參數(shù)計(jì)算流程

2.1 元數(shù)據(jù)解析

元數(shù)據(jù)文件通常以可擴(kuò)展標(biāo)記語(yǔ)言(EXtensible Markup Language，XML)格式隨衛(wèi)星影像產(chǎn)品下發(fā)，XML是一種用于標(biāo)記電子文件使其具有結(jié)構(gòu)性的標(biāo)記語(yǔ)言，以純文本格式存儲(chǔ)，便于讀取、記錄、調(diào)試，使不同系統(tǒng)、不同程序之間的數(shù)據(jù)共享變得更加簡(jiǎn)單。Python自帶的Minidom庫(kù)可以方便地對(duì)XML文件進(jìn)行解析，提取出所需要的數(shù)據(jù)[9]，主要過(guò)程如圖4所示。至此即可獲取左上角坐標(biāo)(ulx、uly)、右上角坐標(biāo)(urx、ury)、左下角坐標(biāo)(llx、lly)、右下角坐標(biāo)(lrx、lry)共計(jì)四對(duì)值。

圖4 元數(shù)據(jù)解析主要代碼

2.2 地理變換參數(shù)解算

地理變換參數(shù)涉及到坐標(biāo)的縮放、旋轉(zhuǎn)和平移，由于地理校正后的影像存在扭曲、變形等原因，并非常規(guī)矩形，利用角點(diǎn)坐標(biāo)不同組合計(jì)算按理論方法計(jì)算得到的值會(huì)存在差異，存在多組已知點(diǎn)時(shí)，需要用最小二乘法解算最佳變換參數(shù)，計(jì)算過(guò)程較為繁瑣。GDAL是一個(gè)在X/MIT許可協(xié)議下的開(kāi)源柵格空間數(shù)據(jù)轉(zhuǎn)換庫(kù)，使用抽象數(shù)據(jù)模型來(lái)解析它所支持的數(shù)據(jù)格式，抽象數(shù)據(jù)模型包括數(shù)據(jù)集、坐標(biāo)系統(tǒng)、仿射地理坐標(biāo)轉(zhuǎn)換、大地控制點(diǎn)等，具有強(qiáng)大的地理數(shù)據(jù)處理能力[10]，利用其提供的函數(shù)可以通過(guò)已知控制點(diǎn)計(jì)算地理變換參數(shù)，具體方法如圖5所示。

圖5 地理變換參數(shù)解析主要代碼

3 結(jié)果評(píng)價(jià)

計(jì)算地理仿射變換參數(shù)并生成坐標(biāo)定位文件后，快視圖就可以在地理信息系統(tǒng)或遙感平臺(tái)軟件里進(jìn)行地理定位，在真實(shí)坐標(biāo)系統(tǒng)下顯示。本研究以產(chǎn)品號(hào)為GF1_WFV1_E121.6_N31.2_20211126_L1A0006081870的高分一號(hào)WFV影像產(chǎn)品為例，對(duì)快視圖地理校正的效果進(jìn)行評(píng)價(jià)。從目視對(duì)比看，恢復(fù)地理信息后的快視圖(左半幅)與實(shí)際影像數(shù)據(jù)(右半幅)的空間位置基本一致(圖6)。由于產(chǎn)品快視圖的分辨率在160 m左右，為影像分辨率的十分之一，無(wú)法準(zhǔn)確定位同名地物的位置，本文以影像四角像素的外側(cè)角點(diǎn)為檢查點(diǎn)，計(jì)算校正偏差(表1)，計(jì)算得到X、Y方向上偏差的中誤差分別為0.010 5、0.000 5(單位/°)，偏差主要體現(xiàn)在經(jīng)度方向上，約6個(gè)像素，該精度已能夠滿(mǎn)足粗略定位和范圍顯示的需要。

圖6 地理校正后的快視圖與影像文件疊加對(duì)比

表1 檢查點(diǎn)校正偏差(單位/°)

4 結(jié) 語(yǔ)

“十四五”期間，我國(guó)將有計(jì)劃地陸續(xù)發(fā)射多顆陸地觀(guān)測(cè)衛(wèi)星，進(jìn)一步推動(dòng)遙感業(yè)務(wù)化應(yīng)用向市場(chǎng)化、產(chǎn)業(yè)化發(fā)展。遙感影像產(chǎn)品直接面向終端用戶(hù)，包括影像文件、元數(shù)據(jù)、快視圖等內(nèi)容，其中快視圖是影像文件的縮略圖，為普通的圖片格式，便于用戶(hù)直接瀏覽，了解影像文件的質(zhì)量概況。本研究利用影像產(chǎn)品附帶元數(shù)據(jù)的四角坐標(biāo)信息構(gòu)建控制點(diǎn)，通過(guò)GDAL空間數(shù)據(jù)處理模塊計(jì)算坐標(biāo)定位參數(shù)，生成坐標(biāo)定位文件，實(shí)現(xiàn)對(duì)快視圖的地理校正，恢復(fù)快視圖的地理信息，為用戶(hù)進(jìn)行影像產(chǎn)品存檔、篩選、預(yù)覽、建庫(kù)等工作提供有效的幫助。