基于主成分分析和Melkman 算法的選星策略

段吉蔚，俞杭華，劉會杰

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

0 引言

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

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

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

1 幾何精度因子

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

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

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

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

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

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

則GDOP 可以表示為［11］

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

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

2 基本原理

2.1 主成分分析的原理

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

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

式中：μ為特征x的均值；σ為特征x的標(biāo)準(zhǔn)差。

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

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

式中：vi為第i個特征向量；λi為第i個特征向量對應(yīng)的特征值。

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

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

以上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 算法原理

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

X的凸包CH(X)可以用X內(nèi)所有點(diǎn)的線性組合來構(gòu)造：

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

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

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

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

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

4）重復(fù)步驟3，直到所有點(diǎn)都被處理。

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

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

3 基于PCA 和Melkman 的選星算法

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

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

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

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

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

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

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

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

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

3.2 選星算法仿真分析

為了驗(yàn)證本文提出選星算法的性能，本文用投影到站心坐標(biāo)系XOY平面的二維凸包Melkman 選星算法與基于PCA 和Jarvis 步進(jìn)法的選星算法作為對比實(shí)驗(yàn)，對上海市（121.47°E，31.23°N）位置進(jìn)行選星實(shí)驗(yàn)，選取導(dǎo)航星座中的BDS、GPS 和GLONASS 系統(tǒng)的衛(wèi)星，為保證實(shí)驗(yàn)準(zhǔn)確性，衛(wèi)星位置由星歷計算得出，同時將衛(wèi)星的高度截止角設(shè)定為5°，仿真實(shí)驗(yàn)的起始時間為UTC 時間2023 年2 月26 日00：00：00，持續(xù)時間為24 h，每1 h 進(jìn)行一次采樣，共計24 個采樣點(diǎn)。選星算法的流程如圖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 平面散點(diǎn)Fig.6 Scatter plots in the two-dimensional PCA plane

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

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

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

圖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 選星算法運(yùn)行時間對比Fig.10 Comparison of the running time obtained by the XOY-Melkman and PCA-Melkman satellite selection algorithms

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

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

4 結(jié)束語

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