HY-1C/D 衛星中國海洋水色水溫掃描儀幾何定位方法

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

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

1 引言

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

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

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

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

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

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

2 COCTS 觀測特性

2.1 COCTS 介紹

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

表1 COCTS 地球區圖像數據傳輸格式Table 1 COCTS earth area image data transmission format

2.2 地面景物與COCTS 成像幾何關系

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

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

3 幾何定位方法

本文提出的COCTS 幾何定位方法算法流程如圖2所示。首先，構建ORB 到ECR 的轉換矩陣TECR/ORB。從0 級數據中提取星下點計數值，計算出星下點時間，根據相鄰采樣時間關系計算出各點采樣時間。在提取的星歷數據中采用插值法獲取各采樣時間對應的ECR 衛星位置和速度，通過在ECR 坐標系中建立ORB 坐標系來構造ORB 到ECR 的轉換矩陣TECR/ORB。

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

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

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

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

3.1 坐標系及轉換矩陣 TECR/ORB的 構建

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

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

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

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

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

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

則ORB 到ECR 的轉換矩陣為TECR/ORB=[XORBYORBZORB]。

3.2 插值算法

星下點位置是地球采樣區的中間位置，從0 級輔助數據中提取星下點時間計數值，計算出星下點時間，相鄰采樣點間隔124 μs，所以根據相鄰采樣時間關系可以推導出各點采樣時間，本節插值算法即是針對采樣時間進行的。從0 級GPS 定位廣播數據中提取并計算出GPS 絕對定位時間以及對應的衛星ECR位置和速度，但是由于重復的數據較多，無法直接用于后續的插值算法，所以需要先進行重復數據的剔除處理。

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

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

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

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

最后采樣時間對應的衛星ECR 位置、速度計算公式如下：

3.3 各采樣點ECR 視矢量的構建

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

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

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

利用3.1 節構建的ORB 到ECR 轉換矩陣TECR/ORB，則ECR 坐標系下各采樣點視矢量計算如下：

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

利用地球橢球體長半軸a、短半軸b重新調節ECR 視矢量uECR和 衛星位置矢量pECR，下列等式未考慮光傳輸時間和因衛星移動或相對效應造成的異常，這將導致一點系統性的偏差。

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

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

3.5 ECR 到GEO 坐標系轉換

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

則大地經度（L on ）、緯度（Lat）可表示為

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

4 結果驗證

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

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

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

4.1 HY-1C 衛星COCTS 與MODIS 海岸線定性比較

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

1）NDVI 處理

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

2）Canny 邊緣檢測

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

3）人工優化

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

4.1.1 區域1—阿拉伯半島

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

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

4.1.2 區域2—渤海灣

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

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

4.2 HY-1C 衛星COCTS 海岸線經緯度與GSHHG 數據庫定量比較

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

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

圖9 阿拉伯半島遙感圖參考海岸線疊加（a），星下點附近中間放大區域（b），星下點兩側邊緣放大區域（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 衛星COCTS 星下點經緯度與衛星工具軟件結果定量比較

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

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

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

表2 計算值與參考值誤差統計Table 2 Error statistics between calculated and reference values

4.4 HY-1D 衛星與HY-1C 衛星COCTS 經緯度定量比較

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

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

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

表3 HY-1C/D 衛星特征區域采樣點比較Table 3 Sampling points comparison of HY-1C/D satellite feature areas

5 像元尺度效應對幾何定位誤差影響分析

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

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

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

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

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

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

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

續表 6

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

圖13 赤道及兩側特征點采樣（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 結論

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

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

（2）針對海岸線定量驗證，COCTS 幾何定位方法計算的海岸線經緯度與GSHHG 參考海岸線坐標相比，誤差在2 個像元內，即2 km 內。

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

（4）COCTS 幾何定位方法對HY-1D 衛星與HY-1C 衛星COCTS 的定位結果一致。

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

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