大型高精度GNSS基線(xiàn)向量網(wǎng)并行抗差估計(jì)

崔 陽(yáng)，呂志平，張友陽(yáng)，李林陽(yáng).信息工程大學(xué)地理空間信息學(xué)院，河南鄭州45000；.鄭州輕工業(yè)學(xué)院機(jī)電工程學(xué)院，河南鄭州45000

大型高精度GNSS基線(xiàn)向量網(wǎng)并行抗差估計(jì)

崔 陽(yáng)1，呂志平1，張友陽(yáng)2，李林陽(yáng)1
1.信息工程大學(xué)地理空間信息學(xué)院，河南鄭州450001；2.鄭州輕工業(yè)學(xué)院機(jī)電工程學(xué)院，河南鄭州450002

針對(duì)大型GNSS基線(xiàn)向量網(wǎng)的特點(diǎn)，在改進(jìn)的相關(guān)觀(guān)測(cè)抗差估計(jì)RECO方案的基礎(chǔ)上采用并行計(jì)算技術(shù)進(jìn)行相關(guān)抗差估計(jì)的并行計(jì)算（簡(jiǎn)稱(chēng)“并行抗差估計(jì)”），并給出了可行的解算步驟。通過(guò)算例分析，驗(yàn)證了大型GNSS基線(xiàn)向量網(wǎng)的并行相關(guān)抗差估計(jì)，不僅有效抑制了觀(guān)測(cè)異常對(duì)參數(shù)估值的影響，而且顯著提高了計(jì)算效率。

相關(guān)抗差估計(jì)；GNSS；基線(xiàn)向量網(wǎng)；等價(jià)權(quán)矩陣；并行估計(jì)

1　引 言

隨著越來(lái)越多的國(guó)家級(jí)或區(qū)域級(jí)CORS網(wǎng)絡(luò)被部署，高精度GNSS基線(xiàn)向量網(wǎng)規(guī)模也越來(lái)越大［1］。受觀(guān)測(cè)條件的影響，GNSS觀(guān)測(cè)值中不可避免地帶有粗差，而高精度GPS網(wǎng)粗差比例比常規(guī)大地測(cè)量網(wǎng)相對(duì)較大［2-3］，使得大型高精度GNSS基線(xiàn)向量網(wǎng)粗差的處理更加繁瑣。目前對(duì)于觀(guān)測(cè)值帶有相關(guān)性和異常影響的高精度基線(xiàn)向量網(wǎng)平差處理，通常采用相關(guān)觀(guān)測(cè)抗差估計(jì)方法，如相關(guān)分析的方法［2］、基于標(biāo)準(zhǔn)化殘差法［4］、IGGⅢ方案［3］及其改進(jìn)［5］和RECO方案［6-7］及其改進(jìn)［8］等，而這類(lèi)相關(guān)抗差估計(jì)方法的實(shí)現(xiàn)均要進(jìn)行選權(quán)迭代計(jì)算，使得涉及大型矩陣運(yùn)算的相關(guān)抗差估計(jì)計(jì)算開(kāi)銷(xiāo)也大大增加。對(duì)于大型GNSS基線(xiàn)向量網(wǎng)的處理，由于粗差的復(fù)雜性和抗差能力的有限性，需要頻繁試算以選擇合適的抗差估計(jì)模型與參數(shù)；而且由于科學(xué)研究的需要和應(yīng)用周期的縮短，經(jīng)常面臨采用不同的方案進(jìn)行多次的平差計(jì)算，長(zhǎng)時(shí)間的平差試算與分析降低了科學(xué)研究的效率。隨著用戶(hù)對(duì)大規(guī)模GNSS數(shù)據(jù)處理實(shí)時(shí)性的要求越來(lái)越高，研究提高抗差計(jì)算效率，將“事后計(jì)算”提高到“準(zhǔn)實(shí)時(shí)計(jì)算”在定軌、機(jī)動(dòng)測(cè)量保障、網(wǎng)絡(luò)定位服務(wù)等應(yīng)用領(lǐng)域具有重要意義。

對(duì)于運(yùn)算復(fù)雜的大規(guī)模網(wǎng)平差問(wèn)題，過(guò)去由于計(jì)算機(jī)內(nèi)存和計(jì)算能力的限制，文獻(xiàn)［2］推出了一種多維二次型矩陣壓縮存儲(chǔ)的快速直觀(guān)算法以提高運(yùn)算速度；文獻(xiàn)［9］利用解算約化代替求逆約化，提出了一種分區(qū)平差快速解算新算法。隨著個(gè)人計(jì)算機(jī)硬件平臺(tái)的高速發(fā)展，提高硬件平臺(tái)的性能利用效率和大規(guī)模數(shù)據(jù)處理的效率，對(duì)現(xiàn)有測(cè)量數(shù)據(jù)處理軟件代碼并行化或獨(dú)立開(kāi)發(fā)并行分布式軟件系統(tǒng)成為解決測(cè)量數(shù)據(jù)處理的新方法［1015］。針對(duì)大型高精度GNSS向量網(wǎng)的抗差估計(jì)問(wèn)題，本文在改進(jìn)RECO方案基礎(chǔ)上利用并行計(jì)算技術(shù)進(jìn)行并行相關(guān)抗差估計(jì)的研究。針對(duì)基于消息傳遞接口（message passing interface，MPI）的并行算法設(shè)計(jì)復(fù)雜、消息傳遞和數(shù)據(jù)移動(dòng)開(kāi)銷(xiāo)大問(wèn)題，著重考慮當(dāng)前主流的共享內(nèi)存的多核系統(tǒng)優(yōu)勢(shì)，結(jié)合使用Microsoft在2010年發(fā)布的.NET4并行擴(kuò)展開(kāi)發(fā)工具并行任務(wù)庫(kù)（Task Parallel Library，TPL），提出了一種基于任務(wù)劃分的相關(guān)抗差估計(jì)多核并行算法，通過(guò)算例驗(yàn)證了大型GNSS基線(xiàn)向量網(wǎng)相關(guān)抗差估計(jì)并行算法的有效性，提高了硬件資源利用效率和抗差計(jì)算效率。

2　相關(guān)抗差估計(jì)模型的并行化

2.1 相關(guān)抗差估計(jì)模型

大型高精度GNSS網(wǎng)通常采用劃分子網(wǎng)模型進(jìn)行處理［2，16-17］。假設(shè)一個(gè)GNSS網(wǎng)劃分為p個(gè)子網(wǎng)，得到m組基線(xiàn)向量組，各向量組的向量數(shù)分別為n1、n2、…、nm（n1+n2+…+nm=n），第i（i=1，2，…，m）個(gè)基線(xiàn)向量組的相關(guān)基線(xiàn)向量觀(guān)測(cè)值Li建立的誤差方程為

式中，常數(shù)項(xiàng)li=Li-AiX0，相關(guān)觀(guān)測(cè)值Li的先驗(yàn)權(quán)矩陣為Pi。對(duì)m個(gè)基線(xiàn)向量組進(jìn)行整網(wǎng)平差，則整網(wǎng)平差的誤差方程為

式中由于向量組之間的觀(guān)測(cè)值是獨(dú)立的，所以全網(wǎng)的權(quán)矩陣P是由各個(gè)組向量權(quán)矩陣Pi組成的一準(zhǔn)對(duì)角陣，則基于抗差M估計(jì)的選權(quán)迭代解式為

式中，ˉP稱(chēng)為相關(guān)等價(jià)權(quán)陣。在雙因子等價(jià)權(quán)模型中

式中γii和γjj稱(chēng)為自適應(yīng)降權(quán)因子或收縮因子。在第k次迭代時(shí)，γii可取

式中

k0和k1稱(chēng)為閾值，根據(jù)實(shí)際情況選用不同的數(shù)值，一般取k0=1.0～1.5，k1=2.5～5.5。

式（7）單位權(quán)中誤差計(jì)算公式中，n′等于觀(guān)測(cè)量個(gè)數(shù)n減去淘汰觀(guān)測(cè)量（γii為0）個(gè)數(shù)；median是中位數(shù)計(jì)算函數(shù)，通常是選擇抗差性更高的中位數(shù)法計(jì)算單位權(quán)中誤差［18-19］。為了在精度穩(wěn)定的情況下獲得最快求解效率，不同于采用IGGⅢ權(quán)函數(shù)，在文獻(xiàn)［8］RECO改進(jìn)方法的計(jì)算步驟基礎(chǔ)上，將降權(quán)因子γ（k）ij求解公式中的平方改為可變參數(shù)t（t=1，2，…），并進(jìn)一步改正了協(xié)方差陣Qv的求解公式。由于Qv的求解過(guò)程涉及等價(jià)權(quán)的求逆，而迭代過(guò)程中求逆非常耗時(shí)，通常采用不嚴(yán)密的方法進(jìn)行簡(jiǎn)單處理，而在并行計(jì)算技術(shù)提供的高效率支撐下，可以采用更嚴(yán)密的方法進(jìn)行抗差計(jì)算。剔除含粗差的情況會(huì)導(dǎo)致等價(jià)權(quán)求逆失敗，本文采用降維求逆或者對(duì)剔除粗差的等價(jià)權(quán)對(duì)角元素賦固定值（如10-6）的方法。圖1是完善后經(jīng)實(shí)踐證明可行的步驟。

圖1　相關(guān)抗差估計(jì)流程Fig.1 Flow of dependent robust estimation

按照?qǐng)D1流程實(shí)現(xiàn)時(shí)，在相關(guān)抗差的迭代循環(huán)過(guò)程中判斷兩次計(jì)算的所有參數(shù)改正數(shù)之差絕對(duì)值均小于限差eps（如0.001m），則計(jì)算結(jié)束。由圖1可知，相關(guān)抗差估計(jì)流程中通過(guò)選權(quán)迭代進(jìn)行抗差計(jì)算，抗差過(guò)程主要涉及設(shè)計(jì)矩陣的構(gòu)建、法方程組的構(gòu)建與求解和等價(jià)權(quán)的更新等流程。

2.2 相關(guān)抗差估計(jì)并行算法設(shè)計(jì)

由上述分析可知，相關(guān)抗差估計(jì)的并行設(shè)計(jì)歸結(jié)為相關(guān)矩陣的構(gòu)建和矩陣運(yùn)算，而矩陣并行算法在MPI或OpenMp（Open Multi-Processin，共享存儲(chǔ)編程接口）或MPI+OpenMp環(huán)境下的研究比較多［20-22］。MPI主要關(guān)注的是在集群上運(yùn)行的應(yīng)用程序，在共享內(nèi)存的多核系統(tǒng)中會(huì)引起額外內(nèi)存開(kāi)銷(xiāo)，導(dǎo)致并行效率較低。OpenMp雖然是針對(duì)單臺(tái)計(jì)算機(jī)多核/多CPU并行計(jì)算而設(shè)計(jì)的工具，但不適合需要復(fù)雜線(xiàn)程間同步和互斥的場(chǎng)合，需要相應(yīng)的編譯器支持。微軟發(fā)布的.NET4并行擴(kuò)展將并行程序的“構(gòu)造單元”從“線(xiàn)程（thread）”提升到“任務(wù)（task）”，將任務(wù)分派給線(xiàn)程的工作由TPL負(fù)責(zé)完成，大大降低了并行計(jì)算程序的開(kāi)發(fā)難度，提升了并行開(kāi)發(fā)的效率［23］。

任何設(shè)計(jì)的并行程序都是通過(guò)并行代碼來(lái)表達(dá)的，利用并行編譯器將部分串行執(zhí)行部分轉(zhuǎn)換為并行執(zhí)行的模式，是對(duì)原有串行代碼改動(dòng)最小的并行化方法，但通常其并行效率不高。通過(guò)重新設(shè)計(jì)并行算法，建立基于任務(wù)并行的并行程序，避免了傳統(tǒng)的多線(xiàn)程重量級(jí)編程模型，而且能夠獲得較高并行效率。結(jié)合.NET的輕量級(jí)并發(fā)模型，以下對(duì)抗差估計(jì)從數(shù)據(jù)并行和任務(wù)并行兩方面進(jìn)行設(shè)計(jì)以實(shí)現(xiàn)抗差估計(jì)的并行化。

2.2.1 命令式數(shù)據(jù)并行

在.NET Framwork4中，TPL引入了基于任務(wù)并行的新名稱(chēng)空間System.Threadiing.Tasks，利用其中的Parallel靜態(tài)類(lèi)所提供的方法即可實(shí)現(xiàn)負(fù)載均衡的并行執(zhí)行，不需要從底層建立管理、控制和同步復(fù)雜的多個(gè)線(xiàn)程，也支持通過(guò)自定義分區(qū)器掌控?cái)?shù)據(jù)的分發(fā)。對(duì)于采用子網(wǎng)劃分模型處理大型GNSS網(wǎng)得到的m個(gè)基線(xiàn)向量組，各個(gè)向量組分別保存在不同的文件中。通常對(duì)m個(gè)觀(guān)測(cè)向量文件的讀取，觀(guān)測(cè)矩陣Ai、常向量li和權(quán)矩陣ˉPi的建立采用循環(huán)的方式進(jìn)行構(gòu)建，簡(jiǎn)要的串行循環(huán)版本（C#語(yǔ)言）如下：

大量基線(xiàn)向量觀(guān)測(cè)數(shù)據(jù)構(gòu)建矩陣變量的過(guò)程中，每一份基線(xiàn)向量觀(guān)測(cè)數(shù)據(jù)均執(zhí)行相同的操作，則可以采用TPL提供的數(shù)據(jù)并行形式實(shí)現(xiàn)并行化。對(duì)上述代循環(huán)的串行代碼進(jìn)行并行化，引入Parallel提供的循環(huán)并行類(lèi)，如下是利用Parallel .For的并行化版本：

由上可知，TPL不需要改變?cè)醒h(huán)內(nèi)部代碼即可實(shí)現(xiàn)對(duì)固定數(shù)目的獨(dú)立循環(huán)迭代的并行執(zhí)行，而且能夠?qū)崿F(xiàn)所有任務(wù)在大部分時(shí)間內(nèi)保持繁忙計(jì)算，實(shí)現(xiàn)了負(fù)載均衡的潛在執(zhí)行。

2.2.2 命令式任務(wù)并行

針對(duì)很多可以并發(fā)運(yùn)行的不同操作，TPL也提供通過(guò)任務(wù)并行發(fā)揮并行化優(yōu)勢(shì)的形式。采用基于任務(wù)的編程模式不需要編寫(xiě)底層的、復(fù)雜的重量級(jí)線(xiàn)程代碼。需要強(qiáng)調(diào)的是，TPL建立的任務(wù)并不是線(xiàn)程，而任務(wù)代碼的運(yùn)行需要使用底層的線(xiàn)程，而且任務(wù)和線(xiàn)程之間并不是一對(duì)一的關(guān)系。如對(duì)于m個(gè)觀(guān)測(cè)矩陣Ai、常向量li和權(quán)矩陣建立相應(yīng)的法矩陣和法向量，若每一個(gè)建立過(guò)程定義為一個(gè)任務(wù)，由于每個(gè)觀(guān)測(cè)向量的數(shù)目不盡相當(dāng)，導(dǎo)致各個(gè)任務(wù)的計(jì)算開(kāi)銷(xiāo)也各不相同，則通過(guò)運(yùn)行在多個(gè)線(xiàn)程中的同步代碼實(shí)現(xiàn)是非常復(fù)雜的?；谌蝿?wù)的編程方法將工作調(diào)度等一系列的同步問(wèn)題隱藏在幕后，提供了一種新的能夠利用多核的輕量級(jí)并行機(jī)制。假定有8個(gè)基線(xiàn)向量觀(guān)測(cè)文件，定義8個(gè)建立法矩陣和法向量的任務(wù)（task（i），i=1，2，…，8），采用任務(wù)模式的并行和并發(fā)執(zhí)行完全由底層線(xiàn)程執(zhí)行，相應(yīng)的運(yùn)行過(guò)程如圖2所示。

圖2 任務(wù)并行執(zhí)行過(guò)程圖Fig.2 Process of task parallel execution

圖2中，橫向?yàn)闀r(shí)間軸，8個(gè)計(jì)算任務(wù)通過(guò)4個(gè)線(xiàn)程并行執(zhí)行完成。最終的并行時(shí)間等于線(xiàn)程1的task1和task6兩個(gè)任務(wù)執(zhí)行時(shí)間之和，而串行時(shí)間等于8個(gè)任務(wù)分別執(zhí)行的時(shí)間之和。圖2中任務(wù)所包含的代碼會(huì)在一個(gè)線(xiàn)程中運(yùn)行，任務(wù)調(diào)度器會(huì)自動(dòng)根據(jù)所有可用的邏輯內(nèi)核嘗試對(duì)工作進(jìn)行負(fù)載均衡，在幕后使用工作竊取隊(duì)列找到一個(gè)最合算的線(xiàn)程，并將任務(wù)加入隊(duì)列中，如圖2中的task5和task3，降低了并行開(kāi)發(fā)的難度。通過(guò)使用任務(wù)對(duì)8個(gè)基線(xiàn)向量觀(guān)測(cè)文件建立8個(gè)法矩陣和法向量的任務(wù)進(jìn)行并行化代碼如下：

對(duì)于大型GNSS網(wǎng)的m個(gè)不相關(guān)的基線(xiàn)向量組，可采用并行累加的方法實(shí)現(xiàn)式（3）的計(jì)算。因此，結(jié)合TPL提供的數(shù)據(jù)并行與任務(wù)并行的形式，本文設(shè)計(jì)如圖3的基線(xiàn)向量整網(wǎng)相關(guān)抗差估計(jì)粗粒度的并行計(jì)算流程圖。

圖3　大型GNSS網(wǎng)相關(guān)抗差并行估計(jì)流程Fig.3 Process of parallel robust estimation for large GNSS network

圖3中，基于數(shù)據(jù)并行形式，讀取各個(gè)基線(xiàn)向量組的觀(guān)測(cè)文件信息，并行建立基線(xiàn)向量組誤差方程式；基于任務(wù)并行形式，對(duì)各個(gè)向量組進(jìn)行法方程并行累加計(jì)算和等價(jià)權(quán)并行更新。由于抗差估計(jì)過(guò)程中需要求取Qδ^X，因此對(duì)于上圖中法方程的求解過(guò)程需直接對(duì)法矩陣N=∑ATiˉPiAi求逆。為了進(jìn)一步提高求逆的并行效率，對(duì)法矩陣N采用分塊方法分成n×n塊（圖中n取為4為例），采用文獻(xiàn)［14］提出的分塊法矩陣Cholesky并行分解方法進(jìn)行并行Cholesky分解，將法矩陣N的逆陣的每個(gè)列向量根據(jù)Cholesky分塊分解結(jié)果進(jìn)行并行方程解算。

3　試驗(yàn)算例及結(jié)果分析

3.1 算例1

取IGS發(fā)布的1017周—1717周之間每隔100周的8個(gè)周解SINEX文件，各個(gè)周解文件IGS站數(shù)目如圖4。

圖4　IGS測(cè)站統(tǒng)計(jì)情況Fig.4 Statistics of IGS stations

由圖4可看出，隨著全球可用CORS站的增多，IGS選用的測(cè)站數(shù)也在不斷的遞增，從第1017周的134個(gè)測(cè)站遞增到第1717周的470個(gè)測(cè)站。IGS公布的周解坐標(biāo)形成的基線(xiàn)向量精度較高，可看作沒(méi)有粗差的觀(guān)測(cè)向量。以此模擬每個(gè)周的7個(gè)日解文件，并隨機(jī)對(duì)其中2%的觀(guān)測(cè)向量隨機(jī)添加3～6倍觀(guān)測(cè)值中誤差的粗差，通過(guò)7個(gè)日解文件計(jì)算出周解。用IGGⅢ權(quán)函數(shù)的雙因子等價(jià)權(quán)原理的相關(guān)抗差估計(jì)流程進(jìn)行抗差估計(jì)（抗差因子k0=1，k1=4.5，限差eps=0.005），按照本文提出的并行抗差流程進(jìn)行并行抗差估計(jì)。采用VS平臺(tái)C#語(yǔ)言實(shí)現(xiàn)抗差流程，其中TPL技術(shù)實(shí)現(xiàn)并行計(jì)算流程。試驗(yàn)的硬件平臺(tái)為神州K590S-i7電腦，CPU主頻2.40GHz，內(nèi)存4GB，物理核為4核。統(tǒng)計(jì)各周的串行抗差估計(jì)和并行抗差估計(jì)時(shí)間開(kāi)銷(xiāo)，如圖5。

圖5 串行與并行抗差估計(jì)計(jì)算時(shí)間Fig.5 Runtime of serial and parallel robust estimation

圖5中，“RS-S”代表串行抗差估計(jì)消耗時(shí)間，“RS-P”代表并行抗差估計(jì)消耗時(shí)間，縱坐標(biāo)表示計(jì)算消耗時(shí)間，單位為秒，橫坐標(biāo)表示GPS周數(shù)。由圖5可知，隨著測(cè)站規(guī)模的增加，抗差估計(jì)的消耗時(shí)間呈幾何倍遞增，第1717周傳統(tǒng)串行抗差估計(jì)需要519.32s完成，采用并行計(jì)算后僅需要168.37s，加速比達(dá)到3.08，即4核環(huán)境下并行抗差估計(jì)效率比串行抗差估計(jì)效率提高了3.08倍。這是計(jì)算一次的時(shí)效，而通常的抗差計(jì)算，經(jīng)常需要選擇不同的參數(shù)進(jìn)行數(shù)十次甚至百次的計(jì)算。

3.2 算例2

取IGS公布的第1721周的7個(gè)日解SINEX文件，每個(gè)文件作為一個(gè)同步網(wǎng)基線(xiàn)文件，對(duì)其進(jìn)行聯(lián)合平差。IGS公布的日解坐標(biāo)形成的基線(xiàn)向量的精度較高，可將原始觀(guān)測(cè)向量看作沒(méi)有粗差的觀(guān)測(cè)向量。由于IGS每天選取的站點(diǎn)是不固定的，第1721周7個(gè)日解文件的公共點(diǎn)共365個(gè)，為方便比較，僅對(duì)公共點(diǎn)進(jìn)行處理，取其中1個(gè)點(diǎn)的周解作為已知點(diǎn)，其余點(diǎn)為未知點(diǎn)。隨機(jī)對(duì)其中1%的觀(guān)測(cè)向量隨機(jī)加上3～6倍觀(guān)測(cè)值中誤差的粗差，然后用最小二乘估計(jì)和IGGⅢ權(quán)函數(shù)的雙因子等價(jià)權(quán)原理建立的相關(guān)抗差估計(jì)進(jìn)行平差計(jì)算，其中抗差估計(jì)分別進(jìn)行串行估計(jì)和并行估計(jì)。將IGS發(fā)布的周解坐標(biāo)值作為真值進(jìn)行比較。試驗(yàn)的硬件平臺(tái)與算例1相同。

最小二乘估計(jì)與抗差估計(jì)（參數(shù)t分別取1和2）的計(jì)算時(shí)間、驗(yàn)后單位權(quán)中誤差^σ0統(tǒng)計(jì)如表1所示。

表1　最小二乘估計(jì)與抗差估計(jì)的計(jì)算時(shí)間與驗(yàn)后單位中誤差統(tǒng)計(jì)表Tab.1 Statistics of runtime and0for the least squares estimation and robust estimation

表1　最小二乘估計(jì)與抗差估計(jì)的計(jì)算時(shí)間與驗(yàn)后單位中誤差統(tǒng)計(jì)表Tab.1 Statistics of runtime and0for the least squares estimation and robust estimation

^σ0串行時(shí)間/s并行時(shí)間/s加速比LS 4.09 　39.25 　14.01 　2.80 　Robust（t=1）1.29 　339.14 　121.44 　2.79 　Robust（t=2）1.30 　208.40 　73.79 　2.82

由表1可知，有粗差時(shí)，抗差估計(jì)的驗(yàn)后單位權(quán)中誤差均小于最小二乘平差的驗(yàn)后單位權(quán)中誤差，反映出采用等價(jià)權(quán)法抗差估計(jì)的平差精度明顯優(yōu)于最小二乘平差，但抗差計(jì)算時(shí)間大大增加。當(dāng)參數(shù)t選取不同的值時(shí)，抗差效率也互不相同。當(dāng)t=1時(shí)抗差串行計(jì)算時(shí)間為339.14s，而當(dāng)t= 2時(shí)抗差串行計(jì)算時(shí)間為208.40s。采用本文設(shè)計(jì)的并行策略后，并行后的計(jì)算效率都比串行計(jì)算效率得到明顯提高，t=1時(shí)，并行時(shí)間為121.44s，t=2時(shí)，并行時(shí)間僅為73.79s。針對(duì)當(dāng)前算例，t=2時(shí)比t=1時(shí)的抗差效率更高，而抗差效果基本一致，因此，針對(duì)不同的抗差問(wèn)題，應(yīng)選擇具體的抗差參數(shù)以確保抗差效率與精度的最優(yōu)。不同估計(jì)方法與周解的坐標(biāo)差比較如圖6所示，其中，縱坐標(biāo)表示不同估計(jì)方法得出的坐標(biāo)值與周解坐標(biāo)值之間的差，單位為米。

圖6　最小二乘、不同參數(shù)的抗差估計(jì)與周解的坐標(biāo)分量差值比較Fig.6 Coordinate difference comparison for least square and robust solutions

從圖6中可看出，在有粗差情況下，最小二乘平差的坐標(biāo)分量與周解的坐標(biāo)分量的差值在部分點(diǎn)的分量有明顯差異，最大達(dá)到1.8cm。當(dāng)選取恰當(dāng)?shù)呐R界值，采用等價(jià)權(quán)法抗差估計(jì)方案有效地抑制了粗差的影響，得到較優(yōu)的坐標(biāo)分量，與周解的坐標(biāo)分量差值基本都在5mm以?xún)?nèi)。

由于大規(guī)模GNSS數(shù)據(jù)處理本身的復(fù)雜性，以及對(duì)粗差問(wèn)題認(rèn)識(shí)的局限性，任何一種抗差估計(jì)方法的實(shí)際抗差能力都是有一定條件的［6］，如粗差的分布情況、粗差的數(shù)值大小和數(shù)目多少等。限于篇幅，本文不詳細(xì)討論此問(wèn)題。本文進(jìn)行了大量抗差估計(jì)的試驗(yàn)，試驗(yàn)時(shí)并行抗差效率都比串行抗差效率得到提高。

4　結(jié) 論

由于科學(xué)研究及應(yīng)用的需要，全球CORS站網(wǎng)規(guī)模不斷增加，大型GNSS網(wǎng)絡(luò)數(shù)據(jù)處理越來(lái)越頻繁，對(duì)常規(guī)的GNSS數(shù)據(jù)存儲(chǔ)、計(jì)算等技術(shù)帶來(lái)一系列新的挑戰(zhàn)。針對(duì)含有粗差的基線(xiàn)向量網(wǎng)抗差計(jì)算，頻繁的抗差計(jì)算時(shí)間降低了科學(xué)研究的效率，而采用并行計(jì)算技術(shù)實(shí)現(xiàn)相關(guān)抗差估計(jì)的并行計(jì)算，不僅能顯著提高抗差效率，而且提高了現(xiàn)有計(jì)算資源的利用效率。本文的主要工作有：

（1）將相關(guān)抗差雙因子等價(jià)權(quán)模型中的降權(quán)因子的平方改為可變參數(shù)t，同時(shí)完善了殘差協(xié)方差矩陣的求解公式，在并行計(jì)算技術(shù)的高效率支撐下，采用更嚴(yán)密的抗差計(jì)算步驟進(jìn)行抗差估計(jì)。

（2）從數(shù)據(jù)并行和任務(wù)并行兩方面實(shí)現(xiàn)了抗差估計(jì)的并行化，基于數(shù)據(jù)并行形式并行讀取觀(guān)測(cè)文件信息和建立基線(xiàn)向量組誤差方程，基于任務(wù)并行形式對(duì)各個(gè)向量組進(jìn)行法方程并行累加計(jì)算和等價(jià)權(quán)并行更新，設(shè)計(jì)了一種粗粒度的相關(guān)抗差估計(jì)并行計(jì)算流程。

（3）利用TPL并行計(jì)算技術(shù)建立了基于任務(wù)劃分的相關(guān)抗差估計(jì)并行算法，不僅避免了傳統(tǒng)多線(xiàn)程重量級(jí)編程模型，降低了并行計(jì)算程序的開(kāi)發(fā)難度，而且能獲得較高的抗差并行效率，通過(guò)試驗(yàn)算例驗(yàn)證了并行抗差估計(jì)算法的可靠性和高效性。

（4）算例表明并行計(jì)算技術(shù)在不改變抗差結(jié)果的前提下能有效提高抗差計(jì)算效率，而抗差模型參數(shù)的選擇嚴(yán)重影響抗差計(jì)算時(shí)間，在保證抗差穩(wěn)定和精度的前提下，應(yīng)將并行計(jì)算技術(shù)和高效的抗差模型與步驟相結(jié)合，以保證大型GNSS基線(xiàn)網(wǎng)抗差效率與精度的最優(yōu)。

［1］ GEOFFERY B.Fixed Point Theorems of GPS Carrier Phase Ambiguity Resolution and Their Application to Massive Network Processing：Ambizap［J］.Journal of Geophysical Research，2008.doi：10.1029/JB005736.

［2］ SHI Chuang.Large-scale High-precision GPS Network Adjustment and Analysis Theory and Its Application［M］. Beijing：Surveying and Mapping Press，2002.（施闖.大規(guī)模高精度GPS網(wǎng)平差與分析理論及其應(yīng)用［M］.北京：測(cè)繪出版社，2002.）

［3］ YANG Yuanxi.Robust Estimation Theory and Its Applications［M］.Beijing：Publishing House of PLA，1993.（楊元喜.抗差估計(jì)理論及其應(yīng)用［M］.北京：解放軍出版社，1993.）

［4］ YU Xuexiang，LüWeicai.Robust Estimation Model for Correlated Observations Based on Standardized Residuals［J］.Journal of Wuhan Technical University of Surveying and Mapping，1999，24（10）：75-78.（余學(xué)祥，呂偉才.基于標(biāo)準(zhǔn)化殘差的相關(guān)觀(guān)測(cè)抗差估計(jì)模型［J］.武漢測(cè)繪科技大學(xué)學(xué)報(bào)，1999，24（10）：75-78.）

［5］ WU Jiangfei，YANG Yuanxi.Robust Estimation for Correlated GPS Baseline Vector Network［J］.Acta Geodaetica et Cartographica Sinica，2001，30（3）：247-251.（吳江飛，楊元喜.相關(guān)GPS基線(xiàn)向量網(wǎng)的抗差估計(jì)［J］.測(cè)繪學(xué)報(bào)，2001，30（3）：247-251.）

［6］ YANG Yuanxi，SONG Lijie，XU Tianhe.Robust Parameter Estimation for Geodetic Correlated Observations［J］.Acta Geodaetica et Cartographica Sinica，2002，31（2）：95-100.（楊元喜，宋力杰，徐天河.大地測(cè)量相關(guān)觀(guān)測(cè)抗差估計(jì)理論［J］.測(cè)繪學(xué)報(bào)，2002，31（2）：95-100.）

［7］ YANG Yuanxi，SONG Lijie，XU Tianhe.New Robust Estimator for the Adjustment of Correlated GPS Networks［C］∥Proceedings of First International Symposium on Robust Statistics and Fuzzy Techniques in Geodesy and GIS：295.Berlin：IGP-Bericht，2001：199-208.

［8］ ZHAO Qinghai.Robust Estimation for High Precision GPS Vector Network［J］.Acta Geodaetica et Cartographica Sinica，2004，33（1）：43-46.（趙慶海.高精度GPS向量網(wǎng)抗差估計(jì)［J］.測(cè)繪學(xué)報(bào)，2004，33（1）：43-46.）

［9］ SONG Lijie，OUYANG Guichong.A Fast Method of Solving Partitioned Adjustment for Super Large-scale Geodetic Network［J］.Acta Geodaetica et Cartographica Sinica，2003，32（3）：204-207.（宋力杰，歐陽(yáng)桂崇.超大規(guī)模大地網(wǎng)分區(qū)平差快速解算方法［J］.測(cè)繪學(xué)報(bào)，2003，32（3）：204-207.）

［10］ HENNO B.Distributed Processing for Large Geodetic Solutions.Reference Frames for Application in Geosciences［C］∥Proceedings of the International Association of Geodesy Symposia：138.Berlin：Springer-Verlag，2013：102-114.

［11］ HIROSH T，TETSURO K，YASUHIRO K，et al.A VLBI Correlator with Internet-based Distributed Computing［EB/OL］.［2014-04-16］.http：∥www.ursi.org/proceedings/procga05/pdf/j06.10（0886）.pdf.

［12］ MARC B，LUC G，SERGE G，et al.Parallel Tools for Solving Incremental Dense Least Squares Problem：Application to Space Geodesy［J］.Journal of Algorithms &Computational Technology，2009，3（1）：117-133.

［13］ ZOU Xiancai，LI Jiancheng，WANG Haihong，et al. Application of Parallel Computing with OpenMP in Data Processing for Satellite Gravity［J］.Acta Geodaetica et Cartographica Sinica，2010，39（6）：636-641.（鄒賢才，李建成，汪海洪，等.OpenMP并行計(jì)算在衛(wèi)星重力數(shù)據(jù)處理中的應(yīng)用［J］.測(cè)繪學(xué)報(bào)，2010，39（6）：636-641.）

［14］ CUI Yang，LüZhiping，CHEN Zhengsheng，et al.Research of Parallel Data Processing for GNSS Network Adjustment under Multi-core Environment［J］.Acta Geodaetica et Cartographica Sinica，2013，42（5）：661-667.（崔陽(yáng)，呂志平，陳正生，等.多核環(huán)境下的GNSS網(wǎng)平差數(shù)據(jù)并行處理研究［J］.測(cè)繪學(xué)報(bào)，2013，42（5）：661-667.）

［15］ LüZhiping，CHEN Zhengsheng，CUI Yang，et al.DistributedProcessing of Large-scale CORS Network’s Baseline Vectors［J］.Journal of Geomatics Science and Technology，2013，30（4）：433-438.（呂志平，陳正生，崔陽(yáng)，等.大型CORS網(wǎng)基線(xiàn)向量的分布式處理［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2013，30（4）：433-438.）

［16］ DONG D，F(xiàn)ANG P，BOCK Y，et al.Anatomy of Apparent Seasonal Variations from GPS-derived Site Position Time Series［J］.Journal of Geophysical Research，2002，107（B4）：2075.doi：10.1029/2001jb000573.

［17］ HERRING T.Aspects of Large Station Networks for GPS Orbits and Clocks［C］∥Proceedings of IGS Workshop and Symposium.Berne：IGS，2004.

［18］ YANG Ling，SHEN Yunzhong，LOU Lizhi.Equivalent Weight Robust Estimation Method Based on Median Parameter Estimates［J］.Acta Geodaetica et Cartographica Sinica，2011，40（1）：28-32.（楊玲，沈云中，樓立志.基于中位參數(shù)數(shù)初值的等價(jià)權(quán)抗差估計(jì)方法［J］.測(cè)繪學(xué)報(bào)，2011，40（1）：28-32.）

［19］ LI Yan，HU Xiangyun，JIN Gangxie，et al.Research of 1-D Magnetotelluric Parallel Computation Based on MPI［J］.Progress in Geophysics，2010，25（5）：1612-1616.（李焱，胡祥云，金剛燮，等.基于MPI的一維大地電磁并行計(jì)算研究［J］.地球物理學(xué)進(jìn)展，2010，25（5）：1612-1616.）

［20］ NIE Linjuan，SHEN Wenbin，WANG Zhengta，et al. Application of Parallel Solutions Technology in Satellite Gravity Measurement［J］.Journal of Geodesy and Geodynamics，2012，32（2）：64-73.（聶琳娟，申文斌，王正濤，等.基于超級(jí)計(jì)算機(jī)平臺(tái)的并行解技術(shù)在衛(wèi)星重力測(cè)量中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012，32（2）：64-73.）

［21］ YANG Kai.Research on High Performance GNSS Normal Equation Processing Method Based on OpenMP［J］. Journal of Geodesy and Geodynamics，2013，33（1）：141-144.（楊凱.基于OpenMp的高性能GNSS法方程解算方法研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2013，33（1）：141-144.）

［22］ HILLAR G C.Professional Parallel Programming with C#：Master Parallel Extensions with.NET4［M］. Beijing：Tsinghua University Press，2012.（HILLAR G C.C#并行編程高級(jí)教程：精通.NET4Parallel Extensions［M］.北京：清華大學(xué)出版社，2012.）

（責(zé)任編輯：叢樹(shù)平）

E-maiI：ssscenter@126.com

ParaIIeI Robust Estimation for Large-scaIe High-precision GNSS BaseIine Vector Network

CUI Yang1，LüZhiping1，ZHANG Youyang2，LI Linyang1
1.Institute of GeospatiaI Information，Information Engineering University，Zhengzhou 450001，China；2.Institute of MechanicaI and EIectricaI Engineering，Zhengzhou University of Light Industry，Zhengzhou 450002，China

According to the characteristics of Iarge-scaIe GNSS baseIine vector network，the reIevant observations robust estimation modeI which is based on RECO scheme is improved by paraIIeI computing.A feasibIe process steps are aIso presented.Using the improved modeI，a numericaI exampIe is given.It is demonstrated that the Iarge-scaIe GNSS baseIine vector reIated robust paraIIeI estimation not onIy can inhibit the infIuence of the outIiers to the parameter estimates，but aIso significantIy improve the computationaI efficiency.

dependent robust estimation；GNSS；baseIine vector network；equivaIent weight matrix；paraIIeI estimation

The NationaI NaturaI Science Foundation of China（No.41274015）；The NationaI Hightech Research and DeveIopment Program of China（863 Program）（No.2013AA122501）

CUI Yang（1988—），maIe，PhD candidate，majors in the theory and method of Iarge-scaIe survey data processing.

LüZhiping

P228

A

1001-1595（2015）05-0495-08

國(guó)家自然科學(xué)基金（41274015）；國(guó)家863計(jì)劃（2013AA122501）

CUI Yang，LüZhiping，ZHANG Youyang，et al.Parallel Robust Estimation for Large-scale High-precision GNSS Baseline Vector Network［J］.Acta Geodaetica et Cartographica Sinica，2015，44（5）：495-502.（崔陽(yáng)，呂志平，張友陽(yáng)，等.大型高精度GNSS基線(xiàn)向量網(wǎng)并行抗差估計(jì)［J］.測(cè)繪學(xué)報(bào)，2015，44（5）：495-502.）

10.11947/j.AGCS.2015.20140154

2014-03-28

崔陽(yáng)（1988—），男，博士生，研究方向?yàn)榇笠?guī)模測(cè)量數(shù)據(jù)處理方法與理論。

呂志平

修回日期：2014-07-09

E-maiI：cuiyangtj@126.com