|
李建成2, 孔祥雪3 1. 東北大學(xué)資源與土木工程學(xué)院, 遼寧 沈陽 110819; 收稿日期：2019-04-05；修回日期：2019-09-28 基金項目：國家自然科學(xué)基金(41804002)；中央高?；究蒲袠I(yè)務(wù)專項資金資助項目(170103009)；武漢大學(xué)地球空間環(huán)境與大地測量教育部重點(diǎn)實驗室開放基金(17-01-04) 第一作者簡介：張勝軍(1987-), 男, 博士, 講師, 研究方向為衛(wèi)星測高數(shù)據(jù)處理及應(yīng)用研究。E-mail:zhangshengjun@mail.neu.edu.cn 摘要：海洋重力場模型反演的質(zhì)量主要依賴于采用測高數(shù)據(jù)的精度、空間分辨率和數(shù)據(jù)分布密集程度。本文聯(lián)合Geosat GM/ERM、ERS-1 GM/ERM、TOPEX/Poseidon、Envisat、Cryosat-2、Jason-1 ERM/GM和SARAL/AltiKa等多種測高觀測數(shù)據(jù)集，深入比較了多種波形重跟蹤算法的效果，回波數(shù)據(jù)重跟蹤處理后的沿軌海面高標(biāo)準(zhǔn)差。統(tǒng)計表明，Sandwell算法優(yōu)于MLE-4算法、Davis閾值法、改進(jìn)閾值法和β參數(shù)擬合法；基于不同測高數(shù)據(jù)波形重采樣的結(jié)果給出了沿軌海面梯度計算中低通濾波的參數(shù)選擇方法，并采用Sandwell提出的垂線偏差法，反演了全球海域1'×1'的重力場模型。檢核表明，反演結(jié)果與DTU13和SIO V23.1模型檢核的差值均方根分別為3.4、1.8 mGal，與NGDC船測數(shù)據(jù)的檢核精度為4~8 mGal，且本文模型在部分典型海區(qū)內(nèi)精度更優(yōu)。 關(guān)鍵詞：衛(wèi)星測高 波形重跟蹤 垂線偏差 重力異常 Inversion of global marine gravity anomalies with vertical deflection method deduced from Laplace equation ZHANG Shengjun1 1. School of Resources and Civil Engineering, Northeastern University, Shenyang 110819, China; Foundation support: The National Natural Science Foundation of China (No. 41804002); The Fundamental Research Funds for the Central Universities (No. 170103009); The Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University(No. 17-01-04) First author: ZHANG Shengjun(1987—), male, PhD, lecturer, majors in satellite altimeter data processing and application.E-mail:zhangshengjun@mail.neu.edu.cn. Abstract: The quality of marine gravity field derived from satellite altimetry mainly relies on denser data coverage, enhanced spatial resolution as well as improved range precision. In this study, a comparison of application effect for several waveform retracking methods is firstly executed on the basis of multi-satellite altimetry dataset sincluding Geosat GM/ERM, ERS-1 GM/ERM, TOPEX/Poseidon, Envisat, Cryosat-2, Jason-1 ERM/GM and SARAL/AltiKa. According to the statistics of along-track standard deviation, the two-step retracker proposed by Sandwell is superior to other retrackers such as MLE-4, threshold, advanced threshold, and β-5 parameter fitting method. Tuned parameters for certain low-pass filtersare given within the resampling procedure for each altimeter mission in order to calculate accurate along-track slopes. Then a 1'×1' global marine gravity field model is inversed based on the remove-restore procedure by introducing the EGM2008 as reference model. The root mean squares (rms) of difference are respectively 3.4 mGal and 1.8 mGal by comparing with DTU13 and SIO V23.1 model. Meanwhile, the comparison with ship-measured data results a RMS range of 4~8 mGal, while better coincidence shows up in typical ocean areas. In summary, both the comparisons indicate that the new model has a comparable accuracy level with DTU13 and V23.1 at designed 1'×1' grid interval. Key words: satellite altimetry waveform retracking vertical deflection gravity anomaly 精細(xì)海洋重力場是確定高階全球重力場模型的前提，在海洋大地測量學(xué)、物理海洋學(xué)領(lǐng)域均有廣泛應(yīng)用。衛(wèi)星測高資料是海洋重力場反演的重要數(shù)據(jù)源，海洋重力場模型的持續(xù)改進(jìn)依賴于測高數(shù)據(jù)的積累和數(shù)據(jù)處理水平的提升。國際上長期致力于海洋重力場模型反演的研究機(jī)構(gòu)包括丹麥科技大學(xué)空間研究中心(DTU)和美國加利福尼亞大學(xué)圣迭戈分校斯克利普斯海洋研究所(Scripps Institution of Oceanography，SIO)[1-5]，典型代表如DTU13和SIOV23.1模型等，格網(wǎng)間距均為1′×1′，模型更新的主要工作在于新增衛(wèi)星測高數(shù)據(jù)的融合使用及原有測高數(shù)據(jù)的精細(xì)處理[6-8]。我國研究學(xué)者在海洋重力場模型反演上作了許多工作[9-12]，但模型分辨率和精度水平與國際模型仍有一定差距。 海洋重力場模型的質(zhì)量包括精度和分辨率兩個方面，其精度主要依靠反演方法和相鄰測高數(shù)據(jù)的相對精度，而分辨率取決于采用測高數(shù)據(jù)及反演中相應(yīng)輸入數(shù)據(jù)的覆蓋率。因此，提升海洋重力場模型質(zhì)量的途徑有：①測高數(shù)據(jù)的精細(xì)預(yù)處理提高輸入數(shù)據(jù)精度水平，包括波形重跟蹤、高頻誤差處理等；②低頻率測高數(shù)據(jù)的波形重采樣，提升輸入數(shù)據(jù)的空間分辨率；③引入新的高質(zhì)量測高數(shù)據(jù)，增加輸入數(shù)據(jù)的空間覆蓋。 近年來，測高數(shù)據(jù)波形重跟蹤算法得到長足發(fā)展，且Cryosat-2衛(wèi)星已積累了多個完整的369 d周期觀測數(shù)據(jù)，Jason-1衛(wèi)星在失效前已實施完成一個406.5 d完整周期的漂移軌道任務(wù)，提供了均勻密集覆蓋的觀測數(shù)據(jù)，精度水平相對GeosatGM和ERS-1 GM數(shù)據(jù)有較大提升；SARAL/AltiKa衛(wèi)星采用增加的工作頻率、帶寬、脈沖重復(fù)頻率以及減小的天線寬度，減小足印區(qū)面積并提升了測量精度。這些測高資料不僅提供了高質(zhì)量的觀測值，而且豐富了地面軌跡的分布方式，為更高精度更高分辨率的海洋重力場模型構(gòu)建提供了契機(jī)[13-15]。因此，本文擬聯(lián)合不同任務(wù)階段的多星測高資料，在對比分析多種波形重跟蹤算法改進(jìn)效果的基礎(chǔ)上，采用垂線偏差法構(gòu)建南北緯85°之間全球海域1′×1′的海洋重力場模型。 1 算法原理 文獻(xiàn)[16]指出衛(wèi)星測高資料反演海洋重力場的起算信息分為大地水準(zhǔn)面高和垂線偏差信息兩種，通?；谀鍿tokes公式、逆Vening-Meinesz公式或最小二乘配置法等反演重力異常。利用衛(wèi)星測高數(shù)據(jù)計算垂線偏差的過程，有效抑制了徑向軌道誤差、觀測值校正誤差等長波長項的影響，無須交叉點(diǎn)平差，廣泛應(yīng)用于海洋重力場模型構(gòu)建[17-19]。此外，文獻(xiàn)[4]基于拉普拉斯方程推導(dǎo)的垂線偏差與重力異常之間關(guān)系式，在無須給定復(fù)雜核函數(shù)的情況下可獲取較高的解算精度，該方法目前國內(nèi)學(xué)者應(yīng)用相對較少。本文利用衛(wèi)星測高資料反演海洋重力場將基于該方法展開，下面給出簡化的推導(dǎo)過程。 給定一個地心空間直角坐標(biāo)系(x, y, z)，擾動位T滿足Laplace方程，擾動重力δg可表示為擾動位的徑向梯度。顧及Bruns公式給定的大地水準(zhǔn)面與擾動位之間的關(guān)系，大地水準(zhǔn)面梯度有關(guān)的垂線偏差分量信息ξ和η可分別表示為擾動位的水平梯度。根據(jù)Laplace方程整理出擾動重力與垂線偏差之間的關(guān)系式
引入δg(K, z)、ξ(K)和η(K)分別表示δg(x, y, z)、ξ(x, y)和η(x, y)的二維傅里葉變換，其中K滿足K=(kx, ky), kx=1/λx, ky=1/λy,
基于Laplace方程波數(shù)域內(nèi)的解，可以建立空間任一高程z處的擾動重力與地面擾動重力之間的關(guān)系。沿徑向取導(dǎo)并取高程為零時的值，得
理論上，式(2)與式(3)等號右邊部分等值，聯(lián)立兩式得
引入傅里葉逆變換算子F－1，式(4)變?yōu)?/p>
式(5)即為基于Laplace方程推導(dǎo)出的垂線偏差分量解算擾動重力的傅里葉變換表達(dá)式。根據(jù)重力測量的基本微分方程，文獻(xiàn)[20]給出了重力異常與擾動重力之間的關(guān)系式
式中，擾動位T(x, y)以及大地水準(zhǔn)面N(x, y)信息可采用位模型計算。顧及地球正常重力γ0和地球半徑R的平均值，擾動重力與重力異常轉(zhuǎn)換的改正數(shù)可近似表示為-0.307 6N, N表示單位為米的大地水準(zhǔn)面高。 2 數(shù)據(jù)選取與預(yù)處理 顧及近年新補(bǔ)充的測高資料不僅提供了高質(zhì)量的觀測值，多個漂移軌道任務(wù)及不同的軌道傾角設(shè)計豐富了地面軌跡的分布方式，為海洋重力場模型反演的潛在改進(jìn)提供了契機(jī)。本文首先選取包含波形信息的Jason-1、Cryosat-2和SARAL/AltiKa傳感器地球物理數(shù)據(jù)記錄(SGDR)，其中，Jason-1數(shù)據(jù)包括重復(fù)軌道和大地測量任務(wù)，Cryosat-2數(shù)據(jù)包括低分辨率模式(LRM)、合成孔徑雷達(dá)(SAR)和干涉測量(SIN)3種模式。在此基礎(chǔ)上，聯(lián)合使用SIOV18.1海洋重力場模型反演采用的Geosat、ERS-1、Envisat、TOPEX/Poseidon數(shù)據(jù)集，其中，Geosat和ERS-1數(shù)據(jù)也同時包括了重復(fù)軌道和漂移軌道觀測數(shù)據(jù)，Envisat數(shù)據(jù)包括35 d和30 d兩種重復(fù)周期的觀測數(shù)據(jù)。 在全球海域范圍內(nèi)，存在大量的近岸和島嶼觀測數(shù)據(jù)，受雷達(dá)測高觀測技術(shù)限制，這些區(qū)域內(nèi)測高數(shù)據(jù)本身觀測質(zhì)量不高，目前，較為成熟的方法是利用波形重跟蹤方法重新處理這些區(qū)域內(nèi)的測高觀測數(shù)據(jù)[21-24]。文獻(xiàn)[25]提出了一種顧及回波到達(dá)時間和有效波高之間固有關(guān)聯(lián)性的波形重跟蹤算法，通過沿軌平滑的海面有效波高作為先驗信息和兩次波形擬合估算最終的目標(biāo)跟蹤門位置。該算法能夠最小化海面高梯度的計算誤差，多個算例的研究表明兩步擬合算法針對不同測高任務(wù)的多種回波波形均有效，在近岸和開闊海域取得了顯著改進(jìn)[8, 15, 26]。國內(nèi)眾多學(xué)者也提出了多種不同改進(jìn)理念的重跟蹤算法，均取得了較為顯著的改進(jìn)效果[27-29]。然而由于觀測條件和方案設(shè)計的改進(jìn)原則不同，波形重跟蹤算法通常具有各自的最優(yōu)化應(yīng)用方向。因此，本文的首要工作是對比選取適用于垂線偏差法反演海洋重力場的重跟蹤算法。 本文以Jason-1 GMcycle501的回波數(shù)據(jù)為例，選取中國海及鄰近海域(100°—140°E, 0°—40°N)作為研究區(qū)域，選擇MLE-4算法、Davis閾值法、多波形前緣改進(jìn)閾值法和β參數(shù)擬合法與Sandwell算法進(jìn)行了比較。參照EGM2008模型，沿地面軌跡統(tǒng)計重跟蹤距離改正后的海面高觀測值標(biāo)準(zhǔn)差，同時提取反映海況條件的有效波高信息，對應(yīng)關(guān)系如圖 1所示。
圖 1中紅色點(diǎn)對應(yīng)MLE-4重跟蹤方法改正結(jié)果，綠色和藍(lán)色點(diǎn)分別表示Sandwell算法的兩步擬合結(jié)果，β參數(shù)擬合法、Davis閾值法和改進(jìn)閾值法的統(tǒng)計信息依次由黑色、青色和紫色點(diǎn)表示。圖中平滑曲線表示按照0.5 m區(qū)間的有效波高值對應(yīng)的標(biāo)準(zhǔn)差統(tǒng)計結(jié)果的算術(shù)平均值，數(shù)值越小表明噪聲水平越低，其中Sandwell算法改正的結(jié)果優(yōu)于其他算法。此外，為防止刪除兩步擬合失效的回波波形，Sandwell算法還提供閾值法重跟蹤改正值作為備用。因此，本文統(tǒng)一選擇Sandwell算法重跟蹤處理不同測高任務(wù)的回波，兩步擬合過程中的噪聲改進(jìn)水平如圖 2所示，除了SAR和SIN模式的Cryosat-2數(shù)據(jù)，其他測高數(shù)據(jù)在第2次擬合過程中均顯著降低了噪聲水平。此外，SARAL/AltiKa數(shù)據(jù)的噪聲水平最低，Cryosat-2 SIN數(shù)據(jù)的精度最差。
為提高海洋重力場反演輸入數(shù)據(jù)的分辨率，對重跟蹤處理后的上述所有測高資料，采用等權(quán)平均方式將沿軌高采樣率數(shù)據(jù)統(tǒng)一重采樣為5 Hz。為了抑制高頻噪聲的影響，重采樣處理結(jié)合了低通濾波過程，采用有限長單位沖激響應(yīng)(FIR)濾波器，設(shè)計準(zhǔn)則為半增益值(0.5 Gain)對應(yīng)波長為6.7 km。根據(jù)不同測高任務(wù)觀測值的時空間隔，對應(yīng)設(shè)計的重采樣過程所用的FIR濾波器如圖 3所示，10、18、20和40 Hz的沿軌數(shù)據(jù)采用的濾波器長度依次為49、99、99和199。
根據(jù)統(tǒng)一重采樣為沿軌5 Hz的多星測高資料[4, 30]，考慮各種距離延遲對應(yīng)的校正項，由相鄰海表面高度觀測值結(jié)合觀測時間信息計算海面梯度，同時根據(jù)EGM2008模型提供的沿軌大地水準(zhǔn)面高計算相應(yīng)的沿軌大地水準(zhǔn)面梯度，根據(jù)兩者差值的直方圖統(tǒng)計結(jié)果和3倍標(biāo)準(zhǔn)差準(zhǔn)則選定用于剔除粗差數(shù)據(jù)的閾值，結(jié)果見表 1?？紤]到沿軌海面梯度計算的差分過程對高頻噪聲的放大影響，對沿軌海面梯度值再次進(jìn)行了低通濾波。根據(jù)10 km波長信號對應(yīng)增益值為0.5的設(shè)計原則，不同測高衛(wèi)星采用的Parks-McClellan濾波器如圖 4所示，濾波器長度均為30。本文選取參與構(gòu)建全球海洋重力異常的多星測高資料的使用與預(yù)處理概況見表 2。 表 1 不同衛(wèi)星測高資料沿軌海面梯度粗差剔除的閾值Tab. 1 Threshold for editing gross errors for multi-satellite altimeter along-track slopes m/s
表選項
表 2 構(gòu)建全球海洋重力場模型時不同測高資料的使用與預(yù)處理概況Tab. 2 Data records of multi-satellite altimeter material for constructing global marine gravity field model
表選項 從表 2中可以看出，測高資料的數(shù)據(jù)量非常龐大，篩選粗差數(shù)據(jù)能夠剔除較大比例的觀測值，高緯度海域(南北緯60°以上)內(nèi)的Envisat、Cryosat-2/SAR和Cryosat-2/SIN數(shù)據(jù)尤為明顯。此外，Envisat/ERM、Jason-1/ERM、SARAL/ AltiKa，以及軌道調(diào)整后的T/P數(shù)據(jù)未經(jīng)共線處理而直接計算沿軌殘余垂線偏差，重復(fù)周期的沿軌計算結(jié)果在格網(wǎng)點(diǎn)最小二乘解算的過程中視作相互獨(dú)立的輸入數(shù)據(jù)。 3 模型反演與檢核 基于多星測高海表面高度觀測值，扣除DOT2008A穩(wěn)態(tài)海面地形和EGM2008大地水準(zhǔn)面高計算殘余大地水準(zhǔn)面高，根據(jù)相鄰數(shù)值差分和觀測時間信息計算沿軌殘余大地水準(zhǔn)面梯度，結(jié)合測高衛(wèi)星地面軌跡處的速度近似公式計算沿軌殘余垂線偏差。在每個格網(wǎng)點(diǎn)處建立沿軌殘余垂線偏差和殘余垂線偏差方向分量的關(guān)聯(lián)，根據(jù)最小二乘準(zhǔn)則計算殘余垂線偏差北向分量和東向分量，進(jìn)而推求殘余重力異常值，如圖 5所示。最后恢復(fù)參考模型EGM2008計算的重力異常，得到全球大部分海域1′×1′的重力場模型(圖 6)。
為了評價反演模型的精度，首先采用DTU13和SIO V23.1海洋重力場模型進(jìn)行了檢核，差值統(tǒng)計信息見表 3?？梢钥闯?，反演模型與不同檢核模型之間的數(shù)值差異不大，均方根分別為1.8、3.4 mGal。其中，反演模型與V23.1模型符合度更好，原因在于兩者均采用基于垂線偏差的反演方法。 表 3 反演結(jié)果與已發(fā)布模型的檢核信息Tab. 3 Validation information of inversed result with published gravity models mGal
表選項 此外，還采用NGDC提供的船測數(shù)據(jù)進(jìn)行了外部檢核，檢核數(shù)據(jù)集包含683條測線，具體分布如圖 7所示。通過與EGM2008模型重力異常值對比逐測線進(jìn)行系統(tǒng)偏差改正，根據(jù)穩(wěn)健性異常值探測算法(RODA)剔除了約3.17%的粗差觀測值，預(yù)處理后的NGDC船測數(shù)據(jù)集按照太平洋、大西洋、西北大西洋、印度洋和北冰洋大致區(qū)域進(jìn)行劃分，最終統(tǒng)計的檢核結(jié)果見表 4?？梢钥闯?，本文反演模型與DTU13和V23.1模型的整體精度相當(dāng)，已接近或達(dá)到目前國際上基于多星測高資料反演海洋重力場的先進(jìn)水平。
表 4 反演結(jié)果與NGDC船測數(shù)據(jù)的檢核信息Tab. 4 Validation information of inversed result with NGDC shipborne gravity data mGal
表選項
本文聯(lián)合采用多種衛(wèi)星測高觀測任務(wù)高采樣率的波形資料，基于兩步擬合波形重跟蹤算法提升相鄰測高觀測值的相對精度，顧及觀測精度和空間覆蓋密集度兩方面因素沿軌重采樣為5 Hz觀測數(shù)據(jù)集，并在此基礎(chǔ)上基于Laplace方程導(dǎo)出的垂線偏差法，反演了南北緯85°之間全球海域1′×1′格網(wǎng)分辨率的海洋重力場模型。 以Jason-1 GM數(shù)據(jù)為例，比較了MLE-4算法、Davis閾值法、多前緣改進(jìn)閾值法、β-5參數(shù)擬合法和兩步擬合波形重跟蹤算法等方法的改進(jìn)效果。結(jié)果表明Sandwell提出的兩步擬合波形重跟蹤算法改正精度更優(yōu)。通過調(diào)整波形后緣參數(shù)，成功擬合了不同測高任務(wù)的回波波形，除了Cryosat-2 SAR和SIN模式的數(shù)據(jù)兩步擬合未現(xiàn)改進(jìn)效果，其余測高任務(wù)在采用有效波高值作為先驗信息的第2步波形擬合過程均顯著提升約1.5倍的海面高精度。 針對研究區(qū)域內(nèi)反演得到的1′×1′海洋重力場模型，采用已權(quán)威發(fā)布的海洋重力場模型和NGDC船測數(shù)據(jù)分別進(jìn)行了檢核。結(jié)果表明，本文反演模型與DTU13和V23.1模型的精度相當(dāng)，已達(dá)到目前基于測高資料反演海洋重力場的國際先進(jìn)水平。在四大洋海域分別選取的船測數(shù)據(jù)結(jié)果表明，本文反演模型與船測數(shù)據(jù)符合更好，表明在融合Jason-1 GM、Cryosat-2和SARAL/AltiKa Ka波段衛(wèi)星測高資料，有效提升了海洋重力場模型構(gòu)建的水平。 致謝: 衷心感謝AVISO提供的測高數(shù)據(jù)，感謝DTU提供的重力場模型，感謝NGDC提供的船測數(shù)據(jù)，感謝Sandwell提供的測高數(shù)據(jù)和V23.1模型。 【引文格式】張勝軍, 李建成, 孔祥雪. 基于Laplace方程的垂線偏差法反演全球海域重力異常. 測繪學(xué)報，2020，49(4)：452-460. DOI: 10.11947/j.AGCS.2020.20190108 |
|
|
, LI Jiancheng2, KONG Xiangxue3 
(1)
和ky表征波數(shù)域內(nèi)x軸和y軸方向的分量，λx與λy表示對應(yīng)坐標(biāo)軸方向上的波長。根據(jù)導(dǎo)數(shù)傅里葉變換定理，將微分方程式轉(zhuǎn)換為復(fù)數(shù)域代數(shù)方程式
(2)
(3)
(4)
(5)
(6)
