三軸穩(wěn)定型靜止軌道遙感衛(wèi)星指向確定的地標(biāo)匹配方法

丁 璐，童曉沖，秦志遠(yuǎn)，賴(lài)廣陵

1. 信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450001； 2. 河南省城建學(xué)院，河南 平頂山 467036

靜止軌道遙感衛(wèi)星與地球處于相對(duì)靜止?fàn)顟B(tài)，可以實(shí)現(xiàn)對(duì)地面目標(biāo)的持續(xù)觀測(cè)，同時(shí)具備動(dòng)態(tài)目標(biāo)探測(cè)能力和動(dòng)態(tài)目標(biāo)指示潛力，特別是近年來(lái)三軸穩(wěn)定型靜止軌道遙感衛(wèi)星的發(fā)展，對(duì)進(jìn)一步提升靜止軌道上持續(xù)觀測(cè)的時(shí)間和空間分辨率都具有重要的意義[1-4]。當(dāng)前，靜止軌道遙感衛(wèi)星已經(jīng)成為國(guó)家的基礎(chǔ)性、戰(zhàn)略性資源，美、歐、日、印、韓等國(guó)均投入了大量資源進(jìn)行研發(fā)。美國(guó)、日本于2016年11月分別將新一代三軸穩(wěn)定的靜止軌道氣象衛(wèi)星GOES-R[5]和“葵花9號(hào)”[6]發(fā)射升空；我國(guó)于2015年12月和2016年12月分別發(fā)射了具有自主產(chǎn)權(quán)的三軸穩(wěn)定型靜止軌道遙感衛(wèi)星“高分四號(hào)”[7]和“風(fēng)云四號(hào)”氣象衛(wèi)星[8]；歐洲下一代地球靜止軌道氣象衛(wèi)星MTG-I系列成像衛(wèi)星01星計(jì)劃于2021年發(fā)射升空[9]；印度也于2016年9月發(fā)射了一顆先進(jìn)的氣象衛(wèi)星INSAT-3DR[10]。

遙感衛(wèi)星在軌運(yùn)行時(shí)，受到內(nèi)部系統(tǒng)和外部環(huán)境因素的影響，會(huì)導(dǎo)致遙感儀器掃描鏡指向發(fā)生變化，使遙感圖像產(chǎn)生幾何失真，即指向偏差[11-13]。靜止軌道遙感衛(wèi)星屬于高軌衛(wèi)星，運(yùn)行在距離地球3.6萬(wàn)千米的地球靜止軌道上，衛(wèi)星的微小偏差在地面上會(huì)產(chǎn)生較大的地面誤差[14]，該問(wèn)題成為所有靜止軌道遙感衛(wèi)星不可規(guī)避的問(wèn)題[15]。為保證觀測(cè)結(jié)果的準(zhǔn)確性和可靠性，需要對(duì)衛(wèi)星的指向偏差進(jìn)行修正。目前，國(guó)內(nèi)外多采用地標(biāo)匹配的方法來(lái)獲取地標(biāo)偏移量，并由地標(biāo)偏移量確定導(dǎo)航調(diào)整量，以實(shí)現(xiàn)對(duì)遙感衛(wèi)星姿態(tài)的糾正[16-20]。其中存在的問(wèn)題在于現(xiàn)有的算法大多針對(duì)的是自旋穩(wěn)定式靜止軌道遙感衛(wèi)星和極軌衛(wèi)星獲取的數(shù)據(jù)，這些數(shù)據(jù)的空間分辨率普遍偏低，只能得到像元級(jí)的匹配結(jié)果，適用于整體剛性圖像的指向計(jì)算，其粗差剔除方法僅針對(duì)誤差分布均勻的情況，而且上述兩類(lèi)衛(wèi)星均不受熱形變的影響，不需要每次都進(jìn)行指向計(jì)算，只需要在間隔較長(zhǎng)的時(shí)間進(jìn)行一個(gè)系統(tǒng)差計(jì)算即可。而最近一兩年發(fā)射的三軸穩(wěn)定型高分辨率靜止軌道遙感衛(wèi)星對(duì)于地標(biāo)匹配能力提出了更高的要求，一方面要求地標(biāo)匹配能夠精度更高，穩(wěn)健性更強(qiáng)；另一方面為了達(dá)到更高定位精度，還需要具備局部剛性圖像的指向計(jì)算能力，特別是掃描成像的地球圓盤(pán)數(shù)據(jù)。除此之外，靜止軌道遙感衛(wèi)星在拍攝地球圓盤(pán)的過(guò)程中，由于不同圓盤(pán)區(qū)域受到云層、光照或者地標(biāo)位置區(qū)影像特征、地標(biāo)本身幾何形態(tài)等問(wèn)題的影響，導(dǎo)致不同地標(biāo)會(huì)產(chǎn)生匹配質(zhì)量的差異，如果直接使用這些地標(biāo)匹配結(jié)果，將會(huì)給指向計(jì)算帶來(lái)較大的偏差，因此需要在指向偏差計(jì)算中考慮觀測(cè)量的粗差剔除、不等權(quán)處理等問(wèn)題，保證算法的穩(wěn)健性，使其能夠適應(yīng)一天內(nèi)多種類(lèi)型成像條件下的圖像計(jì)算問(wèn)題。本文圍繞三軸穩(wěn)定型靜止軌道遙感衛(wèi)星拍攝地球圓盤(pán)時(shí)的指向確定，針對(duì)現(xiàn)有方法存在的問(wèn)題以及新型傳感器的成像特點(diǎn)，提出了一套基于地標(biāo)匹配的指向計(jì)算方法，并已經(jīng)在我國(guó)相關(guān)衛(wèi)星數(shù)據(jù)的在線處理中得到應(yīng)用。

1 靜止軌道遙感衛(wèi)星指向確定的地標(biāo)匹配方法

1.1 基本思路

實(shí)現(xiàn)三軸穩(wěn)定型靜止軌道遙感衛(wèi)星的指向偏差計(jì)算，主要是利用地球表面固定、明顯的地標(biāo)進(jìn)行匹配，由于靜止軌道遙感衛(wèi)星的分辨率普遍在百米量級(jí)(例如：美國(guó)GOES-R最高分辨率為星下點(diǎn)500 m[5])，而一般性的地標(biāo)物中，能較均勻地覆蓋整個(gè)地球圓盤(pán)并且能夠被大多數(shù)可見(jiàn)光、紅外波段觀測(cè)到的地標(biāo)物只有海陸分界線，因此靜止軌道遙感圖像的地標(biāo)參照應(yīng)選擇地球圓盤(pán)上相對(duì)固定的海陸分界線上的特征區(qū)域。結(jié)合以上分析，靜止軌道遙感衛(wèi)星指向確定的地標(biāo)匹配方法主要包括：數(shù)據(jù)準(zhǔn)備、地標(biāo)匹配及指向計(jì)算3個(gè)部分。具體為：

(1) 數(shù)據(jù)準(zhǔn)備：按照靜止軌道遙感衛(wèi)星成像方式的不同，制作與成像方式類(lèi)似的地標(biāo)模板數(shù)據(jù)，其中涉及地標(biāo)數(shù)據(jù)的選擇、成像投影網(wǎng)格的計(jì)算、地標(biāo)數(shù)據(jù)的投影與重采樣、根據(jù)地標(biāo)形態(tài)和分布形成地標(biāo)特征區(qū)域和地標(biāo)庫(kù)。

(2) 地標(biāo)匹配：根據(jù)實(shí)時(shí)獲取的靜止軌道遙感數(shù)據(jù)，對(duì)整個(gè)圓盤(pán)區(qū)域分布的地標(biāo)特征區(qū)域進(jìn)行云篩選，選擇無(wú)云或者少云的區(qū)域作為待匹配區(qū)，利用海陸模板形成的地標(biāo)區(qū)域和遙感圖像處理后得到的邊界信息進(jìn)行匹配處理，選取合適的整體特征(梯度和)完成粗匹配，并得到像平面的幾何變形參數(shù)，在此基礎(chǔ)上，針對(duì)局部特征點(diǎn)的復(fù)合特征(灰度相關(guān)、形態(tài)相關(guān)、矩相關(guān))完成最小二乘精匹配。

(3) 指向計(jì)算：將地標(biāo)匹配的結(jié)果用于指向偏差角的計(jì)算，利用粗差剔除和抗差估計(jì)的方法，得到指向角計(jì)算結(jié)果，并針對(duì)晝夜各半的特殊圓盤(pán)數(shù)據(jù)進(jìn)行多波段合成處理等工作。

總體技術(shù)路線如圖1所示。

1.2 地標(biāo)數(shù)據(jù)準(zhǔn)備

按照三軸穩(wěn)定型靜止軌道遙感衛(wèi)星的成像規(guī)律，拍攝地球表面的區(qū)域數(shù)據(jù)主要有兩種方式：一是畫(huà)幅式成像，即采用面陣CCD對(duì)某區(qū)域進(jìn)行一次性成像，這種情況往往面陣較大，成像投影方式采用的是面陣中心投影的方法，例如高分四號(hào)的對(duì)地成像[4，21]；二是掃描鏡成像，為了快速獲取大面積的遙感數(shù)據(jù)，特別是地球完整圓盤(pán)數(shù)據(jù)，常采用線陣掃描成像方式，利用東西、南北兩個(gè)方向的快速掃描與步進(jìn)完成大幅面區(qū)域的覆蓋[22]，這種成像方式常采用的是規(guī)范化地球靜止投影(normalized geostationary projection，NGP)[23]。因此，在地標(biāo)數(shù)據(jù)準(zhǔn)備的第1步，需要根據(jù)衛(wèi)星的成像特點(diǎn)，選擇恰當(dāng)?shù)耐队胺绞剑⒃谔囟ǖ耐队胺绞较?，形成作為參考基?zhǔn)的標(biāo)稱(chēng)網(wǎng)格數(shù)據(jù)。

圖1 總體技術(shù)路線流程Fig.1 Flowchart on general technical roadmap

第2步，根據(jù)矢量的海陸邊界數(shù)據(jù)，結(jié)合成像投影，生成柵格化的海陸邊界掩碼數(shù)據(jù)。海陸掩碼需要結(jié)合全球海岸線矢量數(shù)據(jù)和標(biāo)稱(chēng)網(wǎng)格數(shù)據(jù)計(jì)算獲得。其基本流程如圖2所示。

圖2 海陸掩碼文件生成流程Fig.2 Flowchart on the generation of sea-land mask file

第3步，建立地標(biāo)數(shù)據(jù)集。自動(dòng)由海陸邊界模板生成地標(biāo)控制點(diǎn)，以地標(biāo)控制點(diǎn)為中心的一定大小的區(qū)域就是地標(biāo)，所有地標(biāo)構(gòu)成地標(biāo)數(shù)據(jù)集。該部分的技術(shù)難點(diǎn)在于地標(biāo)控制點(diǎn)檢測(cè)算法以及地標(biāo)的分布等問(wèn)題。檢測(cè)海陸邊界模板中的地標(biāo)控制點(diǎn)，地標(biāo)控制點(diǎn)一般選取海陸邊界上的角點(diǎn)、“檢”型連接點(diǎn)或者其他高曲率點(diǎn)。

1.3 圖像地標(biāo)匹配

本節(jié)重點(diǎn)討論圖像與地標(biāo)的匹配問(wèn)題，具體的步驟如下。

1.3.1 匹配范圍的確定

首先是匹配初值的問(wèn)題，目前三軸穩(wěn)定型靜止軌道遙感衛(wèi)星，利用星上傳感器測(cè)量數(shù)據(jù)進(jìn)行直接定位會(huì)帶來(lái)大約數(shù)十公里的偏差，也就是說(shuō)，折算成特定分辨率的圖像(以500 m分辨率為例)，通過(guò)定位計(jì)算得到的圖像上像素的幾何位置與該像素的理想幾何位置(標(biāo)稱(chēng)網(wǎng)格位置)之間會(huì)有十幾個(gè)到幾十個(gè)像元的差異，那么可以以該誤差的上限作為匹配初值，設(shè)定搜索的半徑R。

1.3.2 對(duì)待匹配區(qū)域進(jìn)行云判處理

衛(wèi)星觀測(cè)的遙感數(shù)據(jù)中有大量的云覆蓋于地表信息之上，對(duì)后續(xù)的海陸邊緣檢測(cè)及地標(biāo)匹配帶來(lái)影響，因此必須進(jìn)行云檢測(cè)，對(duì)云覆蓋區(qū)域不進(jìn)行地標(biāo)匹配處理。在處理過(guò)程中，引入了最小交叉熵的概念，采用基于最小交叉熵的云區(qū)域檢測(cè)方法[24]，獲取一幅二值的云掩膜圖像M(x,y)。

1.3.3 圖像的粗匹配與變形參數(shù)解算

在靜止軌道遙感數(shù)據(jù)和地標(biāo)的匹配過(guò)程中，由于圖像和地標(biāo)數(shù)據(jù)一個(gè)是灰度圖，一個(gè)是黑白柵格圖，二者的相似度較低，直接進(jìn)行灰度特征相關(guān)很難達(dá)到可靠的相關(guān)結(jié)果，因此需要尋找其他的測(cè)度。對(duì)于地標(biāo)數(shù)據(jù)，圖像上與其最相似的實(shí)際上就是灰度的梯度信息，也就是圖像上的邊界。

需要說(shuō)明的是，靜止軌道遙感圖像f(x,y)常見(jiàn)的方式是采用單中心小面陣拼接的方式，使用局部較大區(qū)域進(jìn)行匹配存在幾何上的形變，該形變直接影響匹配的精度與可靠性。另外，即使在地標(biāo)數(shù)據(jù)和圖像數(shù)據(jù)保持同樣投影的過(guò)程中，仍然存在較小的空間幾何變換關(guān)系導(dǎo)致差異的產(chǎn)生，因此在匹配的過(guò)程中需要同步解算像面幾何形變與空間的變換，才能使地標(biāo)和圖像可靠地匹配上。在匹配測(cè)度上，采用梯度信息最大為目標(biāo)，與大部分匹配以圖像與基礎(chǔ)，地標(biāo)為目標(biāo)進(jìn)行匹配相反，該方法以地標(biāo)D(x,y)為基礎(chǔ)，圖像為目標(biāo)進(jìn)行匹配，其根本原因是最終需要在圖像上形成梯度。該方法的步驟包括：

(1)

需要說(shuō)明的是，(i′,j′)如果對(duì)應(yīng)的圖像是有云區(qū)域，需要使用云掩膜M(x,y)處理，D對(duì)應(yīng)的區(qū)域是選擇的待匹配無(wú)云地標(biāo)區(qū)域的全集。

(2) 給出幾何變換F的變換形式。由于該問(wèn)題是以地標(biāo)D(x，y)為基礎(chǔ)，圖像為目標(biāo)的匹配，因此可以建立一個(gè)虛擬成像的幾何過(guò)程，即從地標(biāo)格網(wǎng)到真實(shí)圖像的虛擬投影。下面給出的是基本思路：

糾正后格網(wǎng)坐標(biāo)通過(guò)若干次旋轉(zhuǎn)后，綜合旋轉(zhuǎn)矩陣是R，由式(2)

(2)

得到真實(shí)像空間坐標(biāo)系下的坐標(biāo)。其中，綜合旋轉(zhuǎn)矩陣R是遙感圖像定位過(guò)程中，從真實(shí)像空間到衛(wèi)星本體坐標(biāo)系綜合旋轉(zhuǎn)矩陣R′的逆矩陣(R′)-1，在每次成像過(guò)程中是已知量。

(3)

(3) 對(duì)式(1)進(jìn)行優(yōu)化，采用LM算法(Levenberg-Marquard)[25]進(jìn)行最大值的求解與計(jì)算，求解參數(shù)[αβγa0a1a2a3a4a5b0b1

b2b3b4b5]′(共15個(gè)參數(shù))和綜合旋轉(zhuǎn)矩陣R。其中，初值α0=β0=γ0=0，[a0a1a2a3a4a5

b0b1b2b3b4b5]′=[0 1 0 0 0 0 0 0 0 1 0

0 0]′。

該過(guò)程的核心是計(jì)算待求解變量的偏導(dǎo)數(shù)問(wèn)題，首先需要將[x′y′]T與[XYZ]T建立等式關(guān)系，由式(3)得到式(4)

(4)

式中，D′(·)運(yùn)算表示從真實(shí)的像空間坐標(biāo)系到真實(shí)像平面坐標(biāo)系的轉(zhuǎn)換，再對(duì)式(4)進(jìn)行線性化，過(guò)程略。利用LM算法進(jìn)行迭代優(yōu)化，其中初值

α0=β0=γ0=0

[a0a1a2a3a4a5b0b1b2b3b4b5]′=

[0 1 0 0 0 0 0 0 1 0 0 0]′

(4) 根據(jù)整體匹配的結(jié)果以及計(jì)算出的幾何變換參數(shù)，將每一個(gè)地標(biāo)的網(wǎng)格坐標(biāo)(x,y)代入式(3)可以得到對(duì)應(yīng)的圖像坐標(biāo)(i′,i′)，形成粗匹配的匹配對(duì)。

1.3.4 遙感圖像邊界提取與精匹配

遙感圖像中存在各種地物的紋理特征，邊界特征不明顯，無(wú)法直接與地標(biāo)模板進(jìn)行直接匹配，因此，需要對(duì)遙感圖像進(jìn)行邊緣提取。對(duì)灰度圖像的邊緣提取較常采用梯度算子，運(yùn)算效率高且適用于全圓盤(pán)大范圍的遙感圖像，但梯度算子容易出現(xiàn)漏檢邊緣，即間斷的地方。還有方法采用邊界跟蹤來(lái)提取邊緣，但得到的邊界點(diǎn)就會(huì)在角點(diǎn)附近出現(xiàn)大段缺失，因而丟失許多有用的邊界信息[26]。結(jié)合兩種方法，研究中采用圖像相加的方式，將兩種方法處理得到的邊界圖進(jìn)行疊加，從而得到一幅效果更加良好的邊界圖。

在邊界圖像和粗匹配的基礎(chǔ)上，本文設(shè)計(jì)了遙感圖像與地標(biāo)的精匹配方法，精匹配的思路采用最小二乘匹配，以梯度與海陸地標(biāo)相似性及梯度方向的相關(guān)性作為匹配測(cè)度。下面給出精確匹配的步驟：

(1) 針對(duì)選擇的多塊地標(biāo)特征區(qū)域的海陸掩碼圖像，采用高斯模糊的方式，將海陸掩碼的0-1圖像進(jìn)行模糊處理，得到0～1之間的浮點(diǎn)數(shù)組成的浮點(diǎn)矩陣。

該步驟是最小二乘方法精匹配能夠收斂的關(guān)鍵性問(wèn)題之一，因?yàn)楹ｊ懷诖a圖像本身是0-1圖像，像素之間的數(shù)值并不連續(xù)，而遙感圖像梯度所形成的圖像是連續(xù)變化的，兩者進(jìn)行相關(guān)得到的相關(guān)系數(shù)也是一個(gè)不連續(xù)的曲面，而最小二乘圖像相關(guān)的關(guān)鍵是相關(guān)函數(shù)是連續(xù)的，不連續(xù)的相關(guān)函數(shù)將直接影響到最小二乘算法的收斂程度。采用高斯模糊的方式將0-1圖像轉(zhuǎn)變成相對(duì)連續(xù)的浮點(diǎn)圖像，可以保證相關(guān)函數(shù)的連續(xù)性，使得最小二乘精匹配可以收斂。

(2) 由于遙感圖像在作邊緣檢測(cè)的過(guò)程中，容易受到各種噪聲和地物變化的影響，采用單一的測(cè)度存在可靠性較低的情況，因此為了增加匹配的可靠性，本項(xiàng)研究引入了多重判據(jù)匹配的思路，其目的是希望利用多個(gè)測(cè)度得到一致性的結(jié)果，這樣匹配結(jié)果的可靠性大大增強(qiáng)。除了灰度相關(guān)系數(shù)(適合圖像灰度分布特征)匹配方法外，此處還引入Zernike矩差異測(cè)度[27](適合圖像統(tǒng)計(jì)的匹配測(cè)度)和誤差橢圓之差[28](適合形態(tài)的匹配測(cè)度)作為匹配測(cè)度，其目的是由于這種測(cè)度與相關(guān)系數(shù)測(cè)度有所不同，并且地標(biāo)特征點(diǎn)與邊緣圖像點(diǎn)的誤差橢圓應(yīng)該是接近的，因此選擇這3種測(cè)度作為匹配測(cè)度。

考慮到相關(guān)系數(shù)測(cè)度的優(yōu)良特性，本文將其作為主要的測(cè)度，其他兩個(gè)測(cè)度作為一種有力的補(bǔ)充，在相關(guān)系數(shù)測(cè)度得到的結(jié)果可靠性不高的情況下尋找更為可靠的結(jié)果。相關(guān)系數(shù)測(cè)度結(jié)果的可靠性依據(jù)相關(guān)系數(shù)值來(lái)確定，需要給定相關(guān)系數(shù)一個(gè)最大閾值Tmax和一個(gè)最小閾值Tmin，對(duì)任意點(diǎn)的匹配過(guò)程如下：

步驟1：在給定大小窗口下計(jì)算相關(guān)系數(shù)γ；

步驟2：對(duì)γ進(jìn)行判斷，若γ≥Tmax，則認(rèn)為匹配結(jié)果可靠，接受該結(jié)果；若γ≤Tmin，則認(rèn)為匹配結(jié)果錯(cuò)誤，將其舍棄；若Tmin<γ

步驟3：加入Zernike矩差異與誤差橢圓之差兩個(gè)測(cè)度，并將三者的結(jié)果進(jìn)行比較，如果三者結(jié)果一致，則接受它；如果三者結(jié)果不一致，則逐漸擴(kuò)大目標(biāo)窗口并重復(fù)步驟1—步驟3；

步驟4：如果到目標(biāo)窗口擴(kuò)大終止仍得不到可接受的結(jié)果，則將其標(biāo)示為可疑點(diǎn)，進(jìn)入下一步的匹配。

設(shè)計(jì)這樣的匹配過(guò)程，是因?yàn)椴捎枚嘀嘏袚?jù)匹配會(huì)帶來(lái)時(shí)間的消耗，而對(duì)于那些相關(guān)系數(shù)較大或較小的點(diǎn)，無(wú)需進(jìn)一步判斷和計(jì)算，僅需要對(duì)那些相關(guān)系數(shù)處于中間的點(diǎn)進(jìn)行多重判據(jù)匹配，這樣在提高匹配成功率的基礎(chǔ)上有效地控制了計(jì)算效率。

1.4 指向偏差計(jì)算方法

經(jīng)過(guò)上述地標(biāo)匹配工作，將得到一系列均勻分布在圖像上的地標(biāo)匹配點(diǎn)對(duì)，使用這些點(diǎn)對(duì)可以直接解算光軸的指向偏差。文獻(xiàn)[17]中給出了風(fēng)云二號(hào)衛(wèi)星的指向計(jì)算方法，在上一代自旋穩(wěn)定的平臺(tái)上，對(duì)地指向精度較高，利用地標(biāo)指向修正并不是一項(xiàng)經(jīng)常性的工作，常采用人工檢查地標(biāo)匹配點(diǎn)后再計(jì)算。而新一代三軸穩(wěn)定的靜止遙感衛(wèi)星由于指向精度受到更多因素影響[29]，這項(xiàng)處理變得更加頻繁。由于靜止軌道遙感圖像的匹配本質(zhì)上是圖像與矢量邊界的匹配，受到圖像中線狀地物的影響較大，特別是受到復(fù)雜的云跡、薄云等現(xiàn)象的影響，而這些云極難被云檢測(cè)算法準(zhǔn)確定位，因此地標(biāo)匹配很難保證嚴(yán)格的準(zhǔn)確性。另外，由于圖像和地標(biāo)格網(wǎng)本身并不是簡(jiǎn)單的平移旋轉(zhuǎn)關(guān)系，還存在空間的變換，因此直接使用簡(jiǎn)單的地標(biāo)偏差統(tǒng)計(jì)方法進(jìn)行剔除(例如2σ的方式)很難得到理想的效果。

本文工作的試驗(yàn)表明，由于地標(biāo)匹配精度受圖像中復(fù)雜云場(chǎng)景的變化影響，無(wú)法做到一天內(nèi)每次觀察都能穩(wěn)定給出高質(zhì)量的匹配結(jié)果，可能出現(xiàn)許多不同等級(jí)誤差的匹配點(diǎn)，而圖像匹配中常使用的RANSAC方法[30]剔除匹配誤差的方式，在直接使用效果并不理想。為了得到穩(wěn)健的指向偏差計(jì)算方法，本文結(jié)合RANSAC方法、抗差估計(jì)[31]與擬合外推方法，設(shè)計(jì)了下面的思路：

(5)

(2) 使用Ransac算法對(duì)匹配點(diǎn)對(duì)的粗差進(jìn)行第1次剔除，其中約束的幾何模式采用式(5)的方程。根據(jù)幾何關(guān)系與試驗(yàn)給出設(shè)定參數(shù)：適用于模型的最小數(shù)據(jù)集數(shù)據(jù)個(gè)數(shù)為4個(gè)，判定新的數(shù)據(jù)是否適用于數(shù)據(jù)集模型的標(biāo)準(zhǔn)是代入后偏差(均方差)小于數(shù)據(jù)集殘差的2倍(2σ)，新的數(shù)據(jù)適用于數(shù)據(jù)集模型的數(shù)據(jù)數(shù)目占總匹配數(shù)目的60%，迭代次數(shù)設(shè)定為100次。

(3) 由于海陸模板與靜止軌道圖像匹配的復(fù)雜性，會(huì)出現(xiàn)許多不同等級(jí)誤差的匹配點(diǎn)，針對(duì)該問(wèn)題本文在第2層次誤差處理的過(guò)程中，采用將粗差歸入隨機(jī)模型的粗差定位法[32](基于選權(quán)迭代的抗差估計(jì)方法)，試驗(yàn)表明，使用IGGIII方案[33-35]給出的等價(jià)權(quán)函數(shù)，如式(6)，可以達(dá)到較好的粗差剔除效果

(6)

(4) 第3層次的工作是為了保證長(zhǎng)時(shí)間多次指向計(jì)算的結(jié)果的穩(wěn)定，避免出現(xiàn)由于地標(biāo)匹配誤差而導(dǎo)致的指向計(jì)算不穩(wěn)定的現(xiàn)象。在該過(guò)程中，考慮到指向偏差變化與軌道位置、溫差變化的物理因素相關(guān)，變化本身具有一定周期特性，因此采用Fourior擬合的方式對(duì)指向偏差角的三軸(α，β，γ)分別進(jìn)行擬合處理，擬合的數(shù)據(jù)來(lái)自本次解算之前的若干次數(shù)據(jù)，然后外推本次的三軸，再與直接解算的結(jié)果進(jìn)行對(duì)比，設(shè)置一定閾值，大于閾值范圍認(rèn)為本次計(jì)算失敗，反之則成功，并將本次解算三軸的結(jié)果存入數(shù)據(jù)庫(kù)中，為下一次解算結(jié)果的驗(yàn)證提供依據(jù)。

在實(shí)現(xiàn)中，以軌道周期(24H)為基頻，進(jìn)行1、2、4、8、16、32倍頻共6階傅里葉函數(shù)的擬合，如式(7)，其中共需要解算a0—a6，b1—b6，c0—c6，d1—d6，e0—e6，f1—f6，ω共40個(gè)參數(shù)，因此至少需要之前14次觀測(cè)的結(jié)果，才能擬合外推本次的結(jié)果。方程式(7)是非線性的方程，本文采用LM算法進(jìn)行優(yōu)化，得到擬合參數(shù)

(7)

2 試驗(yàn)與分析

試驗(yàn)以某靜止軌道遙感衛(wèi)星數(shù)據(jù)為例來(lái)驗(yàn)證算法。試驗(yàn)選擇該衛(wèi)星多波段掃描輻射計(jì)獲取的第2波段(500 m分辨率，550 nm波段)數(shù)據(jù)進(jìn)行試驗(yàn)處理，選擇某日4幅不同時(shí)段(UTC時(shí)間)的全圓盤(pán)成像數(shù)據(jù)為例，如圖3所示。

圖3 一天內(nèi)不同時(shí)段的全圓盤(pán)數(shù)據(jù)(UTC時(shí)間)Fig.3 Full-disc data in different time periods within one day (UTC)

主要完成下面兩個(gè)方面的試驗(yàn)：

(1) 標(biāo)稱(chēng)數(shù)據(jù)及標(biāo)稱(chēng)模板生成：針對(duì)該衛(wèi)星多波段掃描輻射計(jì)的成像掃描方式和定點(diǎn)位置，生成全球標(biāo)稱(chēng)格網(wǎng)數(shù)據(jù)，利用全球海陸矢量邊界數(shù)據(jù)，生成地標(biāo)模板數(shù)據(jù)。

(2) 地標(biāo)匹配與指向計(jì)算：通過(guò)獲取的4幅全圓盤(pán)數(shù)據(jù)，開(kāi)展地標(biāo)匹配試驗(yàn)，并利用地標(biāo)匹配結(jié)果解算指向偏差，最后利用解算后的指向偏差，對(duì)影像進(jìn)行重新糾正并重采，檢查與標(biāo)準(zhǔn)地標(biāo)數(shù)據(jù)的偏差。

2.1 標(biāo)稱(chēng)數(shù)據(jù)及標(biāo)稱(chēng)模板生成

試驗(yàn)采用適合于該型衛(wèi)星成像規(guī)律的標(biāo)稱(chēng)格網(wǎng)數(shù)據(jù)(NPG投影的一種)，其成像規(guī)律是利用南北、東西兩個(gè)方向的掃描鏡的機(jī)械掃描來(lái)實(shí)現(xiàn)地面區(qū)域的覆蓋，而標(biāo)稱(chēng)格網(wǎng)的計(jì)算核心是其入射光線的計(jì)算，式(8)給出了入射光線的形式

(8)

式中，α是南北鏡轉(zhuǎn)角；β是東西鏡轉(zhuǎn)角。

以衛(wèi)星定點(diǎn)位置104.7°E為例，結(jié)合式(8)的光線方程，形成標(biāo)稱(chēng)格網(wǎng)的數(shù)據(jù)(格網(wǎng)點(diǎn)對(duì)應(yīng)的地面點(diǎn)大地坐標(biāo))，其經(jīng)緯度范圍為：緯度-81.253 585°S—81.253 585°N，經(jīng)度23.421 845°E—174.021 845°W，尺寸21 984×21 984(500 m分辨率)，該格網(wǎng)為地標(biāo)模板生成提供空間參考基礎(chǔ)。

試驗(yàn)中所用的矢量海岸線數(shù)據(jù)是GSHHG數(shù)據(jù)(Version 2.3.1，2014)，GSHHG數(shù)據(jù)包括世界海岸線、河流、湖泊、島嶼、冰凍圈等水文要素內(nèi)容[36]。GSHHG數(shù)據(jù)坐標(biāo)系為WGS-84坐標(biāo)系，矢量節(jié)點(diǎn)的坐標(biāo)為大地經(jīng)緯度坐標(biāo)，數(shù)據(jù)包括5種精度的數(shù)據(jù)，分別以f、h、i、l、c表示，其含義為：全分辨率(full resolution)、高分辨率(high resolution)、中分辨率(intermediate resolution)、低分辨率(low resolution)、粗分辨率(rude resolution)。

在數(shù)據(jù)相關(guān)說(shuō)明中并沒(méi)有對(duì)具體5種數(shù)據(jù)的具體比例尺及精度進(jìn)行說(shuō)明，以該衛(wèi)星的多波段掃描輻射計(jì)為例，其中包括：500 m、1 km、2 km、4 km 4種不同的分辨率，在生成標(biāo)稱(chēng)格網(wǎng)時(shí)，同時(shí)也生成對(duì)應(yīng)4種分辨率的格網(wǎng)。對(duì)不同分辨率標(biāo)稱(chēng)格網(wǎng)對(duì)應(yīng)的海陸邊界模板采用不同精度的矢量海岸線數(shù)據(jù)，其中500 m、1000 m分辨率的海陸邊界用f精度數(shù)據(jù)，2000 m分辨率的海陸邊界采用h精度數(shù)據(jù)，4000 m、8000 m分辨率的海陸邊界采用i精度數(shù)據(jù)。同時(shí)各精度數(shù)據(jù)只應(yīng)用觀測(cè)范圍內(nèi)幾大洲的海岸線，以及面積較大的湖泊、島嶼的數(shù)據(jù)用于地標(biāo)生成。圖4是通過(guò)重采樣方法計(jì)算的不同分辨率標(biāo)稱(chēng)格網(wǎng)下的海陸掩碼圖像數(shù)據(jù)，這些數(shù)據(jù)是地標(biāo)匹配的數(shù)據(jù)基礎(chǔ)。圖5是500 m標(biāo)稱(chēng)格海陸掩碼圖像上的特征提取結(jié)果(“+”)，并根據(jù)全圓盤(pán)范圍內(nèi)的海岸線分布以及常年的云層分布特點(diǎn)，形成了全球多塊不同的地標(biāo)區(qū)域。

2.2 地標(biāo)匹配與指向計(jì)算

采用地標(biāo)匹配的方法，進(jìn)行靜止軌道遙感數(shù)據(jù)與地標(biāo)區(qū)域的匹配，得到相關(guān)的匹配結(jié)果。以圖3(c)中的數(shù)據(jù)為例，主要經(jīng)過(guò)下面幾個(gè)步驟：

(1) 圖像的粗匹配與變形參數(shù)解算。進(jìn)行圖像按區(qū)域的整體匹配與變形參數(shù)解算，得到圖6的結(jié)果，其中圖6(a)中是不解算變形參數(shù)，直接按照式(1)進(jìn)行匹配的結(jié)果，可以發(fā)現(xiàn)圖中下半部分的海陸邊界基本匹配較好，但是“○”標(biāo)注的湖泊，發(fā)生了偏差；(b)中是匹配過(guò)程中，直接解算變形參數(shù)的方式，可以發(fā)現(xiàn)無(wú)論是下半部海陸邊界還是湖泊，都匹配得很好。

圖4 不同分辨率標(biāo)稱(chēng)格數(shù)據(jù)形成的海陸掩碼圖像Fig.4 Sea-land mask images based on the nominal grid with a varying resolution

圖5 500 m標(biāo)稱(chēng)格海陸掩碼圖像上的特征提取結(jié)果與地標(biāo)區(qū)分布Fig.5 Results of characteristic extraction and distribution of landmark regions in the sea-land mask image based on the nominal grid with the resolution of 500 m

圖6 圖像的粗匹配與變形參數(shù)解算過(guò)程中的對(duì)比(局部：青藏高原)Fig.6 Rough image matching and comparison in the calculation of deformation parameters (local：tibet plateau)

在粗匹配的基礎(chǔ)上，完成了地標(biāo)精確匹配的過(guò)程，精匹配是按照?qǐng)D5(b)的全球地標(biāo)區(qū)域進(jìn)行匹配的，全球范圍內(nèi)共18塊區(qū)域，圖7是其中兩塊區(qū)域的精確匹配的結(jié)果。

圖8是局部區(qū)域利用地標(biāo)糾正指向后的圖像結(jié)果與糾正前的圖像定位結(jié)果的對(duì)比；圖9是一天內(nèi)不同時(shí)次(橫坐標(biāo)是UTC時(shí)間)圖像(78次)利用地標(biāo)進(jìn)行糾正后再進(jìn)行檢測(cè)的結(jié)果，一天內(nèi)有78次數(shù)據(jù)觀測(cè)任務(wù)，每次都需要批量進(jìn)行地標(biāo)處理與糾正，由于數(shù)據(jù)變化的周期是以天為單位進(jìn)行變化的，因此分析一天內(nèi)的糾正情況是合理的，統(tǒng)計(jì)的內(nèi)容包括不同時(shí)次的平均偏差和最大偏差(每一個(gè)時(shí)次觀測(cè)圖像都會(huì)有若干個(gè)地標(biāo)匹配點(diǎn)對(duì))，其中一天內(nèi)平均偏差的均值：X方向0.690 325 874像元，Y方向0.677 241 358個(gè)像元；一天內(nèi)最大偏差的均值X方向3.395 537 354個(gè)像元，Y方向2.626 338 072個(gè)像元；圖10是一天內(nèi)利用地標(biāo)進(jìn)行圖像糾正得到的指向偏差的統(tǒng)計(jì)分析結(jié)果，是在圖9的基礎(chǔ)上，對(duì)所有任務(wù)中，每次任務(wù)所有地標(biāo)點(diǎn)檢測(cè)偏差均值的直方圖(X/Y方向)和偏差最大值的直方圖(X/Y方向)。

需要說(shuō)明的是，圖8(a)和(b)是沒(méi)有經(jīng)過(guò)地標(biāo)指向糾正的原始圖像，通過(guò)地標(biāo)匹配(1.3.3節(jié)、1.3.4節(jié))，并經(jīng)過(guò)選權(quán)迭代的指向偏差計(jì)算(1.4節(jié))，一方面對(duì)匹配地標(biāo)的進(jìn)行粗差剔除，另一方面同步得到指向偏差的三軸關(guān)系(如式(5)中的指向偏差角(α,β,γ))；圖8(c)和(d)是在指向角的基礎(chǔ)上，再將圖像逐點(diǎn)按照指向角進(jìn)行空間相似變換(利用式(5))，完成圖像的重采樣工作后，再與標(biāo)準(zhǔn)地標(biāo)數(shù)據(jù)進(jìn)行套合，得到的對(duì)比結(jié)果。另外，圖9中一天內(nèi)偏差的統(tǒng)計(jì)，是對(duì)每一次任務(wù)都完成上述的地標(biāo)匹配與指向計(jì)算，形成外推模式(式(7))，每一次任務(wù)使用前面多次指向角擬合后外推的指向角再進(jìn)行糾正，然后本次任務(wù)再進(jìn)行本次地標(biāo)匹配并得到本次指向角，加入到指向角序列中，為下一次任務(wù)外推服務(wù)，其中本次地標(biāo)匹配的結(jié)果即為本次的偏差統(tǒng)計(jì)量，整個(gè)流程是一個(gè)完整迭代的過(guò)程。圖10展示的是偏差的分布情況，從中可以看出，一天內(nèi)的每一次任務(wù)的偏差均值/最大值基本呈現(xiàn)正態(tài)分布的特征，其中，均值表現(xiàn)出更好的正態(tài)特性，最大值相對(duì)更加發(fā)散，具體的統(tǒng)計(jì)值(均值/均方差)如圖10所示。

從上述試驗(yàn)可以看出，本文提出的基于地標(biāo)匹配的指向偏差計(jì)算方法，能夠較好地對(duì)靜止軌道遙感數(shù)據(jù)進(jìn)行糾正，糾正后平均差異在0.6個(gè)像元左右，每次最大誤差平均在3個(gè)像元左右。通過(guò)對(duì)一天內(nèi)78次全球圓盤(pán)觀測(cè)的試驗(yàn)結(jié)果進(jìn)行分析發(fā)現(xiàn)，本文提出的算法能夠較穩(wěn)定地保證圖像的糾正精度，算法整體的穩(wěn)定度較好，但由于靜止軌道遙感數(shù)據(jù)的全圓盤(pán)圖像是采用拼圖方式完成，受到機(jī)械控制精度等各方面因素的影響，很難保證圖像內(nèi)部的嚴(yán)格幾何關(guān)系，因此在匹配過(guò)程中會(huì)出現(xiàn)全圖中不同區(qū)域受到同一組全局指向偏差糾正后，糾正的精度存在差異的情況；一天內(nèi)太陽(yáng)光照變化也會(huì)給圖像匹配精度和指向偏差計(jì)算的效果帶來(lái)影響，這兩部分都將是后序工作需要繼續(xù)改進(jìn)的內(nèi)容。

圖9 一天內(nèi)利用地標(biāo)進(jìn)行圖像糾正得到的指向偏差結(jié)果對(duì)比Fig.9 Comparison of pointing error results after image rectification through landmark matching within one day

圖10 一天內(nèi)利用地標(biāo)進(jìn)行圖像糾正得到的指向偏差的統(tǒng)計(jì)分析結(jié)果Fig.10 Pointing error statistical analysis results after image rectification through landmark matching within one day

3 結(jié) 論

本文圍繞三軸穩(wěn)定型靜止軌道遙感衛(wèi)星的指向偏差計(jì)算問(wèn)題開(kāi)展了研究，利用分布在地球表面的海陸邊界數(shù)據(jù)作為匹配模板與地標(biāo)，研究了一套精度較高且能夠穩(wěn)健運(yùn)行的地標(biāo)匹配方法，重點(diǎn)突破了靜止地標(biāo)數(shù)據(jù)的制作、圖像與地標(biāo)的高精度匹配及基于粗差剔除的指向計(jì)算等方法。試驗(yàn)表明，該方法在解決靜止軌道遙感圖像指向偏差方面具有較好效果，并已經(jīng)在我國(guó)相關(guān)衛(wèi)星數(shù)據(jù)的在線處理中得到應(yīng)用。該方法目前還存在不少問(wèn)題，主要集中在圖像與地標(biāo)匹配的效率上，目前的全圓盤(pán)地標(biāo)匹配算法的效率在幾分鐘到十幾分鐘之間(業(yè)務(wù)系統(tǒng)中)，匹配精度的穩(wěn)定性還不能完全適應(yīng)全天24 h的各時(shí)段、各類(lèi)型數(shù)據(jù)，仍然有待提升。下一步工作中，由于靜止軌道遙感衛(wèi)星的時(shí)間分辨率(從十幾分鐘變成幾分鐘)、空間分辨率(數(shù)百米變成數(shù)十米)越來(lái)越高，對(duì)地標(biāo)匹配的效率也要求越來(lái)越高，需要進(jìn)一步突破精度更高、更加穩(wěn)健、效率更快(準(zhǔn)實(shí)時(shí))的匹配與指向計(jì)算方法，且能夠適應(yīng)更加復(fù)雜的圖像場(chǎng)景。