聶建亮 張雪萍 郭鑫偉 王莉莉 趙文普

1 陜西測繪地理信息局，西安市友誼東路334號，710054 2 自然資源部大地測量數(shù)據(jù)處理中心，西安市友誼東路334號，710054 3 自然資源部第二地形測量隊，西安市測繪路6號，710054

為滿足社會發(fā)展需要，我國建立了新一代中國似大地水準(zhǔn)面模型CQG2000。21世紀(jì)初，陳俊勇等[1-3]和晁定波[4]基于地面重力、衛(wèi)星測高、船測重力等數(shù)據(jù)，采用EGM96全球重力場，建立重力似大地水準(zhǔn)面模型。利用高程異?？刂凭W(wǎng)對其進(jìn)行擬合糾正，獲取15′×15′分辨率的似大地水準(zhǔn)面模型，其中102°E以東地區(qū)中誤差小于±0.3 m，102°E以西、36°N以北地區(qū)中誤差小于±0.4 m，102°E以西、36°N以南地區(qū)中誤差小于±0.6 m。近年來，各省市陸續(xù)建立了高精度、高分辨率的區(qū)域似大地水準(zhǔn)面模型[5-11]。雖然高精度似大地水準(zhǔn)面模型能夠滿足精度要求，但高程轉(zhuǎn)換計算量巨大。以1∶10 000圖幅為例，單個圖幅的高程基準(zhǔn)轉(zhuǎn)換通常需要幾十萬到上百萬個要素點，然而圖幅高程基準(zhǔn)轉(zhuǎn)換所需的高程精度卻較低，在遙感行業(yè)，精度一般為dm級甚至m級[12-13]。區(qū)域似大地水準(zhǔn)面在地球表面空間域內(nèi)變化較緩，但在地形劇烈變化區(qū)域或地類變化過渡區(qū)域的變化較為明顯。梯度作為數(shù)值變化的表征量，可以反映數(shù)值在一定區(qū)域內(nèi)變化的大小、方向。CQG2000的梯度用格網(wǎng)點在經(jīng)、緯度方向上的高程異?？臻g變化率表示，該數(shù)值可以量化不同比例尺圖幅內(nèi)高程異常的變化幅度，可用于分析圖幅內(nèi)平均高程異常對不同比例尺高程精度的適用性。此外，CQG2000的梯度也能反映我國大陸范圍內(nèi)重力場的空間變化特性，有助于解譯地球內(nèi)部質(zhì)量分布等相關(guān)問題。

基于此，本文利用梯度定量描繪CQG2000在全國范圍空間域內(nèi)的變化趨勢，分析CQG2000在空間域中的相關(guān)性，為建立實用化區(qū)域似大地水準(zhǔn)面模型提供基礎(chǔ)數(shù)據(jù)，同時也為深入解釋似大地水準(zhǔn)面變化機(jī)理提供科學(xué)依據(jù)。

1 梯度模型及計算方法

1.1 梯度模型

區(qū)域似大地水準(zhǔn)面受地球質(zhì)量分布的影響，在平原、丘陵及高山區(qū)的變化趨勢不同，變化劇烈區(qū)域的結(jié)果可能未達(dá)到行業(yè)應(yīng)用的精度標(biāo)準(zhǔn)。為更好地表述區(qū)域似大地水準(zhǔn)面格網(wǎng)模型變化，引入梯度函數(shù)，即格網(wǎng)點(B,L)處的梯度：

(1)

式中，f(B,L)為區(qū)域似大地水準(zhǔn)面模型，B、L為緯度和經(jīng)度，ΔB、ΔL為緯度、經(jīng)度方向上的分辨率，m、n為格網(wǎng)模型的行、列數(shù)，i、j為梯度矢量的分量方向。

對比分析多種方法發(fā)現(xiàn)，三階反距離平方權(quán)差分方法計算效果較好[14]。因此，本文采用坡度的三階反距離權(quán)差分方法計算CQG2000格網(wǎng)點的梯度數(shù)值。

1.2 三階反距離平方權(quán)差分方法

采用三階反距離平方權(quán)差分方法選取3×3的移動窗口，計算似大地水準(zhǔn)面規(guī)則格網(wǎng)點的梯度(圖1)，即

(2)

(3)

式中，fB、fL為緯度、經(jīng)度方向上的梯度分量，ΔB、ΔL為緯度、經(jīng)度方向上的分辨率，ζ1～ζ9分別為9個格網(wǎng)點的高程異常值。

圖1 梯度計算示意圖Fig.1 Schematic diagram of gradient calculation

梯度方位角α為：

(4)

1.3 梯度精度推導(dǎo)

假設(shè)在3×3的窗口內(nèi)，區(qū)域似大地水準(zhǔn)面可以表示為曲面ζk=f(Bk,Lk)，且三階連續(xù)可微，按泰勒級數(shù)展開并取至三次項，即

f(B+ΔB,L)=f(B,L)+fB(B,L)ΔB+

(5)

f(B-ΔB,L)=f(B,L)-fB(B,L)ΔB+

(6)

在8個格網(wǎng)點處按泰勒級數(shù)展開，整理后的fB為：

(f?B(ζB,L-ΔL)-f?B(γB,L-ΔL))+

(f?B(ζB,L+ΔL)-f?B(γB,L+ΔL))]}

(7)

令MB、ML為f?B、f?L關(guān)于B、L的三階導(dǎo)數(shù)最大范圍，其誤差具有隨機(jī)性，則有：

(8)

式中，fB一般取式(8)第1項，第2、3項為數(shù)據(jù)誤差及其他誤差。

同理可得fL為：

(9)

假設(shè)區(qū)域似大地水準(zhǔn)面格網(wǎng)點中誤差為σζ，MB、ML中誤差相等且均為M，則σfB為：

(10)

同理可得σfL為：

(11)

由于M數(shù)值較小，因此式(10)、(11)中的第2項可以忽略。

根據(jù)誤差傳播定律，由式(4)、(10)、(11)推導(dǎo)出梯度方位角誤差σα為：

(12)

2 CQG2000的梯度分析

為更加精細(xì)地描繪CQG2000的空間域變化，在計算梯度時，對原有15′×15′分辨率的模型進(jìn)行雙線性內(nèi)插處理成5′×5′分辨率的格網(wǎng)模型(忽略內(nèi)插誤差影響)。按照§1.3中公式計算格網(wǎng)點梯度數(shù)值及方向，獲取全國似大地水準(zhǔn)面的梯度，如圖2所示。同時，根據(jù)表1、2中的全國梯度信息，按照CQG2000精度劃分的3個區(qū)域范圍統(tǒng)計102°E以東(區(qū)域1)，102°E以西、36°N以北(區(qū)域2)和102°E以西、36 °N以南(區(qū)域3)3個區(qū)域的梯度分量精度和梯度方位角平均精度，結(jié)果見表3、圖2。

表1 CQG2000的梯度數(shù)值

表2 CQG2000的梯度方位角

表3 CQG2000的梯度精度

圖2 CQG2000的梯度分布Fig.2 The distribution of gradient for CQG2000

綜合分析表1～3和圖2可知：

1)102°E以東地區(qū)CQG2000的梯度方位角方向主要為東，量級小于0.2 m/(′)，且變化平緩。其中，東北地區(qū)、華南地區(qū)梯度方向主要為東南，角度約為東向南45°；30°～40°N區(qū)間內(nèi)的梯度方位角方向為正東。102°E以西地區(qū)CQG2000梯度在不同地類區(qū)域的方向不同，且梯度數(shù)值變化幅度較大。

2)我國東部、西北、西南地區(qū)的梯度分量精度σfB、σfL隨似大地水準(zhǔn)面格網(wǎng)點精度σζ的增大而增大，隨分辨率的增大而減小。西南地區(qū)σα誤差最大，為50.1°；西北地區(qū)σα誤差最小，為16.4°。造成上述差異的主要原因為梯度方位角精度σα主要受梯度分量及其精度的影響，較小的梯度數(shù)值變化會導(dǎo)致較大的梯度方位角精度變化。進(jìn)一步說明，相比于梯度分量精度σfB、σfL，梯度方位角誤差σα對梯度方位角α更加敏感。

3)新疆北部與南部、青海、成都平原等地存在多個梯度擴(kuò)散中心，區(qū)域似大地水準(zhǔn)面格網(wǎng)點高程異常值由中心向四周逐漸增大；喜馬拉雅山脈、新疆中部地區(qū)出現(xiàn)區(qū)域似大地水準(zhǔn)面梯度會聚現(xiàn)象；河西走廊也出現(xiàn)區(qū)域似大地水準(zhǔn)面梯度會聚現(xiàn)象，甘南西部梯度方向為正南，同時順時針旋轉(zhuǎn)；寧夏南部、甘肅東北部地區(qū)區(qū)域似大地水準(zhǔn)面梯度方向為正東，同時逆時針旋轉(zhuǎn)。喜馬拉雅山脈、新疆南部昆侖山脈等地的區(qū)域似大地水準(zhǔn)面梯度數(shù)值較大，最大值位于西藏墨脫縣(0.625 m/(′))；青藏高原絕大部分區(qū)域的區(qū)域似大地水準(zhǔn)面梯度較小，變化較平緩。

4)CQG2000梯度變化劇烈區(qū)域與我國大山脈走向一致，說明區(qū)域似大地水準(zhǔn)面與地球內(nèi)部質(zhì)量分布密切相關(guān)。在喜馬拉雅山、昆侖山、天山、秦嶺等山脈區(qū)域，區(qū)域似大地水準(zhǔn)面變化速度較大，梯度數(shù)值變化為小-大-小。

5)CQG2000自身精度為dm級，實際分辨率為15′×15′。雖然區(qū)域似大地水準(zhǔn)面梯度精度不高，但能整體反映我國范圍內(nèi)似大地水準(zhǔn)面空間域的變化特征。此外，CQG2000梯度與全國范圍內(nèi)的垂線偏差分布[15]一致，僅在西藏東南部地區(qū)存在一定差異，可能是該區(qū)域缺乏重力等基礎(chǔ)資料、似大地水準(zhǔn)面精度較低所致。

6)由于CQG2000是利用EGM96參考重力場模型獲取高精度重力似大地水準(zhǔn)面，再通過GNSS水準(zhǔn)點糾正融合得到，因此重力似大地水準(zhǔn)面梯度與似大地水準(zhǔn)面梯度的總體趨勢一致，二者系統(tǒng)偏差較小，梯度大小與方向差異也較小。CQG2000梯度與高精度EGM2008重力場模型梯度在我國中東地區(qū)的差異較小，在喜馬拉雅山、昆侖山、天山、秦嶺等山脈區(qū)域存在一定差異，最大差值位于喜馬拉雅山地區(qū)，約為0.18 m/(′)，可能是CQG2000在該區(qū)域缺乏基礎(chǔ)資料所致。

3 結(jié) 語

梯度信息可反映我國似大地水準(zhǔn)面CQG2000在空間域內(nèi)的整體變化情況，描述不同范圍內(nèi)高程異常變化的量級和方向，可為建立行業(yè)所需的實用化模型提供基礎(chǔ)數(shù)據(jù)。整體上看，東北地區(qū)CQG2000梯度方位角方向為正東且運動規(guī)律與地殼水平運動趨勢一致；西部地區(qū)梯度方位角方向隨地形分布逐漸變化，梯度劇烈變化區(qū)域與昆侖山等大山脈走向一致。本文可為進(jìn)一步開展地球動力學(xué)機(jī)理研究提供科學(xué)依據(jù)。此外，高分七號等遙感衛(wèi)星能夠開展大范圍1∶10 000立體測圖工作，滿足各行業(yè)對高精度產(chǎn)品的需求，但如何綜合利用區(qū)域似大地水準(zhǔn)面及梯度信息快速、高效地實現(xiàn)衛(wèi)星遙感影像數(shù)據(jù)處理的高程自動化轉(zhuǎn)換，是下一步的研究重點。