基于PS-InSAR監(jiān)測的石油開采地表下沉模型反演

王遠堅，姜 岳，MISA R，SROKA A，姜 巖

(1.山東科技大學測繪與空間信息學院，山東 青島 266590；2.中國礦業(yè)大學環(huán)境與測繪學院，江蘇 徐州 221116；3.波蘭科學院巖層力學研究所，克拉科夫 30-059)

石油是重要的基礎能源和工業(yè)原材料，關系著國民經(jīng)濟命脈和國家安全。石油開采導致的緩慢地表移動變形現(xiàn)象在全球廣泛分布，這不僅會破壞油井的生產(chǎn)設施，給石油開采造成巨大的損失，同時會對建筑物和周邊生態(tài)環(huán)境產(chǎn)生損害性影響[1-2]。某油田是一個具有多套生油層系、多種儲集類型、多種油氣藏的復式油氣區(qū)，地表移動變形現(xiàn)象在持續(xù)發(fā)展，因此必須采取相應手段對該油田進行監(jiān)測。傳統(tǒng)的監(jiān)測手段需要耗費大量人力、物力、財力，而合成孔徑雷達干涉測量技術(interferometric synthetic aperture radar，InSAR)是一種空間對地觀測技術，具有高空間分辨率、高精度、監(jiān)測范圍廣等優(yōu)勢，已經(jīng)被廣泛應用于地表形變監(jiān)測、災害監(jiān)測與預警、地球物理參數(shù)反演等工作中[3-4]。針對石油開采導致的地表下沉具有緩慢累積的特點，采用PS-InSAR(permanent scatters InSAR)技術進行監(jiān)測，該技術能夠?qū)﹂L時間以固定重訪周期穩(wěn)定、連續(xù)獲得的SAR影像進行時間序列處理，利用在長時間和空間基線下能夠保持高相干的點目標，克服時間、空間去相干對干涉信號的影響，并利用形變相位和大氣效應相位在時間序列上的不同時空特性將兩者加以分離，得到地表形變演變參數(shù)，這使該技術對于監(jiān)測緩慢地表形變具有較大的優(yōu)勢[5]。

目前關于油田開采導致的地表形變監(jiān)測研究較多[6-7]，但對于地表形變的進一步研究則相對較少。在油田地表形變模型研究方面有部分學者將火山形變反演中的Mogi模型、地震形變反演中的Okada

模型引入[8-9]，這些模型能夠大致反映油田區(qū)域的地表下沉趨勢，但對于地表累積下沉量較大區(qū)域及下沉邊界區(qū)域模型擬合結果存在較大誤差。因此，根據(jù)地表實際下沉情況，確定一個符合目標區(qū)域的地表移動和變形模型能夠為評估開采引起的地表移動變形的損害風險提供依據(jù)。

1 研究區(qū)概況

該油田是一個以重油開發(fā)為主的復雜斷塊油田，面積約200 km2，采油廠位置及數(shù)據(jù)處理范圍如圖1所示。

本文采用覆蓋研究區(qū)域的2017年3月至2019年11月的33景C波段升軌Sentinel-1A數(shù)據(jù)進行下沉監(jiān)測，數(shù)據(jù)具體情況見表1。數(shù)據(jù)的距離向分辨率為5 m，方位向分辨率為20 m，數(shù)據(jù)處理范圍如圖1所示，面積約為1 736 km2。

圖1 石油開采區(qū)域及數(shù)據(jù)處理區(qū)域Fig.1 Oil production area and data processing area

表1 Sentinel-1A數(shù)據(jù)參數(shù)表Table 1 Sentinel-1A data parameters

續(xù)表1

2 PS-InSAR技術基本原理及地表移動變形模型

2.1 PS-InSAR技術基本原理

PS-InSAR技術能夠克服時間、空間去相關對干涉信號的影響，并利用大氣效應相位和形變相位在時間序列上的不同時空特性將兩者加以分離，得到地表形變演變參數(shù)[10-11]。PS-InSAR技術的基本原理是假設有K+1幅SAR單視復數(shù)影像，經(jīng)配準、輻射定標、PS探測和干涉處理，并借助已知DEM進行差分干涉處理后得到K幅干涉和差分干涉圖、H個PS點以及各PS點在各差分干涉圖中的差分干涉相位集。在考慮地表形變、高程誤差、大氣影響及失相關的情況下，每個PS點在每幅差分干涉圖上的差分干涉相位組成，見式(1)[12-13]。

φdiff=φtopo+φdef+φorb+φatm+φnoise

(1)

通過設定振幅離差指數(shù)、相干系數(shù)、旁瓣閾值提取PS候選點，后對這些PS候選點進行噪聲計算，去除僅在個別干涉圖中相位穩(wěn)定或受到鄰近PS點影響而表現(xiàn)出PS點特征的象元，剩余的即為PS點。估計并剔除PS點相位的空間不相關誤差，并更新PS點相位；對更新后的PS點相位進行時空三維解纏；分析各干擾誤差相位在時間和空間上的特征差異，對解纏相位進行時、空濾波，分離并剔除各誤差相位，最后提取形變相位。PS-InSAR的解算結果是視線向的形變值，為了準確反映研究區(qū)域的下沉變化，需將視線向的解算成果轉(zhuǎn)換為垂直向的結果。該轉(zhuǎn)換可通過式(2)來進行[14]。

(2)

式中：d為下沉值；Δr為視線向形變值。

2.2 地表移動變形模型

為了研究開采對地表油井、建筑物及環(huán)境的影響，需要根據(jù)PS-InSAR監(jiān)測結果，模擬計算出地表下沉盆地數(shù)學模型，為計算地表移動與變形提供基礎。關于石油開采地表下沉計算有許多方法[15-17]，但是，部分方法計算比較復雜，本文結合研究區(qū)域的地表下沉分布特征，將地表下沉視作橢圓形彈性薄板撓曲，構建基于橢圓薄板模型的地表下沉模型[18-20]。如圖2所示，圖中q為薄板上方的橫向載荷；w為地表下沉量；a為橢圓下沉盆地長軸半長，即橢圓板長半軸；b為橢圓下沉盆地短軸半長，即橢圓板短半軸。

圖2 彈性橢圓薄板的坐標和位移示意圖Fig.2 Coordinates and displacements ofelastic elliptic thin plates

彈性橢圓薄板的撓度計算模型邊界條件表達式見式(3)[21]。

(3)

基于彈性橢圓薄板的撓度計算模型邊界條件及下沉區(qū)域空間形態(tài)，可建立監(jiān)測區(qū)域地表下沉模型表達式見式(4)。

(4)

式中，k為常數(shù)，即為需要反演的值。

根據(jù)地表移動變形理論[22-23]，由式(4)可以求出，地表傾斜i(x,y,φ),曲率K(x,y,φ),水平移動U(x,y,φ),水平變形ε(x,y,φ)。

沿φ方向的傾斜，見式(5)。

(5)

沿φ方向的曲率，見式(6)。

(6)

沿φ方向的水平移動，見式(7)。

U(x,y,φ)=Bri(x,y,φ)

(7)

沿φ方向的水平變形，見式(8)。

ε(x,y,φ)=BrK(x,y,φ)

(8)

式中：B為水平移動系數(shù)；r為主要影響半徑；φ為沿x軸正向逆時針計算到指定方向的角值。

3 監(jiān)測結果及時空分析

以2018年6月17日獲取的影像為主影像，其他32景影像為輔影像，運用PS-InSAR技術對獲取的影像進行處理，獲取了2017年3月至2019年11月的監(jiān)測成果，在數(shù)據(jù)處理范圍內(nèi)共提取到401 728個PS點，采油廠監(jiān)測范圍內(nèi)共提取到17 796個PS點，如圖3所示。由于該區(qū)域內(nèi)存在大面積的蘆葦蕩，造成眾多區(qū)域失相干，未提取到PS點。從圖3中

可以看到，該采油田范圍內(nèi)存在明顯的下沉漏斗，下沉現(xiàn)象較嚴重，下沉區(qū)域面積達29.5 km2，最大下沉速率為-210.9 mm/a，最大累積下沉量達-585.6 mm，這會對地表的建構筑物和公路等設施產(chǎn)生影響，甚至造成破壞。

圖3 2017年3月至2019年11月下沉速率圖Fig.3 Subsidence rate map from March 2017to November 2019

為了更深入準確研究該區(qū)域的下沉狀況，按季度刻畫出該區(qū)域下沉時間序列圖，如圖4所示。由圖4可知，從宏觀上直觀反映出該采油廠在時間和空間上的下沉演變，該區(qū)域下沉是一個緩慢發(fā)展的過程，隨著時間的推移，下沉區(qū)域累積下沉量不斷增大，下沉范圍也在持續(xù)擴大。

圖4 下沉時間序列圖Fig.4 Subsidence time series

根據(jù)監(jiān)測結果可知該下沉區(qū)域形狀近似橢圓形(圖5)，對PS-InSAR監(jiān)測的矢量成果插值獲取下沉區(qū)域沿圖5中A線和B線的方向的下沉剖面線，分別如圖6和圖7所示，A線和B線分別視作橢圓形下沉區(qū)域的長軸和短軸。根據(jù)A線和B兩條剖面線數(shù)據(jù)進一步反映的下沉情況可以看到A線和B線的下沉量關于剖面線上的最大下沉量處對稱，可以判斷下沉盆地的中心位于A線和B線的交點周圍。

研究區(qū)域位于遼河盆地，自更新世以來，整個盆地以下沉為主，呈北東向展布，與凹陷區(qū)的構造較為一致，在遼陽、鞍山一帶下沉速度為-3 mm/a，在臨近渤海的盤山，最大下沉速度為-5 mm/a[24]。根據(jù)國家重點基礎理論研究(攀登)項目研究報告[25]，中國學者對中國東南沿海地區(qū)地殼與陸地垂直運動的研究，其中遼寧營口地區(qū)遼河平原地殼與陸地平均下沉速度為3 mm/a左右。根據(jù)上述文獻，研究區(qū)域陸地自然下沉量為-3～-5 mm/a，所以3年累計最大下沉為-15 mm，因此在PS-InSAR監(jiān)測的下沉數(shù)據(jù)中，包含著石油開采、地下水開采及陸地自然下沉等復雜影響因素。但在約3年的監(jiān)測時間跨度內(nèi)，PS-InSAR技術獲得的地表最大下沉超過-580 mm，自然下沉量僅占2.6%，自然下沉相對開采引起的下沉量較小，且由圖6和圖7可以看出，地表下沉是連續(xù)的，沒有出現(xiàn)非連續(xù)臺階狀下沉，說明構造對地表下沉的影響沒有顯現(xiàn)出來，因無法區(qū)分開采下沉與自然下沉，因此，把地表下沉監(jiān)測結果視作開采引起的，忽略了自然下沉等次要因素影響。

圖5 下沉盆地區(qū)域及剖面線位置示意圖Fig.5 Schematic diagram of subsidence basin areaand section line position

圖6 下沉盆地長軸剖面線Fig.6 Long axis section line of subsidence basin

4 地面移動變形模型反演及下沉

本文中以下沉量為10 mm的外部邊界區(qū)域作為下沉盆地邊界以求取橢圓下沉盆地長軸和短軸的長度，以A、B兩條剖面線的數(shù)據(jù)為已知量對地表下沉模型進行反演以確定模型參數(shù)。根據(jù)下沉盆地剖面線的數(shù)據(jù)可以確定橢圓下沉盆地的長半軸a和短半軸b的長度分別為3 350 m、2 050 m，最終模型反演，見式(9)。

(9)

式中，w(x,y)為累積下沉量，m。

通過利用反演得到的模型對下沉盆地兩條剖面線位置的累積下沉量進行正演，比較模型擬合精度，結果如圖8和圖9所示。

圖7 下沉盆地短軸剖面線Fig.7 Short axis section line of subsidence basin

圖8 下沉盆地長軸監(jiān)測結果與正演結果比較Fig.8 Comparison of long axis monitoring resultsand forward results in subsidence basin

圖9 下沉盆地短軸監(jiān)測結果與正演結果比較Fig.9 Comparison of short axis monitoring resultsand forward results in subsidence basin

該模型的相關系數(shù)為0.996 5，正演結果的均方根誤差為23.4 mm，該模型對該下沉盆地累積下沉量的正演結果具有較高的精度。

基于該下沉模型及地表移動變形相關理論可分別計算得到該下沉區(qū)域的傾斜i(x,y,φ)、曲率K(x,y,φ)、水平移動U(x,y,φ)、水平變形ε(x,y,φ)表達式分別為式(10)、式(11)、式(12)和式(13)。

(10)

(11)

結合該區(qū)域地質(zhì)情況和監(jiān)測數(shù)據(jù)，取B=0.25，根據(jù)主要影響半徑定義計算得到r=1 985 m。

U(x,y,φ)=Bri(x,y,φ)=

(12)

ε(x,y,φ)=BrK(x,y,φ)=

(13)

基于上述得到的該采油廠區(qū)域地表移動變形模型可計算得到橢圓下沉盆地各點的傾斜、曲率、水平移動和水平變形值。下沉盆地長軸A上各點沿φ=0°方向和短軸B上各點沿φ=90°的傾斜和水平變形分別如圖10和圖11所示。

圖10 A軸傾斜和水平變形Fig.10 A-axis tilt and horizontal deformation

圖11 B軸傾斜和水平變形Fig.11 B-axis tilt and horizontal deformation

5 結 論

1) 本文采用PS-InSAR技術對某油田進行了下沉監(jiān)測，通過對監(jiān)測結果分析發(fā)現(xiàn)：該油田的下沉演變是一個由小到大，不斷增加的過程。在監(jiān)測時間段內(nèi)，該采油廠監(jiān)測范圍內(nèi)共提取到17 796個PS點，下沉區(qū)域面積達29.5 km2，PS-InSAR技術為對該區(qū)域下沉過程監(jiān)測提供了高效的技術方法。

2) 應用橢圓彈性薄板撓度模型，建立該區(qū)域地表移動和變形模型，并結合實測結果反演得到模型參數(shù)，模型計算均方根誤差為23.4 mm，模型的相關系數(shù)為0.996 5，具有較高的擬合精度，所建立的計算模型可以作為地表下沉變形評價的基準。

3) 截至2019年11月，監(jiān)測區(qū)域地表最大下沉585.6 mm，最大下沉速率為-210.9 mm/a。根據(jù)反演模型計算，地表最大傾斜0.47 mm/m，最大水平拉伸為0.42 mm/m，最大壓縮變形分別為-0.34 mm/m。計算結果為評價該區(qū)域開采損害風險評估提供依據(jù)。

4) 目前關于石油開采引起的地表下沉規(guī)律研究的還不充分，應通過對地表下沉監(jiān)測獲取更多的地表下沉信息與儲層變化的相互關系，建立更加精確的地表下沉與變形計算模型，為評價開采對地面設施與環(huán)境的影響提供依據(jù)。