段吉蔚，俞杭華，劉會杰

（1.中國科學院 微小衛(wèi)星創(chuàng)新研究院，上海 201203；2.中國科學院大學，北京 100049；3.上?？萍即髮W 信息科學與技術學院，上海 201210）

0 引言

隨著北斗衛(wèi)星導航系統(tǒng)（Beidou Navigation Satellite System，BDS）的全面運行，北斗、全球定位系統(tǒng)（Global Positioning System，GPS）、Galileo和Glonass 系統(tǒng)等4 大全球衛(wèi)星導航系統(tǒng)（Global Navigation Satellite System，GNSS）在軌衛(wèi)星數(shù)量共計已超過百顆［1］。與利用單系統(tǒng)導航定位的方式相比，多系統(tǒng)融合定位可以充分利用多系統(tǒng)的可見衛(wèi)星，實時可見衛(wèi)星數(shù)量更多、分布更好，有利于提高導航定位服務的精度和可用性。然而，如果將多系統(tǒng)所有可見衛(wèi)星均用于導航定位計算，雖然理論上定位精度最優(yōu)，但計算效率降低，對應用終端的資源需求大，且嚴重影響導航定位的實時性。合適的選星策略可以在保證定位精度的同時，降低計算復雜度，提高導航定位的實時性，增強大多數(shù)GNSS 應用。因此，對GNSS 多系統(tǒng)融合導航定位選星策略的研究具有重要的理論和實際意義。

目前，選星算法的研究方向主要有以下3 個方面。1）利用智能優(yōu)化算法進行選星。文獻［2］提出了一種基于灰狼優(yōu)化（Grey Wolf Optimizer，GWO）算法的快速選星方法，該算法利用自適應收斂因子和信息反饋機制，在局部尋優(yōu)與全局搜索之間實現(xiàn)平衡。文獻［3］提出了一種基于NSGA-Ⅱ算法的多目標快速選星方法，該算法通過不斷的選擇不同的遺傳算子和效用函數(shù)做實驗，最終選取一組合適的遺傳算子和效用函數(shù)做出最優(yōu)的選星決策，并且不依賴于衛(wèi)星的幾何位置分布，可適用于有障礙或者遮擋的復雜情況。文獻［4］提出了一種基于K-means++聚類算法的快速選星算法，該方法優(yōu)先選擇仰角最大的一顆衛(wèi)星后，對剩余衛(wèi)星進行K-means++聚類算法，然后在每個簇中選擇一顆衛(wèi)星，得到最終的選星組合。以上算法會導致較高的計算復雜度，增加選星所需要的時間，對實時性要求較高的應用場景可能不太適用。2）根據(jù)可見衛(wèi)星的幾何空間分布進行選星。文獻［5］通過實驗發(fā)現(xiàn)了可見星的空間分布規(guī)律，提出了一種基于衛(wèi)星高度角和方位角進行分區(qū)篩選的快速選星算法。文獻［6］根據(jù)最佳選星方案的衛(wèi)星分布特性，通過利用衛(wèi)星高度角和方位角信息進行區(qū)域劃分，并引入代價函數(shù)篩選中仰角區(qū)域的衛(wèi)星，從而得到最終的選星方案。文獻［7］分析了幾何精度因子（Geometric Dilution Precision，GDOP），與可見衛(wèi)星布局的關系，確定各個衛(wèi)星對GDOP 的貢獻值，選擇貢獻值大的衛(wèi)星作為最終用于定位解算的衛(wèi)星組合。3）利用計算幾何學中的凸包算法進行選星。文獻［8］提出用三維凸包算法進行選星定位，雖然該算法考慮了衛(wèi)星在三維空間中的分布，能夠在選星過程中找到更有利于定位的衛(wèi)星子集，但是三維凸包算法在計算時的復雜度較高，無法滿足實時性的要求。因此為了提高實時性，文獻［9］提出將三維衛(wèi)星數(shù)據(jù)降為二維，然后利用二維凸包Graham 算法進行快速選星，但是該選星算法將三維衛(wèi)星數(shù)據(jù)直接投影到站心坐標系下的水平平面，但可能會忽略高度方向上的影響，導致選擇的衛(wèi)星幾何分布不夠理想，從而影響定位精度。

綜上所述，基于二維凸包算法的選星策略在處理三維衛(wèi)星數(shù)據(jù)降維時，常忽略垂直方向的高度位置信息，為解決該問題，本文提出了一種融合主成分分析（Principal Component Analysis，PCA）與二維凸包Melkman 算法的選星策略。該策略利用PCA 技術在降維過程中最大限度地減小信息損失，確保生成的二維數(shù)據(jù)不僅保留了原始三維數(shù)據(jù)的水平信息，還包含了高度方向的信息。在后續(xù)選星計算中，能更準確地描述衛(wèi)星的位置信息，其對于選星策略至關重要。此外，本文采用時間復雜度較低的Melkman 算法進行選星計算，以提高選星速度。

1 幾何精度因子

GDOP 是衛(wèi)星導航系統(tǒng)中一個重要指標，其代表衛(wèi)星導航系統(tǒng)測距誤差造成的接收機與空間衛(wèi)星間的距離矢量放大因子［10］。導航衛(wèi)星的定位誤差的標準偏差δ定義如下：

式中：σUERE為偽距測量誤差。

從式（1）可知，在偽距測量誤差相同的情況下，導航定位精度完全由GDOP 決定，衛(wèi)星幾何分布越優(yōu)，GDOP 越小，導航精度越高。GDOP 的推導過程如下：

式中：H為BDS/GPS/GLONASS 中某單一衛(wèi)星系統(tǒng)觀測矩 陣的前3 列；α為導航衛(wèi)星的 仰角；β為導航衛(wèi)星的方位角；N為單一系統(tǒng)參與定位解算的衛(wèi)星的數(shù)量。

由于本文選取導航星座中的BDS、GPS 和GLONASS 這3 個系統(tǒng)衛(wèi)星進行選星，因此這里只考慮BDS/GPS/GLONASS 組合導航系統(tǒng)的觀測矩陣：

式中：1x、0x（x為BDS、GPS、GLONASS）分別為k×1 維的全1、全0 向量；k為當前系統(tǒng)x用于定位解算的衛(wèi)星數(shù)量。

則GDOP 可以表示為［11］

GDOP 數(shù)值與定位精度的關系見表1［12］。

表1 GDOP 與定位精度的評價關系Tab.1 Evaluation relationship between the GDOP and positioning accuracy

2 基本原理

2.1 主成分分析的原理

PCA 方法的基本思想是通過將數(shù)據(jù)投影到最大方差的方向上，找到數(shù)據(jù)中最主要的模式和結構，并且保留盡可能多的信息［13］。利用PCA 算法將三維衛(wèi)星數(shù)據(jù)降為二維時，通常需要經過以下步驟［14-16］：

1）對原始數(shù)據(jù)進行標準化。假設有n個衛(wèi)星數(shù)據(jù)樣本，每個樣本有3 個特征(x，y，z)，對于每個特征，計算其均值μ和標準差σ，然后將每個特征的值減去其均值，再除以其標準差，以使每個特征具有相同的權重，將標準化后的數(shù)據(jù)整理為一個n×3的矩陣Xstd。以特征x為例，標準化的過程可以用以下公式描述：

式中：μ為特征x的均值；σ為特征x的標準差。

2）計算標準化后的數(shù)據(jù)矩陣Xstd的協(xié)方差矩陣S。協(xié)方差矩陣可以用以下公式計算：

3）對協(xié)方差矩陣進行特征值分解，得到特征值和特征向量。特征值和特征向量可以用以下公式計算：

式中：vi為第i個特征向量；λi為第i個特征向量對應的特征值。

4）選擇主成分。將特征值按降序排列，選取前k個最大特征值對應的特征向量作為數(shù)據(jù)降維后的新坐標系的基向量，組成新的特征向量矩陣Vk。在本文中，由于要將三維數(shù)據(jù)降為二維，因此k=2。

5）將原始數(shù)據(jù)投影到新的特征向量上，得到降維后的數(shù)據(jù)。新矩陣Y的每一行表示原始數(shù)據(jù)在新特征空間的坐標。

以上5 個步驟的流程如圖1 所示。

圖1 PCA 算法流程Fig.1 Flowchart of the PCA algorithm

通過遵循上述5 個步驟，可以使用PCA 算法將原始三維衛(wèi)星數(shù)據(jù)有效地降至二維空間。由于該新的二維空間中的基向量是原始數(shù)據(jù)空間的線性組合，且相互正交，因此降維后的二維數(shù)據(jù)不僅保留了原始三維衛(wèi)星數(shù)據(jù)的水平位置信息，還保留了高度位置信息，為后續(xù)選星計算提供了一個更精確地描述衛(wèi)星位置信息的有效方式。

2.2 Melkman 算法原理

凸包的公式表達如下：在一個實數(shù)向量空間V中，對于給定集合X，所有包含X的凸集的交集CH(X)被稱為X的凸包［17］，即：

X的凸包CH(X)可以用X內所有點的線性組合來構造：

二維凸包是將平面上給定的點集用一個最小的凸多邊形包圍的形狀。對于二維凸包的計算，可以使用不同的算法來實現(xiàn)，包括暴力、Graham 掃描、Jarvis 步 進、QuickHull 算 法、Melkman 算法等［18-21］。其中，Melkman 算法的時間復雜度為O（n），其比大多數(shù)其他算法（如Graham 掃描法、分治法和Andrew′s Monotone Chain Algorithm）的O（nlogn）時間復雜度更低［22］，這意味著在處理二維凸包問題時，Melkman 算法可以更快地找到解決方案。因此本文采用Melkman 算法進行凸包計算。Melkman算法的核心思想是利用Deque 來維護凸包邊界上的點。算法的具體步驟如下［23］：

1）對輸入數(shù)據(jù)進行預判斷，確保至少有3 個點。如果輸入數(shù)據(jù)少于3 個點，則無法構成多邊形，無需計算凸包。

2）初始化雙端隊列（Double-ended Queue，Deque）。在所有點中找出具有最小x值的點A，然后確定一個處理順序，例如逆時針方向，計算剩余點與點A形成的向量與y軸負方向的夾角，再根據(jù)這些夾角將點集按從小到大的順序進行排序，最后，將排序后的前兩個點（點A和點B）加入Deque。

3）按順序處理剩余的點，記當前處理的點為C。進行以下操作：①檢查隊列的頭部。如果C與隊列頭部的前2 個點構成一個右拐（順時針方向），則隊列頭部的第一個點出隊。重復此操作，直到隊列頭部的2 個點與C構成一個左拐（逆時針方向）為止。② 檢查隊列的尾部。如果C與隊列尾部的后2個點構成一個左拐（逆時針方向），則隊列尾部的第1 個點出隊。重復此操作，直到隊列尾部的2 個點與C構成一個右拐（順時針方向）為止。③將點C添加到隊列的尾部。

值得注意的是：判斷3 點之間的相對方向（左拐或右拐）可以通過計算向量積（Cross Product）來實現(xiàn)。例如，有3 個點A、B和C，那么如果向量AB與向量BC的向量積大于0，則表示點C在AB的逆時針方向（左拐）；反之，如果向量積小于0，則表示點C在AB的順時針方向（右拐）。

4）重復步驟3，直到所有點都被處理。

5）當所有點處理完畢后，Deque 中的點即為凸包上的頂點。從隊列頭部到尾部，按順序連接這些頂點，就可以得到凸包的邊界。Melkman 算法的流程如圖2 所示。

圖2 Melkman 算法流程Fig.2 Flowchart of the Melkman algorithm

3 基于PCA 和Melkman 的選星算法

3.1 選星算法執(zhí)行步驟

本文針對BDS/GPS/GLONASS 三模導航定位系統(tǒng)提出了一種基于PCA 和Melkman 的選星算法，具體實驗步驟如下所述：

1）根據(jù)衛(wèi)星廣播星歷信息，分別計算衛(wèi)星在BDS、GPS、GLONASS 坐標系統(tǒng)下的坐標、仰角和方位角。雖然3 個導航系統(tǒng)的坐標系統(tǒng)不相同，但是其均為笛卡爾坐標系，可以根據(jù)文獻［24］和文獻［25］將3 個坐標系統(tǒng)一轉換到GPS 系統(tǒng)的WGS-84坐標系下進行定位解算。

2）將衛(wèi)星原始三維數(shù)據(jù)通過2.1 章分析的PCA 技術投影到新的二維坐標系下。

3）對二維平面的數(shù)據(jù)進行預處理，篩選掉不在凸包邊界上的衛(wèi)星數(shù)據(jù)。具體步驟如下：①找到x坐標、y坐標最小、最大的4 個點，分別記為Px-min（xmin，y1），Px-max（xmax，y2），Py-min（x1，ymin），Py-max（x2，ymax），顯然，這4 個點都在凸包的邊界上。② 將橫坐標介于x1和x2之間，縱坐標介于y1和y2之間的點（圖3 中的D 區(qū)域）從點集中剔除掉，因為D 區(qū)域的點不可能在凸包的邊界上。

圖3 數(shù)據(jù)預處理的衛(wèi)星分類Fig.3 Satellite classification for data preprocessing

為了提高Melkman 算法在二維凸包掃描中的優(yōu)化效果，首先分析了已有凸包掃描結果，發(fā)現(xiàn)當衛(wèi)星位于圖3 中的區(qū)域D 時，假設該點在凸包邊界上，則該點與圖3 中已知的4 個凸包邊界上的點圍成的“凸包”違反了凸包的性質：凸包上的任意2 點之間的線段都完全位于凸包內部或邊界上，因此，該假設不成立，可以確定區(qū)域D 內的衛(wèi)星一定不在二維凸包邊界上。基于該發(fā)現(xiàn)，可以針對衛(wèi)星數(shù)據(jù)進行預處理，減少后續(xù)Melkman 算法需要掃描的衛(wèi)星數(shù)量，從而提高整體選星效率。

4）對新的二維數(shù)據(jù)預處理之后，基于2.2 章分析的Melkman 算法找到篩選后的點集的凸包。

5）凸包點即為最終選擇用于用戶定位解算的衛(wèi)星，找到凸包點對應的衛(wèi)星三維數(shù)據(jù)，計算選星后的GDOP 值。

3.2 選星算法仿真分析

為了驗證本文提出選星算法的性能，本文用投影到站心坐標系XOY平面的二維凸包Melkman 選星算法與基于PCA 和Jarvis 步進法的選星算法作為對比實驗，對上海市（121.47°E，31.23°N）位置進行選星實驗，選取導航星座中的BDS、GPS 和GLONASS 系統(tǒng)的衛(wèi)星，為保證實驗準確性，衛(wèi)星位置由星歷計算得出，同時將衛(wèi)星的高度截止角設定為5°，仿真實驗的起始時間為UTC 時間2023 年2 月26 日00：00：00，持續(xù)時間為24 h，每1 h 進行一次采樣，共計24 個采樣點。選星算法的流程如圖4所示。

圖4 選星算法的流程Fig.4 Flowchart of the satellite selection algorithm

UTC 時間14：00：00 使用PCA 算法生成的新的二維平面如圖5 所示。UTC 時間14：00：00 使 用PCA 算法將原始三維衛(wèi)星數(shù)據(jù)投影至新的二維平面的分布如圖6 所示。

圖5 PCA 算法生成的新的二維平面Fig.5 New two-dimensional plane generated by the PCA algorithm

圖6 二維PCA 平面散點Fig.6 Scatter plots in the two-dimensional PCA plane

預處理剔除的無關衛(wèi)星數(shù)量變化如圖7 所示?？梢园l(fā)現(xiàn)，無論是投影到XOY坐標系，還是PCA 投影生成的新坐標系下，在大多數(shù)時刻都會剔除較多的無關衛(wèi)星數(shù)，在19：00 時刻甚至可以減少18 顆無關衛(wèi)星數(shù)據(jù)，可為后續(xù)Melkman算法減少了大量的掃描衛(wèi)星數(shù)，提高了整體算法掃描二維凸包的選星速度。

投影到站心坐標系XOY平面的二維凸包Melkman 選星算法和本文提出選星算法這2 種算法的選星數(shù)量如圖8 所示?？梢园l(fā)現(xiàn)，在大部分時間里，本文所提算法的選星數(shù)量要多于投影到站心坐標系XOY平面的二維凸包Melkman 選星算法的選星數(shù)量，為后續(xù)定位解算提供了更多的觀測數(shù)據(jù)，從而進一步提高定位精度。

圖8 XOY-Melkman 與PCA-Melkman 算法的選星數(shù)量Fig.8 Numbers of selected satellites obtained by the XOYMelkman and PCA-Melkman algorithms

XOY-Melkman 與PCA-Melkman 這2 種選星算法GDOP 曲線和運行時間曲線，如圖9 和圖10 所示。PCA-Jarvis 與PCA-Melkman 選星算法運行時間曲線，如圖11 所示。3 種選星算法的選星數(shù)、GDOP 和運行時間進行列表分析，見表2。如圖9 所示，本文所提選星算法在60%的時間里實現(xiàn)了GDOP 小于2，而在其余40%的時間里GDOP 值介于2～3.5 之間，根據(jù)表1 中GDOP 與定位精度的關系，說明本文所提選星算法在一天中的60%時間能夠達到優(yōu)秀的定位精度，而在剩下的40%時間能夠獲得良好的定位精度。相對而言，投影到站心坐標系XOY平面的二維凸包Melkman選星算法僅有25%的時間達到優(yōu)秀的定位精度。結合圖10 和表2，可以發(fā)現(xiàn)本文所提選星算法的平均運行時間為0.002 0 s，投影到站心坐標系XOY平面的二維凸包Melkman 選星算法的平均運行時間為0.001 1 s，雖然運行效率降低，但是本文所提選星算法能夠實現(xiàn)更優(yōu)的GDOP 性能，即平均GDOP 性能可以提升0.384 0。結合圖11 和表2，本文所提選星算法與PCA-Jarvis 選星算法獲取的GDOP 性能一樣，但是本文所提算法的平均運行時間卻只有PCA-Jarvis 選星算法的三分之一。因此，從選星結果的角度來看，本文提出的選星算法是一種可行的選星策略，可以為實時定位、精確測量等領域提供高效可靠的解決方案。

圖9 XOY-Melkman 與PCA-Melkman 選 星算 法GDOP 性能對比Fig.9 Comparison of the GDOP obtained by the XOYMelkman and PCA-Melkman satellite selection algorithms

圖10 XOY-Melkman 與PCA-Melkman 選星算法運行時間對比Fig.10 Comparison of the running time obtained by the XOY-Melkman and PCA-Melkman satellite selection algorithms

圖11 PCA-Jarvis與PCA-Melkman選星算法運行時間對比Fig.11 Comparison of the running time obtained by the PCA-Jarvis and PCA-Melkman satellite selection algorithms

表2 3 種選星算法衛(wèi)星數(shù)、GDOP 與運行時間對比Tab.2 Comparison of the results of satellite number，GDOP，and running time obtained by three satellite selection algorithms

4 結束語

本文針對二維凸包選星算法在處理三維衛(wèi)星數(shù)據(jù)降維過程中對高度位置信息的忽略，提出了一種基于PCA 和Melkman 算法的選星策略。該策略運用PCA 技術在降維時最大程度地減小信息損失，生成的二維數(shù)據(jù)不僅保留了原始三維數(shù)據(jù)的水平信息，還包含了高度方向的信息，從而在后續(xù)選星計算中能更準確地反映衛(wèi)星位置。同時，該策略采用時間復雜度較低的Melkman 算法計算凸包，以提高選星速度。實驗結果顯示，相較于將數(shù)據(jù)投影到站心坐標系XOY平面的二維凸包 Melkman 選星算法，本文提出的選星算法實現(xiàn)了更優(yōu)的GDOP 性能，平均GDOP性能提升了0.384 0；相較于PCA-Jarvis 選星算法，本文所提算法的平均運行時間僅有PCA-Jarvis 選星算法的三分之一。因此，本文所提出的選星策略在運行速度和定位精度方面具有明顯優(yōu)勢。