一種利用InSAR反演地下流體體積變化和三維地表形變的方法
【技術領域】
[0001] 本發明屬于基于遙感影像的大地測量領域，尤其涉及一種利用InSAR反演地下流 體體積變化和三維地表形變的方法。
【背景技術】
[0002] 差分 InSAR (Differential InSAR，D-InSAR)技術是目前國際上在 InSAR 應用上 最為成熟的技術，它最主要的目的就是監測地球表面厘米級甚至毫米級的形變。多時相 InSAR (Multi-Temporal InSAR, MT-InSAR)技術則是近20多年來在D-InSAR技術的基礎上 發展起來的一種地表形變監測方法，如PS、SBAS和TCP等。通過對多景SAR影像的聯合分 析，MT-InSAR技術可以更好地抑制干涉圖中時空失相關和大氣噪聲的影響，并能提供地表 在SAR影像獲取時刻的形變序列。然而，到目前為止，無論是D-InSAR還是MT-InSAR技術都 只能利用單一平臺、單一軌道的SAR數據，因此只能監測地表在雷達視線方向（即斜距向） 上的一維形變結果。但是在現實中，地表形變是發生在三維空間框架下的，即所謂的三維形 變場。因此InSAR的斜距向形變測量結果不能反映出真實的地表形變。
[0003] 如何將InSAR技術監測的一維形變拓展為三維，目前國際上已經有一些學者進行 了探索性的研究，大致可以分為以下幾類：（1)多方向InSAR斜距向形變測量值融合法，即 利用三個或以上的InSAR斜距向形變測量值反演三維地表形變場，但由于目前的SAR衛星 都是極軌飛行軌道，這種方法只在高煒度地區適用；（2)升降軌InSAR斜距向和方位向形變 測量值融合法，即融合D-InSAR提供的升軌、降軌斜距向形變測量值和Offset-Tracking或 MAI提供的升軌、降軌方位向形變測量值反演三維地表形變場，但受限于Offset-Tracking 或MAI技術的分米到米級的測量精度，這種方法目前只適用于監測地震、火山噴發和冰川 移動等大型地表形變；（3) InSAR和GPS融合法，即融合D-InSAR或MT-InSAR提供的斜距向 測形變量值和GPS提供的離散點上的三維形變測量值反演三維地表形變場，但這種方法的 精度強烈依賴于GPS臺站的密度和分布，在GPS臺站較少的地區不適用。
[0004] 而在現實中，往往存在著一些緩慢的地表形變需要進行大范圍、全方位和長期的 監測，如地下水抽取、油氣田開采和巖漿活動等地下流體運動引起的地表形變。這些地表形 變已經成為影響區域經濟和社會可持續發展的重要因素。例如，由于長期過度集中開采地 下水資源，我國迄今為止已有70多個城市發生了不同程度的地面沉降，其中上海和天津市 區的最大沉降已超過兩米，嚴重制約了城市社會經濟的發展，同時對人們的生活也構成了 極大的威脅。然而根據上述分析，目前InSAR技術難以廣泛地用于監測由地下流體運動引 起的三維地表形變。
【發明內容】
[0005] 本發明的目的在于，克服現有InSAR技術在監測地下流體運動引起的地表形變的 不足和局限性，提供一種利用InSAR斜距向形變測量值監測三維地表形變，并同步反演出 地下流體體積變化量的方法。
[0006] 一種利用InSAR反演地下流體體積變化和三維地表形變的方法，包括以下步驟：
[0007] 步驟1 :利用D-InSAR或多時相InSAR技術獲取待監測的受地下流體運動影響的 地區的地表在雷達視線方向的形變測量值，并對其進行地理編碼；
[0008] 所述雷達視線方向即為斜距向；
[0009] 步驟2 :根據SAR衛星成像幾何，按照以下公式構建InSAR斜距向形變測量值I (X1) 和三維地表形變之間的函數模型：
[0010]
[0011] 其中，X1表示地表觀測點，i = 1，2,…，M，共有M個地表觀測點；
[0012] S1(X1)、S2(X1)及S3(X 1)分另Ij為地表觀測點X1的東西、南北和垂直向地表形變在 InSAR 斜距向上的投影系數，S1(Xi) = -cos (α -3 π /2) sin Θ，S2(Xi) = -sin (α -3 π /2) sin θ，S3(Xi) = cos θ ; θ和α分別為雷達局部入射角和衛星飛行方向角；
[0013] Cl1(X1)為地表觀測點位置X1的形變量，1 = 1，2, 3分別代表東西、南北和垂直向上 的三個分量；
[0014] n (X1)為 InSAR 觀測誤差；
[0015] 【該值很小，通過最小二乘平差可以使其對InSAR斜距向形變測量值的影響最 小；】
[0016] 步驟3 :根據彈性半空間理論建立每個觀測點上的三維地表形變與地下流體的體 積變化之間的函數模型：
[0017]
[0018] 其中，G1 (X，y)為格林函數，
V為泊松比，S為地下塊源y 到地表觀測點X之間的距離，共有N個塊源(X)和P1 (y)分別為地表觀測點X和塊源y 的三維空間位置，
為地下塊源y到地 面觀測點X之間的距離；Dv(y)為地下半空間體積V中塊源y的流體體積變化量；Vy表示單 位塊源的體積，由單位塊源的尺寸和地下流體的厚度計算獲得；ε Jx1)表示模型與真實形 變之間的殘差；
[0019] 步驟4 :利用SAR衛星成像幾何將步驟3得到的模型轉化為地表每個觀測點上的 InSAR斜距向形變測量值與地下流體的體積變化之間的函數關系：
[0020]
[0021] 步驟5 :令地下流體中的所有塊源的體積變化AV(y)相同，利用地下流體場的先 驗信息確定泊松比V和地下流體的層數，利用M個地表觀測點上的InSAR斜距向形變測量 值，通過非線性方法反演出塊源的體積變化AV(y)、地下流體的深度h和厚度t ;
[0022] 步驟6 :構建以所有地表觀測點的三維形變量和所有地下塊源的體積變化量為未 知參數的聯合模型：
[0023]
[0024] 其中Ω為4MX 1的觀測矩陣，由M個地理編碼后的InSAR斜距向形變測量值和3M 個虛擬觀測量組成：Ω = [I (X1)…I (χΜ) O O O ......... O O 0]τ;
[0025] Δ為4ΜΧ 1的殘差矩陣，由M個InSAR觀測誤差和3Μ個模型誤差組成：
[0026] A = [ n (X1) ..· η (χΜ) ε ! (X1) ε 2(Χι) ε 3(Xl) ......... E1(Xm) ε2(χΜ) ε 3 (χμ)]
[0027] Γ為（3Μ+Ν) X 1的待求參數矩陣，由M個地面觀測點上的三維形變量和N個地下 塊源的體積變化量組成：
[0028]
[0029] B為4ΜΧ (3Μ+Ν)的設計矩陣：
[0030]
[0031] 步驟7 :利用稀疏最小二乘算法求解聯合模型，計算所有觀測點上的三維形變量 和所有塊源的體積變化量，即得到整個監測地區的三維地表形變量Cl 1 (Xl)、d2(Xl)和d3(Xl) 以及所有地下流體塊源的體積變化量D v (yj)。
[0032] 所述地下塊源的數量少于或等于所述地表觀測點的數量。
[0033] 【以保證解算聯合模型時不出現秩虧情況。】
[0034] 所述泊松比V取值為0. 25,地下流體的層數為1-3層。
[0035] 有益效果
[0036] 本發明提供了一種利用InSAR反演地下流體體積變化和三維地表形變的方法，1) 利用InSAR技術獲取該地區的地理編碼后的斜距向地表形變場的測量值；2)利用SAR衛星 成像幾何建立InSAR斜距向形變測量值與三維地表形變之間的函數模型，并基于彈性半空 間理論建立三維地表形變與地下流體體積變化之間的函數模型；3)設定地下流體塊源體 積均一變化，利用InSAR斜距向形變測量值非線性反演出地下流體塊源的深度和厚度；4) 利用上述的兩個函數模型建立針對所有地面觀測點和地下流體塊源的聯合解算模型；最后 利用最小二乘方法對聯合模型進行解算，得到所有地面觀測點的三維形變量和所有地下流 體塊源的體積變化量。該方法實現簡單，不受區域的限制，且不依賴于任何其它大地測量數 據（如GPS)，對地下流體運動引起的三維地表形變而言是一種高精度、大范圍、低成本和切 實可行的監測方法。不僅突破目前InSAR難以獲得高精度三維地表形變場的技術瓶頸，積 極推動InSAR大地測量技術向實用化和市場化發展，而且還能同時獲得地下流體的體積變 化結果，對研究地球內部的地球物理過程具有重要的科學價值和指導意義。
【附圖說明】
[0037] 圖1是本發明所述方法的流程圖；
[0038] 圖2是地下流體體積變化引起地表形變的示意圖；
[0039] 圖3是SAR衛星的成像幾何圖；
[0040] 圖4是模擬的地下流體體積變化；
[0041] 圖5是由模擬的地下流體體積變化造成的地表形變，其中，（a)為東西向形變；(b) 為南北向形變；（c)為垂直向形變；（d)為含噪的