基于SPAC法探測松科二井深層地?zé)醿?chǔ)水構(gòu)造

應(yīng)恒成, 李洪強(qiáng), 張玉敏, 靳中原, 張家鐸, 符 偉, 高 磊,安棟召, 潘宗棟, 汪 偉, 侯賀晟*

1)中國地質(zhì)科學(xué)院, 北京 100037; 2)黑龍江省自然資源調(diào)查院, 黑龍江 哈爾濱 150036

中國作為世界上最大的碳排放國家與工業(yè)門類最為齊全的國家, 正處在城鎮(zhèn)化快速發(fā)展階段,面臨著經(jīng)濟(jì)轉(zhuǎn)型、環(huán)境保護(hù)、應(yīng)對氣候變化等多重挑戰(zhàn)(張希良等, 2022)。習(xí)近平總書記于 2020年 9月22日宣布, 中國將力爭于2030年前完成碳達(dá)峰,爭取在2060年前實(shí)現(xiàn)碳中和。而地?zé)崮苁且环N綠色低碳、可循環(huán)利用的可再生能源。在清潔能源(太陽能、風(fēng)能、生物質(zhì)能、地?zé)崮艿?大家族中, 地?zé)豳Y源作為可再生能源, 是清潔能源家族中的重要成員,具有豐富的資源潛力, 利用系數(shù)高, 生命周期長,二氧化碳排放低, 潛在的低成本優(yōu)勢。相比于水力發(fā)電, 地?zé)崮馨l(fā)電能經(jīng)濟(jì)地建立相對較小的發(fā)電機(jī)組。并且, 我國擁有豐富的地?zé)豳Y源, 初步估算在3000 m的深度范圍內(nèi), 地?zé)豳Y源可利用的熱能相當(dāng)于25億t標(biāo)準(zhǔn)煤。它在未來能源供應(yīng)與節(jié)能減排方面有著巨大潛力(王貴玲等, 2000; 周總瑛等, 2015;田寶卿等, 2020)。隨著國家提出碳中和目標(biāo), 日后地?zé)崮茉丛谖覈偰茉吹恼急葧?huì)不斷增加。目前,鉆探和地球物理勘探技術(shù)是常用的地?zé)峥碧椒椒?其中地球物理勘探的主要作用是圈定含水破碎帶和熱儲(chǔ)水層的區(qū)域分布, 是勘查地?zé)峥瞻讌^(qū)的有效手段(李安寧和呂金波, 2001; 鄭東明等, 2005)。

在城市地區(qū), 由于環(huán)境條件復(fù)雜, 常規(guī)的地質(zhì)鉆探和傳統(tǒng)的物探手段在勘查中受到很大局限和約束, 如鉆探工程, 施工成本較大, 并受地下的各類管線與人防工程及現(xiàn)有建筑物等限制, 有些鉆孔無法實(shí)施, 會(huì)出現(xiàn)地質(zhì)信息盲點(diǎn); 傳統(tǒng)物探手段則會(huì)受到噪音、電磁等各種干擾, 導(dǎo)致數(shù)據(jù)失真, 影響探測準(zhǔn)確性。同時(shí), 含水破碎帶和熱儲(chǔ)水層有時(shí)埋深較大, 也加大了工作難度, 僅依靠鉆探及傳統(tǒng)物探手段很難達(dá)到勘查目的。因此, 嘗試采用新的物探手段與多種傳統(tǒng)物探方法聯(lián)合勘探, 以提高地?zé)峥辈樘綔y精度這一觀點(diǎn)具有重要意義(付微等,2012)。

微動(dòng)勘探因?yàn)槠渚哂袩o需人工源、方便快捷、對環(huán)境影響小的優(yōu)勢, 非常合適于城市或農(nóng)村進(jìn)行地?zé)峄驕\層速度異?？碧?。Xu(2012)通過在微動(dòng)速度剖面的低速異常體上打地?zé)峋姆绞? 驗(yàn)證了低速異常體的構(gòu)造意義, 得到了微動(dòng)探測得到的低速異常體很大概率為斷層/破碎帶、含水巖層或石灰?guī)r溶洞這一結(jié)論。本文將在這一結(jié)論與微動(dòng)勘探方法的基礎(chǔ)上, 結(jié)合在黑龍江省安達(dá)市松科二井周邊的地?zé)豳Y源調(diào)查實(shí)踐, 展示微動(dòng)探測的實(shí)際應(yīng)用效果。

1 微動(dòng)原理與探測方法

1.1 微動(dòng)定義

“微動(dòng)”是地球表面無時(shí)無刻都存在的微弱震動(dòng)(Okada et al., 2004)。微動(dòng)是一種在時(shí)間、空間域都極不規(guī)則的震動(dòng)現(xiàn)象, 其振動(dòng)幅值大約為10-4～10-2mm, 加速度只有 1～100 μm/s2, 頻率范圍主要集中在0.1～50 Hz范圍。既可以由諸如氣壓、風(fēng)速、海浪、潮汐變化等自然現(xiàn)象產(chǎn)生, 也源自車輛行駛、機(jī)器運(yùn)行等過程, 以及人們的日常生活、生產(chǎn)活動(dòng)等行為(徐佩芬等, 2013a)。微動(dòng)是一種由體波(P波和S波)和面波(瑞雷波和勒夫波)組成的復(fù)雜震動(dòng), 且面波能量占信號總能量的 70%以上(Toksoz and Lacoss, 1968)。

1.2 微動(dòng)探測方法

要獲取瑞雷波中的信息, 首先要從微動(dòng)信號中提取出面波(瑞雷波)信號, 在提取面波環(huán)節(jié), 一般是以微動(dòng)信號的垂直分量作為資料, 因?yàn)槔辗虿]有垂直分量, 由觀測垂直分量得到的面波僅為瑞雷波, 不會(huì)混入勒夫波, 并且對于垂直分量的觀測也比較容易做到(王振東, 1990)。

目前, 常用的提取面波的方法有兩種, 一是Aki(1957)提出的空間自相關(guān)法(Spatial Autocorrelation Method, 簡稱 SPAC 法); 二是頻率-波數(shù)法(Frequency-Wavenumber Spectral Method, 簡稱F-K法)(Ohori et al., 2002)。本文使用的就是空間自相關(guān)法。

1.3 基于SPAC法的S波速度結(jié)構(gòu)的估算

SPAC法基于兩點(diǎn)基本假設(shè): (1)微動(dòng)在時(shí)空上符合平穩(wěn)隨機(jī)過程; (2)微動(dòng)所包含的各種成分的波中, 基階面波占優(yōu)。其基本原理如下:

圓心臺(tái)陣圓心O(0, 0)和半徑為r的圓周上A(r,θ)點(diǎn)的微動(dòng)信號垂直分量表達(dá)式為(Okada, 2006):

該式表明, 波的相速度 c(ω)可以通過第一類零階貝塞爾函數(shù)和空間自相關(guān)系數(shù)求得。

下一步將圓心和圓周上兩點(diǎn)的微動(dòng)信號通過傅里葉變換轉(zhuǎn)化到頻率域的表示, 得到空間自相關(guān)系數(shù)在頻率域下的表示, 將其代入式(6)可得相速度c(ω), 獲得相速度頻散曲線(徐佩芬等, 2013b)。

獲得頻散曲線之后, 利用全局尋優(yōu)算法——分歧型遺傳算法(fGA)(Cho et al., 1999)反演地下S波速度結(jié)構(gòu)。

2 應(yīng)用實(shí)例

2.1 測區(qū)概況

本次測區(qū)位于黑龍江安達(dá)市南部松科二井及周邊地區(qū)(如圖1)。松科二井完鉆井深7018 m, 是亞洲國家實(shí)施的最深的大陸科學(xué)鉆井, 也是國際大陸科學(xué)鉆探計(jì)劃(ICDP)成立 22年來實(shí)施的最深鉆井(侯賀晟等, 2018)。松科二井位于中國東北的松遼盆地, 松遼盆地是一個(gè)整體形狀呈現(xiàn)為菱形的中—新生代沉積盆地, 是沉積盆地型地?zé)豳Y源在我國的主要分布區(qū)之一(王貴玲等, 2017)。松遼盆地沉積了7000 m厚的中—新生代陸源碎屑巖, 是良好的熱儲(chǔ)層。具體到地層上, 姚家組、青二段、青三段、泉四段及泉三段都是可能的儲(chǔ)層。而青山口組的底部青一段、嫩江組的底部嫩一、二段為盆地分別處于最大沉降期時(shí)所形成的半深湖-深湖相沉積, 沉積環(huán)境穩(wěn)定, 泥巖純且厚度大, 分布面積廣, 封蓋性能好, 是非常好的區(qū)域性隔水層(朱煥來, 2011)。在熱源方面, 根據(jù)中國大陸地區(qū)大地?zé)崃鲾?shù)據(jù)匯編(第四版)顯示, 松遼盆地在統(tǒng)計(jì)上屬于我國的高熱流帶。地表熱流數(shù)據(jù)反映了巖石圈淺部熱狀態(tài)和能量均衡信息(羅凡等, 2022)。在具體的熱流數(shù)據(jù)上,中國大陸地區(qū)(含渤海海域)平均熱流值為(61.5±13.9) mW·m-2, 而松遼盆地的平均熱流值高達(dá)(70.9±14.4) mW·m-2(姜光政等, 2016)。同時(shí), 松遼盆地中部位于我國地溫梯度最高值分布區(qū)域(張薇等, 2019), 表明該區(qū)域具有很高的地溫背景。松科二井的多次井溫測井結(jié)果顯示, 盆地具有較大的地溫梯度, 大地?zé)崃髦颠_(dá)到84 mW·m-2, 2000 m以深地溫梯度達(dá) 3.7 ℃/100 m, 具有良好的深部熱傳導(dǎo)和淺部熱聚集構(gòu)造機(jī)制, 地?zé)豳Y源賦存條件較好(侯賀晟等, 2018)。

圖1 微動(dòng)測區(qū)以及探測點(diǎn)位置圖Fig. 1 Survey location and layout of the microtremor station point site

2.2 儀器一致性測試

正式進(jìn)行數(shù)據(jù)采集之前, 為保證實(shí)驗(yàn)結(jié)果的可靠性, 需要保證每臺(tái)儀器間具有良好的一致性。因此, 本次實(shí)驗(yàn)使用的10臺(tái)SmartSolo IGU-16HR在正式采集數(shù)據(jù)前, 由深圳面元智能科技有限公司統(tǒng)一進(jìn)行一致性檢測, 使儀器一致性滿足微動(dòng)探測的需求。深圳面元智能科技有限公司提供的儀器的一致性表現(xiàn)圖在文章最后的附錄中展示。

2.3 微動(dòng)臺(tái)陣布設(shè)

空間自相關(guān)法和頻率-波數(shù)法對臺(tái)站布設(shè)的要求不同, 本次實(shí)驗(yàn)采用空間自相關(guān)法, 組成臺(tái)陣一般需要十臺(tái)拾震器, 至少為四臺(tái)。將一臺(tái)拾震器布置在圓心, 剩下的儀器等間距布置在半徑分別為 r,2r, 4r的圓周上(半徑r由測深需求和現(xiàn)場環(huán)境決定),在本次實(shí)驗(yàn)中三個(gè)同心圓的半徑分別為 200 m、400 m和800 m, 臺(tái)陣布設(shè)如圖2。儀器布設(shè)完成后對每個(gè)觀測點(diǎn)進(jìn)行6小時(shí)的連續(xù)觀測, 共10個(gè)觀測點(diǎn), 每個(gè)觀測點(diǎn)上使用10臺(tái)節(jié)點(diǎn)式地震儀, 最終獲得100份三分量地震數(shù)據(jù)。微動(dòng)臺(tái)陣觀測對布局的要求很嚴(yán)格, 需要對各臺(tái)站進(jìn)行精確定位, 且臺(tái)陣中的各臺(tái)站要同時(shí)觀測, 才能保證數(shù)據(jù)的準(zhǔn)確。

臺(tái)陣的半徑對探測深度起主要決定作用。一般來說, 在使用圖 2中的圓形臺(tái)陣觀測時(shí), 勘探深度是圓形臺(tái)陣半徑r的3～5倍。

圖2 微動(dòng)臺(tái)陣布設(shè)以及本次實(shí)驗(yàn)采用的布設(shè)參數(shù)Fig. 2 Layout of the microtremor station points and the designed parameters

2.4 數(shù)據(jù)處理流程

一致性測試達(dá)標(biāo)后, 就要布置臺(tái)陣并采集數(shù)據(jù)。本次實(shí)驗(yàn)中, 臺(tái)陣半徑分別設(shè)置為200 m、400 m、800 m, 觀測系統(tǒng)中共10個(gè)探測點(diǎn), 探測點(diǎn)間距為1000 m。利用RTK(Real-time kinematic)在實(shí)地找到預(yù)設(shè)的點(diǎn)位, 由于本次實(shí)驗(yàn)預(yù)設(shè)的點(diǎn)位有少數(shù)在田間小道旁, 甚至有個(gè)別接近鄉(xiāng)道, 因此實(shí)驗(yàn)中在滿足接近預(yù)設(shè)點(diǎn)位的同時(shí), 對部分點(diǎn)位的儀器布設(shè)進(jìn)行點(diǎn)位偏移, 以確保盡量能夠避免不必要的干擾,比如道路上的行車, 村民的家畜擅自觸碰儀器等情況, 導(dǎo)致數(shù)據(jù)采集不準(zhǔn)確。采集時(shí)間設(shè)定為 6個(gè)小時(shí)。使用SmartSolo IGU-16HR進(jìn)行采集。獲得原始數(shù)據(jù)后的數(shù)據(jù)處理流程(何正勤等, 2007)如圖 3所示。①對原始數(shù)據(jù)進(jìn)行預(yù)處理, 并分析其頻譜, 檢查是否存在明顯的干擾或數(shù)據(jù)缺失; ②對預(yù)處理后的數(shù)據(jù)進(jìn)行數(shù)字濾波, 計(jì)算其空間自相關(guān)系數(shù), 并進(jìn)行方位平均; ③經(jīng)過以上處理后可以獲得瑞雷波的相速度頻散, 手動(dòng)或自動(dòng)選取頻散譜的峰值以獲得頻散曲線; ④對頻散曲線進(jìn)行反演得到S波速度結(jié)構(gòu)并繪制成圖。

圖3 微動(dòng)臺(tái)陣探測數(shù)據(jù)處理流程圖Fig. 3 Flow chart of data processing for microtremor array exploration

圖4展示了2021年7月10日探測點(diǎn)記錄的野外記錄?？梢钥闯鍪_(tái)拾震儀接受到的垂直分量信號在連續(xù)性上表現(xiàn)較好, 表明儀器在觀測期間一直處于正常記錄狀態(tài)。但是原始記錄上存在較多毛刺,即地震記錄中的突兀的高峰, 這是地震儀器偶爾出現(xiàn)問題導(dǎo)致的, 可以在數(shù)據(jù)預(yù)處理部分改正。此外,從圖中還能看到 st.04記錄的數(shù)據(jù)存在一個(gè)較小的傾斜趨勢, 這也可以在預(yù)處理步驟中進(jìn)行改正。

圖4 2021年7月10日探測點(diǎn)原始微動(dòng)數(shù)據(jù)Fig. 4 Microtremor recording acquired on July 10, 2021

獲得原始數(shù)據(jù)后, 首先需要對數(shù)據(jù)進(jìn)行預(yù)處理。利用SAC(Seismic Analysis Code)讀入原始數(shù)據(jù),成圖并檢查每天獲得的地震波波形, 排除受人為因素而嚴(yán)重干擾的數(shù)據(jù)或存在儀器沒能正常記錄的數(shù)據(jù)。對質(zhì)量較好的原始數(shù)據(jù)進(jìn)行預(yù)處理, 進(jìn)行如下操作: 一是去均值、去線性趨勢和波形尖滅, 因?yàn)椴ㄐ螖?shù)據(jù)總會(huì)存在一個(gè)非零的均值或者存在一個(gè)長周期的線性趨勢, 這會(huì)影響到數(shù)據(jù)的分析, 必須在數(shù)據(jù)分析前去除。另一方面, 在對數(shù)據(jù)進(jìn)行譜域操作(如FFT、濾波等)時(shí), 若數(shù)據(jù)的兩端不為零, 則會(huì)出現(xiàn)譜域假象, 因而實(shí)際數(shù)據(jù)經(jīng)常需要做尖滅處理,使得數(shù)據(jù)兩端在短時(shí)間窗內(nèi)逐漸變成零值。二是去毛刺, 地震儀器偶爾會(huì)出現(xiàn)問題, 導(dǎo)致連續(xù)地震數(shù)據(jù)流中出現(xiàn)尖鋒或者數(shù)據(jù)丟失。這些所謂的毛刺,肉眼很容易識別, 但是在使用程序自動(dòng)處理數(shù)據(jù)時(shí)卻很容易被誤認(rèn)為是地震信號, 因而需要在數(shù)據(jù)分析之前將毛刺去除。三是對原始數(shù)據(jù)進(jìn)行濾波, 本文使用的是10 Hz的低通濾波器。最后要對數(shù)據(jù)進(jìn)行分段, 目的是在提取頻散曲線時(shí)將不同時(shí)間段的頻散數(shù)據(jù)進(jìn)行疊加, 獲得更為準(zhǔn)確的結(jié)果。本次的數(shù)據(jù)原長6 h, 按45 min為標(biāo)準(zhǔn)分成若干段。

完成數(shù)據(jù)預(yù)處理后, 需要①定義 10個(gè)檢波器的排列參數(shù), 排列參數(shù)按照儀器的實(shí)際排列情況給出; ②進(jìn)行頻散曲線的提取, 以2021年7月 10日采集的數(shù)據(jù)為例, 在獲得瑞雷波頻散譜后, 手動(dòng)選取頻散譜的峰值, 此時(shí)提取的頻散點(diǎn)不宜太過密集,因?yàn)槊芗念l散點(diǎn)對最終結(jié)果并無太大意義, 且會(huì)極大的降低頻散曲線的反演效率。下面是經(jīng)過疊加后得到的瑞利波相速度頻散曲線(圖5)。③使用分歧型遺傳算法(fGA)反演獲得觀測點(diǎn)地下 S波速度結(jié)構(gòu), 如圖6。圖中的綠線為初始模型, 初始模型可以自動(dòng)創(chuàng)建或者手動(dòng)給定, 這里選擇利用松科二井的地層分層資料給定。④將10個(gè)觀測點(diǎn)的地下S波速度結(jié)構(gòu)整合, 獲得測區(qū)內(nèi)一條橫跨 9 km,深度到地下4000 m的地下二維S波速度結(jié)構(gòu)剖面,結(jié)果如圖7b。

圖5 瑞利波頻散曲線Fig. 5 Rayleigh wave dispersion curve

圖6 橫波速度模型Fig. 6 Shear wave velocity mode

3 結(jié)果與討論

本次微動(dòng)實(shí)驗(yàn), 獲得了一條橫跨9 km, 深度到地下4000 m的地下二維S波速度結(jié)構(gòu)剖面。根據(jù)剖面圖可以發(fā)現(xiàn)地下S波速度總體呈現(xiàn)西側(cè)速度高,東側(cè)速度低的趨勢; 1000 m以上速度基本成層分布,但是在松科二井(圖7b中星形所在位置)西側(cè)100 m、東側(cè)2300 m、東側(cè)3700 m處, 深度在1400～2800 m范圍內(nèi)存在三個(gè)明顯的低速異常體。根據(jù)Xu(2012)的研究與實(shí)驗(yàn)結(jié)果表明, 隱伏斷裂破碎帶在微動(dòng)視S波速度剖面上顯示出明顯的低速異常, 這成為微動(dòng)剖面解釋隱伏地?zé)針?gòu)造的重要標(biāo)志, 為地?zé)峋倪x址提供重要依據(jù)。除此之外, 此類速度異常還可能是含水巖層或石灰?guī)r溶洞。根據(jù)這一結(jié)論, 說明圖中顯示的低速異常體很可能存在斷層或隱伏斷裂構(gòu)造。結(jié)合松科二井地層分層資料(侯賀晟等, 2018),如圖 7a所示, 在 1675～2 964.92 m的深度范圍內(nèi),地層以泉頭組和登婁庫組為主, 泉頭組部分(1675～2540 m), 泉三段巖性為砂巖和泥巖的互層;泉四段, 泉二段和泉一段巖性以泥巖為主, 泉一段底部存在粉砂質(zhì)泥巖層。登婁庫組(2540～2 964.92 m)部分, 除了登二段的底部存在部分砂礫巖層外, 其余部分均為不同粒徑的砂巖與泥巖的互層。因此從整體上看, S波速度剖面存在低速異常的深度范圍內(nèi), 巖性整體以砂巖和泥巖互層為主, 而砂巖孔隙度大, 滲透率高, 適合水熱的儲(chǔ)存, 說明此處具有良好的熱儲(chǔ)存條件。其次, 結(jié)合測區(qū)概況部分提及的松遼盆地的地?zé)岜尘? 相關(guān)研究顯示松遼盆地?zé)崃髦颠h(yuǎn)高于全球平均熱流值, 且其下地層如姚二、三段、姚一段、青二、三段等地?zé)崴匀划a(chǎn)能豐富。因此這些位置可能存在豐富的地?zé)豳Y源。

圖7 松科二井地下1500～3100 m范圍的巖性剖面(a)與微動(dòng)橫波速度剖面(b)Fig. 7Lithologic profile of the SK-2 area between1500and 3100 mbelow the surface (a) andmicrotremor Vs cross section(b)

相較于傳統(tǒng)的地震勘探方法, 本次實(shí)驗(yàn)使用的微動(dòng)探測方法,該方法不受城市中強(qiáng)電磁環(huán)境、人文工業(yè)活動(dòng)的干擾; 其無需主動(dòng)源, 因此對城市環(huán)境沒有損害(田寶卿和丁志峰, 2021)。本次研究區(qū)在東北廣闊的玉米地中, 對周圍環(huán)境的影響極小, 即使測線附近均為當(dāng)?shù)卮迕穹N植的玉米, 儀器依舊能按照預(yù)期位置布置而不影響居民的日常生活及其農(nóng)業(yè)生產(chǎn)。同時(shí), 儀器的布置以及數(shù)據(jù)的采集效率高,每天能采集一到兩個(gè)測點(diǎn)的數(shù)據(jù), 并能及時(shí)對數(shù)據(jù)質(zhì)量做出檢查, 若是儀器及人員充足, 每天能采集的數(shù)據(jù)數(shù)量也能進(jìn)一步增加。

數(shù)據(jù)處理方面, 目前相關(guān)軟件已經(jīng)較為成熟,在保證數(shù)據(jù)質(zhì)量的前提下, 可以通過及時(shí)的數(shù)據(jù)處理結(jié)果顯示, 對相關(guān)參數(shù)進(jìn)行改正, 以獲取更準(zhǔn)確、直觀的結(jié)果。

4 結(jié)論

本文通過微動(dòng)探測方法, 獲得了一條長 9 km,探測深度可達(dá)4000 m的二維S波地下速度剖面。本次實(shí)驗(yàn)很好地發(fā)揮了微動(dòng)探測效率高, 抗干擾能力強(qiáng)的優(yōu)點(diǎn), 野外施工能夠高效、且沒有對測區(qū)內(nèi)農(nóng)民的生產(chǎn)生活造成影響的情況下完成。應(yīng)用目前已經(jīng)成熟的各種軟件, 在保證數(shù)據(jù)質(zhì)量的前提下,可以及時(shí)的通過顯示的數(shù)據(jù)處理結(jié)果, 對相關(guān)參數(shù)進(jìn)行改正, 以獲取更準(zhǔn)確、直觀的結(jié)果。在獲得的S波反演速度剖面上識別出了三個(gè)低速異常體, 依據(jù)前人進(jìn)行相關(guān)工作得到的結(jié)論(Xu et al., 2012),以及結(jié)合測區(qū)的地?zé)岜尘?姜光政等, 2016; 王貴玲等, 2017; 張薇等, 2019), 結(jié)合松科二井獲得的地下分層數(shù)據(jù)(侯賀晟等, 2018), 判斷這三個(gè)低速異常體很有可能對應(yīng)了深部的地?zé)醿?chǔ)水構(gòu)造。

在現(xiàn)今世界, 碳中和已經(jīng)是世界關(guān)注的熱點(diǎn),地球科學(xué)可以在其中發(fā)揮巨大作用(馬冰等, 2021)。微動(dòng)方法憑借其應(yīng)用區(qū)域廣、成本低、探測效率高的優(yōu)勢, 相信能在地?zé)崽綔y方面做出貢獻(xiàn), 為東北經(jīng)濟(jì)振興助力。

致謝: 感謝深圳面元智能科技有限公司提供的儀器與相關(guān)技術(shù)支持。

Acknowledgements:

This study was supported by China Geological Survey (Nos. DD20190010 and DD20221643-2), and Institute of Natural Resources of Heilongjiang Province (No. HX2020-23).