龐正輝,張品超,李金濤,吳小串,李群喜,鄧廷銀
(1. 中工武大設計集團有限公司,湖北 武漢 430079;2. 南昌市自然資源和規劃局,江西 南昌 330038)
測繪基準體系是國民經濟、社會發展、國家安全和信息化建設的重要基礎,主要包括大地基準、高程基準和重力基準。2000年已有學者開始研究我國現代大地測量基準建設工作,2008年我國發布使用CGCS2000 坐標系。國家現代測繪基準體系基礎設施建設工程于2012年啟動,2017年5 月通過驗收,建成了高精度、三維、動態、幾何與物理基準融合統一的我國現代測繪基準體系[1-6]。我國全部省級區域、部分城市群和多數地市已進行區域現代測繪基準體系建設,形成了一些研究成果[7-10]。2015年江西省啟動了江西省現代大地基準完善項目,并于2017年通過驗收,建立了與國家基準協調一致的高精度、三維地心、動態實用的江西省現代測繪基準體系[11]。2017年南昌市建成了南昌市連續運行衛星定位服務系統(NCCORS),建立并啟用推廣南昌2000 坐標系[12];同步啟動南昌市市域現代測繪基準體系建設工程,主要包括南昌市市域平面和高程基礎控制網的復測與聯測項目以及南昌市似大地水準面精化模型建立項目,并分別于2018年、2020年通過驗收,建立了既與國家基準體系協調一致,又能滿足南昌市市域城鄉經濟建設、重大工程施工需求的高精度、三維地心、動態實用的市域現代測繪基準體系,極大地提高了南昌市測繪地理信息保障服務能力[13-14]。
南昌市市域已有測繪基準不統一、現勢性較差,覆蓋范圍有限,平面、高程成果分離,高精度CORS定位三維成果與二維坐標框架不匹配,繁瑣的低等級水準測量制約了測繪生產效率,無法滿足高質量經濟社會發展的需求。已有基準主要為CGCS2000坐標系、南昌2000 坐標系、1954年北京坐標系、1980 西安坐標系、南昌國土坐標系和南昌城建坐標系。江西省連續運行衛星定位服務系統(JXCORS)和NCCORS分別于2009年和2017年建成并投入試運行,在南昌市市域及其周邊設有國家和江西省共用基準站1 座、江西省基準站2座、江西省與南昌市共用基準站8座、南昌市基準站1座。現有國家布設的GNSS B級點、一等水準點,江西省布設的GNSS B級點、C級點和二等水準點。
在充分利用已有測繪基礎設施與成果的條件下,本文采用精密水準測量、衛星定位測量等現代大地測量技術,通過新建、整合、改造和完善等方式,布設加密全市域GNSS B級、C級點和二等水準點[15];聯測確定全市各平面坐標系間的有效聯系;結合重力數據、地形資料和重力場模型,建立覆蓋南昌市域的高精度似大地水準面模型;推廣應用南昌2000 坐標系;基于NCCORS平臺,實現全市域高精度實時三維衛星定位服務。南昌市市域現代測繪基準體系見圖1。

圖1 南昌市市域現代測繪基準體系圖
1)控制網布測。控制網采用南昌2000 坐標系,1985國家高程基準。現有9 座NCCORS 基準站采用厘米級實時動態定位新技術,建成了以北斗為主、兼容GPS、GLONASS 和Galileo 等系統的多模CORS 系統,并利用JXCORS 基準站實現了站網服務覆蓋全市域及其鄰域一定范圍。
基于CORS網,利用市域GNSS B級點,加密布設B 級控制網點,形成的骨干網覆蓋全市域及其鄰域一定范圍;基于骨干網,利用市域GNSS C 級點,加密布設C級控制網點,利用市域國家一等和二等水準點成果,布設二等水準網并進行聯測,形成的基本網覆蓋全市域及其鄰域一定范圍。
2)南昌2000 坐標系的應用推廣。通過聯測已有GNSS控制點,求取南昌2000 坐標系與1954年北京坐標系、1980西安坐標系、南昌國土和城建坐標系的轉換關系,并研制了坐標轉換軟件。基于此,開展市域有關單位現有測繪成果和相關專題數據庫的數據轉換工作,促進全市域各類測繪地理信息成果共享,避免重復測繪、節約財政資金。
3)市域高精度似大地水準面模型的建立。以基本控制網的GNSS/水準共用點成果為基礎,根據收集的重力數據,融合兩類似大地水準面高,建立南昌市域2′×2′格網高精度大地水準面模型,模型精度優于±5 mm,開發完成NCGEOID系統應用軟件,插值有效范圍為南昌市市域外擴5 km區域。市域似大地水準面等值圖見圖2。

圖2 南昌市市域2′×2′似大地水準面等值圖
穩定的控制點仍可能存在極其緩慢微小的基礎沉降或反彈,雖不對控制點精度造成影響,但因GNSS和水準外業觀測存在時間差,若能有效控制GNSS 與水準觀測的同期性,將使這種微小的豎向變化的影響降到最低。本文的GNSS 和水準外業觀測方案經嚴格設計后實際作業時間控制在2 個月以內,可為似大地水準面內業計算提供變化量極小的大地高和水準高同期數據,有效提高似大地水準面模型精度。
1)基線處理策略。B、C 級網基線處理采用GAMIT 軟件,考慮了多頻多模GNSS 數據的綜合處理,包括BDS、GPS和GLONASS的相關觀測數據,優化了各觀測值的權比,最終得到融合一致的高可靠性基線向量。基線解算主要考慮的因素為:①衛星鐘差模型改正用廣播星歷中的鐘差參數;②接收機鐘差模型改正用根據偽距觀測值計算得到的鐘差;③電離層折射影響用LC 觀測值消除;④對流層折射根據標準大氣模型用薩斯坦莫寧模型改正,采用分段線形的方法估算折射量偏差數;⑤衛星相位中心偏差采用IGS_08模型,接收機天線相位中心偏差采用GPS:IGS_08模型;⑥測站位置的潮汐改正;⑦觀測值采用雙頻消電離層組合觀測值;⑧數據解算模式為周跳自動修復技術;⑨B 級網考慮衛星軌道誤差,C 級網不考慮,即松馳IGS軌道和固定IGS軌道。
2)重復基線、同步環和異步環精度統計。B、C級網分別有536 組和924 組重復基線,整網的基線向量重復性見表1,可以看出,B、C級網基線各分量整體重復性的固定誤差和比例誤差較小,基線整體解算精度較高。由于GAMIT軟件采用的是網解(即全組合解),其同步環閉合差在基線解算時已進行了分配,可以解的NRMS值作為檢驗同步環質量好壞的一個指標,一般要求B、C 級網的NRMS 值分別小于0.3 和0.5。本文B、C級網同步環NRMS最大值分別為0.194 7、0.220 8,說明GNSS 網的整體外業觀測質量好,基線解的精度高。異步環閉合差反映的是整個GNSS 網的外業觀測質量和基線解算質量的可靠性,相對于同步環閉合差,異步環閉合差對GNSS 成果質量更重要。本文B、C 級網中分別檢核最簡異步環201 個和274 個,結果見表2,可以看出,B、C 級網所有異步環精度均遠優于規范要求,基線解算精度優異。

表1 基線向量重復性統計表/mm+10-8

表2 異步環精度統計表
4.3.1 平差軟件和觀測量
三維平差采用武漢大學編制的POWERNET 科研版軟件,以各同步觀測網的獨立基線向量及其全協方差矩陣作為觀測量。
4.3.2 平差策略
以兼容性好的江西和南昌CORS站點的CGCS2000成果為基準,對B 級網進行平差,再以B 級網平差結果為基準,對C級網點進行平差。考慮基于現勢框架歷元的GNSS網觀測值與CGCS2000坐標系下的約束點可能存在系統性差異,本文創新性地引入空間平差的思想,尤其是在約束平差數據處理中重點考慮了不同基準間的尺度和旋轉等系統誤差影響,將系統誤差分為尺度基準系統誤差和方位基準系統誤差,并估計各種系統誤差參數。
1)顧及基線向量尺度基準系統誤差的模型。基線向量尺度基準系統誤差是指GNSS 網基線向量整體尺度上的系統性誤差。在整體平差時,采用附加未知尺度參數的函數模型可消除該誤差。假設GNSS 網相對于平差基準的系統誤差為m,則基線向量的誤差方程為:
式中,VΔX、VΔY、VΔZ為觀測值改正數;dXi、dYi、dZi和dXj、dYj、dZj分別為i、j兩點坐標從近似值到平差值的改正數;ΔXij、ΔYij、ΔZij為i、j兩點的近似坐標差值;為基線向量觀測值。
通過整體平差即可得到消除GNSS 網整體尺度誤差的平差結果,并同時求出觀測量相對于平差基準的系統尺度誤差。
2)顧及基線向量方位基準系統誤差的模型。基線向量方位基準系統誤差是指整個GNSS 網在空間直角坐標系中整體旋轉的系統誤差,可表述為繞X、Y、Z軸的旋轉角。整體平差時,在函數模型中加入這3個旋轉參數可消除該誤差。假設GNSS網相對于平差基準的旋轉參數為wx、wy、wz,則基線向量的誤差方程為:
3)同時顧及基線向量尺度基準和方位基準系統誤差的模型。其誤差方程為:
4.3.3 平差精度統計
本文GNSS 網約束平差最終獲得了高精度的平差成果,點位精度統計見表3、4。

表3 點位精度中誤差統計/m

表4 點位精度中誤差區間個數統計/個
4.4.1 確定高分辨率高精度地面局部重力場的總體思路
1)構建高分辨率(2′×2′)地面重力異常格網數字模型。采用地形均衡重力歸算,以高分辨率地形均衡重力異常為基礎進行擬合內插,擬合插值函數采用連續曲率張力樣條函數。該方法對實測重力數據中包含的地形影響敏感,在山區均衡重力異常變化反映地形變化,張力樣條函數對局部地形起伏有較好的自適應特性,特別適宜在重力數據稀疏地區的擬合內插。
2)精化地形影響。采用Helmert 第二類(地形)凝集法,將地形質量垂向壓縮凝集在大地水準面上成薄層,這種地形移去—恢復模式優于其他模式,且物理概念簡明,能較方便地計算地形間接影響。
3)精化各類地形影響的球面積分。傳統積分采用平面近似,隨后發展為球面近似,但均忽略了其中球面曲率的影響項,研究推導嚴格顧及該項的精密球面積分,可使精度提高到亞厘米級。
4)精化局部重力場短波和甚短波分量。將確定局部重力場擾動位的Stokes積分公式等價分解為多分辨率表達形式,即S(Δg+δA)=S(Δg)+S(δA),S(?)為Stokes積分算子,Δg為較低分辨率的格網值,δA為地形直接影響。計算S(δA)時Stokes 核函數采用高階或甚高階球諧展開級數,利用高分辨率DEM計算高分辨率或甚高分辨率δA格網數據,可使S(δA)的分辨率達到10 km,甚至更高,可最大限度地恢復局部地形變化的貢獻。
4.4.2 球面近似和橢球改正
在局部重力場的計算中,由于積分計算的范圍較小,當只要求米級或分米級精度時,為簡化計算,通常又將球面積分作平面近似,忽略了球面彎曲,即地球曲率影響;當要求厘米級或更優的精度時,必須采用更嚴密的球面積分。Stokes-Helmert 邊值問題涉及地球表面和大地水準面(或似大地水準面)兩個邊界面,都是非常復雜的曲面,數學上不可能用解析式精確表達,因此在這種高粗糙度的表面上實施Stokes 積分,無法給出其積分核函數的封閉解析表達式,必須將其簡化為一個規則光滑的曲面。大地水準面代表地球總體形狀,粗略一點接近球體,準確一點接近一個旋轉橢球體,但要構建橢球面上Stokes 積分核函數也很困難,因此現有的計算界面上擾動位T的積分公式以及計算地形影響dV 和dA 的積分都是球面近似公式。同時,各種計算公式中涉及的空間點位、方向和觀測量都是在以參考橢球為基準面的大地坐標系中(基準方向是橢球面法線方向)定義的,而采用球面近似的大地水準面相對以橢球面為基準的大地水準面將受到歪曲,觀測的重力方向也偏離了橢球法線方向,產生垂線偏差,使在大地坐標系中計算的點位和方位角受到影響。利用球面近似計算的擾動位去計算垂線偏差又受到間接影響。理論上參考橢球不能取球面代替,在高精度全球和局部重力場的逼近中所涉及的量值,特別是作為基礎數據的重力異常Δg都應加橢球改正,本文著重討論Δg力的橢球改正[16-18]。Δg的橢球改正包括垂線偏差影響(εp)和橢球法向與地心矢徑方向之間夾角的影響(εe)。
式中,ME、NE分別為子午圈和卯酉圈的曲率半徑。
圖3是在局部直角坐標系中表示垂線偏差,X軸指向北,Y軸指向東,Z軸與橢球法線重合;eφ、eλ、eh分別為對應坐標軸的單位向量;u為總垂線偏差;eg為重力方向的單位向量,在該坐標系中表示為:

圖3 垂線偏差
將式(5)代入式(4)可得:
Δg是按其定義顧及垂線偏差影響的理論值,但作為Stokes 積分邊值的重力異常觀測值Δg'未顧及垂線偏差影響,直接采用重力異常向量(gP-γQ)在橢球法向上的分量,即
顧及cosu×ξ≈ξ,cosu×η≈η,則有:
圖4中eh為橢球面法線外方向上的單位向量;er為地心矢徑方向上的單位向量;eθ為在子午面內與er正交的單位向量,指向地心余緯θ增加的方向,則有:

圖4 橢球面法線方向與地心矢徑方向的偏差
利用式(9)可將屬于橢球法線的向量轉換為屬于橢球矢徑的向量(er,eθ),但推求εe難以利用推導εp類似的幾何關系,因此利用Molodensky問題嚴格線性化邊值條件,將邊界正常的地球表面視作正常橢球面(中心重合,扁率相同),且等天頂線方向τ與橢球法線方向h近似重合,并顧及,則法線方向偏導數與地心矢徑方向偏導數的關系為:
梯度向量?T在徑向標架(er,eθ,eλ)中可表示為:
由橢球大地測量學可知, sinδ≈e2sinθcosθ,cosδ≈1 ,則有:
由正常位U(r,θ)可導出和,略去推導細節,則有:
式中,J2為二階帶諧系數(動力形狀因子);P2(cosθ)為二階Legendre多項式;ω為地球自轉角速度;實用上可將εh和εγ做成數字格網,形成格網平均值。
4.4.3 利用高分辨率地形數據恢復甚短波擾動重力場
在局部重力場的確定中,目前可利用5′×5′ (相當于10 km×10 km)地面格網平均重力異常和2 160階次的參考重力場EGM08恢復近乎全波段的局部擾動重力場,包括擾動位、大地水準面、高程異常和垂線偏差。這里可將高于5′ ×5′ 分辨率(對應階次n>2 160)的重力場譜分量認為是甚短波重力場,若需獲得更高分辨率的局部重力場參數,理論上要求有更高分辨率的重力數據,可能要大幅提高成本投入;而中等和高山區正是重力數據稀疏、重力測量困難的地區,如何在這些地區在不增加重力測量的條件下,恢復更高分辨率的擾動位,是本文研究的問題。
1)理論依據。本文以地形對地面點引力的直接影響δgt為邊值數據實施廣義Stokes 積分,計算地形對地面點擾動位和垂線偏差的貢獻。地形高采用超高分辨率(如7.5″×7.5″)格網化數據,Stokes核函數則采用相應超高階球諧級數展開。
該方法的基本理論依據為地面甚短波局部擾動重力場完全由局部地形質量產生。地面上一點的擾動重力場參數、重力異常、擾動位和垂線偏差均可在頻譜域分解為長波、中波、短波和甚短波4 個波段譜分量的迭加。地球外部(包括地球表面)擾動重力場是地球自然重力場與正常地球重力場之差,后者是由正常橢球產生的正常重力場。假定正常橢球包含的質量與地球總質量相等、橢球中心與地球質心重合、橢球短軸與地球自轉軸重合、旋轉角速度相同以及正常橢球表面是一個封閉的水準面,其位常數與大地水準面的位常數相等,則可認為擾動重力場由地球體中擾動質量分布產生。擾動質量分布又歸結為地球密度異常分布,是地球密度分布相對于正常橢球某種理論密度分布之差;正常橢球外部無物質分布,因此大地水準面外部的地形質量全部構成一種擾動質量分布或密度異常分布,在空間尺度上,產生中波、短波和甚短波擾動重力場;大地水準面下面至地球深部擾動質量分布,產生長波和中波擾動重力場。因此,可以認為短波和甚短波擾動重力場,特別是甚短波,完全由地形質量產生,與大地水準面以下異常物質分布無關。
目前擾動重力場的長波和中波段可由最新高精度全球擾動位模型以優于1 cm的精度確定,100 km分辨率的短波部分也可達到厘米級精度(2~3 cm),隨著GOCE 衛星重力梯度數據的積累和精化處理,新一代全球擾動位模型提供的短波擾動重力場的精度也可望達到1 cm精度,但垂線偏差不可能達到1″的精度,超出了衛星重力探測敏感地形影響的極限能力,必須在局部重力場確定時利用超高分辨率的地形數據恢復短波和甚短波擾動重力場,使垂線偏差的精度達到或優于1″的水平。
2)由地形數據恢復甚短波地面擾動位。地面Helmert空間異常可表示為:
右邊第一項相當于Stokes 方法中大地水準面上的普通空間異常,用Δg表示,則有:
由Stokes公式,地面擾動位可表示為:
將式(21)分項積分,用兩項積分之和表示,即
其中第一項積分是利用地面重力異常計算地面擾動位的普通Stokes 積分,其分辨率取決于Δg格網數據的分辨率,這里可以是5′ ×5′ 或2.5′ ×2.5′ ;第二項積分中的邊值數據δgt是完全由地形數據計算的地形對地面點引力的直接影響,與重力觀測數據無關,這里通過7.5″ ×7.5″ 地形高格網數據計算。S(r,ψ)的球諧級數展開式為:
式中,r為計算點的地心距;Pn(cosψ)為n階Legendre多項式。
設t=cosψ,則有以下簡單遞推公式:
式(23)中S(r,ψ)級數的最大截斷階nmax應與7.5″ ×7.5″分辨率一致,則nmax=86 400。Pn(cosψ)遞推計算簡單快速高效,雖然總體增加了計算機耗時,但可以很高的精度恢復甚短波地面擾動位。
4.4.4 似大地水準面計算實施情況
南昌市市域似大地水準面模型覆蓋面積近8 000 km2,在全市域似大地水準面計算中,7.5″ ×7.5″ SRTM 用于Airy-Haiskanen 地形均衡歸算、地形位和地形引力影響的計算,采用4 250個點重力數據和80個GNSS水準資料,以EIGEN6C4 地球重力場模型為參考重力場,由第二類Helmert 凝集法計算重力似大地水準面,與80 個GNSS 水準資料的比較精度為±7 mm,均勻分布的22 個未參與大地水準面計算的獨立GNSS 水準點外部檢核精度為±7 mm(表5)[18-20]。利用球冠諧[21-22]調和分析方法將GNSS水準與重力似大地水準面聯合求解似大地水準面,其精度可達±4.5 mm,若采用四舍五入的小數進位法后得出似大地水準面精度為±5 mm;若采用奇進偶不進的小數進位法可得似大地水準面精度為±4 mm,是目前國內精度最高的地級市市域似大地水準面之一。

表5 GNSS水準與GNSS似大地水準面高程殘差統計/m
南昌市市域現代測繪基準體系的建設,實現了市域覆蓋、動靜結合、鄰域關聯,完善了南昌市域現代化、高精度的測繪基準框架,可為社會各界提供三維、動態、高精度的空間定位保障,為數字南昌、國土空間規劃、城鄉建設、自然資源調查和災害監測等各項經濟建設提供廣泛的測繪服務,具有極其重要的科學意義和社會效益。NCCORS 與亞厘米級市域高精度似大地水準面的結合應用,真正實現了GNSS 技術在幾何和物理意義上的三維定位功能,極大地改善了傳統測量作業模式,大幅提高了生產效率,具有巨大的經濟效益。通過研究GNSS 和水準觀測的同期性、GNSS 數據處理以及大地水準面計算等關鍵技術,提升了測繪技術水平,促進了科技進步,可為國內城市建立市域現代測繪基準體系提供借鑒。。