劉建陽，毛志華, ，陶邦一，馬力，朱乾坤，黃海清，劉建強，丁靜

(1.上海交通大學(xué) 海洋學(xué)院，上海 200240；2.自然資源部第二海洋研究所 衛(wèi)星海洋環(huán)境動力學(xué)國家重點實驗室，浙江杭州 310012；3.國家衛(wèi)星海洋應(yīng)用中心，北京 100081)

1 引言

20 世紀(jì)70 年代美國發(fā)射了第一代海色探測器海岸帶水色掃描儀（CZCS），用于監(jiān)測海洋水色和浮游生物情況。搭載中分辨率成像光譜儀（MODIS）的兩顆衛(wèi)星Terra 和Aqua 分別于1999 年12 月、2002 年5 月發(fā)射成功，至今仍然在軌運行。MODIS 共有36 個光譜波段，其中分辨率為250 m 的波段有2 個，分辨率為500 m 的波段有5 個，分辨率為1 km 的波段有29 個，其在110°的視場角內(nèi)實現(xiàn)2 330 km 的掃描寬度[1]。Nishihama 等[2]在1997 年提出MODIS 幾何定位算法，從焦平面位置開始，經(jīng)過轉(zhuǎn)換與計算得到探測器坐標(biāo)系下視矢量，再經(jīng)過衛(wèi)星（SC）坐標(biāo)系、軌道（ORB）坐標(biāo)系、地心慣性（ECI）坐標(biāo)系、地心旋轉(zhuǎn)（ECR）坐標(biāo)系轉(zhuǎn)換，得到ECR 視矢量，結(jié)合衛(wèi)星ECR 位置，計算ECR 視矢量與地球的交點，最后轉(zhuǎn)換交點為地理經(jīng)緯度。

HY-1A/B 衛(wèi)星水色儀的幾何定位方法利用下行數(shù)據(jù)中的GPS 數(shù)據(jù)和軌道擬合技術(shù)得到衛(wèi)星軌道參數(shù)，計算衛(wèi)星在掃描時刻的ECI 位置和速度，結(jié)合水色掃描儀工作原理和球面幾何關(guān)系，求解掃描點在ECI 下的位置，過程中涉及SC 坐標(biāo)系、ORB 坐標(biāo)系、ECI 坐標(biāo)系的轉(zhuǎn)換，然后再轉(zhuǎn)到ECR 坐標(biāo)系，最后計算出掃描點的地理經(jīng)緯度[3]。HY-1A/B 衛(wèi)星水色儀與MODIS 幾何定位的區(qū)別在于，一是利用GPS 數(shù)據(jù)和軌道擬合技術(shù)來計算衛(wèi)星的位置和速度，二是未采用視矢量。

風(fēng)云三號（FY-3）衛(wèi)星微波成像儀采用天線繞軸旋轉(zhuǎn)形成圓錐形跨軌的掃描方式，GPS 接收機提供衛(wèi)星三維位置，利用兩組以上衛(wèi)星位置計算出衛(wèi)星的ECI 實時速度。根據(jù)探測器觀測幾何、探測器空間位置和指向，建立觀測像元與地面位置的關(guān)系模型。先將計算的天線視矢量經(jīng)過坐標(biāo)系轉(zhuǎn)換，生成ECR 視矢量，然后結(jié)合衛(wèi)星ECI 位置到ECR 位置的轉(zhuǎn)換，即可計算出ECR 視矢量與地球的交點，最后再把交點轉(zhuǎn)換成地理經(jīng)緯度[4]。不同于MODIS，F(xiàn)Y-3 微波成像儀根據(jù)GPS 數(shù)據(jù)中衛(wèi)星位置計算速度。

以上幾種探測器的幾何定位方法一致，關(guān)鍵過程是計算出探測器視線與地球橢球體的交叉點，過程涉及坐標(biāo)系轉(zhuǎn)換與計算，最后轉(zhuǎn)化為地理經(jīng)緯度。該幾何定位方法同樣應(yīng)用于陸地遙感衛(wèi)星多光譜掃描儀（MSS）、可見光紅外成像輻射儀（VIIRS）、甚高分辨率掃描輻射計（AVHRR）等探測器中[5]，應(yīng)用過程中存在各種變形和改進，如傳感器視線呈現(xiàn)形式不同，獲得衛(wèi)星位置速度的方式不同，初始視矢量所在坐標(biāo)系不同，坐標(biāo)系變換步驟不同，經(jīng)緯度計算公式不同等。

潘德爐等[6]在1997 年提出了基于衛(wèi)星軌道參數(shù)和探測器觀測角的幾何定位方法，用于衛(wèi)星仿真，既不采用地面控制點，也不采用插值法近似逼近。根據(jù)先驗知識獲得的已知點作為地理參考點，利用衛(wèi)星軌道參數(shù)和探測器觀測角等參數(shù)把圖像平面坐標(biāo)依次轉(zhuǎn)化為以星下點為球心的球面坐標(biāo)、以星下點為極點的極坐標(biāo)、以北極為極點的極坐標(biāo)，最后轉(zhuǎn)化為像元地理坐標(biāo)，未采用常規(guī)坐標(biāo)系。研究過程中通過逐步逼近法，推算出所有行的中心像元點地理坐標(biāo)和衛(wèi)星飛行方向，以此減少誤差。

2018 年9 月7 日，我國在太原衛(wèi)星發(fā)射中心成功發(fā)射HY-1C 衛(wèi)星，2020 年6 月11 日又成功發(fā)射HY-1D衛(wèi)星，兩星組網(wǎng)，進一步提高了全球海洋水色、海洋帶資源與生態(tài)環(huán)境的有效觀測能力。作為HY-1A/HY-1B 衛(wèi)星的后繼者，HY-1C/D 衛(wèi)星的正常運行標(biāo)志著海洋一號系列衛(wèi)星業(yè)務(wù)化運行的成功。HY-1C/D 衛(wèi)星搭載中國海洋水色水溫掃描儀（COCTS），COCTS幾何定位是遙感資料預(yù)處理的核心要求，需要通過方法研究來提高定位精度，精確的地理定位數(shù)據(jù)可以提高COCTS 后續(xù)產(chǎn)品的質(zhì)量。

2 COCTS 觀測特性

2.1 COCTS 介紹

HY-1C/D 衛(wèi)星搭載5 個有效載荷，本文重點研究其中的水色儀COCTS。COCTS 重達50 kg，功率為70 W，數(shù)據(jù)傳輸速率為670 kb/s。COCTS 以1 km 的空間分辨率、114°的全視域角、2 900 km 的地面幅寬掃描監(jiān)測中國鄰海以及全球海洋，提供了海洋水色、植被產(chǎn)品以及海表溫度等數(shù)據(jù)。沿著HY-1C/D 衛(wèi)星飛行方向，COCTS 從右向左掃描，四元同步逐點擺掃，每元有1 664 個采樣點，每個采樣點有10 個波段（B1-B10），包括8 個可見波段和近紅外波段、2 個熱紅外波段。水色儀地球區(qū)圖像數(shù)據(jù)傳輸格式如表1 所示。

表1 COCTS 地球區(qū)圖像數(shù)據(jù)傳輸格式Table 1 COCTS earth area image data transmission format

2.2 地面景物與COCTS 成像幾何關(guān)系

COCTS 采用45°掃描鏡和四元探測器掃描成像形式，其掃描方向表示為+Y，指向地心方向表示為+Z，飛行方向表示為+X，垂直于Y軸、Z軸所在平面指向外側(cè)，建立的坐標(biāo)系如圖1 所示。COCTS 每360°掃描周期內(nèi)都有一段對地球區(qū)進行掃描成像，圖1 中地球區(qū)以星下點為中心，左右各轉(zhuǎn)57°，即從成像開始逆時針旋轉(zhuǎn)114°的區(qū)域。COCTS 安裝在衛(wèi)星載荷艙對地面，依靠45°掃描鏡穿軌方向的掃描和衛(wèi)星的飛行，獲取地球區(qū)圖像數(shù)據(jù)，10 波段焦平面與地面景物成像幾何關(guān)系如圖1，獲取的波段數(shù)據(jù)將用于后續(xù)遙感圖的繪制。

圖1 水色儀COCTS 周期掃描成像幾何關(guān)系示意Fig.1 Diagram of the geometric relation of COCTS periodic scanning imaging

3 幾何定位方法

本文提出的COCTS 幾何定位方法算法流程如圖2所示。首先，構(gòu)建ORB 到ECR 的轉(zhuǎn)換矩陣TECR/ORB。從0 級數(shù)據(jù)中提取星下點計數(shù)值，計算出星下點時間，根據(jù)相鄰采樣時間關(guān)系計算出各點采樣時間。在提取的星歷數(shù)據(jù)中采用插值法獲取各采樣時間對應(yīng)的ECR 衛(wèi)星位置和速度，通過在ECR 坐標(biāo)系中建立ORB 坐標(biāo)系來構(gòu)造ORB 到ECR 的轉(zhuǎn)換矩陣TECR/ORB。

圖2 COCTS 幾何定位方法流程Fig.2 Flow chart of COCTS geometric positioning method

其次，計算掃描行各采樣點的ORB 視矢量。探測器（INST）坐標(biāo)系下原點視矢量為OZ=[0,0,1]，忽略儀器安裝誤差，理論上認為安裝矩陣是單位矩陣。INST 到SC 坐標(biāo)系轉(zhuǎn)換矩陣TSC/INST為單位矩陣，熱脹冷縮是影響該矩陣動態(tài)誤差的因素，誤差相對較小且難被分析，可忽略。SC 到ORB 坐標(biāo)系轉(zhuǎn)換矩陣TORB/SC由衛(wèi)星姿態(tài)決定，使用HY-1C/D 衛(wèi)星雙敏感器與陀螺的搭配來確定衛(wèi)星姿態(tài)，精度更高，穩(wěn)定性更好，故本文采用基于擴展卡爾曼濾波（EKF）姿態(tài)確定算法，分成兩個子系統(tǒng)進行研究[7]。由于此方法難度較大，本文暫時將TORB/SC簡單視為單位矩陣，所以O(shè)RB 坐標(biāo)系原點視矢量為OZ=[0,0,1]，掃描行的各采樣點ORB 視矢量可以通過將原點視矢量繞X、Y軸旋轉(zhuǎn)相應(yīng)角度的方式來獲得。

然后，計算采樣點ECR 視矢量與地球橢球體的交叉點，并轉(zhuǎn)換為地理經(jīng)緯度。結(jié)合轉(zhuǎn)換矩陣TECR/ORB，進一步計算出各采樣點的ECR 視矢量，建立ECR 視矢量與地球橢球體交叉模型，從而求解出交叉點ECR 坐標(biāo)，利用英國學(xué)者鮑林的研究思路導(dǎo)出的公式即可轉(zhuǎn)換交叉點ECR 坐標(biāo)為地理經(jīng)緯度[8]。

最后，結(jié)合波段數(shù)據(jù)繪制遙感圖。根據(jù)上述過程計算的各采樣點經(jīng)緯度以及提取的對應(yīng)波段數(shù)據(jù)，設(shè)置像元尺度，投影各采樣點，各投影像元匹配對應(yīng)采樣點的平均波段數(shù)據(jù)，據(jù)此可以繪制遙感圖。

3.1 坐標(biāo)系及轉(zhuǎn)換矩陣 TECR/ORB的 構(gòu)建

在水色儀COCTS 幾何定位算法中，涉及到3 個主要坐標(biāo)系，分別是ORB 坐標(biāo)系、ECR 坐標(biāo)系、GEO坐標(biāo)系，坐標(biāo)系及幾何定位示意圖如圖3 所示。ORB坐標(biāo)系基于慣性空間的衛(wèi)星位置，以衛(wèi)星質(zhì)心為中心，Z軸從衛(wèi)星質(zhì)心指向地球質(zhì)心，Y軸是Z軸與瞬時速度矢量的標(biāo)準(zhǔn)叉乘，X軸為Y軸與Z軸的標(biāo)準(zhǔn)叉乘。

圖3 坐標(biāo)系及幾何定位示意Fig.3 Coordinate system and geometric positioning diagram

ECR 坐標(biāo)系以地球質(zhì)心為起點，X軸從地球質(zhì)心指向格林威治子午線與赤道的交點，Z軸從地球質(zhì)心指向北半球極點，Y軸為Z軸與X軸的叉乘，ECR 坐標(biāo)系是一個空間直角坐標(biāo)系，所以某點空間位置可以用該點在坐標(biāo)軸上的投影值來表示。

GEO 坐標(biāo)系是基于WGS-84，用經(jīng)度、緯度表示坐標(biāo)的坐標(biāo)系。經(jīng)度定義為本地子午線與格林威治子午線的夾角，格林威治子午線向東為正，稱為東經(jīng)，向西為負，稱為西經(jīng)，變化范圍是環(huán)全球經(jīng)度；緯度定義為橢球體法線與赤道的夾角，向北為正，稱為北緯，向南為負，稱為南緯，變化范圍是90°S～90°N[8]。

傳統(tǒng)的ORB 到ECR 坐標(biāo)系的轉(zhuǎn)換需要兩個步驟，一是從ORB 轉(zhuǎn)換到ECI，該過程需要用到衛(wèi)星的ECI位置和速度；二是從ECI 轉(zhuǎn)換到ECR，該過程需要考慮地球的自轉(zhuǎn)、極移、章動、進動等因素的影響[9]。由于地球方位信息動態(tài)變化，所以需要實時從國際地球自轉(zhuǎn)和參考系服務(wù)網(wǎng)站獲取更新的信息。

傳統(tǒng)方法較為復(fù)雜，考慮因素較多，并且需要實時更新信息，因此采用一種簡捷、不需實時獲取更新信息的直接轉(zhuǎn)換方法顯得尤為必要。本文根據(jù)衛(wèi)星的ECR 位置和速度，建立了ORB 與ECR 的聯(lián)系，本質(zhì)上是ECR 坐標(biāo)系分別繞著3 個坐標(biāo)軸旋轉(zhuǎn)一定角度，最終轉(zhuǎn)換到ORB 坐標(biāo)系。假設(shè)某一時刻衛(wèi)星的ECR 位置和速度分別為PECR、VECR，則ORB 坐標(biāo)系的單位矢量計算如下[10]：

則ORB 到ECR 的轉(zhuǎn)換矩陣為TECR/ORB=[XORBYORBZORB]。

3.2 插值算法

星下點位置是地球采樣區(qū)的中間位置，從0 級輔助數(shù)據(jù)中提取星下點時間計數(shù)值，計算出星下點時間，相鄰采樣點間隔124 μs，所以根據(jù)相鄰采樣時間關(guān)系可以推導(dǎo)出各點采樣時間，本節(jié)插值算法即是針對采樣時間進行的。從0 級GPS 定位廣播數(shù)據(jù)中提取并計算出GPS 絕對定位時間以及對應(yīng)的衛(wèi)星ECR位置和速度，但是由于重復(fù)的數(shù)據(jù)較多，無法直接用于后續(xù)的插值算法，所以需要先進行重復(fù)數(shù)據(jù)的剔除處理。

衛(wèi)星位置的每一個元素被當(dāng)作獨立的一維函數(shù)，速度作為其函數(shù)的一階導(dǎo)數(shù)；每組位置和速度組合起來進行插值，從而保證插值位置和速度的一致性。從0 級輔助數(shù)據(jù)中提取的GPS 時間都是整數(shù)時間，相鄰間隔為1 s，每個時間對應(yīng)一組三維ECR 位置和速度，使用插值法計算各采樣時間對應(yīng)的衛(wèi)星ECR 位置和速度，在GPS 時間范圍外的采樣時間無法進行插值計算。插值法示意如圖4 所示。

圖4 插值法示意Fig.4 Diagram of insertion method

為計算采樣時間Ts對應(yīng)的衛(wèi)星ECR 位置和速度，在參考時間中找到采樣時間Ts的兩個最近相鄰整數(shù)時間T1、T2，dT是相鄰時間T1、T2之差，根據(jù)時間T1、T2及其對應(yīng)的位置P1、P2和速度V1、V2，計算多項式系數(shù)a0、a1、a2、a3，公式如下[11]：

兩個位置P1、P2是 連續(xù)相鄰的，兩個速度V1、V2也是連續(xù)相鄰的。為了插入到合適的點，對采樣時間Ts進行歸一化處理，轉(zhuǎn)化為0 與1 之間的數(shù)值：

最后采樣時間對應(yīng)的衛(wèi)星ECR 位置、速度計算公式如下：

3.3 各采樣點ECR 視矢量的構(gòu)建

針對ORB 坐標(biāo)系下視矢量的求解，以下降過程為例加以說明，如圖5 所示，交點表示水色儀掃描地面所采樣的4×1 664 個點。沿著衛(wèi)星運行方向，COCTS探頭從右向左掃描，圖5 中最左邊是第1 列，最右邊是最后一列；每掃描行有四元即4 個探頭，從下向上分別是第一元到第四元。坐標(biāo)原點設(shè)置在正中間位置，衛(wèi)星運行方向表示X軸正方向，掃描方向表示Y軸正方向，構(gòu)建圖5 所示的右正左負、下正上負坐標(biāo)系。假設(shè)Z軸垂直X、Y軸所在平面指向地面，則水色儀在原點處ORB 視矢量可視為OZ=[0,0,1]。

圖5 水色儀降COCTS 軌掃描示意Fig.5 Diagram of COCTS downward scanning

掃描行每元都有1 664 個采樣點，各列采樣通過COCTS 探頭在太空中旋轉(zhuǎn)一定角度來實現(xiàn)。以星下點為中心，COCTS 探頭向左、向右最大旋轉(zhuǎn)角度設(shè)計值都是57°，在空中的整個視域角為114°，實際在軌運行時視域角約為116°。COCTS 旋轉(zhuǎn)掃描一圈的周期為640 ms，每個周期對應(yīng)360°，相鄰采樣點掃描間隔為124 μs，那么相鄰列之間旋轉(zhuǎn)角roll=124/(640 000)×2π，原點ORB 視矢量需要繞X軸旋轉(zhuǎn)角度XR=roll×[831.5,-1,-831.5]來獲得各列視矢量；X軸方向上的最小旋轉(zhuǎn)角度pitch=0.001 38，原點ORB 視矢量需要繞Y軸旋轉(zhuǎn)角度YR=pitch×[1.5,0.5,-0.5,-1.5]來獲得各元視矢量。所以各采樣點ORB 視矢量根據(jù)原點視矢量OZ 分別繞X軸、Y軸旋轉(zhuǎn)一定角度來獲得[12]，每掃描行共計4×1 664 個視矢量，即為ORB 坐標(biāo)系下的視矢量uORB。

利用3.1 節(jié)構(gòu)建的ORB 到ECR 轉(zhuǎn)換矩陣TECR/ORB，則ECR 坐標(biāo)系下各采樣點視矢量計算如下：

3.4 ECR 視矢量與地球交叉算法

利用地球橢球體長半軸a、短半軸b重新調(diào)節(jié)ECR 視矢量uECR和 衛(wèi)星位置矢量pECR，下列等式未考慮光傳輸時間和因衛(wèi)星移動或相對效應(yīng)造成的異常，這將導(dǎo)致一點系統(tǒng)性的偏差。

使用二次方程式求解與單位球面交叉的視矢量u′的 縮放系數(shù)d，下式d是兩個交點中靠近衛(wèi)星的一個解，即較小的那個解，再使用d計算橢球體交叉矢量x。

所以，根據(jù)ECR 視矢量uECR、衛(wèi)星位置矢量pECR以及推算出的縮放系數(shù)d，可以計算出ECR 視矢量與地球橢球體交叉點坐標(biāo)pECR+duECR，即在ECR 空間直角坐標(biāo)系下的坐標(biāo)。

3.5 ECR 到GEO 坐標(biāo)系轉(zhuǎn)換

ECR 到GEO 坐標(biāo)系轉(zhuǎn)換是根據(jù)英國科學(xué)家鮑林研究思路推演而成的計算公式，把ECR 坐標(biāo)系下的坐標(biāo)（X，Y，Z）轉(zhuǎn)化為大地坐標(biāo)，即地理經(jīng)緯度[8]。a、b分別為地球橢球體的長半軸、短半軸；f為橢球扁率；據(jù)此計算橢球第一偏心率平方e2、第二偏心率平方e′2：

則大地經(jīng)度（L on ）、緯度（Lat）可表示為

本文未考慮地形對于幾何定位精度的影響，根據(jù)后文對幾何定位結(jié)果的驗證，可以得知地形校正前該幾何定位方法對經(jīng)緯度誤差的影響程度，并且在后續(xù)的研究中使用數(shù)字高程模型DEM 對地形加以校正，從而減少誤差。

4 結(jié)果驗證

在環(huán)全球經(jīng)度范圍，緯度范圍90°S～90°N 的地圖上，進行網(wǎng)格化處理，橫向經(jīng)度分成36 000 個格子，縱向緯度分成18 000 個格子，最終形成方格邊長為0.01°的36 000×18 000 網(wǎng)格，表征經(jīng)度范圍為360°、緯度范圍為180°的投影地圖。上述幾何定位方法計算出的經(jīng)緯度，可以轉(zhuǎn)化為在網(wǎng)格中的行列號，再轉(zhuǎn)化行列號為在網(wǎng)格中的序號，網(wǎng)格按從上向下，從左向右順序標(biāo)序號。把投影在同一網(wǎng)格序號中的采樣點的R、G、B 波段數(shù)據(jù)分別取平均，表示該網(wǎng)格序號對應(yīng)的三波段數(shù)據(jù)，最后顯示得到赤道區(qū)域分辨率約為1 km 的彩色遙感圖。

選 擇COCTS 從2020 年10 月27 日 02:33:00 UTC開始掃描到12:35:00 UTC 結(jié)束的0 級數(shù)據(jù)，衛(wèi)星在該段時間內(nèi)共繞地運行6 軌，掃描范圍可以覆蓋歐亞非區(qū)域。按上述方法計算該區(qū)域的經(jīng)緯度以及提取波段數(shù)據(jù)，作圖6a 所示歐亞非大陸遙感圖，圖6a 為1 km分辨率遙感圖，1 km 高分辨率體現(xiàn)地圖的精細化信息，根據(jù)遙感圖估算的經(jīng)緯度數(shù)據(jù)也更加準(zhǔn)確。MODIS 1 級產(chǎn)品光譜分辨率有250 m、500 m、1 000 m 3 種模式，圖6b 為相同區(qū)域的MODIS 遙感圖。從圖6 可以看出，COCTS 與MODIS 遙感圖輪廓大體一致。

圖6 歐亞非大陸遙感圖Fig.6 Remote sensing image of Europe,Asia and Africa

4.1 HY-1C 衛(wèi)星COCTS 與MODIS 海岸線定性比較

選擇阿拉伯半島和渤海灣兩區(qū)域進行海岸線定性驗證。阿拉伯半島海陸分界線輪廓清晰，經(jīng)過Canny 邊緣檢測、人工優(yōu)化處理即可提取海岸線，然后進行后續(xù)的輪廓匹配；而渤海灣區(qū)域具有光譜反射率低、上空多云多霧等特征，造成陸海分界線在視覺上不明顯，需要先對波段數(shù)據(jù)進行歸一化植被指數(shù)（NDVI）處理，以此增強海陸色度差，提高區(qū)分度，再實施Canny 邊緣檢測、人工優(yōu)化來提取海岸線[13]，最后對提取的海岸線進行輪廓匹配。

1）NDVI 處理

NDVI 是衡量植被覆蓋率的一個參數(shù)，根據(jù)NIR和R 兩個波段反射值之差與之和的比來計算，范圍處于-1 與1 之間。根據(jù)R、NIR 波段范圍，從COCTS的10 個波段中分別選擇第6、8 波段來計算NDVI，正值表征植被，植被覆蓋率越高，NDVI 越大。通過對NDVI 進行適合閾值設(shè)定，從而剔除絕大部分白云、水、巖石、裸土以及綠色植被覆蓋較少的地方，剩下的深綠區(qū)域植被覆蓋率較高[14-15]。

2）Canny 邊緣檢測

首先，高斯濾波確定高斯核的尺寸為5×5，標(biāo)準(zhǔn)差為0.5，與圖像進行離散卷積；其次，考慮到Sobel 算子邊緣強弱分明，抗噪聲好，選定該算子計算X方向和Y方向像素梯度；然后，將當(dāng)前像素梯度強度與沿正負梯度方向上的相鄰像素梯度強度進行比較，抑制非極大值像素梯度，消除邊緣檢測帶來的雜散響應(yīng)；最后，設(shè)置高低閾值，定義真實邊緣點[16-18]。

3）人工優(yōu)化

經(jīng)過前面兩步處理，海陸交界處像素梯度顯著，所以海岸線被提取得很好。然而海岸圖仍有尚未剔除的干擾線段，所以最后需要對海岸圖進行人工優(yōu)化，用干擾線段附近的圖像數(shù)據(jù)取代干擾區(qū)域，從而屏蔽干擾線段，得到最終清晰的海岸線。

4.1.1 區(qū)域1—阿拉伯半島

根據(jù)MODIS 在2020 年10 月27 日的1 級產(chǎn)品畫出1 km 空間分辨率的阿拉伯半島地區(qū)遙感圖，如圖7a所示。對其進行Canny 邊緣檢測，調(diào)整閾值參數(shù)，得到圖7b，該海岸線還存在干擾線段，因此還需進行人工優(yōu)化，剔除噪聲干擾，優(yōu)化如圖7c 或圖7d所示，該圖即為根據(jù)MODIS 數(shù)據(jù)最終提取的阿拉伯半島地區(qū)海岸線。從2020 年10 月27 日 05:54:23 UTC 至09:15:17 UTC，COCTS 繞地運行兩軌，覆蓋阿拉伯半島，該時間段的0 級數(shù)據(jù)按上述方法計算出1 級數(shù)據(jù)，據(jù)此作遙感圖7e，與圖7a 同樣經(jīng)緯度范圍。利用CorelDRAW 軟件把圖7d 海岸線疊加到圖7e 遙感圖中，通過調(diào)整透明度進而驗證海岸線吻合情況，因此從圖7f 中可以定性判斷MODIS 白色海岸線與HY-1C 衛(wèi)星COCTS 遙感圖海岸線匹配一致，二者變化趨勢吻合。

圖7 阿拉伯半島海岸線提取與匹配Fig.7 Extraction and match of coastline in the Arabian Peninsula

4.1.2 區(qū)域2—渤海灣

COCTS 從2020 年10 月19 日 02:12:25 UTC 至02:22:37 UTC 掃描探測渤海灣區(qū)域，該時間段對應(yīng)的0 級數(shù)據(jù)按上述幾何定位方法計算出經(jīng)緯度，提取波段信息，據(jù)此作1 km 分辨率遙感圖，依次進行圖8 所示處理。圖8e 為根據(jù)同一天的MODIS 1 級產(chǎn)品作同一經(jīng)緯區(qū)域的遙感圖，圖8f 為COCTS 提取的海岸線與MODIS 遙感圖海岸線輪廓匹配結(jié)果，因此可以定性判斷COCTS 白色海岸線與MODIS 遙感圖海岸線匹配一致，二者變化趨勢吻合。

圖8 渤海灣海岸線提取與匹配Fig.8 Extraction and match of coastline in the Bohai Gulf

4.2 HY-1C 衛(wèi)星COCTS 海岸線經(jīng)緯度與GSHHG 數(shù)據(jù)庫定量比較

GSHHG 是全球一致、分層、高分辨率地理數(shù)據(jù)庫，由3 個數(shù)據(jù)庫構(gòu)成：世界矢量海岸線、世界數(shù)據(jù)庫II 和冰凍圈地圖集，分別是海岸線、湖泊和南極海岸線的基礎(chǔ)。GSHHG 地理數(shù)據(jù)分辨率有5 個級別，選取最高分辨率“full”模式。因為“high”模式空間分辨率為200 m，“full”模式所占磁盤空間是“high”模式的4 倍多，所以“full”分辨率是優(yōu)于50 m 的[19]。為了檢驗幾何定位的精度，根據(jù)幾何定位方法計算的結(jié)果作1 km 分辨率遙感圖，從GSHHG 數(shù)據(jù)庫中提取相同區(qū)域的海岸線坐標(biāo)，在該遙感圖上疊加相應(yīng)的海陸分界線，該海路分界線即為參考，通過其與實際遙感圖海岸線的吻合程度來判斷幾何定位的精度。

選擇上述COCTS 在2020 年10 月27 日掃描覆蓋阿拉伯半島的數(shù)據(jù)，作1 km 分辨率遙感圖（圖9a）。選定阿拉伯半島區(qū)域重點研究，用海岸線提取工具，選定始末坐標(biāo)，開始于35.91°N，30.78°E，結(jié)束于11.34°N，60.84°E。該始末坐標(biāo)范圍內(nèi)的海岸線坐標(biāo)都可以提取出來，疊加到圖9a 中，紅色即為參考海岸線。選取中間區(qū)域的波斯灣北部，進一步放大如圖9b 所示；選取兩側(cè)區(qū)域的波斯灣東南部，進一步放大如圖9c 所示，由此二圖可以清晰看到像元級別。由于COCTS的掃描特征，越靠近掃描區(qū)域的兩側(cè)，采樣越稀疏，采樣不足導(dǎo)致存在無效數(shù)據(jù)點，在遙感圖上體現(xiàn)為黑色像元點。紅色參考海岸線與實際作的遙感圖海岸線吻合程度反映了幾何定位精度，如圖9 所示，COCTS幾何定位方法計算的經(jīng)緯度誤差在2 個像元內(nèi)，即2 km 內(nèi)。

圖9 阿拉伯半島遙感圖參考海岸線疊加（a），星下點附近中間放大區(qū)域（b），星下點兩側(cè)邊緣放大區(qū)域（c）Fig.9 Reference coastline overlaying Arabian remote sensing image (a),middle enlargement area near Nadir (b),marginal enlargement area away from Nadir (c)

4.3 HY-1C 衛(wèi)星COCTS 星下點經(jīng)緯度與衛(wèi)星工具軟件結(jié)果定量比較

成熟的衛(wèi)星工具軟件（Satellite Tool Kit，STK）是航天領(lǐng)域先進的系統(tǒng)分析軟件，用于分析復(fù)雜的陸地、海洋等任務(wù)以及提供精確的分析報告。高精度軌道生成模塊考慮點重力模型、日月重力影響、大氣阻力、章動、自旋、質(zhì)心變化等很多因素，以此確保生成高精度軌道。

選擇某一區(qū)域數(shù)據(jù)，從0 級文件中提取GPS 時間及對應(yīng)的位置和速度，利用STK 軟件的高精度軌道生成模塊模擬衛(wèi)星軌道，然后根據(jù)輸入的星下點時間，輸出對應(yīng)的衛(wèi)星位置和速度。利用獲得的星下點時間、位置、速度在STK 軟件上再次模擬出衛(wèi)星軌道，星下點作為特征點，一個星下點時間對應(yīng)一組星下點經(jīng)緯度。然后導(dǎo)出LLA 文件，文件中包括星下點經(jīng)緯度[20-21]，該數(shù)據(jù)作為參考標(biāo)準(zhǔn)，與本文幾何定位方法計算的星下點經(jīng)緯度進行比較和分析。

首先選擇COCTS 從2020 年10 月19 日 02:12:25 UTC 至02:22:37 UTC 掃描覆蓋渤海灣區(qū)域的0 級數(shù)據(jù)，利用上述定位方法計算957 組星下點經(jīng)緯度，STK輸出對應(yīng)的參考值。計算值與參考值對比結(jié)果為，誤差在0.001°～0.01°之間的星下點占22.46%，誤差在0.01°～0.02°之間的星下點占77.54%，誤差在2 個像元內(nèi)。再選取長時間大范圍的數(shù)據(jù)，選擇COCTS 從2020 年10 月27 日 07:34:50 UTC 至09:15:17 UTC 掃描經(jīng)過阿拉伯半島的一整軌0 級數(shù)據(jù)，同樣方法計算出9 418 組星下點經(jīng)緯度以及輸出STK 對應(yīng)參考值。為了進一步對比分析，把數(shù)據(jù)分成南北緯0°～30°、南北緯30°～60°、南北緯60°～90°共低、中、高緯3 部分，計算值與參考值對比，誤差占比統(tǒng)計結(jié)果如表2 所示。分析星下點誤差可知，本文幾何定位方法計算的星下點誤差都在0.02°內(nèi)，即2 個像元內(nèi)；從低緯到高緯誤差變化趨勢可以判斷，赤道誤差最小，越往兩極，誤差越大。

表2 計算值與參考值誤差統(tǒng)計Table 2 Error statistics between calculated and reference values

4.4 HY-1D 衛(wèi)星與HY-1C 衛(wèi)星COCTS 經(jīng)緯度定量比較

HY-1D 衛(wèi)星與HY-1C 衛(wèi)星功能一致，共同監(jiān)測全球水色和海表溫度，兩顆衛(wèi)星組網(wǎng)，提高全球掃描覆蓋能力，共同構(gòu)成中國首個海洋民用衛(wèi)星星座。兩星搭載同樣的水色儀COCTS，COCTS 的硬件結(jié)構(gòu)、運行模式和工作原理等完全一致。HY-1D 衛(wèi)星與HY-1C 衛(wèi)星的區(qū)別是，在遙感圖上HY-1D 衛(wèi)星從右下往左上運行成像，而HY-1C 是從右上往左下運行成像。

將2020 年9 月23 日 HY-1D COCTS 的遙感圖與上述2020 年10 月27 日 HY-1C COCTS 的遙感圖進行比較。選擇波斯灣西部一個特征小區(qū)域，遙感圖放大6 倍，如圖10 所示，圖10a 為HY-1C 特征區(qū)域，圖10b 為HY-1D 衛(wèi)星特征區(qū)域，兩圖中紅色像元點為用于定量比較的選取點，把選取點在全球范圍遙感圖中的行列號記錄在表3。經(jīng)統(tǒng)計，在1 km 分辨率前提下，HY-1C 衛(wèi)星與HY-1D 衛(wèi)星的緯度、經(jīng)度之差都為0 像元，所以COCTS 幾何定位方法對HY-1D 衛(wèi)星與HY-1C 衛(wèi)星都適用，二者定位結(jié)果一致。

圖10 HY-1C 衛(wèi)星（a）與HY-1D 衛(wèi)星（b）對應(yīng)特征點采樣Fig.10 Corresponding feature points sampled from HY-1C satellite (a) and HY-1D satellite (b)

表3 HY-1C/D 衛(wèi)星特征區(qū)域采樣點比較Table 3 Sampling points comparison of HY-1C/D satellite feature areas

5 像元尺度效應(yīng)對幾何定位誤差影響分析

COCTS 幾何定位方法的地理定位結(jié)果存在誤差，本節(jié)分析像元尺度效應(yīng)對幾何定位誤差的影響。地球是個不規(guī)則的三維橢球體，兩極略尖，赤道略鼓，從而半軌成像具有兩極寬、赤道窄的特點。選取COCTS 在2020 年10 月27 日掃描覆蓋阿拉伯半島的半軌經(jīng)緯度數(shù)據(jù)，對于1 664 列數(shù)據(jù)，從中選擇7 對相鄰樣本列，分別是1 和2 列、208 和209 列、416 和417 列、832 和833 列、1 247 和1 248 列、1 454 和1 455 列、1 663 和1 664 列。為了詳細分析半軌經(jīng)緯度數(shù)據(jù)，把半軌掃描區(qū)域分成3 部分：靠近北極的前1/3 區(qū)域、中間1/3 區(qū)域、靠近南極的后1/3 區(qū)域。分別計算各區(qū)域每對樣本列的平均距離，記錄如表4。根據(jù)表4中數(shù)據(jù)作7 對樣本列距離變化趨勢圖，便于直觀看出經(jīng)緯度數(shù)據(jù)的特征趨勢，如圖11 所示，黑色、綠色、藍色虛線分別表示前1/3 區(qū)域、中間1/3 區(qū)域、后1/3 區(qū)域，紅色實線表示半軌區(qū)域，同時用柱形圖表示該半軌區(qū)域的7 對樣本列距離。結(jié)合水色儀COCTS的掃描原理，分析表4 數(shù)據(jù)和圖11 的變化趨勢，可以得出，從星下點到兩側(cè)邊緣，從赤道到兩極，采樣像元尺度呈增大趨勢。

圖11 7 對樣本列距離變化趨勢Fig.11 Distance change trends of seven pairs of adjacent sample arrays

表4 相鄰樣本列距離（單位：（°））Table 4 Distance of adjacent sample arrays (unit:(°))

半軌掃描區(qū)域如圖12a 中黑色線框范圍，為了具體分析兩側(cè)和星下點附近中間區(qū)域的誤差特征，從兩側(cè)和中間區(qū)域各選取40 個海岸特征點，分別用紅色、黃色標(biāo)志，進一步放大如圖12b 所示。兩側(cè)和中間區(qū)域的特征點在1 km 分辨率遙感圖中的行列號以及與GSHHG 數(shù)據(jù)庫海岸線坐標(biāo)的像元誤差分別記錄在表5 和表6 中，經(jīng)誤差統(tǒng)計，兩側(cè)區(qū)域的緯度方向、經(jīng)度方向平均誤差分別是 0.95、1.25 個像元，即誤差為1.57 個像元；中間區(qū)域的緯度方向、經(jīng)度方向平均誤差分別是 0.30、0.65 個像元，即誤差為0.72 個像元，所以掃描區(qū)域兩側(cè)的幾何定位誤差比中間大。理論分析結(jié)合實際結(jié)果驗證可知，星下點誤差最小，越接近兩側(cè)邊緣，采樣像元畸變越大，定位誤差隨之增大。因此，從星下點到兩側(cè)邊緣，定位誤差隨采樣像元尺度增大而呈增大趨勢，二者呈正相關(guān)。

圖12 中間及兩側(cè)邊緣特征點采樣（a）及放大（b）Fig.12 Feature points sampled in the middle and marginal sides (a) and enlargement (b)

表5 靠近邊緣兩側(cè)區(qū)域特征點誤差Table 5 Errors of feature points on both edge-nearing sides of the track

表6 中間區(qū)域特征點誤差Table 6 Errors of feature points in the middle area of the track

續(xù)表 6

為了具體分析從赤道到兩極的誤差特征，在掃描范圍內(nèi)的赤道附近區(qū)域和南北緯35°附近區(qū)域分別選取20 個海岸特征點，并且特征點都分布在星下點附近區(qū)域，如圖13a 所示，進一步放大如圖13b 所示。赤道附近、35°N 附近和35°S 附近的特征點在遙感圖中的行列號，以及特征點與GSHHG 數(shù)據(jù)庫海岸線坐標(biāo)的像元誤差分別記錄在表7、表8、表9 中，經(jīng)誤差統(tǒng)計，赤道附近特征點的緯度方向、經(jīng)度方向平均誤差分別是0.20、0.40 個像元，即誤差為0.45 個像元；35°N 附近特征點的緯度方向、經(jīng)度方向平均誤差分別是0.55、0.80 個像元，即誤差為0.97 個像元；35°S附近特征點的緯度方向、經(jīng)度方向平均誤差分別是0.50、0.70 個像元，即誤差為0.86 個像元，因此掃描區(qū)域赤道附近定位誤差比兩側(cè)小。理論分析結(jié)合實際結(jié)果驗證可知，赤道定位誤差最小，越靠近兩極，采樣像元畸變越大，定位誤差隨之增大。因此，從赤道到兩極，定位誤差隨采樣像元尺度增大而逐漸增大，二者呈正相關(guān)。

圖13 赤道及兩側(cè)特征點采樣（a）及放大（b）Fig.13 Feature points sampled near the equator and its both sides (a) and enlargement (b)

表7 赤道附近特征點誤差Table 7 Errors of feature points near the equator

表8 35°N 附近特征點誤差Table 8 Errors of feature points near 35°N

表9 35°S 附近特征點誤差Table 9 Errors of feature points near 35°S

6 結(jié)論

HY-1C/D 衛(wèi)星水色儀COCTS 幾何定位方法是遙感資料預(yù)處理系統(tǒng)的重要組成部分，精確的地理定位數(shù)據(jù)即1 級產(chǎn)品可以提高HY-1C/D 衛(wèi)星后續(xù)產(chǎn)品的質(zhì)量，因此建立一套將水色儀0 級產(chǎn)品處理成1 級產(chǎn)品的系統(tǒng)，包括數(shù)據(jù)解包、幾何定位，顯得尤為重要。本文幾何定位方法取代了傳統(tǒng)需要根據(jù)6 個軌道參數(shù)來計算衛(wèi)星位置的方法，也簡化了傳統(tǒng)ORB到ECR 的復(fù)雜轉(zhuǎn)換過程，在衛(wèi)星星歷中采用插值法獲得采樣時間對應(yīng)的衛(wèi)星位置和速度，進而直接計算ORB 到ECR 的轉(zhuǎn)換矩陣?；贑OCTS 逐點擺掃、四元并掃的掃描方式以及相關(guān)參數(shù)，計算采樣點ORB 視矢量，建立采樣點視矢量與地球交叉點關(guān)系模型，從而對COCTS 采樣像元進行地理定位。通過對HY-1C/D 衛(wèi)星水色儀COCTS 幾何定位結(jié)果的驗證與分析，得出以下結(jié)論：

（1）針對海岸線定性驗證，COCTS 與MODIS 遙感圖海岸線匹配一致，二者變化趨勢吻合。

（2）針對海岸線定量驗證，COCTS 幾何定位方法計算的海岸線經(jīng)緯度與GSHHG 參考海岸線坐標(biāo)相比，誤差在2 個像元內(nèi)，即2 km 內(nèi)。

（3）針對星下點定量驗證，星下點坐標(biāo)計算值與STK 生成的參考值相比，誤差在2 個像元內(nèi)且從赤道至兩極，星下點誤差呈遞增趨勢。

（4）COCTS 幾何定位方法對HY-1D 衛(wèi)星與HY-1C 衛(wèi)星COCTS 的定位結(jié)果一致。

（5）像元尺度效應(yīng)與幾何定位誤差呈正相關(guān)。從星下點到兩側(cè)邊緣，定位誤差隨采樣像元尺度增大而呈增大趨勢，星下點誤差最小；從赤道到兩極，定位誤差也隨采樣像元尺度增大而逐漸增大，赤道定位誤差最小。

綜上所述，經(jīng)過定性、定量驗證與分析，以衛(wèi)星軌道參數(shù)和水色儀COCTS 參數(shù)為基礎(chǔ)的幾何定位方法是可行的，滿足一定的定位精度要求，可以用于COCTS 遙感圖像預(yù)處理的幾何定位。