郝 蓉，張菊清

(長安大學地質工程與測繪學院，陜西西安710054)

一、引 言

InSAR由于具有高精度、高分辨率等特點，是目前地表變化監(jiān)測的一種重要手段。然而，由于受時間、空間去相干等因素的影響，提取的形變信息中經常存有缺失現(xiàn)象，需要后期對其進行填補處理。對于缺失數(shù)據(jù)的擬合，常采用空間數(shù)據(jù)插值方法進行擬合填補。根據(jù)模型不同，空間數(shù)據(jù)插值法可分為函數(shù)模型法和隨機模型法兩種［1］。常用的函數(shù)模型法有曲線內插法、曲面內插法、多面函數(shù)內插法、多項式內插法及移動曲面擬合法等，隨機模型法主要有反距離加權法、Kriging插值法和擬合推估法等。

在眾多的插值和擬合方法中，擬合推估由于其既能較好地反映趨勢性變化成分，又能很好地刻畫帶有局部相關性的隨機特征，已廣泛應用于重力場推估［2-4］、高程轉換［5-7］、地殼形變分析［8-10］、地圖數(shù)字化糾正及坐標變換［11-14］之中。然而在擬合推估求解過程中，隨機信號協(xié)方差函數(shù)的確定是難點和關鍵，已引起眾多學者的關注。劉念等針對異常誤差對觀測值的影響提出了協(xié)方差函數(shù)的抗差擬合方法［15-16］;劉漢江等結合擬合推估實際應用，進行了局部協(xié)方差函數(shù)的擬合研究［17-18］;金保平等則利用遺傳規(guī)劃理論研究了協(xié)方差函數(shù)的確定問題［19］?？紤]到經擬合得到的信號向量協(xié)方差與觀測向量的先驗協(xié)方差很難一致，進而影響推估精度的問題，楊元喜等提出了應用方差分量估計方法來重新調整信號向量與觀測向量的先驗權比［20-21］，并在此基礎上進一步提出了自適應擬合推估方法，即利用方差分量估計構建觀測誤差與隨機信號之間的自適應因子，以此來調整觀測向量與信號向量對模型參數(shù)估計的貢獻［22-23］，并在高程轉換、地殼形變分析等領域得到了很好的應用。

然而，現(xiàn)有的協(xié)方差函數(shù)大都是基于隨機信號具有各向同性的二階寬平穩(wěn)隨機過程的假設條件完成的，而實際上各向異性現(xiàn)象更加普遍。Steffen H G針對大氣噪聲的各向異性進行了協(xié)方差函數(shù)的擬合研究，并有效地改善了形變監(jiān)測精度［24］。考慮到變異函數(shù)具有刻畫區(qū)域化變量的結構性和隨機性的能力，它通過本身的結構及各項參數(shù)可以從不同角度反映空間的變異特性［25］，進而能更好地展現(xiàn)空間變量的異常特征，突出異常的局部變化性和空間結構信息?；诖?，本文提出利用幾何各向異性的特性，借助變異函數(shù)并通過線性變換將隨機信號由各向異性轉換為各向同性，在此基礎上擬合協(xié)方差函數(shù)，并通過方差分量估計構建自適應因子，以此調節(jié)觀測噪聲與隨機信號對擬合推估解的影響，最后將其應用于InSAR缺失數(shù)據(jù)的擬合后處理中。

二、變異函數(shù)與協(xié)方差函數(shù)

變異函數(shù)是地統(tǒng)計學特有的基本工具，它既可以描述區(qū)域化變量的空間結構性，又可以解釋其隨機性。變異函數(shù)可表示為

式中，γ為變異函數(shù);Z為區(qū)域化變量;E()表示數(shù)學期望。若隨機變量Z具有二階寬平穩(wěn)隨機特性，則其變異函數(shù)γ僅與距離h有關，與位置x無關。于是式(1)可簡寫為

顯然，變異函數(shù)反映的是數(shù)據(jù)空間不相關程度與其距離的關系。從式(2)可以看出，空間不相關程度隨距離的增大而增大。

協(xié)方差函數(shù)是另一種表達空間相關性的量，可定義為:隨機場Z(x)在空間點x與x+h處的兩個隨機變量Z(x)與Z(x+h)的二階混合中心距，即

式中，μ為隨機向量的期望值。

可以證明，在二階寬平穩(wěn)假設條件下，隨機向量Z(x)的協(xié)方差函數(shù)C(h)和變異函數(shù)γ(h)之間具有如下的函數(shù)關系

式(4)可理解為先驗方差C(0)等于協(xié)方差與變異函數(shù)之和。其關系曲線如圖1所示。

圖1 協(xié)方差函數(shù)和變異函數(shù)關系

顯然，協(xié)方差函數(shù)反映的是區(qū)域化變量的相關程度，而變異函數(shù)反映的則是區(qū)域化變量的變異程度，它們從兩個不同側面反映了區(qū)域化變量間的關系。

三、基于各向異性協(xié)方差函數(shù)的擬合

1.各向同性協(xié)方差函數(shù)的擬合

隨機信號協(xié)方差函數(shù)的確定是擬合推估求解的關鍵。由于協(xié)方差函數(shù)的嚴密表達式往往很難獲得，通常利用觀測所獲得的地表信息對其進行擬合得到。即統(tǒng)計不同距離間的協(xié)方差，選擇一種合適的協(xié)方差函數(shù)，利用最小二乘估計擬合相應的參數(shù)值，最終確定協(xié)方差函數(shù)表達式。

常用的協(xié)方差函數(shù)有如下幾種，其中C表示協(xié)方差，它僅與數(shù)據(jù)點之間的距離d有關，C0為信號方差。

1)高斯(Gauss)函數(shù)為

式中，A2為待定參數(shù)。

2)希爾方納(Hirvonen)協(xié)方差函數(shù)為

式中，B2為待定參數(shù)。

3)指數(shù)函數(shù)為

式中，K為待定參數(shù)。

2.各向異性協(xié)方差函數(shù)的擬合

在常規(guī)的協(xié)方差函數(shù)擬合過程中，常將隨機信號視為各向同性的二階寬平穩(wěn)隨機過程，然而各向異性特征更具有普遍性。通過對各個方向上的變異函數(shù)進行研究，可以顯示出隨機函數(shù)各種各向異性的性質。各向異性是指區(qū)域化隨機變量在各個不同方向上的性質各不相同。根據(jù)性質的不同，各向異性又可分為幾何各向異性和帶狀各向異性兩類［26］。為簡單起見，本文僅討論幾何各向異性。幾何各向異性是指區(qū)域化變量在不同方向上具有相同的變異程度但具有不同的連續(xù)性，即不同方向的變異函數(shù)基臺值相同而變程不同，如圖2所示。

圖2 幾何各向異性

通常不同方向上變異性之差可以通過變程之比表示

式(8)意味著a2方向上距離為h的兩點間的平均變異程度與a1方向上距離為kh的兩點間的平均變異程度相同，即γa2(h)=γa1(kh)。

顯然，若a1=a2，則表明隨機向量各向同性，其對應的方向變程圖為圓;若a1≠a2，則隨機向量各向異性，其方向變程圖可視為橢圓(如圖3所示)。

為了將具有各向異性的隨機向量轉換成各向同性來處理，一個簡單的辦法是先通過旋轉，使其方向與橢圓軸向重合;然后通過伸縮處理，即可獲取具有相同統(tǒng)計特性的隨機向量，如圖4、圖5所示。

圖3 方向—變程圖

圖4 旋轉變換

圖5 伸縮變換

因此，在實際計算過程中，首先須根據(jù)已測信息，統(tǒng)計計算隨機向量的橢圓三參數(shù)(即旋轉角度φ、主軸長度amax與次軸長度amin)，以及主次軸長度比值 k=amax/amin。

其次通過旋轉矩陣將坐標軸旋轉角度φ，使其與橢圓軸相重合，旋轉變換矩陣為

然后進行伸縮變換，將橢圓變換為半徑等于橢圓長半軸的圓，其伸縮變換矩陣為

最后，逆旋轉－φ角，恢復原來坐標方向，其旋轉變換矩陣為

經過3次坐標變換，最后的坐標為

其中

經式(12)的線性變換，即可將隨機向量由幾何各向異性轉化為各向同性。于是可按各向同性的隨機向量進一步計算處理。

四、基于各向異性的自適應擬合推估

基于各向異性的協(xié)方差函數(shù)能夠較好地顧及隨機向量的統(tǒng)計特性，由此擬合得到的協(xié)方差函數(shù)能夠更加合理地表達向量間的相關性，但難以保證與觀測噪聲具有相同的單位權方差因子。因此為了協(xié)調隨機信號與觀測噪聲對擬合推估解的影響，楊元喜［22-23］等提出了自適應擬合推估。其基本思想是通過引入自適應因子來平衡兩者的貢獻。

假設擬合推估函數(shù)模型為

式中，A、B為系數(shù)矩陣;L為觀測向量;X為非隨機參數(shù)向量;Y為隨機信號向量，包括已測點信號S和未測點信號S'，即Y=[S S']T;Δ為觀測噪聲。

按照自適應擬合推估準則

即可獲得

式中，V為誤差向量;S為信號向量;PΔ和PS分別為V和S的權陣;α(0≤α≤1)為自適應因子;DΔ為觀測值協(xié)因數(shù)陣;DS為隨機信號協(xié)方差(BDSBT+αDΔ)－1。自適應因子 α 可通過方差比進行構建，即，其中隨機信號協(xié)方差因子估計值和觀測噪聲協(xié)方差因子估計值可應用常用的方差分量估計(如Helmert方差分量估計、最大似然方差分量估計等)進行計算得到。

五、實例驗證

為了驗證基于各向異性自適應擬合推估的有效性，本文以某區(qū)域的部分地表形變監(jiān)測信息(如圖6所示)為試驗數(shù)據(jù)，隨機設置一系列缺失數(shù)據(jù)(如圖7所示)作為檢核數(shù)據(jù)。擬合推估計算中傾向性參數(shù)A、X采用二次多項式模型f=a0+a1x+a2y+a3xy+a4x2+a5y2，隨機信號的協(xié)方差函數(shù)采用Gauss函數(shù);并在不同的先驗觀測噪聲(σ0=0.1，1，10)下，分別采用基于各向同性與各向異性的常規(guī)擬合推估方法與自適應擬合推估方法對試驗數(shù)據(jù)進行處理。計算結果通過已測信號點的殘差均方根誤差和檢核點不符值均方根誤差來體現(xiàn)。

圖6 試驗數(shù)據(jù)

圖7 缺失數(shù)據(jù)

計算結果見表1。

表1 各方案均方根誤差計算結果

分析上述計算結果，可以發(fā)現(xiàn):

1)基于各向同性與各向異性的常規(guī)擬合推估方法與自適應擬合推估方法都取得了較好的擬合效果，表明擬合推估方法在缺失數(shù)據(jù)擬合應用中具有可行性。

2)相比于各向同性的擬合推估，基于各向異性的擬合推估方法內部符合精度和外部檢核精度均有不同程度的提高，表明基于各向異性的協(xié)方差矩陣能夠更加合理地刻畫隨機向量的空間特性，從而具有更高的擬合填補精度。

3)正常擬合推估方法已測點與檢核點推估解受觀測噪聲初始方差影響嚴重，而自適應擬合推估法則不受觀測噪聲初始方差的影響，表明通過自適應因子的引入，合理地調節(jié)了觀測噪聲與信號之間的權比，進而保證兩者協(xié)方差矩陣的協(xié)調一致。

六、結束語

由于受各種因素限制，InSAR數(shù)據(jù)往往存在缺失，因此需要對缺失數(shù)據(jù)進行填補擬合。擬合推估方法既能較好地反映其空間分布規(guī)律，又能充分利用地表信息的空間相關特性，因此比較適合于缺失數(shù)據(jù)的擬合與推估。自適應擬合方法通過自適應因子的引入，保證了信號向量先驗協(xié)方差矩陣與觀測向量協(xié)方差矩陣的協(xié)調一致，進而平衡了兩者的貢獻。

此外，InSAR作為地表形變監(jiān)測的重要手段，其數(shù)據(jù)信息不僅存在一定的空間分布規(guī)律，往往還存在局部變異性，呈現(xiàn)出各向異性的特點?；诟飨虍愋缘臄M合推估方法利用變異函數(shù)可透過隨機性反映區(qū)域化變量結構性的特點，先采用變異函數(shù)對形變信息進行各向異性研究，更好地對空間形變信息進行分析，將各向異性轉化為各向同性后，再利用擬合推估方法進行擬合。這種方法充分考慮到地表形變監(jiān)測信息的結構性與隨機性，較常規(guī)擬合推估可以更好地對數(shù)據(jù)進行擬合，提高擬合精度。

此方法僅針對較為簡單的幾何各向異性進行研究，對于更為復雜的各向異性將在后續(xù)進行更為深入的研究。此外，擬合推估中協(xié)方差函數(shù)表達式的選擇也十分關鍵，其選擇也將直接影響擬合結果的可靠性。

［1］魏玉明.抗差最小二乘配置在重復重力測量數(shù)據(jù)分析中的應用［D］.西安:長安大學，2007.

［2］夏哲仁，林麗.局部重力異常協(xié)方差函數(shù)逼近［J］.測繪學報，1995，24(1):23-27.

［3］楊元喜，劉長建.重力異常的擬合推估迭代解算模型及算法［J］.地震學報，1996，18(4):475-479.

［4］章傳銀，丁劍，晁定波.局部重力場最小二乘配置通用表示技術［J］.武漢大學學報:信息科學版，2007，32(5):431-434.

［5］楊鳳蕓，張旭東.采用抗差推估法剔除GPS高程數(shù)據(jù)粗差［J］.測繪通報，2005(10):9-11.

［6］王增利，黃騰，鄧標.基于二次曲面的擬合推估法在GPS高程測量中的應用［J］.測繪工程，2009，18(1):50-52.

［7］許雙安.最小二乘配置和普通Kriging法在GPS高程轉換中的應用［J］.測繪地理信息，2012，37(6):64-67.

［8］李沖，李建成，瞿偉.基于移動原理的擬合推估模型建立區(qū)域地殼運動速率場［J］.大地測量與地球動力學，2012，32(5):33-36.

［9］江在森，劉經南.應用最小二乘配置建立地殼運動速度場與應變場的方法［J］.地球物理學報，2010，53(5):1109-1117.

［10］曾安敏，劉光明，秦顯平，等.應用擬合推估法建立我國大陸地殼水平運動速度場的研究［J］.測繪通報，2012(S0):11-15.

［11］張菊清，楊元喜，張亮.附有限制條件的擬合推估模型及其在地圖數(shù)字化誤差糾正中的應用［J］.大地測量與地球動力學，2008，28(5):91-95.

［12］張亮.抗差最小二乘配置在數(shù)字化地圖幾何糾正中的應用討論［J］.測繪與空間地理信息，2008，31(2):145-148.

［13］曾安敏.基于擬合推估的1980西安坐標系到 2000國家坐標系的變換［J］.大地測量與地球動力學，2008，28(5):125-128.

［14］曾安敏，張菊清.基于擬合推估兩步極小解法的地圖坐標變換［J］.大地測量與地球動力學，2008，28(1):81-84.

［15］劉念，胡榮明.協(xié)方差函數(shù)的待定參數(shù)對最小二乘擬合推估精度的影響［J］.西安科技學院學報，2001，21(1):61-64.

［16］劉念.協(xié)方差函數(shù)的抗差擬合［J］.測繪科學，2001，26(3):25-28.

［17］文漢江.最小二乘配置法中局部協(xié)方差函數(shù)的計算［J］.測繪科學，2000，25(3):37-39.

［18］劉敏，王昆，鄧凱亮，等.最小二乘配置中兩種局部協(xié)方差函數(shù)的比較［J］.海洋測繪，2013，33(2):16-19.

［19］金保平，岳東杰，李成仁.基于遺傳規(guī)劃的協(xié)方差函數(shù)計算及其在GPS高程擬合中的應用［J］.測繪科學技術學報，2012，29(1):16-19.

［20］楊元喜，張菊清，張亮.基于方差分量估計的擬合推估及其在GIS誤差糾正的應用［J］.測繪學報，2008，37(2):152-157.

［21］張菊清，張亮.基于Helmert方差分量估計的擬合推估法及其在地圖數(shù)字化誤差糾正中的應用［J］.測繪通報，2008(2):35-37，51.

［22］YANG Y X，XU TH.An Adaptive Kalman Filter Based on Sage Windowing Weights and Variance Components［J］.The Journal of Navigation，2003，56(2):231-240.

［23］楊元喜，曾安敏，吳富梅.基于歐拉矢量的中國大陸地殼水平運動自適應擬合推估模型［J］.中國科學:地球科學，2011，41(8):1116-1125.

［24］KNOSPE S，JóNSSóN S.Covariance Estimation for dIn-SAR Surface Deformation Measurements in the Presence of Anisotropic Atmospheric Noise［J］.IEEE Transactions on Geoscience and Remote Sensing，2010，48(4):2057-2062.

［25］程勖，楊毅恒，丁建華，等.變異函數(shù)在異?？臻g插值中的應用［J］.世界地質，2007，26(3):298-303.

［26］孫洪泉.地質統(tǒng)計學及其應用［M］.徐州:中國礦業(yè)大學出版社，1990.

［27］侯景儒.實用地質統(tǒng)計學［M］.北京:地質出版社，1998.

［28］MORITZ H.Statistical Foundations of Collocation［J］.Bollettino di Geodesia e ScienzeAffini，1980(2):131-150.