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

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

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

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

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

1　引 言

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

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

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

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

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

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

式中

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

式（7）單位權(quán)中誤差計(jì)算公式中，n′等于觀測量個(gè)數(shù)n減去淘汰觀測量（γ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的求解過程涉及等價(jià)權(quán)的求逆，而迭代過程中求逆非常耗時(shí)，通常采用不嚴(yán)密的方法進(jì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)過程中判斷兩次計(jì)算的所有參數(shù)改正數(shù)之差絕對(duì)值均小于限差eps（如0.001m），則計(jì)算結(jié)束。由圖1可知，相關(guān)抗差估計(jì)流程中通過選權(quán)迭代進(jìn)行抗差計(jì)算，抗差過程主要涉及設(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)存開銷，導(dǎo)致并行效率較低。OpenMp雖然是針對(duì)單臺(tái)計(jì)算機(jī)多核/多CPU并行計(jì)算而設(shè)計(jì)的工具，但不適合需要復(fù)雜線程間同步和互斥的場合，需要相應(yīng)的編譯器支持。微軟發(fā)布的.NET4并行擴(kuò)展將并行程序的“構(gòu)造單元”從“線程（thread）”提升到“任務(wù)（task）”，將任務(wù)分派給線程的工作由TPL負(fù)責(zé)完成，大大降低了并行計(jì)算程序的開發(fā)難度，提升了并行開發(fā)的效率［23］。

任何設(shè)計(jì)的并行程序都是通過并行代碼來表達(dá)的，利用并行編譯器將部分串行執(zhí)行部分轉(zhuǎn)換為并行執(zhí)行的模式，是對(duì)原有串行代碼改動(dòng)最小的并行化方法，但通常其并行效率不高。通過重新設(shè)計(jì)并行算法，建立基于任務(wù)并行的并行程序，避免了傳統(tǒng)的多線程重量級(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ù)并行的新名稱空間System.Threadiing.Tasks，利用其中的Parallel靜態(tài)類所提供的方法即可實(shí)現(xiàn)負(fù)載均衡的并行執(zhí)行，不需要從底層建立管理、控制和同步復(fù)雜的多個(gè)線程，也支持通過自定義分區(qū)器掌控?cái)?shù)據(jù)的分發(fā)。對(duì)于采用子網(wǎng)劃分模型處理大型GNSS網(wǎng)得到的m個(gè)基線向量組，各個(gè)向量組分別保存在不同的文件中。通常對(duì)m個(gè)觀測向量文件的讀取，觀測矩陣Ai、常向量li和權(quán)矩陣ˉPi的建立采用循環(huán)的方式進(jìn)行構(gòu)建，簡要的串行循環(huán)版本（C#語言）如下：

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

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

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

對(duì)于大型GNSS網(wǎng)的m個(gè)不相關(guān)的基線向量組，可采用并行累加的方法實(shí)現(xiàn)式（3）的計(jì)算。因此，結(jié)合TPL提供的數(shù)據(jù)并行與任務(wù)并行的形式，本文設(shè)計(jì)如圖3的基線向量整網(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è)基線向量組的觀測文件信息，并行建立基線向量組誤差方程式；基于任務(wù)并行形式，對(duì)各個(gè)向量組進(jìn)行法方程并行累加計(jì)算和等價(jià)權(quán)并行更新。由于抗差估計(jì)過程中需要求取Qδ^X，因此對(duì)于上圖中法方程的求解過程需直接對(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測站統(tǒng)計(jì)情況Fig.4 Statistics of IGS stations

由圖4可看出，隨著全球可用CORS站的增多，IGS選用的測站數(shù)也在不斷的遞增，從第1017周的134個(gè)測站遞增到第1717周的470個(gè)測站。IGS公布的周解坐標(biāo)形成的基線向量精度較高，可看作沒有粗差的觀測向量。以此模擬每個(gè)周的7個(gè)日解文件，并隨機(jī)對(duì)其中2%的觀測向量隨機(jī)添加3～6倍觀測值中誤差的粗差，通過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#語言實(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í)間開銷，如圖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可知，隨著測站規(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)基線文件，對(duì)其進(jìn)行聯(lián)合平差。IGS公布的日解坐標(biāo)形成的基線向量的精度較高，可將原始觀測向量看作沒有粗差的觀測向量。由于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%的觀測向量隨機(jī)加上3～6倍觀測值中誤差的粗差，然后用最小二乘估計(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ì)不同的抗差問題，應(yīng)選擇具體的抗差參數(shù)以確?？共钚逝c精度的最優(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以內(nèi)。

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

4　結(jié) 論

由于科學(xué)研究及應(yīng)用的需要，全球CORS站網(wǎng)規(guī)模不斷增加，大型GNSS網(wǎng)絡(luò)數(shù)據(jù)處理越來越頻繁，對(duì)常規(guī)的GNSS數(shù)據(jù)存儲(chǔ)、計(jì)算等技術(shù)帶來一系列新的挑戰(zhàn)。針對(duì)含有粗差的基線向量網(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ù)并行形式并行讀取觀測文件信息和建立基線向量組誤差方程，基于任務(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)多線程重量級(jí)編程模型，降低了并行計(jì)算程序的開發(fā)難度，而且能獲得較高的抗差并行效率，通過試驗(yàn)算例驗(yàn)證了并行抗差估計(jì)算法的可靠性和高效性。

（4）算例表明并行計(jì)算技術(shù)在不改變抗差結(jié)果的前提下能有效提高抗差計(jì)算效率，而抗差模型參數(shù)的選擇嚴(yán)重影響抗差計(jì)算時(shí)間，在保證抗差穩(wěn)定和精度的前提下，應(yīng)將并行計(jì)算技術(shù)和高效的抗差模型與步驟相結(jié)合，以保證大型GNSS基線網(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］.北京：測繪出版社，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)觀測抗差估計(jì)模型［J］.武漢測繪科技大學(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基線向量網(wǎng)的抗差估計(jì)［J］.測繪學(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.（楊元喜，宋力杰，徐天河.大地測量相關(guān)觀測抗差估計(jì)理論［J］.測繪學(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］.測繪學(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.（宋力杰，歐陽桂崇.超大規(guī)模大地網(wǎng)分區(qū)平差快速解算方法［J］.測繪學(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］.測繪學(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.（崔陽，呂志平，陳正生，等.多核環(huán)境下的GNSS網(wǎng)平差數(shù)據(jù)并行處理研究［J］.測繪學(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.（呂志平，陳正生，崔陽，等.大型CORS網(wǎng)基線向量的分布式處理［J］.測繪科學(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］.測繪學(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)星重力測量中的應(yīng)用［J］.大地測量與地球動(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］.大地測量與地球動(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é)任編輯：叢樹平）

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

國家自然科學(xué)基金（41274015）；國家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.（崔陽，呂志平，張友陽，等.大型高精度GNSS基線向量網(wǎng)并行抗差估計(jì)［J］.測繪學(xué)報(bào)，2015，44（5）：495-502.）

10.11947/j.AGCS.2015.20140154

2014-03-28

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

呂志平

修回日期：2014-07-09

E-maiI：cuiyangtj@126.com