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

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

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

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

1 梯度模型及計(jì)算方法

1.1 梯度模型

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

(1)

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

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

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

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

(2)

(3)

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

圖1 梯度計(jì)算示意圖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ù)可微，按泰勒級(jí)數(shù)展開(kāi)并取至三次項(xiàng)，即

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個(gè)格網(wǎng)點(diǎn)處按泰勒級(jí)數(shù)展開(kāi)，整理后的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項(xiàng)，第2、3項(xiàng)為數(shù)據(jù)誤差及其他誤差。

同理可得fL為：

(9)

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

(10)

同理可得σfL為：

(11)

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

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

(12)

2 CQG2000的梯度分析

為更加精細(xì)地描繪CQG2000的空間域變化，在計(jì)算梯度時(shí)，對(duì)原有15′×15′分辨率的模型進(jìn)行雙線性?xún)?nèi)插處理成5′×5′分辨率的格網(wǎng)模型(忽略?xún)?nèi)插誤差影響)。按照§1.3中公式計(jì)算格網(wǎng)點(diǎn)梯度數(shù)值及方向，獲取全國(guó)似大地水準(zhǔn)面的梯度，如圖2所示。同時(shí)，根據(jù)表1、2中的全國(guó)梯度信息，按照CQG2000精度劃分的3個(gè)區(qū)域范圍統(tǒng)計(jì)102°E以東(區(qū)域1)，102°E以西、36°N以北(區(qū)域2)和102°E以西、36 °N以南(區(qū)域3)3個(gè)區(qū)域的梯度分量精度和梯度方位角平均精度，結(jié)果見(jiàn)表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的梯度方位角方向主要為東，量級(jí)小于0.2 m/(′)，且變化平緩。其中，東北地區(qū)、華南地區(qū)梯度方向主要為東南，角度約為東向南45°；30°～40°N區(qū)間內(nèi)的梯度方位角方向?yàn)檎龞|。102°E以西地區(qū)CQG2000梯度在不同地類(lèi)區(qū)域的方向不同，且梯度數(shù)值變化幅度較大。

2)我國(guó)東部、西北、西南地區(qū)的梯度分量精度σfB、σfL隨似大地水準(zhǔn)面格網(wǎng)點(diǎn)精度σζ的增大而增大，隨分辨率的增大而減小。西南地區(qū)σα誤差最大，為50.1°；西北地區(qū)σα誤差最小，為16.4°。造成上述差異的主要原因?yàn)樘荻确轿唤蔷圈姚林饕芴荻确至考捌渚鹊挠绊?，較小的梯度數(shù)值變化會(huì)導(dǎo)致較大的梯度方位角精度變化。進(jìn)一步說(shuō)明，相比于梯度分量精度σfB、σfL，梯度方位角誤差σα對(duì)梯度方位角α更加敏感。

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

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

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

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

3 結(jié) 語(yǔ)

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