最小二乘配置與REHSM 求解GPS應變場的方法

管守奎 瞿 偉 蔣 軍

1 武漢大學測繪學院，武漢市珞瑜路129號，430079

2 長安大學地質工程與測繪學院，西安市雁塔路126號，710054

3 武漢大學GNSS中心，武漢市珞瑜路129號，430079

應變場是描述區域變形的一項重要指標。各國學者基于GPS觀測資料對應變場開展研究，取得許多成果［1－5］。但是當研究區域僅有少量觀測數據時，對應變場的精確求解就會變得困難，而最小二乘配置能夠較好地解決觀測數據較少的問題。本文提出使用剛體旋轉模型（RRM），采用最小二乘配置解法獲得較為充分的虛擬觀測速度值，之后結合REHSM 方法，解決觀測數據量較少時應變場的解算問題。

1 最小二乘配置與REHSM 解算應變場

1.1 REHSM 解算應變場

歐拉定理是現代板塊運動定量描述的基本定理。如果把板塊看成剛性，地球近似為球體，球心看成強制在地球表面運動的剛體板塊運動的固定點，則地殼運動滿足歐拉定理。設一點的經緯度為（λ，φ），則有：

式中，r為地球半徑，Ve、Vn為站心參考系的經向和緯向速度，ωe、ωen、ωn為板塊的旋轉參數分量［6］。該式稱為塊體的剛體旋轉模型（RRM）。

一個塊體在周圍塊體作用下發生整體旋轉的同時，內部也將發生應變。塊體內部的應變可分解為兩部分［7］：一是塊體的純應變，二是由應變引起的塊體旋轉。表示為：

x、y分別為觀測點至研究區域中心經線與緯線間的平行圈弧長與子午線弧長，北東為正。εe、εen、εn分別是塊體的東西向線應變、東西向和南北向之間的剪應變、南北向線應變。該式稱為塊體的整體旋轉與均勻應變模型（REHSM）。

由上式可知，對于一個觀測點，ε矩陣不能直接求得。本文對整個研究區進行柵格劃分，以每個柵格4個頂點的速率信息求出的應變參數來代替柵格中心的應變參數。求得應變參數后，即可進一步計算應變場的特征值［8］。最大主應變率：

最大剪應變率：

面膨脹率：

1.2 最小二乘配置原理

最小二乘平差將所有待估參數當作非隨機量，或不考慮參數的隨機性質。最小二乘配置除包括經典最小二乘的非隨機參數外，還包含隨機參數，按照極大驗后估計、最小方差估計或廣義最小二乘原理求定參數的最佳估值［9］。函數模型為：

式中，Y′是非觀測向量，在式（6）中其系數為0，并未包含在l內，而是依賴于協方差相聯系，在對濾波模型作參數估計的同時對Y′進行推估。

配置模型的最優線性無偏估計為：

在最小二乘配置理論中，關鍵是求取已知點信號的方差－協方差矩陣、已知點信號與未知點信號的協方差矩陣。這兩個矩陣都是由同一個協方差函數計算得到的，協方差函數選擇恰當與否直接影響估計值的精度。本文分別采用以下6個協方差擬合函數進行計算。

1）多項式函數：

2）高斯函數：

3）似高斯函數：

4）希爾沃寧函數：

5）空間調和函數：

Z＊＝Hi＋Hj，Hi、Hj為數據點的海拔高，考慮到d＝0的特殊情況，取Z＊＝｜Hi－Hj｜，但并不影響函數的對稱、正定和協調性。該函數可以作為經驗協方差函數：

式中，D為協方差，d為兩測站點間的距離，D（0）為所有觀測點方差的平均值，k為待定系數。通過不同的協方差擬合函數比較，獲得最優的速度場擬合信息，為應變場的解算提供數據支持［10］。

2 計算實例

2.1 速度場構建

將位于運城盆地的基于“數字地震網絡工程”的19個監測點分為兩組，13個為已測點，6個為檢核點。通過6種協方差擬合函數，比較獲得最優的擬合方式。根據表1，采用希爾沃寧函數為協方差擬合函數時，E方向上為1.046 4mm／a，N方向上為2.009 6mm／a，兩個方向的殘差最大值均小于其他協方差擬合函數的結果。采用希爾沃寧函數為速率場構建時的協方差擬合函數，對所有的觀測速度進行擬合，實測與擬合的速度場如圖1，具體速度殘差見表2。

由圖1可直觀看出，19個監測點的實測結果與擬合結果大小與方向較為符合。由表2 可看出，E方向上的殘差絕對值最大值為0.775 6 mm／a，N方向為1.218 7mm／a，擬合結果良好。

表1 速率檢核統計Tab.1 The statistics of velocity check

圖1 實測速率與擬合速率比較Fig.1 Compared the measured velocity and fitted velocity

表2 速率擬合結果Tab.2 The result of velocity fitting

協方差擬合函數確定后，根據最小二乘配置方法構建的速率場信息如圖2。由圖2 可知，速度場擬合結果比較緊密。研究區域SN 方向的速率較EW 方向小，整體為水平偏南走向，同武艷強［3］、李延興［6］的結果基本一致。

圖2 擬合速率場Fig.2 Fitted velocity field

2.2 應變場構建

有了充足的速度信息，結合REHSM 與柵格化的解算方法，即可得到研究區域應變場特征參數的等值線圖。

由圖3～5 可以看出，高值區域主要分布在WS、EN 區域，低值區域主要分布在WN－ES 區域。高值區域為張性或應變率較為強烈的地區，數量級最大達1×10－8，反映這些區域是應變能高積聚區，發生中強地震的可能性較大。結合運城盆地的地質背景可知，WS、EN 區域為中條山斷裂帶，驗證了應變參數求解的有效性。圖6中黑色箭頭為張性劇烈程度增加走向。可以看出，中條山斷裂兩側向中條山區域兩側膨脹加劇，中條山兩側均受擠壓，在碰撞處向WS、EN 兩個方向釋放能量，導致WS、EN 即中條山斷裂走向區域張性提升。膨脹區域在膨脹交匯處，因各自膨脹方向均不相同，地殼運動活躍。據山西省地震局通報，2010－12－14 11：45在山西省運城市聞喜縣、夏縣交界發生3.3級地震，震中位于35.3°N、111.1°E。此區域正是圖6中的橢圓區域，故而驗證了該區域地殼運動的活躍性。

圖3 面膨脹率等值線圖Fig.3 The isogram of superficial expansion rate

圖4 最大剪應變等值線圖Fig.4 The isogram of shear strain rate

圖5 最大主應變等值線圖Fig.5 The isogram of maximum principal strain rate

3 結 語

本文使用的地殼形變參數求解方法，能夠有效解決監測點數據較少時應變場的解算問題。通過塊體的剛體旋轉模型，采用最小二乘配置以不同的協方差擬合函數獲得較優的速率場，并同之前的全國性速度場結果進行比較。再將研究區域進行柵格化，通過塊體的均勻旋轉與應變模型，解決區域地殼應變場的求解問題。

圖6 面膨脹率分析Fig.6 The superficial expansion rate analysis

［1］Dumak B S，Kotlia B S，Kumar K，et al.Crustal Deformation Revealed by GPS in Kumaun Himalaya India［J］.Journal of Mountain Science，2014（1）：41－50

［2］江在森，馬宗晉，張希，等.GPS初步結果揭示的中國大陸水平應變場與構造變形［J］.地球物理學報，2003（3）：352－358（Jiang Zaisen，Ma Zongjin，Zhang Xi，et al.Horizontal Strain Field and Tectonic Deformation of China Mainland Revealed by Preliminary GPS Result［J］.Chinese Journal of Geophysics，2003（3）：352－358）

［3］武艷強，江在森，楊國華，等.利用最小二乘配置在球面上整體解算GPS應變場的方法及應用［J］.地球物理學報，2009（7）：1 707－1 714（Wu Yanqiang，Jiang Zaisen，Yang Guohua，et al.The Application and Method of GPS Strain Calculation in Whole Mode Using Least Square Collocation in Sphere Surface［J］.Chinese Journal of Geophysics，2009（7）：1 707－1 714）

［4］黃立人，王敏.中國大陸構造塊體的現今活動和變形［J］.地震地質，2003（1）：23－32（Huang Liren，Wang Min.Present－Day Activity and Deformation of Tectonic Blocks in Chinese Continent［J］.Seismology and Geology，2003（1）：23－32）

［5］瞿偉，張勤，王慶良，等.渭河盆地現今地殼水平形變特征及區域構造活動性［J］.武漢大學學報：信息科學版，2011（7）：830－834（Qu Wei，Zhang Qin，Wang Qingliang，et al.Research on Present Crustal Horizontal Deformation Feature of Weihe Basin and Its Tectonic Activity［J］.Geomatics and Information Science of Wuhan University，2011（7）：830－834）

［6］李延興，李智，張靜華，等.中國大陸及周邊地區的水平應變場［J］.地球物理學報，2004，47（2）：222－231（Li Yanxing，Li Zhi，Zhang Jinghua，et al.Horizontal Strain Field in the Chinese Mainland and Its Surrounding Areas［J］.Chinese Journal of Geophysics，2004，47（2）：222－231）

［7］許才軍，張朝玉.地殼形變測量與數據處理［M］.武漢：武漢大學出版社，2009（Xu Caijun，Zhang Chaoyu.Crustal Deformation Measurement and Data Processing［M］.Wuhan：Wuhan University Press，2009）

［8］瞿偉.基于GPS速度場分析研究青藏塊體東北緣的形變特征［D］.西安：長安大學，2008（Qu Wei.The Deformation Characteristic of the Northeastern Margin of Qinghai－Tibet Plateau Constrained by GPS Observations［D］.Xi’an：Chang’an University，2008）

［9］崔希璋，陶本藻.廣義測量平差［M］.北京：測繪出版社，1992（Cui Xizhang，Tao Benzao.Generalized Survey Adjustment［M］.Beijing：Surveying and Mapping Press，1992）

［10］李沖，李建成，瞿偉.基于移動原理的擬合推估模型建立區域地殼運動速率場［J］.大地測量與地球動力學，2012（5）：33－36（Li Chong，Li Jiancheng，Qu Wei.Establishing Regional Crustal Movement Velocity Field with Collocation Based Displacement Principle［J］.Journal of Geodesy and Geodynamics，2012（5）：33－36）