全球/區(qū)域地形和重力異常間相干性分析

范 雕，李姍姍，趙東明，楊軍軍，孟書宇，范昊鵬

（1.信息工程大學(xué)，鄭州 450001；2.中國自然資源航空物探遙感中心，北京 100083；3.32022 部隊，廣州 510630）

海洋觀測是研究海洋、開發(fā)海洋和利用海洋的基本手段，是加強海洋國防建設(shè)、預(yù)警海洋災(zāi)害、開發(fā)海洋資源、發(fā)展海洋經(jīng)濟、謀求新發(fā)展空間等的重要支撐。海底地形作為海洋觀測的重要組成內(nèi)容，其探測技術(shù)歷來是各海洋國家相關(guān)組織機構(gòu)的研究熱點。目前海底地形測量技術(shù)主要包括船基海底地形測量技術(shù)[1]、潛基海底地形測量技術(shù)[2]、機載激光雷達(dá)測深（Airborne Lidar bathymetry，ALB）技術(shù)[3]、星基海底地形測量技術(shù)[4,5]等。船基海底地形測量技術(shù)測量（使用單波束、多波束設(shè)備測量）精度高，但就全球海域而言，存在測量效率低下且測量結(jié)果分布極不均勻等缺點。ALB 技術(shù)測量效率高，主要用于測量海區(qū)水質(zhì)較好的淺海環(huán)境，應(yīng)對深海海域或者嚴(yán)重影響光束傳播海域地形測繪能力有限。星基海底地形測量技術(shù)中基于衛(wèi)星測高技術(shù)構(gòu)建海底地形是目前主要依賴的快速構(gòu)建全球海底地形方法手段，可有效滿足海底地形全球性、連續(xù)性、智能化和動態(tài)化觀測要求。

衛(wèi)星測高技術(shù)構(gòu)建海底地形基本原理是利用覆蓋海面的衛(wèi)星測高重力數(shù)據(jù)，采用相應(yīng)的重力反演海底地形方法間接恢復(fù)海域海底地形數(shù)據(jù)[6-9]。海面重力異常主要受海水質(zhì)量虧損、地殼密度分布不均勻、地殼以下質(zhì)量的均衡補償?shù)纫蛩赜绊慬10]，從而利用衛(wèi)星測高重力數(shù)據(jù)僅可恢復(fù)有限頻段海底地形信息[11]。據(jù)此可知，確定重力反演海底地形有效截止頻段范圍是依據(jù)衛(wèi)星測高技術(shù)構(gòu)建海底地形的方法前提[12,13]?；诖耍疚囊环矫嬉訣GM2008、EIGEN-6C4 全球重力場位系數(shù)模型和DTM2006、EARTH2012 全球地形球諧系數(shù)模型為數(shù)據(jù)輸入，使用全球地形-重力相干性處理方法，解算了全球范圍內(nèi)地形和重力異常頻譜相干性特征，給出了全球范圍地形-重力表現(xiàn)強相干性的頻段參考范圍；另一方面，以實際海域局部海底地形和海面重力異常為真實數(shù)據(jù)參考，使用區(qū)域地形-重力相干性計算方法，獲取了真實海域地形-重力表現(xiàn)強相干性的頻段參考。研究結(jié)果以期為利用衛(wèi)星測高重力數(shù)據(jù)恢復(fù)海底地形信息截止頻段的確定提供參考與借鑒。

1 理論與方法

相干性（Coherency）也稱為幅值平方相干（Magnitude Squared Coherence,MSC），表征輸入信號在頻域空間的線性相關(guān)程度。假設(shè)時域輸入信號分別為x(t) 和y(t)，那么信號頻域相干性可表示為信號互功率譜密度（Cross-Power Spectral Density，cross-PSD）的正?；Y(jié)果：

式中，r2(f)表示相干性；PSD(x,y)表示信號x(t)和y(t)的互功率譜密度；PSD(x)和PSD(y)分別表示信號x(t) 和y(t) 的自功率譜密度（Auto-Power Spectral Density,auto-PSD）。

1.1 球地形與重力場元間相干性

由重力場知識可知，重力異?？杀硎緸椋?/p>

式中，PSD表示功率譜密度；其余參數(shù)意義同上。

利用地形球諧系數(shù)模型求解地形h（高度起算面為平均海平面）公式如下[15,16]：

利用式(2)和式(4)可得地形和重力異常的互功率譜密度：

式中參數(shù)意義同上。

依據(jù)相干性分析可得重力異常和地形的相干性為：

依據(jù)文獻(xiàn)[10]可知空間分辨率（ψ為兩點的地心角距）與模型階數(shù)（n）的關(guān)系為：

1.2 區(qū)域地形與重力場元間相干性

設(shè)海面某類重力場元符號為Γ（重力異常、重力異常垂直梯度等），海底地形為b，則依據(jù)導(dǎo)納函數(shù)可得重力場元Γ 與海底地形b的頻率域表達(dá)式為[17]：

式中，F(xiàn)(Γ)和F(b)分別表示重力場元Γ 與海底地形b的傅里葉變換；ZΓ為重力場元Γ 對應(yīng)的導(dǎo)納函數(shù)，如當(dāng)Γ 表示重力異常，考慮撓曲均衡補償時，導(dǎo)納函數(shù)ZΓ函數(shù)形式為：

式中，G為地球引力常數(shù)；ρc和ρw分別表示地殼密度和海水密度；d為海域平均海深（取負(fù)）；Φe(f)為撓曲響應(yīng)函數(shù)；Tc為地殼厚度；f表示徑向頻率。同理，當(dāng)不考慮均衡補償、考慮Airy 均衡補償以及顧及多層地殼結(jié)構(gòu)時，導(dǎo)納函數(shù)ZΓ均存在類似的表達(dá)形式，具體參見文獻(xiàn)[18]。對式(9)中導(dǎo)納函數(shù)施以傅里葉逆變換操作，同時依據(jù)卷積定理可得：

式中，F(xiàn)-1() 表示傅里葉逆變換；“*”表示卷積運算。從而重力場元Γ 可表示為：

式(12)為依據(jù)海底地形恢復(fù)海面重力場信息的表達(dá)式。實際測量中，重力場元Γ 將不可避免存在噪聲影響，設(shè)噪聲與重力場元Γ 和海底地形b相互獨立，從而式(12)可表示為：

式中，ω表示測量噪聲。依據(jù)Wiener-Khintchine 定理和相關(guān)定理，重力場元Γ 的自功率譜密度可表示為：

式中，R(Γ,Γ) 表示重力場元Γ 的自相關(guān)函數(shù)；F*(x)表示F(x) 的共軛。將式(13)帶入式(14)，據(jù)傅里葉變換線性性質(zhì)、共軛性質(zhì)、時域卷積定理、相關(guān)定理和Wiener-Khintchine 定理，同時顧及噪聲與重力場元和海底地形相互獨立條件可得：

同理可得，海底地形的自功率譜密度為：

重力場元與海底地形互功率譜密度為：

從而，根據(jù)信號間相干性關(guān)系可得重力場元和海底地形相干性為：

令重力場元信號和測量噪聲的信噪比表達(dá)式為：

從而依據(jù)式(18)可知，1）相干性等于0.5 可以理解為：海底地形信號b不含噪聲的情況下，重力場元Γ 的信噪比為1：1；2）當(dāng)重力場元信號Γ 和海底地形信號b均不含噪聲情況下，二者的頻域相干性恒為1。也就是說，當(dāng)重力場元信號與海底地形信號相干性在0～1之間時，原因可能有：

（1）重力場元噪聲影響。若重力場元與海底地形服從線性關(guān)系，理論上二者頻域相干性數(shù)值大小應(yīng)為1，然而由于重力場元信號存在噪聲導(dǎo)致二者頻域線性相干性不為1；如圖1 顯示了重力場元在不同信噪比環(huán)境下信號間相干性變化情況。通過圖1 可以看出，隨著重力場元信號信噪比不斷增大，相干性也隨著不斷增大；即噪聲的影響導(dǎo)致重力場元與海底地形線性相關(guān)程度降低。

圖1 相干性隨信噪比變化曲線Fig.1 Coherency curve with SNR

（2）重力場元信號Γ 與海底地形信號b實際線性相關(guān)程度較低。①重力場元信號和海底地形信號相干性表示二者的線性相關(guān)程度，實際上基于Parker 公式[19]推導(dǎo)的重力導(dǎo)納函數(shù)僅僅是Parker 公式的線性近似，忽略了高次項海底地形對海面重力場元影響，二者實際并非完全線性相關(guān)，從而相干性不等于1；②海面重力場元并非完全受海底地形影響，海底洋殼的均衡狀態(tài)、海底地形的密度差異以及海水結(jié)構(gòu)等因素均會不同程度地改變海面重力場元信息，從而海面重力場元與海底地形相干性表現(xiàn)為0～1 間變化。

重力場元信號和海底地形信號均為功率有限信號時，依據(jù)Wiener-Khintchine 定理和相關(guān)定理，信號間相干性可表示為：

根據(jù)式(20)可知，若輸入信號為功率有限信號，則按照式(20)計算信號間相干性，不論信號是否相關(guān)以及線性相關(guān)強弱，二者頻率域相干性恒為1。針對這一問題，通常應(yīng)對策略為將式(20)解算的功率譜密度進(jìn)行平滑處理：式中，<>表示方位角α范圍內(nèi)的平均值。方位角的取平均值實際操作流程如下[20-22]：

（1）坐標(biāo)轉(zhuǎn)換。將頻率笛卡爾坐標(biāo)系（fx、fy）轉(zhuǎn)換為極坐標(biāo)系（f，α），兩坐標(biāo)系間坐標(biāo)關(guān)系為：

式中，fx和fy分別表示x和y方向上的頻率；其余參數(shù)意義同上。

（2）方位角范圍取平均。徑向頻率取值范圍為[0 min（f x,fy）]，以dq為徑向頻率間隔，分別計算環(huán)徑為dq范圍內(nèi)的各信號功率譜密度平均值，然后依照式(16)解算信號相干性。

2 數(shù)值分析試驗

2.1 算例1：全球地形和重力場元相干性

當(dāng)前國際廣泛使用的超高階（大于1000 階）重力場模型有EGM2008、GECO和EIGEN系列等模型等。本文以EGM2008（2160 階）和EIGEN-6C4（2190 階）地球重力場模型作為計算重力異常譜的重力場輸入，以EARTH2012（2160 階）和DTM2006（2190 階）地形球諧系數(shù)模型作為解算地形譜的數(shù)據(jù)輸入，依據(jù)式(3)和式(5)計算的重力異常功率譜密度和地形功率譜密度結(jié)果如圖2 和圖3 所示，其中紅色曲線分別表示EGM2008 和EARTH2012 位系數(shù)模型計算結(jié)果，藍(lán)色曲線分別表示EIGEN-6C4 和DTM2006 位系數(shù)模型解算結(jié)果。可以看出，重力異常和地形的功率譜密度均隨著模型階數(shù)的增加而減小，說明模型中長波信號占主要部分，高頻信號能量較弱；而EIGEN-6C4 模型功率譜密度在2160 階后出現(xiàn)急劇變化，這與模型構(gòu)建機制有關(guān)，如使用的數(shù)據(jù)分辨率等因素。根據(jù)式(6)計算EGM2008 和EARTH2012 之間的互功率譜密度如圖4 中紅色曲線所示，EIGEN-6C4 和DTM2006 之間的互功率譜密度曲線如圖4 中藍(lán)色曲線所示。利用式(1)獲取的分別以EGM2008 和EARTH2012 模型為輸入，以及以EIGEN-6C4 和DTM2006 模型為輸入的相干性結(jié)果如圖5 中紅色和藍(lán)色曲線所示。

圖2 重力異常自功率譜密度Fig.2 The auto-PSD of gravity anomaly

圖3 地形自功率譜密度Fig.3 The auto-PSD of topography

圖4 重力異常和地形互功率譜密度Fig.4 The cross-PSD of gravity anomaly and topography

圖5 重力異常和地形相干性Fig.5 coherence between gravity and topography

圖5 重力異常和地形的相干性結(jié)果表明，隨著階次的增加，二者相干性呈現(xiàn)先升高后降低的特征。由相干性定義可知，相干性表示一種信號可線性轉(zhuǎn)換成另一種信號的程度：當(dāng)相干性為1 時即兩種信號完全線性相關(guān)；相干性為0 表示兩信號沒有線性關(guān)系[23]。當(dāng)相干性為0.5 可理解為信噪比（signal-to-noise，SNR）為1:1 的信號與另一不含噪聲信號的相干性結(jié)果[24]。從圖5 可以看出，全球范圍內(nèi)重力異常和地形的相干性在0.5 以上對應(yīng)的模型階數(shù)大約為120階～1000階，對應(yīng)的波長范圍為25.45 km～210.61 km，該計算結(jié)果與目前國際國內(nèi)學(xué)者利用重力數(shù)據(jù)恢復(fù)海底地形的波段范圍一致[6,7,17,25-28]。

2.2 算例2：區(qū)域地形和重力場元相干性

選擇某3 °× 3°（12°N～15° N，112 °E～115°E）海域作為研究海區(qū)，遵照以上理論分析研究海區(qū)海底地形和重力異常頻域相干性大小。研究海區(qū)水深數(shù)據(jù)（格網(wǎng)化水深）和衛(wèi)星測高重力異常數(shù)據(jù)（S&SV24.1）空間分辨率為1 角分，空間分布情況如圖6 所示。相應(yīng)數(shù)據(jù)統(tǒng)計結(jié)果見表1，數(shù)據(jù)類別括號內(nèi)為數(shù)據(jù)單位。

表1 研究區(qū)域數(shù)據(jù)概況Tab.1 Information about data in the study area

圖6 研究區(qū)域數(shù)據(jù)Fig.6 Data in the study area

經(jīng)計算，研究區(qū)域東西向總長為324.37 km，南北向總長為333.58 km，其中南北沿緯度方向由于變形在距離計算時施加了部分改正，而東西沿經(jīng)度方向可近似認(rèn)為不發(fā)生形變（不加改正），從而東西向和南北向長度不相等。依據(jù)研究區(qū)域邊長可得沿緯度方向和經(jīng)度方向頻率間隔為0.003/km 和0.0031/km。選擇緯度方向頻率間隔和經(jīng)度方向頻率間隔平方根作為徑向頻率間隔（0.0043/km）,最大徑向頻率為0.2698/km，依據(jù)相干性計算流程，最終得到研究海區(qū)海底地形和重力異常相干性結(jié)果如圖7。從圖7 相干性結(jié)果可以看出，海底地形和重力異常在有限波段表現(xiàn)出較強的相干性。需要說明的是，由于實際數(shù)據(jù)更加復(fù)雜，包含的信息量更多，噪聲模式不盡相同，從而不同區(qū)域地形和重力異常相干性結(jié)果也不盡相同。

圖7 海底地形和重力異常相干性結(jié)果Fig.7 The result about coherency between seafloor topography and gravity anomaly

3 結(jié)論

本文一方面以EGM2008、EIGEN-6C4 全球重力場位系數(shù)模型和DTM2006、EARTH2012 全球地形球諧系數(shù)模型為數(shù)據(jù)輸入，研究了全球范圍地形自功率譜、重力異常自功率譜和地形-重力互功率譜特點，計算分析了全球范圍地形-重力異常相干性頻譜特點。另一方面，使用某3°×3°（12°N～15°N，112 °E～115°E）海域作為真實研究海區(qū)，充分研討了局部區(qū)域海底地形-重力異常間相干性特點，分析了重力場元信噪比、非線性項海底地形和地殼均衡等影響干擾局部地形-重力相干性結(jié)果情況。通過實例計算得出：全球范圍地形-重力異常在120 階～1000 階左右表現(xiàn)強相干性，對應(yīng)波段范圍是25.45 km～210.61 km。研究結(jié)果可為利用衛(wèi)星測高重力數(shù)據(jù)恢復(fù)海底地形信息截止頻段的確定提供參考與借鑒。