動態與靜態平差方法在流動重力數據處理中的對比研究

郝洪濤 隗壽春 韋 進 劉少明 胡敏章

1 中國地震局地震大地測量重點實驗室，武漢市洪山側路40號，430071 2 湖北省地震局，武漢市洪山側路40號，430071 3 中國地震局第二監測中心，西安市西影路316號，710054

地表重力場變化監測是我國開展地震預測研究的重要手段之一[1-2]，目前主要采取絕對重力與相對重力聯測的混合測量模式獲取較大空間范圍內的地表重力場變化[1]。在數據平差處理方面，20世紀80年代，我國地震重力監測人員采用基于最小二乘原理的經典平差法研制出流動重力數據平差處理軟件[3]，并一直沿用至今。但該軟件在算法模型上存在隱含的假設條件，即在重力觀測網單期觀測的整個時段內，測點重力值保持不變，因此該方法也稱為靜態平差方法。然而地球系統動力學過程復雜，地表重力場會隨時間發生變化，因此基于該方法的數據處理結果會存在因重力變化而產生的偏差[4]，進而影響重力變化計算結果以及后續相關應用的研究。早期的重力監測網范圍相對較小且分散獨立[5-6]，單個測網單期觀測時間較短，影響相對較小。近年來，中國地震局建立范圍覆蓋中國大陸的地震重力監測網絡，進行每年1期的整體觀測，并在南北帶、大華北、天山等地區實施1期加密觀測(全年觀測2期)。在觀測組織實施上，重力網全網觀測由近20家單位分測區共同完成，單期觀測通常持續3～4個月，南北帶等地區全年2期觀測之間未出現明顯間斷，因此不可忽視重力網觀測期間重力變化因素對計算結果的影響。

部分學者對基于動態平差模型的流動重力網復測數據處理方法進行研究，其基本思路是在平差觀測方程中將重力值描述為與時間相關的變量，從而在數據處理模型上顧及重力值的時變因素[4,7-8]。Pagiatakis等[7]利用重力隨時間變化的線性速率模型對加拿大重力標準網超過40 a的觀測數據進行動態平差處理，獲得該地區的重力長期變化速率圖；康開軒等[4]利用相同的平差模型對滇西區域重力網資料進行處理，獲得該地區重力場長期變化趨勢；隗壽春等[8]將一種分段線性動態平差模型引入我國地殼運動觀測網絡重力網資料中，認為相比于傳統靜態平差模型，分段線性動態平差模型能更有效地反映真實重力場的變化信息。

上述研究為利用動態平差方法進行流動重力網數據處理提供了良好基礎，但靜態和動態平差方法獲得的重力變化處理結果的差異，尤其是不同時間尺度下重力變化結果的差異對比較為欠缺。基于上述分析，本文以我國南北地震帶南段地區流動重力觀測數據為例，分別采用傳統靜態平差和動態平差2種方法進行處理，對0.5～3 a時間尺度下的重力變化結果進行對比分析，探討重力時變因素對靜態平差方法不同時間尺度下重力變化結果的影響，進而為流動重力網數據處理提供參考。

1 平差模型

1.1 靜態平差模型

傳統靜態平差的基本觀測方程為：

(1)

(2)

根據最小二乘準則，可求得各測點重力值：

(3)

1.2 動態平差模型

(4)

2 數據與處理

南北地震帶南段是我國地震重力監測的重點區域之一。該區域自2014年起已形成常態化的整體觀測模式，其中相對重力觀測由四川省地震局(簡稱“四川局”，承擔四川測區)、云南省地震局(簡稱“云南局”，承擔云南測區)、中國地震局第一監測中心(簡稱“一測中心”，與湖北局共同承擔貴州測區)、湖北省地震局(簡稱“湖北局”，承擔重慶測區，與一測中心共同承擔貴州測區)等4家單位分測區共同完成，絕對重力觀測由湖北局、二測中心、一測中心、中國地震局地球物理研究所(簡稱“地球所”)等4家單位共同完成，圖1為絕對重力測點與各相對重力測區的測線分布。

紅色五角星為絕對重力測點，線段為相對重力測線，各測區以不同顏色區分圖1 南北帶南段地區重力觀測網

在研究重力場變化特征時，通常會分析0.5～3 a時間尺度下的重力場變化圖像[1,9]。本文選取2015年第二期、2016年第二期、2017年第二期和2018年全年2期總計5期的觀測資料，獲取0.5～3 a時間尺度下的重力變化結果。表1為5

表1 相對重力觀測數據概況

期觀測資料概況，由表可見，每期觀測的持續時間均為3個月左右，2018年全年2期觀測中僅存在約1個月的空期。2015年以來，研究區內先后發生3次6級以上地震，分別為2017-08-08九寨溝7.0級地震、2019-06-17長寧6.0級地震和2021-05-21漾濞6.4級地震，上述地震均可作為震例分析研究區重力場變化特征與地震的關系。

采用傳統靜態平差和動態平差2種方法進行數據處理，分別獲取0.5～3 a時間尺度下的重力變化結果，并分析2種方法所得重力變化結果的差異。首先對相對重力觀測數據進行預處理，包括固體潮、氣壓、儀器高和零漂改正等，然后分別采用傳統靜態平差和動態平差方法進行平差計算。

靜態平差計算時，首先采用線性擬合方法將多期絕對重力觀測結果內插到相對重力觀測的平均觀測時間中獲得起算基準值；然后與單期相對重力觀測數據預處理結果進行聯合解算，獲得各期次各測點的重力值；最后采用2期差分結果獲得重力變化。計算結果表明，各期次點值精度均值均優于±10 μGal，各時段重力變化精度均值均優于±15 μGal。

動態平差計算時，只需按照絕對重力觀測的實際測量結果直接給出絕對值和觀測時刻即可得到起算基準值；然后聯合相應2期相對重力觀測數據的預處理結果，計算重力時變速率和設定初始時刻的重力值。

動態平差方法的計算結果為時變速率，而靜態平差方法的計算結果為2期觀測的變化量。為進行對比分析，本文根據不同期次各測點的觀測時間差，將動態平差計算的時變速率結果轉換成重力變化量，再與靜態平差獲得的重力變化結果進行對比。由動態平差方法計算得到的各時段重力變化精度均值均優于±15 μGal，具體結果見表2(單位μGal)，其精度與靜態平差方法基本相當。

表2 靜態和動態平差方法計算的各時段重力變化精度均值

3 結果分析與討論

3.1 2種方法重力變化結果的數值差異

圖2為動態平差和靜態平差2種方法獲得的不同時間尺度下重力變化結果的差異分布累積概率曲線。由圖可見，隨著時間尺度的增加，2種方法獲得的重力變化結果差異逐步減小。在0.5 a時間尺度下，2種方法的最大差異約為19 μGal，平均約為2.9 μGal；在概率分布上，差異為2 μGal以內的占比約50%，差異為5 μGal以內的占比約85%，差異為10 μGal以內的占比約98%。在1 a時間尺度下，2種方法的差異明顯減小，最大約為9.5 μGal，平均值為1.5 μGal；差異為2 μGal以內的占比約80%，差異為5 μGal 以內的占比達95%。在2 a和3 a時間尺度下2種方法的差異進一步減小。2 a時間尺度下的最大差異約為5.5 μGal，均值為0.6 μGal；3 a 時間尺度下的最大差異為4.0 μGal，均值為0.6 μGal，與2 a時間尺度下的結果相同。從概率分布結果來看，2 a和3 a時間尺度下的累積概率分布曲線基本重合，說明對于2 a以上時間尺度下的重力變化結果，單期觀測期間重力測點時變因素的影響已趨于穩定，不會再隨著時間尺度的增加發生快速衰減。

圖2 靜態與動態平差方法重力變化結果的差異

3.2 2種方法重力變化圖像對比

重力變化圖像是分析重力變化特征的主要依據，圖3～6分別為基于2種方法獲得的0.5～3 a時間尺度下的重力變化圖像。

3.2.1 0.5 a尺度

由圖3可見，2種方法的重力變化并不顯著，均在30 μGal以內。2種方法在空間分布態勢上基本一致，主要區別在于部分地區的等值線分布細節和重力變化量級。研究區北部以大范圍的負變化為主，在九寨溝震區西北部等地呈局部正變化。相比于靜態平差結果，動態平差結果的正變化區域相對較大，且重力變化零值線穿過震中地區。在滇西南至川藏交界地區，2種方法的結果均為低值正變化，但等值線分布形態略有差異。在漾濞6.4級地震震區附近，動態平差結果正變化量級比靜態平差結果略小，兩者最大差異約為15 μGal。在大涼山次級塊體東南側的毗鄰地區，動態平差結果正變化量級比靜態平差大約15 μGal。2019-06-17在正變化區域東北緣近零值線地區發生長寧6.0級地震。

圖3 2018-04～08重力變化結果對比

對于約0.5 a時間尺度下的重力變化而言，動態平差和靜態平差2種方法的重力變化圖像在整體空間分布態勢上保持一致，但在具體的空間分布細節和重力變化量級上存在一定的差異。2種方法圖像顯示出的重力變化量級最大差異約為15 μGal，與上述重力變化數值的差異基本一致。

圖4 2017-09～2018-08重力變化結果對比

3.2.2 1 a尺度

由圖4可見，相比于0.5 a尺度，1 a尺度下2種方法的差異明顯減小，僅在重慶北部、大涼山次級塊體西北部等地區等值線的分布細節上略有不同。在重力變化空間分布特征上，相比于0.5 a尺度，1 a尺度下巴顏喀拉塊體內部的正變化區域發生擴展，塊體東緣龍門山斷裂一線地區負變化明顯增強，與鮮水河斷裂交會處的最大負變化可達60 μGal。九寨溝地震仍位于巴顏喀拉塊體東北部正負變化分界的零值線上。蓮峰-昭通斷裂帶南側正變化區域向西收縮，最大量級增大至約45 μGal。滇西南地區仍維持低值正變化，但滇西北和川滇藏交界地區整體轉為負變化，在漾濞 6.4 級地震北側形成近東西向的零值線。整體而言，該時段相對顯著的重力變化區基本沿巴顏喀拉塊體東邊界斷裂和大涼山次級塊體邊界斷裂展布，研究區東部和云南南部地區變化平緩。

3.2.3 2 a尺度

由圖5可見，2種方法的差異進一步減小，僅在少數地區存在等值線分布上的細節差別。在重力變化空間分布態勢上，巴顏喀拉塊體東邊界斷裂、大涼山次級塊體邊界斷裂和小江斷裂帶的控制作用更加明顯，顯著變化區均位于斷裂帶上或以西地區，東部地區變化平緩。相比于1 a尺度，2 a尺度下巴顏喀拉塊體和川滇菱形塊體北部的重力正變化進一步增強，龍門山斷裂帶一線負變化區域向東收縮，形成與斷裂帶大致平行的NNE向梯度帶和零值線，九寨溝7.0級地震位于梯度帶和零值線上。蓮峰-昭通斷裂帶南側至滇西南的NE向條帶地區仍維持正變化，但滇西北地區負變化區域向南突刺，在漾濞6.4級地震震中地區形成近似四象限的中心特征。

圖5 2016-09～2018-08重力變化結果對比

3.2.4 3 a尺度

由圖6可見，2種方法已基本無差異。顯著變化區仍位于巴顏喀拉塊體東邊界斷裂、大涼山次級塊體邊界斷裂和小江斷裂帶上或以西地區，重力變化整體呈西強東弱態勢。相比于2 a尺度，3 a尺度下巴顏喀拉塊體和川滇菱形塊體北部的重力正變化有所減弱，但鮮水河斷裂南段地區的局部高值異常更加突出。龍門山斷裂帶一線地區負變化顯著增強，九寨溝7.0級地震仍位于與斷裂帶大致平行的零值線上。蓮峰-昭通斷裂帶南側地區正變化明顯增強，長寧6.0級地震位于正變化區東北側的零值線上。滇西南地區正變化減弱，漾濞6.4級地震位于其零值線上。

圖6 2015-08～2018-08重力變化結果對比

4 討 論

動態平差與靜態平差的主要區別在于動態平差模型可顧及單期觀測期間測點的重力變化因素，當忽略單期內觀測時間的差異時(即假設單期觀測內各測點的觀測時間相同)，動態平差模型則等同于對多期觀測數據同時進行靜態平差，即同時計算多期觀測期間的測點重力值，此時動態平差與靜態平差2種方法獲得的重力變化結果應該相同。為驗證上述結論，本文將2018年2期觀測數據的觀測時間人為設定為相應期次的平均觀測時間，分別采用動態和靜態平差方法計算重力變化并進行對比分析。圖7為2種方法重力變化結果的數值差異，由圖可見，最大差異僅約0.1 μGal，應為平差計算時觀測方程不同所產生的計算誤差。這也從另一方面證明了本文動態平差方法模型的正確性。

圖7 忽略單期內觀測時間差異時2種方法獲得的約0.5 a時間尺度下的重力變化結果差異

參與動態平差計算的不同期次觀測資料時間跨度越大，動態平差模型越接近于靜態平差。基于2種方法計算的0.5～3 a時間尺度下的重力變化結果也表明，隨著時間尺度的增加，2種方法重力變化結果的差異也逐步減小。對于0.5 a尺度下的重力變化，2種方法在數值上的差異最大約20 μGal，2種方法在整體空間分布態勢上能保持一致，但在具體的空間分布細節和重力變化量級上存在一定的差異。對于1 a尺度下的重力變化，2種方法在數值上的差異最大約10 μGal，部分地區的等值線分布存在細節差異。當重力變化時間尺度在2 a以上時，2種方法獲得的重力變化結果在數值上的差異在5.5 μGal以內，同時重力變化圖像基本不存在差異。上述結論表明，在0.5 a和1 a尺度下的重力變化結果中，單期觀測期間測點的重力變化因素會對基于靜態平差方法的重力圖像造成一定的干擾，而在2 a以上時間尺度下的重力變化結果中單期觀測期間測點重力變化因素的影響較小。

重力變化在空間分布上表現出西強東弱的整體態勢，顯著的重力變化區主要集中在巴顏喀拉塊體、大涼山次級塊體、川滇菱形塊體中北部地區，東部揚子塊體地區重力變化較小。從研究區構造背景來看，西部川滇地區位于青藏高原東南緣，新生代以來在青藏高原地殼物質向東側運移和阿薩姆突刺的共同作用下，該地區新構造變形和地震活動十分強烈，是中國大陸最顯著的強震活動區域[10-11]。研究區東部揚子塊體在新構造運動中屬于較為穩定的地塊，地震活動性相對較弱[12]。由此可知，重力場變化的空間分布可大致反映研究區整體構造活動的強弱。

從重力場變化特征與地震發震地點的關系來看，近期研究區內3次6級以上地震的發震地點均與重力變化零值線具有較好的對應關系。動態平差方法圖像顯示，九寨溝7.0級地震在0.5～3 a 尺度圖像中均位于巴顏喀拉塊體東邊界斷裂帶與重力變化零值線的交會位置；長寧6.0級地震在0.5～3 a尺度圖像中均位于零值線上；漾濞6.4級地震位于川滇菱形塊體南邊界斷裂附近，在3 a尺度圖像中位于零值線上，在2 a尺度圖像中位于四象限中心和近零值線位置。重力場演化規律與構造活動關系的現有研究結果表明，活動地塊的邊界帶在重力場演化過程中往往具有重要的控制作用，大地震一般發生在與主要活動斷裂方向一致的重力場變化正負轉換帶上[13]。基于上述分析可以認為，本文基于動態平差方法獲得的0.5～3 a尺度下的重力場變化圖像可大致反映3次地震的發震背景。

5 結 語

1)南北地震帶南段地區流動重力網單期觀測持續時間約3個月，觀測期間的重力時變因素會導致靜態平差方法計算的重力變化結果存在誤差，但隨著重力變化時間尺度的增加，其影響會逐步減小。對于0.5 a、1 a、2 a和3 a尺度下的重力變化結果，重力時變因素的最大影響分別約為19 μGal、9.5 μGal、5.5 μGal和4.0 μGal。

2)動態平差和靜態平差2種方法的重力變化圖像在整體空間分布態勢上保持一致，但0.5 a和1 a時間尺度圖像在重力變化量級和等值線分布細節上存在一定的差異，2 a和3 a時間尺度下的重力變化圖像則基本相同。因此對于計算0.5 a和1 a時間尺度下的重力變化宜采用動態平差方法。

3)近期研究區內3次6級以上地震的發震地點均與重力變化零值線具有較好的對應關系，基于動態平差方法獲得的0.5～3 a時間尺度下的重力場變化圖像可反映3次地震的發震背景。