范 雕，李?yuàn)檴?，趙東明，楊軍軍，孟書(shū)宇，范昊鵬

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

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

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

1 理論與方法

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

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

1.1 球地形與重力場(chǎng)元間相干性

由重力場(chǎng)知識(shí)可知，重力異?？杀硎緸椋?/p>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

重力場(chǎng)元與海底地形互功率譜密度為：

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

令重力場(chǎng)元信號(hào)和測(cè)量噪聲的信噪比表達(dá)式為：

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

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

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

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

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

根據(jù)式(20)可知，若輸入信號(hào)為功率有限信號(hào)，則按照式(20)計(jì)算信號(hào)間相干性，不論信號(hào)是否相關(guān)以及線性相關(guān)強(qiáng)弱，二者頻率域相干性恒為1。針對(duì)這一問(wèn)題，通常應(yīng)對(duì)策略為將式(20)解算的功率譜密度進(jìn)行平滑處理：式中，<>表示方位角α范圍內(nèi)的平均值。方位角的取平均值實(shí)際操作流程如下[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為徑向頻率間隔，分別計(jì)算環(huán)徑為dq范圍內(nèi)的各信號(hào)功率譜密度平均值，然后依照式(16)解算信號(hào)相干性。

2 數(shù)值分析試驗(yàn)

2.1 算例1：全球地形和重力場(chǎng)元相干性

當(dāng)前國(guó)際廣泛使用的超高階（大于1000 階）重力場(chǎng)模型有EGM2008、GECO和EIGEN系列等模型等。本文以EGM2008（2160 階）和EIGEN-6C4（2190 階）地球重力場(chǎng)模型作為計(jì)算重力異常譜的重力場(chǎng)輸入，以EARTH2012（2160 階）和DTM2006（2190 階）地形球諧系數(shù)模型作為解算地形譜的數(shù)據(jù)輸入，依據(jù)式(3)和式(5)計(jì)算的重力異常功率譜密度和地形功率譜密度結(jié)果如圖2 和圖3 所示，其中紅色曲線分別表示EGM2008 和EARTH2012 位系數(shù)模型計(jì)算結(jié)果，藍(lán)色曲線分別表示EIGEN-6C4 和DTM2006 位系數(shù)模型解算結(jié)果?？梢钥闯?，重力異常和地形的功率譜密度均隨著模型階數(shù)的增加而減小，說(shuō)明模型中長(zhǎng)波信號(hào)占主要部分，高頻信號(hào)能量較弱；而EIGEN-6C4 模型功率譜密度在2160 階后出現(xiàn)急劇變化，這與模型構(gòu)建機(jī)制有關(guān)，如使用的數(shù)據(jù)分辨率等因素。根據(jù)式(6)計(jì)算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)先升高后降低的特征。由相干性定義可知，相干性表示一種信號(hào)可線性轉(zhuǎn)換成另一種信號(hào)的程度：當(dāng)相干性為1 時(shí)即兩種信號(hào)完全線性相關(guān)；相干性為0 表示兩信號(hào)沒(méi)有線性關(guān)系[23]。當(dāng)相干性為0.5 可理解為信噪比（signal-to-noise，SNR）為1:1 的信號(hào)與另一不含噪聲信號(hào)的相干性結(jié)果[24]。從圖5 可以看出，全球范圍內(nèi)重力異常和地形的相干性在0.5 以上對(duì)應(yīng)的模型階數(shù)大約為120階～1000階，對(duì)應(yīng)的波長(zhǎng)范圍為25.45 km～210.61 km，該計(jì)算結(jié)果與目前國(guó)際國(guó)內(nèi)學(xué)者利用重力數(shù)據(jù)恢復(fù)海底地形的波段范圍一致[6,7,17,25-28]。

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

選擇某3 °× 3°（12°N～15° N，112 °E～115°E）海域作為研究海區(qū)，遵照以上理論分析研究海區(qū)海底地形和重力異常頻域相干性大小。研究海區(qū)水深數(shù)據(jù)（格網(wǎng)化水深）和衛(wèi)星測(cè)高重力異常數(shù)據(jù)（S&SV24.1）空間分辨率為1 角分，空間分布情況如圖6 所示。相應(yīng)數(shù)據(jù)統(tǒng)計(jì)結(jié)果見(jiàn)表1，數(shù)據(jù)類別括號(hào)內(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)計(jì)算，研究區(qū)域東西向總長(zhǎng)為324.37 km，南北向總長(zhǎng)為333.58 km，其中南北沿緯度方向由于變形在距離計(jì)算時(shí)施加了部分改正，而東西沿經(jīng)度方向可近似認(rèn)為不發(fā)生形變（不加改正），從而東西向和南北向長(zhǎng)度不相等。依據(jù)研究區(qū)域邊長(zhǎng)可得沿緯度方向和經(jīng)度方向頻率間隔為0.003/km 和0.0031/km。選擇緯度方向頻率間隔和經(jīng)度方向頻率間隔平方根作為徑向頻率間隔（0.0043/km）,最大徑向頻率為0.2698/km，依據(jù)相干性計(jì)算流程，最終得到研究海區(qū)海底地形和重力異常相干性結(jié)果如圖7。從圖7 相干性結(jié)果可以看出，海底地形和重力異常在有限波段表現(xiàn)出較強(qiáng)的相干性。需要說(shuō)明的是，由于實(shí)際數(shù)據(jù)更加復(fù)雜，包含的信息量更多，噪聲模式不盡相同，從而不同區(qū)域地形和重力異常相干性結(jié)果也不盡相同。

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

3 結(jié)論

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