MATLAB在深層構(gòu)造預(yù)測(cè)中的應(yīng)用
——以房山張坊水源地為例

孫少游，馬靜晨，劉全國(guó)，楊 楠，潘 懋

（1. 北京市地質(zhì)工程勘察院，北京 100048； 2. 北京大學(xué)，北京 100087；3. University of British Columbia, Vancouver, Canada V6P4Y 7）

MATLAB在深層構(gòu)造預(yù)測(cè)中的應(yīng)用
——以房山張坊水源地為例

孫少游1，馬靜晨1，劉全國(guó)1，楊 楠2,3，潘 懋2

（1. 北京市地質(zhì)工程勘察院，北京 100048； 2. 北京大學(xué)，北京 100087；3. University of British Columbia, Vancouver, Canada V6P4Y 7）

利用MATLAB強(qiáng)大的分析能力，基于北京房山巖溶區(qū)淺層節(jié)理裂隙數(shù)據(jù)的進(jìn)行性半自動(dòng)化的統(tǒng)計(jì)與分析，繪制淺層節(jié)理裂隙玫瑰花圖，并提取深層構(gòu)造的走向線，最后與深層斷層、向斜與背斜等構(gòu)造進(jìn)行擬合。擬合結(jié)果顯示，北京房山巖溶地區(qū)淺層節(jié)理裂隙的統(tǒng)計(jì)走向與深層構(gòu)造走向基本一致，北東向與北西向兩條主構(gòu)造方向的擬合誤差分別為6.4%和2.1%。這一誤差遠(yuǎn)小于設(shè)置的15%的誤差容限，表明此代碼在利用預(yù)測(cè)深層構(gòu)造的應(yīng)用中具有一定的正確性與潛力，能夠?yàn)橐巴鈹?shù)據(jù)的分析統(tǒng)計(jì)提供簡(jiǎn)潔有效的幫助。但方法的有效性，仍需要通過更多的實(shí)例進(jìn)行驗(yàn)證。

MATLAB; 巖溶水; 深層構(gòu)造預(yù)測(cè); 房山

0 前言

北京地區(qū)水資源匱乏，其地下水資源在城市供水中占據(jù)十分重要的位置。巖溶水作為北京地下水的重要組成部分因賦存條件復(fù)雜、勘探難度大、開發(fā)利用成本高，勘察評(píng)價(jià)研究起步較晚（Pyne，1995）。根據(jù)《北京山前淺埋區(qū)巖溶水資源勘察評(píng)價(jià)》項(xiàng)目初步成果，北京市巖溶水天然補(bǔ)給資源量可達(dá)5.83×108m3/a（韓再生， 2003；Herman，2002）。北京地區(qū)碳酸鹽巖分布面積約4900km2，其中山區(qū)碳酸鹽巖面積約2900km2，占山區(qū)總面積的28%；平原隱伏區(qū)約2000km2，以?shī)W陶系灰?guī)r、寒武系灰?guī)r和薊縣系白云巖為主（北京市水文地質(zhì)工程地質(zhì)大隊(duì)，2001(a)；北京市水文地質(zhì)工程地質(zhì)大隊(duì)，2001(b)；British Geology Survey, 2015；U.S.Department of the Interior et al，2008）。

北京市房山區(qū)張坊水源地是北京市建成的第一個(gè)大型應(yīng)急供水工程，地處房山西南山前地區(qū)的霞云嶺—龍門臺(tái)復(fù)向斜水文地質(zhì)單元，受下游斷裂構(gòu)造控制，形成了相對(duì)獨(dú)立的巖溶水系統(tǒng)。區(qū)內(nèi)巖溶裂隙發(fā)育、含水層厚度大，地下水易采易補(bǔ)。大石河從水源地穿過，在豐水期河水直接補(bǔ)給地下水，使該地區(qū)具備地表水和地下水聯(lián)合調(diào)蓄供水條件。北京市房山地區(qū)多年勘查已經(jīng)積累了一定的基礎(chǔ)資料（Bear et al，1987），但由于巖溶區(qū)構(gòu)造裂隙網(wǎng)絡(luò)空間分布的復(fù)雜性、巖溶水分布的不均一性以及巖溶區(qū)地下水與地表水之間相互轉(zhuǎn)化的復(fù)雜性，目前對(duì)北京房山地區(qū)仍缺乏綜合系統(tǒng)的研究，且尚未制定巖溶水資源合理開發(fā)保護(hù)方案、法規(guī)或條例。

Matlab是矩陣（Matrix）和實(shí)驗(yàn)室（Laboratory）的組合，不僅自帶有強(qiáng)大的函數(shù)庫(kù)，實(shí)現(xiàn)各種運(yùn)算，而且還允許用戶進(jìn)行自編程，實(shí)現(xiàn)特定的功能。Matlab以其強(qiáng)大的數(shù)值運(yùn)算與處理能力，在水資源規(guī)劃、水質(zhì)評(píng)價(jià)、土壤重金屬污染空間分布（李鵬等，2014）等方面得到了廣泛的應(yīng)用。本次利用MATLAB代碼實(shí)現(xiàn)的功能：包括提取三組深層節(jié)理數(shù)據(jù)；提取節(jié)理裂隙玫瑰花圖指示的淺層節(jié)理裂隙數(shù)據(jù)；擬合淺層與深層節(jié)理裂隙數(shù)據(jù)。本文研究的主要原因在于：程序化的方法能夠加快數(shù)據(jù)分析，尤其是大量數(shù)據(jù)存在的情況下每一個(gè)微小的錯(cuò)誤都會(huì)導(dǎo)致結(jié)果的偏差，而重復(fù)分析十分耗費(fèi)時(shí)間；其次，程序化的方法使分析過程更簡(jiǎn)潔，分析結(jié)果的視覺化效果更直觀更完善；最后，本文的代碼可用于其他區(qū)域的檢測(cè)，最終用于深層構(gòu)造的預(yù)測(cè)。這一代碼的實(shí)現(xiàn)將會(huì)加快分析和統(tǒng)計(jì)節(jié)理裂隙的進(jìn)度并減小繪圖的工作量，同時(shí)能夠完善視覺效果，最終為科學(xué)有效地管理和保護(hù)巖溶水資源提供有效的幫助。

1 研究區(qū)域及背景

房山區(qū)位于北京市西南，其地理坐標(biāo)為東經(jīng)115°25′—116°15′，北緯39°30′—39°55′。西、南部與河北省接壤，北接門頭溝區(qū)，東接大興區(qū)，東北部與豐臺(tái)區(qū)相接。區(qū)內(nèi)主要道路有房易路、周張路，京石高速公路從研究區(qū)的東部通過，成為房山區(qū)與北京城聯(lián)系的重要通道（Coppana et al， 1985；U.S. Geological Survey，2016）。

房山區(qū)地處北京西部山區(qū)和永定河沖積平原交匯處（Bouwer，2002），地貌主要為山地和平原兩種類型，山地面積廣大，位于房山區(qū)的西北部，約占總面積的三分之二，如圖1所示。房山區(qū)屬于大清河水系，境內(nèi)地勢(shì)西北高，東南低。以鳳凰山—上寺嶺—霞云嶺分水嶺為界，北部為大石河，南部為拒馬河。北部大石河切割較淺，由霞云嶺向北經(jīng)薊縣系地層進(jìn)入寒武系、奧陶系灰?guī)r地層之中，沿岸現(xiàn)巖溶溶蝕作用形成的洼地，并受構(gòu)造作用影響形成佛子莊泉、河北泉等多處泉。南部低山丘陵區(qū)廣泛分布薊縣系霧迷山組碳酸鹽巖，受拒馬河切割影響，形成兩岸陡峭的峰叢、峰林等巖溶地貌，且?guī)r溶洞穴較發(fā)育，形成石花洞、云水洞等重點(diǎn)景觀巖溶洞穴。山區(qū)巖溶裂隙發(fā)育，地下水入滲補(bǔ)給條件好。

圖1 北京市房山區(qū)地貌圖Fig.1 Terrain of Fangshan District

2 方法與材料

2012年4月—2013年4月，野外調(diào)查小組通過對(duì)北京市房山區(qū)地表各時(shí)期地層野外出露的節(jié)理進(jìn)行詳細(xì)的觀察與測(cè)量，共采集來自295個(gè)觀測(cè)點(diǎn)的2095條節(jié)理裂隙數(shù)據(jù)，其中忽略延伸小于1m的節(jié)理，共計(jì)有效節(jié)理2036條，裂隙數(shù)據(jù)采樣點(diǎn)分布如圖2所示。

圖2 房山區(qū)節(jié)理裂隙數(shù)據(jù)分布點(diǎn)Fig.2 Fangshan Distribution of data collected points

3 結(jié)果與分析

將2036條有效的節(jié)理裂隙數(shù)據(jù)根據(jù)所屬地層劃分為32個(gè)組，并根據(jù)其走向數(shù)據(jù)在MATLAB中繪制玫瑰花圖，這些玫瑰花圖指示的總體方向如圖3所示。但基于所有的走向分布，我們可以利用MATLAB代碼提取出主要的走向，總體來說，北京房山地區(qū)以低應(yīng)變脆性破裂變形為主，形成X形共軛節(jié)理。節(jié)理面平直、傾角較陡、延伸中等。利用獲得的節(jié)理數(shù)據(jù)繪制節(jié)理走向玫瑰花圖可見，房山區(qū)節(jié)理發(fā)育指示出的優(yōu)勢(shì)方位如圖4所示。

在獲取了淺層節(jié)理裂隙所指示的主要走向之后，利用MATLAB去擬合淺層與深層構(gòu)造的基本走向。首先導(dǎo)入基本構(gòu)造圖（圖5）并提取基本構(gòu)造圖中所顯示的深層構(gòu)造走向。隨后，MATLAB代碼將進(jìn)行淺層節(jié)理裂隙與深層構(gòu)造走向的擬合，擬合的基本準(zhǔn)則為：在15%的誤差范圍之內(nèi)尋找與淺層節(jié)理裂隙提取的主要走向吻合的深層構(gòu)造。結(jié)果顯示，兩條淺層節(jié)理裂隙的主方向與深層構(gòu)造的擬合誤差分別為6.4%和2.1%。

圖3 玫瑰花圖指示的走向分布Fig.3 All strike directions in rose diagram

圖4 兩條主要的走向Fig.4 Two major strike directions

圖5 房山區(qū)節(jié)理裂隙與區(qū)域斷裂構(gòu)造關(guān)系示意圖Fig.5 The base map for deep structure

較小范圍的誤差，表明淺層節(jié)理裂隙所反映的構(gòu)造走向與深層構(gòu)造（如斷層、背斜和向斜）的實(shí)際走向基本吻合。因?yàn)閰^(qū)域地質(zhì)構(gòu)造是多期次不同構(gòu)造運(yùn)動(dòng)的結(jié)果，而這些大型的深層構(gòu)造均為地質(zhì)歷史時(shí)期多起構(gòu)造應(yīng)力場(chǎng)作用的結(jié)果，因此會(huì)疊加反映多期構(gòu)造運(yùn)動(dòng)的應(yīng)力場(chǎng)方向，但是簡(jiǎn)而言之，無論在淺層節(jié)理裂隙還是深層構(gòu)造中，區(qū)域主應(yīng)力與構(gòu)造走向垂直。

構(gòu)造應(yīng)力場(chǎng)，就是在一個(gè)空間范圍內(nèi)構(gòu)造應(yīng)力的分布。一個(gè)地區(qū)的構(gòu)造應(yīng)力可以用主應(yīng)力方向來表示，即剪應(yīng)力等于零的方向。若要能夠產(chǎn)生構(gòu)造變形，一般需要多向應(yīng)力，其中3個(gè)主應(yīng)力方向上的應(yīng)力值一般應(yīng)該是不相等的，其中必有一個(gè)為最大主壓應(yīng)力方向（σ1），一個(gè)為中間主應(yīng)力方向，一個(gè)為最小主應(yīng)力方向。當(dāng)3個(gè)主應(yīng)力值相等時(shí)，就處于靜巖壓力狀態(tài)，此時(shí)則不會(huì)發(fā)生構(gòu)造變形。

4 分析與后續(xù)研究

MATLAB具有強(qiáng)大的數(shù)值運(yùn)算與處理能力，本研究雖然整個(gè)分析是在MATLAB中開發(fā)的，但整個(gè)過程并非真正的自動(dòng)化，在某些步驟仍需要人機(jī)交互進(jìn)行。同時(shí)，獲得主要構(gòu)造走向，15%誤差容限和獲取玫瑰圖過程的平均方法都對(duì)結(jié)果具有重要影響?；谏鲜鲂畔?，來自該文章的結(jié)果需要通過更多應(yīng)用程序和實(shí)踐進(jìn)行驗(yàn)證：

（1）誤差的疊加可能大于我們的預(yù)期，因?yàn)檎`差產(chǎn)生于過程的多個(gè)節(jié)點(diǎn)；

（2）方法的有效性需要通過更多的應(yīng)用和實(shí)踐來驗(yàn)證；

（3）MATLAB中的編碼方法可能不是最簡(jiǎn)潔的；

（4）鑒于作者在MATLAB上的知識(shí)有限，MATLAB編碼可能不如我們預(yù)期的有效；

（5）為MATLAB在深層構(gòu)造預(yù)測(cè)中的應(yīng)用提供了范例，也為后續(xù)的研究進(jìn)行了有益的嘗試。

Bear J, Berkowitz B, 1987. Groundwater flow and pollution in fractured rock aquifers[J]. Developments In hydraulic Engineering, (4)： 175-238.

Bouwer H, 2002. Artificial recharge of groundwater：Hydrogeology and engineering[J]. Hydrogeology Journal，10(1)：121-142.

British Geology Survey. Elucidating the hydrogeological issues associated with aquifer storage and recovery in the UK[EB/OL].(2015-11-23)[2016-10-22]. http：//nora.nerc.ac.uk/512266/.

Campana ME，Mahin DA, 1985. Model-derived estimates of groundwater mean ages，recharge rates，effective porosities and storage in a limestone aquifer[J]．Journal of Hydrology，76(3-4)：247-264．

Herman Bouwer, 2002. Artificial recharge of groundwater：hydrogeology and engineering[J]. Hydrogeology Journal, 10(1)： 121-142.

Pyne RDG,1995. Groundwater recharge and wells： A guide to Aquifer Storage Recovery[M]. Florida：Lewis Publishing.

U.S.Department of the Interior,U.S.Geological Survey. Ground-water recharge in the arid and semiarid southwestern United States∥USGS Professional Paper 1703： 375[EB/OL]. (2008-8-8)[2016-10-22]. https：//pubs.usgs.gov/pp/pp1703/.

U.S. Geological Survey.U.S. Geological Survey artificial recharge workshop proceedings,[USGS Open-File Report 02-89] [EB/OL] (2016-11-28)[2016-12-5]. https：//water.usgs.gov/ogw/pubs/ofr0289/.

北京市水文地質(zhì)工程地質(zhì)大隊(duì)，2001. 北京地區(qū)巖溶水資源研究與開發(fā)利用工程可行性研究報(bào)告[R].

北京市水文地質(zhì)工程地質(zhì)大隊(duì)，2001. 北京地區(qū)巖溶水資源研究與開發(fā)利用工程工作總體設(shè)計(jì)[R].

韓再生，2003. 為可持續(xù)利用而管理含水層補(bǔ)給--第四屆國(guó)際地下水人工補(bǔ)給會(huì)議綜述[J]. 地質(zhì)通報(bào)，22(2)：142-143.

李鵬，陳永當(dāng)，葉宏武，等． 基于Matlab 的城市表層土壤重金屬污染空間分布的研究[J]．環(huán)境保護(hù)科學(xué)，2014，40(2)：51-54.

Application of MATLAB in Deep Structure Prediction——Application in Fangshan, Beijing

SUN Shaoyou1, YANG Nan2,3, PAN Mao2

(1. Beijing Institute of Geological & Prospecting Engineering, Beijing100048; 2. Peking University, Beijing 100871; 3. University of British Columbia, Vancouver, Canada, V6P 4Y7)

Taking advantage of the robust analytic ability of MATLAB, the analysis of shallow fractured joints in Fangshan karst area of Beijing were analyzed by a semi-automatic coding to draw the shallow joints, and extracted the deep structure trending lines. The results showed that the statistical trending of shallow joints in the Fangshan karst area of Beijing is basically the same as that of the deep tectonics, and the fi tting errors of the two main tectonic directions are 6.4% and 2.1%, respectively. This error is far less than the 15% error margin of the set, indicating that the code has some correctness and potential in the application of predicting deep structure, which can provide a simple and effective help for the analysis and statistics of the fi eld data, but the method’s effectiveness still need to be verif i ed through more examples.

MATLAB; Karst Water; Deep Structure Prediction; Fangshan District

P542.3

A

1007-1903（2017）01-0117-04

10.3969/j.issn.1007-1903.2017.01.021

孫少游（1983- ），男，工程師，主要研究方向：地下水資源。E-mail：418171619@qq.com