干旱气象
  • CN 62-1175/P
  • ISSN 1006-7639
  • 双月刊
  • 中国科技核心期刊
  • 中国学术期刊综合评价数据库统计源期刊
  • 中文科技期刊数据库收录期刊
作者投稿 专家审稿 编辑办公

干旱气象, 2025, 43(2): 195-206 DOI: 10.11755/j.issn.1006-7639-2025-02-0195

论文

戈壁沙尘运动中风速与颗粒数及动能的标度律

王振亭,1, 孟小楠2,3, 王雪松2, 李庆4

1.干旱区生态安全与可持续发展重点实验室，中国科学院西北生态环境资源研究院，甘肃 兰州 730000

2.北京师范大学地理科学学部，地表过程与资源生态国家重点实验室/防沙治沙教育部工程研究中心，北京 100875

3.北京市第八十中学睿德分校，北京 100176

4.河北省科学院地理科学研究所/河北省地理信息开发应用工程技术研究中心，河北 石家庄 050011

Scaling laws in wind speed, particle number and kinetic energy in Gobi sand movement

WANG Zhenting,1, MENG Xiaonan2,3, WANG Xuesong2, LI Qing4

1. Key Laboratory of Ecological Safety and Sustainable Development in Arid Lands， Northwest Institute of Eco-Environment and Resources， Chinese Academy of Sciences， Lanzhou 730000， China

2. State Key Laboratory of Earth Surface Processes and Resource Ecology/MOE Engineering Research Center of Desertification and Blown-Sand Control，Faculty of Geographical Science，Beijing Normal University， Beijing 100875， China

3. Beijing No. 80 Middle School Rui De Branch， Beijing 100176， China

4. Institute of Geographical Sciences/Hebei Engineering Research Center for Geographic Information Application， Hebei Academy Science， Shijiazhuang 050011， China

责任编辑: 王涓力；校对：邓祖琴

收稿日期: 2024-06-27   修回日期: 2024-08-21  

基金资助: 国家自然科学基金项目(42371005)
国家自然科学基金项目(U21A2001)

Received: 2024-06-27   Revised: 2024-08-21  

作者简介 About authors

王振亭（1975—），男，研究员，主要从事风沙物理与风沙地貌研究。E-mail：ztwang@lzb.ac.cn。

摘要

中国西北极端干旱区戈壁广布，风大沙多，灾害频发。充分了解戈壁沙尘运动基本规律是灾害预警与科学防治的重要前提。鉴于目前难以精准预测瞬时输沙率，探讨风沙事件中气流与沙尘特征物理量的统计规律，进而开展统计预测，也许是在秒及以下时间尺度上建立风和沙之间定量关系的可行之路。本研究借鉴湍流统计理论的思路和方法，利用Hilbert-Huang变换分析4个戈壁沙尘运动野外观测数据集。结果表明，沙尘事件中的风速、跃移沙粒动能与个数、粉尘浓度等时间序列的Hilbert边际谱均符合幂次标度律，沙尘特征量和风速的标度指数范围分别为0.78~1.51和0.59~1.47。

关键词： 戈壁; 风沙; 粉尘; Hilbert边际谱; 标度律

Abstract

In the extreme arid region of Northwest China, the Gobi desert is widely distributed, with strong winds and abundant sand, leading to frequent disasters. A deep understanding of the basic laws of the aeolian sand-dust movement over Gobi surfaces is an important prerequisite for disaster warning and scientific prevention. Given the current difficulty in accurately predicting the instantaneous aeolian sediment transport rate, exploring the statistical laws of airflow and sand-dust characteristic physical quantities in wind-sand events, and then conducting statistical forecasting, may be a feasible way to establish a quantitative relationship between wind and sand on time scales of seconds or less. This study borrows the ideas and methods from the statistical theory of turbulence, and analyzes four field observation datasets of Gobi sand-dust movement using the Hilbert-Huang transform. The results indicate that the Hilbert marginal spectra of the time series of wind speed, kinetic energy and number of saltation sand particle, and dust concentration in the aeolian events all follow the power scaling law. The scaling exponents of variables, characterizing aeolian sand and dust motions, and wind speed range from 0.78 to 1.51 and 0.59 to 1.47, respectively.

Keywords： Gobi; aeolian sand; dust; Hilbert marginal spectra; scaling law

PDF (10364KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

王振亭, 孟小楠, 王雪松, 李庆. 戈壁沙尘运动中风速与颗粒数及动能的标度律[J]. 干旱气象, 2025, 43(2): 195-206 DOI:10.11755/j.issn.1006-7639-2025-02-0195

WANG Zhenting, MENG Xiaonan, WANG Xuesong, LI Qing. Scaling laws in wind speed, particle number and kinetic energy in Gobi sand movement[J]. Arid Meteorology, 2025, 43(2): 195-206 DOI:10.11755/j.issn.1006-7639-2025-02-0195

0 引言

我国戈壁主要分布在西北极端干旱区，面积约5.695×105 km2（王涛和陈广庭，2008）。戈壁地区往往地势平坦，植被稀疏，强风及其携带的沙尘经常破坏基础设施。以新疆戈壁大风区的铁路系统为例，风沙危害具体表现为降低能见度、磨蚀轮轨、损坏车窗、掩埋线路、吹翻列车等（尹永顺，1990；Cheng et al.，2015；Shu et al.，2022）。充分了解戈壁沙尘运动基本规律是科学开展其相关危害预警与防治的重要前提。在过去30余年内，国内风沙工作者以野外观测和室内风洞为主要研究手段，对这种现象的认识越来越深入。与流沙地表附近的风沙流相同，戈壁风沙流的常规研究内容也是起动机制、跃移轨迹、发展过程、风沙流结构与输沙率等。其中，处于核心地位的输沙率指单位时间单位宽度内顺风向输运的沙物质质量。戈壁风沙流输沙特征的主要控制因素包括风况和沙源丰富度（柳本立等，2010）、地表物质机械组成（Tan et al.，2022）、砾石盖度（Liu and Kimura，2018；Tan et al.，2019；Wang et al.，2022）和风程（张正偲等，2011；Zhang et al.，2022）等。地表紧实且风速较大导致跃移轨迹具有起跳角、高度和长度都很大的突出特点（邹学勇等，1995）。实验结果表明，跃移轨迹的高度和长度可达到4.8和9.0 m（Liu et al.，2021）。因沙源有限，地表的松散沙粒难以满足风的最大搬运能力，风沙流常处于非饱和状态，出现异常结构。早期的集沙仪与风速廓线仪观测表明，输沙通量在一定高度处呈现最大值并随风速的增加而增大（屈建军等，2005；张克存等，2005），沙粒的质量浓度则随高度增加而单调下降（蒋富强等，2010；李凯崇等，2010）。多项式函数、指数函数和幂函数均被用于描述风沙流结构（屈建军等，2012；黄雨晖等，2019）。近年来的碰撞传感器和超声风速仪观测显示，沙粒数随高度呈指数形式递减，跃移沙动能通量是高度的Pearson VII函数（Tan et al.，2021）。实际上，输沙率和风沙流结构实验结果均与观测时间有关。在分钟及更长时间尺度上，平均输沙率与风速之间具有良好对应关系。1 min的平均输沙率符合Owen公式（Tan et al.，2013；Tan et al.，2016）；10 min到数小时的平均输沙率是风速的指数（王志强等，2010）或低阶多项式函数（刘旭阳等，2020）；以Bagnold公式为基础，可建立年际尺度的输沙量算法（王旭，2019）。目前，在秒及以下的短时间尺度上预测输沙率是公认难题（Chamecki and Kok，2023），难以在戈壁风沙流瞬时输沙率与风速之间建立对应关系（Tan et al.，2020）。

瞬时输沙率难以预测的根本原因在于气流与运动沙粒之间同时存在不同时空尺度的复杂耦合作用（Wang，2024）。风沙流在空间上表现为一丝一缕的条带状（Baas，2008），在时间上具有明显的间歇性（Stout and Zobeck，1997；Davidson-Arnott et al.，2012；Wang et al.，2014；Sherman et al.，2018；Tan et al.，2020）和非平稳性（Butterfield，1998；Spies et al.，2000；Wang and Zheng，2015；Wang et al.，2018）。这些特征给瞬时风沙流的数学描述与实验数据分析带来极大挑战。从数学描述角度看，虽然平均风速大于起动风速所占的时间比例（Stout and Zobeck，1997）、跃移活动参数（Davidson-Arnott et al.，2012）、沙粒浓度峰度（Wang et al.，2014）等均被作为间歇性指标，但因缺乏统一的测量和量化标准而使不同实验结果之间难以比较（Sherman et al.，2018）。对于驱动沙尘运动的湍流，流体力学研究者们已经给出成熟的统计描述。对沙尘运动进行统计预测可能是替代瞬时输沙率预测的可行方案。物理量随时空尺度的变化规律称为标度律。在湍流统计理论中，标度律专指高雷诺数流动的速度结构函数随空间两点间隔的幂次律。著名的K41理论给出了各向同性均匀湍流在惯性区的标度指数为‒5/3。在阿拉善高原两处流动沙丘迎风坡顶的野外观测初步显示，流沙地表附近几乎饱和风沙流的跃移输沙密度也存在类似规律（Wang，2024）。戈壁地表的风沙流与流沙地表的差别很大，虽然进行了大量观测，但尚未探讨可能存在的标度律。从实验数据分析角度看，沙尘和气流的运动都是非平稳的，统计特征随时间而变（Wang et al.，2018）。许多基于线性或平稳过程发展起来的数学工具不适合非平稳过程。适用于风沙物理量非平稳时间序列分析的有效方法包括小波变换（Baas，2006；Ellis and Sherman，2023）、Hilbert-Huang变换（Wang et al.，2018；Wang，2024）和动态时间规整等。其中，Hilbert-Huang变换的适应性很强，提取出的信号特征通常具有明确的物理意义（Huang and Wu，2008；de Souza et al.，2022）。本文借鉴湍流的统计描述思想，利用Hilbert-Huang变换分析处理中国西部典型戈壁区的观测数据，探讨近地表气流、沙和尘物理量的统计特征，以期为沙尘运动的统计预测提供方法参考和数据支撑。

1 数据来源

为了深入了解戈壁风沙运动规律，数个风沙研究团队近年来在中国西北干旱区开展了大量野外观测工作。本文数据来自4个野外实测数据集，物理量、观测高度、传感器型号、采样频率和精度等基本信息列于表1。三维超声风速传感器WindMaster Pro和Young-81000均可在较高采样频率和精度下工作，广泛应用于风沙环境中的风速测量。跃移传感器H11-Lin是比较成熟的常用仪器，但个体之间的灵敏度差异较大，测量值与真实值往往有一定偏差。该类传感器测量给定迎风面积内所有跃移沙粒的总平动动能（Kinetic Energy，KE）和个数（Particle Count，PC）随时间的变化。KE和PC的具体数值依赖于传感器的迎风面积。DustTrak-II粉尘仪也适用于恶劣的野外环境，记录了人工和天然两种不同戈壁地表附近粉尘浓度DA（Dust Concentration over Artificial Gobi Surfaces，DA）和DN（Dust Concentration over Natural Gobi Surfaces，DN）随时间的变化。

表1   戈壁沙尘运动数据集及观测条件

Tab.1  Datasets of aeolian sand and dust motions over Gobi and measurement conditions

数据来源高度/m物理量传感器（厂家）采样频率/Hz采样精度
Wang et al.，20220.30、0.83、3.00U、V、WWindMaster Pro（Gill）200.001 m·s-1
Tan et al.，20200.05、0.12、0.38、0.80、1.38KE、PCH11-Lin（Sensit）11×10-7 J·s-1、s-1
0.70、2.00U、φ、θYoung-81000（RM Young）200.01 m·s-1、0.1°
Wang et al.，20230.05、0.12、0.38、0.80、1.41KE、PCH11-Lin（Sensit）11×10-7 J·s-1、s-1
1.20、2.80U、V、WWindMaster Pro（Gill）10.001 m·s-1
本研究团队0.10、0.20、0.50、1.00、1.50、2.00DA、DNDustTrak-II（TSI）0.20.001 mg·m-3
0.50、1.00U、V、WYoung-81000（RM Young）200.01 m·s-1

注：U为球坐标系下的风速大小或直角坐标系下的水平风速分量；V为直角坐标系下的另一个水平风速分量；W为垂向风速；φ为仰角，θ为方位角；KE为动能，PC为个数；DA为人工地表粉尘浓度，DN为天然地表粉尘浓度。下同。

新窗口打开| 下载CSV


4个野外观测地点位置见图1，基本信息如下：

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   4个野外观测地点分布

注:基于国家测绘地理信息局标准地图服务网站标准地图[审图号为GS ( 2016 ) 2923号]制作，底图无修改。

Fig.1   The distribution of 4 field observation sites


1 ）观测地点位于吐鲁番—哈密盆地中部十三间房的天山冲洪积扇上（91.73°E，43.13°N），鹅卵石覆盖度高达87%左右，年平均风速7.13 m·s-1，最大瞬时风速56.60 m·s-1，年均大风（10 m高度处风速超过17.10 m·s-1）日数179 d（Cheng et al.，2015），处于高风能环境中。仪器布设详见Wang等（2022）文献，观测的物理量包括气温、气压、空气湿度、风速、输沙率等。本文仅分析一个由862 000条风速记录组成的子数据集。

2 ）观测地点位于阿尔金山西麓新疆若羌县米兰镇的阿尔金山冲洪积扇上（89.00°E，39.13°N），40%~50%的地表被细小砾石覆盖，年均大风日数为13 d，处于高风能环境中。仪器布设详见Tan等（2020）文献。该数据集共有25 200条跃移沙粒记录和504 000条风速记录，每条记录包含5个高度处的沙粒动能和个数与2个高度处的风速、仰角和方位角，详细描述及下载链接见Mendeley网站（https://doi.org/10.17632/jdjmg8s544.1）。

3 ）观测地点位于吐鲁番—哈密盆地东南部的南湖戈壁（92.77°E，42.25°N），36.59%的地表被8~10 mm大小的砾石所覆盖，处于中等风能环境中。仪器布设详见Wang等（2023）文献。该数据集共有19 799条记录，每条记录包含5个高度处的沙粒动能和个数与2个高度处的3个风速分量，详细描述及下载链接见Mendeley网站（https://doi.org/10.17632/mwyftts329.1）。

4 ）观测地点位于党河南山以北的阿克塞盆地中（94.40°E，39.70°N），分别在去除植被并对上风向地表实施化学固沙措施的人工观测场和约15%植被盖度的自然观测场中进行观测，所测量的PM10质量浓度分别称为人工地表和天然地表粉尘浓度。两个观测场地表物质均以粗砂和中砂为主，质量比83%~91%。受地形影响，观测地点处于低风能环境中。对观测场的具体描述见相关文献（刘旭阳等，2020；Zhang et al.，2022）。该数据集包括14 d人工观测场和5 d自然观测场内5个高度处的PM10质量浓度和2个高度处的3个风速分量，由本研究团队于2017年3月16日—5月2日期间陆续观测获得。

2 研究方法

描述戈壁沙尘运动的特征物理量选可以直接测量的风速U和W、跃移沙粒动能KE与个数PC、粉尘浓度DA和DN。利用Hilbert-Huang变换探讨标度律的具体步骤依次为特征物理量时间序列选取、模态分解、Hilbert边际谱计算、标度指数拟合。

2.1 时间序列选取

只有在大于起动风速时，才会有风沙活动。因此，沙尘输运过程具有非常显著的间歇特点。只有先定义沙尘事件，才能进行精细的数据分析。当没有沙尘活动时，表1中的跃移沙粒动能KE、个数PC、粉尘浓度DA和DN均为零。参照对流沙地表附近风沙事件的定义（Wang，2024），这里仍将任意一个沙尘特征物理量在某段时间内均不为零的序列定义为一个戈壁沙尘事件。据此，沙尘事件持续的最短时间为一个采样间隔，即表1中采样频率的倒数。除了采样频率，沙尘事件还依赖于传感器灵敏度和精度等。在相同环境条件下，不同传感器探测到的沙尘事件时间和强度往往有所不同，从而由观测数据定义的沙尘事件有一定程度的不确定性。如果沙尘运动服从某种统计规律，那么由不同时间序列得到的统计特征及其定量描述中的概率密度函数和相关参数应该一致。

分析沙尘事件时，时间序列越长，得到的特征参数越准确。受间歇性影响，持续时间大于10 min的跃移事件实际上非常稀少。本文根据各数据集的数据量选定进一步分析的沙尘事件。跃移沙粒KE和PC、粉尘浓度DA和DN序列的最短时长分别取60 s、60 min；风速的非零数据连续性较好，由4个数据集中截取水平风速U和垂向风速W的时间序列时，时长分别取25、25、200和60 s。需要说明的是，本文重点关注各物理量的统计特征，未探讨相互关系，各时间序列都按连续非零原则独立获取。

2.2 Hilbert-Huang变换

对沙尘运动特征物理量的时间序列进行Hilbert-Huang变换的目的是计算能谱，进而探讨可能存在的标度律。下面以水平风速U为例说明计算过程。

首先，通过经验模态分解将U写为趋势项和少数几个脉动项之和：

Ut=U0t+U1t=U0t+∑i=1mUi't

式中：t为时间，U0为趋势项，U1为总脉动项，Ui'为有限阶脉动项，m为项数。每个脉动项又称模态，必须满足极值点与过零点个数最多相差一个且上下包络之和处处为零两个条件。最初的经验模态分解算法存在分解不彻底与模态混杂问题（Huang et al.，1998；de Souza et al.，2022）。通过在每个分解步骤中添加自适应白噪声信号，可以保证模态之间的正交性（Colominas et al.，2014）。

然后，分解后的每个模态都有自己的特征频率且与其他模态正交，通过Hilbert变换构成一个复信号：

U˜i't=Ui't+j1π∫-∞+∞Ui'τt-τdτ=aite jϕit

式中：j为虚数单位，ait和ϕit分别为第i个复信号U˜i'的模和幅角。

由公式（2），可定义瞬时频率：

fi=12πdϕidt

脉动速度Ui'和速度U的Hilbert谱为

Hif, t=ait,0,  如果f=fi其他Hf, t=∑i=1mHif, t              

最后，速度U的Hilbert边际谱为

EUf=1T∫t1t2Hf, tdt

式中：t1和t2分别为沙尘事件的起止时刻，时长T=t2-t1。Hilbert边际谱EUf表征特征频率为f的脉动风速的能量密度，是传统Fourier谱的一种替代，二者总体特征类似，但局部细节有所区别（de Souza et al.，2022）。同样地，可以得到垂向风速W、跃移沙粒动能KE和个数PC、粉尘浓度DA和DN的Hilbert边际谱。

在充分发展的各向同性均匀湍流中，脉动速度能量谱与频率之间符合幂次标度律：

EUf~f-βU

式中：指数βU=5/3。受地面和空中沙尘影响，水平风速U和垂向风速W的Hilbert边际谱指数不一定符合5/3的理论值。

类似地，其他物理量的Hilbert边际谱也可能是频率的幂函数：

EKEf~f-βKEEPCf~f-βPCEDAf~f-βDAEDNf~f-βDNEWf~f-βW

式中：EKE、EPC、EDA、EDN、EW分别为跃移沙粒动能KE、个数PC、粉尘浓度DA和DN、垂向风速W的边际谱，βKE、βPC、βDA、βDN、βW为对应的标度指数。

具体计算时，首先按Colominas等（2014）改进的自适应噪声完备集合经验模态分解（Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise， ICEEMDAN）算法实现公式（1）形式的正交分解，然后根据公式（2）~（5）得到Hilbert边际谱，最后探讨各Hilbert边际谱是否符合公式（6）和（7）的形式，若符合则利用最小二乘法拟合得到各指数值。

3 结果分析

3.1 时间序列分解

按照沙尘事件定义和时间序列选取方法，从4个数据集中筛选出1~37个时长不等的沙尘事件，并依次选取1 724、1 006、98和486个风速序列。实际计算表明，脉动项个数m一般为4~10。每个风速脉动项具有不同的特征时间尺度，彼此叠加即可构成复杂的湍流相干结构。除了地表影响，沙尘脉动主要由风速脉动引起。如果表征风与沙的脉动项之间的小波相干系数较大，那么二者之间大致成正比例关系（Wang et al.，2018）。每个脉动项都有自己的特征频率，可通过构造复频率方式定量描述风与沙脉动项之间的相互作用（Wang，2024）。不同于小波变换，Hilbert-Huang变换尚缺乏严格的数学基础（de Souza et al.，2022）。这里仅将模态分解算法作为计算Hilbert谱的一种有效手段。

3.1.1 跃移沙粒动能与个数

第2个数据集中0.05 m高度处一组跃移沙粒动能与个数时间序列的分解结果见图2，模态项数均为4。可以看出，动能变化趋势为先缓慢上升，然后在40 s左右快速下降；沙粒数先快速上升，约30 s后平缓变化。图3为分解后的动能脉动项KEi' （i=1， 2， 3， 4），可见各模态的特征频率存在明显差别，振幅均随时间变化完美地刻画了图2中原始信号的局部变化细节。原始信号的高频部分主要体现在前两阶模态KE1'和KE2'中，低频部分则分解到后两阶模态KE3'和KE4'中。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   跃移沙粒动能（a）与个数（b）时间序列分解

[数据来自Tan等（2020）文献，测量高度0.05 m；物理量下标0为趋势项，1为脉动项，下同]

Fig.2   Decomposition of time series of kinetic energy （a） and count （b） of saltating sand grains

(Data are sourced from the ref. (Tan et al., 2020), measured at 0.05 m height. Subscripts of variables 0 is for trend, 1 is for fluctuation, the same as follows)


图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   跃移沙粒动能分解后的脉动模态

（a）一阶模态,（b）二阶模态,（c）三阶模态,（d）四阶模态

[数据来自Tan等（2020）文献，测量高度0.05 m]

Fig.3   Fluctuation modes of kinetic energy of saltating sand grains after decomposition

(a) the first mode, (b) the second mode, (c) the third mode, (d) the fourth mode

(Data are sourced from the ref. (Tan et al., 2020), measured at 0.05 m height)


第3个数据集中0.05 m高度处一组跃移沙粒动能与个数时间序列的分解结果见图4，模态项数均为6。该事件的持续时间约为图2所示事件的6倍，特征量在数值上也远大于前者。随着时间增加，动能和沙粒数均呈现出缓慢下降趋势；二者脉动项在个别时刻的幅值非常大。需要特别说明的是，动能和沙粒数的变化规律并不完全一致；沙粒数代表风沙流中固体物质的浓度，动能由沙粒的个数和速度共同决定，二者既有关联又有一定的独立性。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   跃移沙粒动能（a）与个数（b）的时间序列分解

[数据来自Wang等（2023）文献，测量高度0.05 m]

Fig.4   Decomposition of time series of kinetic energy （a） and count （b） of saltating sand grains

(Data are sourced from the ref. (Wang et al., 2023), measured at 0.05 m height)


3.1.2 粉尘浓度

0.10 m高度处两组粉尘浓度时间序列的分解结果见图5，模态项数均为6。其中，人工地表粉尘浓度专指在去除植被且对上风向边界处地表做过不可蚀处理的观测场（刘旭阳等，2020；Zhang et al.，2022）内测量的粉尘浓度。因为并非同步观测，图5中人工和天然戈壁观测场的粉尘浓度不具有可比性。在1 h内，人工观测场近地表粉尘浓度大致表现为先增加后稳定再下降的趋势，因趋势项数值较低，脉动现象比较明显；天然观测场的沙尘浓度一直呈缓慢上升趋势，因为趋势项数值较大，脉动项更像在零点附近波动的噪声，实际并非如此，脉动项振幅与人工观测场很接近，亦有一定的动态规律性。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   人工地表（a）与天然地表（b、c）粉尘浓度（测量高度0.10 m）时间序列分解

Fig.5   Decomposition of time series of dust concentration of artificial Gobi surface (a) and natural Gobi surface (b, c) measured at 0.10 m height


3.1.3 风速

图6为第4个数据集中人工戈壁观测场0.50 m高度处的一组风速分解结果。从数值上看，水平风速与垂向风速的趋势项差别较大，最大脉动幅度则比较接近。虽然垂向风速的趋势项接近于零，总脉动项近似等于原始信号，但是各阶模态的正交性极大方便了Hilbert边际谱的计算。风速的模态分解属于湍流研究中的常规工作，不再给出各模态的具体结果。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   水平风速（a）与垂向风速（b）（测量高度0.50 m）时间序列分解

Fig.6   Decomposition of time series of horizontal (a) and vertical (b) wind speeds measured at 0.50 m height


3.2 标度律

图2、图4~6各时间序列对应的Hilbert边际谱见图7~10。可以看出，虽然不同序列的脉动项差异很大，但谱特征非常一致。这是因为通过定义沙尘事件剔除了时间序列中无风沙活动的影响。所有计算结果均表明，无论风、沙还是尘的各种特征量的Hilbert边际谱均随着频率升高先迅速上升后缓慢下降。这说明公式（6）和（7）的幂次标度律猜测比较合理。可能轻视风沙流的间歇性、小波基选取不当或者对极少序列的不恰当平均等原因，Baas（2006）得到的风速和输沙率密度小波谱并没有如此良好的规律性。由Karman经验谱公式能够得到很好的全局性风速能谱拟合结果（刘江和王元，2009）。在壁湍流的不同频率区间内，风速能谱分别符合经典的“-1”和“-5/3”标度律（Marusic and Monty，2019）。可能受沙尘运动或地形起伏影响，我们得到的能谱难以进一步按频率区间进行区分。因此，将Hilbert边际谱最大值时的频率作为幂次标度律区间起点来拟合公式（6）和（7）中的各标度指数。对于同一个沙尘事件，跃移沙粒动能和个数的标度指数可能存在较大差异（图7），也可能非常接近（图8）。跃移沙粒动能和个数是相互独立的物理量，在个别事件中的标度指数接近属于巧合。人工地表和天然地表的粗糙元差别很大，前者以砾石为主，后者存在植被，从而导致二者的近地表气流场结构及相同高度处粉尘浓度标度指数均存在很大差异。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   跃移沙粒动能（a）与个数（b）的Hilbert边际谱

（相应的时间序列见图2）

Fig.7   Hilbert marginal spectra of kinetic energy (a) and count (b) of saltating sand grains

(The corresponding time series is shown in fig. 2)


图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   跃移沙粒动能（a）与个数（b）的Hilbert边际谱

（相应的时间序列见图4）

Fig.8   Hilbert marginal spectra of kinetic energy (a) and count (b) of saltating sand grains

(The corresponding time series is shown in fig. 4)


图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   人工地表（a）与天然地表（b）粉尘浓度的Hilbert边际谱

（相应的时间序列见图5）

Fig.9   Hilbert marginal spectra of artificial (a) and natural (b) surface dust concentration

(The corresponding time series is shown in fig. 5)


图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   水平（a）和垂向（b）风速的Hilbert边际谱

（相应的时间序列见图6）

Fig.10   Hilbert marginal spectra of horizontal (a) and vertical (b) wind speed

（The corresponding time series is shown in fig. 6)


表2列出了不同观测条件下标度指数的平均值、标准差和根据t分布计算的95%置信区间。作为气流场的边界，地面影响着附近的风速和沙粒运动。空中运动的沙尘颗粒也改变了湍动能的传递和耗散过程。因此，沙尘运动各物理量标度指数与净气流速度的“-1”或“-5/3”有所差别。从均值看，跃移沙粒动能KE的指数大致随高度上升而下降；沙粒个数PC的指数随高度上升而下降（第2数据集）或者先降后升（第3数据集）。这是因为前者缺少0.80 m及以上高度的时间序列。人工戈壁观测场粉尘浓度的指数基本随高度上升而下降，天然观测场的则无明显规律。从数值上看，沙尘特征量和风速的各指数差别较大，范围分别为0.78~1.51和0.59~1.47。此前，利用Hilbert-Huang变换得到的流沙地表4 cm高度内跃移沙粒数序列和平坦地面上10 m高度处水平风速序列的Hilbert边际谱指数分别为1（Wang，2024）和1.27（Wei et al.，2017），包含在表2中的能谱指数平均值范围之内。

表2   不同观测条件下的沙尘运动各特征物理量能谱的标度指数拟合结果

Tab.2  The fitting results of scaling exponents of energy spectrum of physical variables describing aeolian sand and dust motions under different measurement conditions

数据来源高度/m物理量序列数均值标准差95%置信区间物理量序列数均值标准差95%置信区间
Wang et al.，20220.30U1 7241.4380.418[1.418， 1.457]W1 7241.3360.491[1.313， 1.360 ]
0.831.4450.345[1.429， 1.461]1.3600.456[1.339， 1.382 ]
3.001.4680.360[1.451， 1.485]1.4040.508[1.380， 1.428]
Tan et al.，20200.05KE371.1410.175[1.083， 1.200]PC341.3930.295[1.290， 1.496]
0.12371.1100.193[1.046， 1.174]341.3800.198[1.311， 1.449]
0.38181.0380.165[0.956， 1.120]241.2240.212[1.135， 1.314]
0.80PC11.451
0.70U1 0061.1920.055[1.188， 1.195]W1 0061.3130.143[1.305， 1.322]
2.001.3130.143[1.305， 1.322]1.1340.184[1.123， 1.146]
Wang et al.， 20230.05KE301.0210.208[0.943， 1.098]PC231.2340.145[1.172， 1.297]
0.12240.9290.215[0.839， 1.020]281.1760.186[1.104，1.248]
0.38150.9920.138[0.916， 1.068]141.0550.185[0.948，1.161]
0.80210.9070.131[0.847， 0.966]151.1210.137[1.045，1.197]
1.41190.7810.183[0.692， 0.869]211.2040.216[1.106，1.302]
1.20U980.8350.106[0.813， 0.856]W980.5900.089[0.572，0.608]
2.800.8910.100[0.871， 0.911]0.6290.104[0.608，0.650]
本研究0.10DA111.5130.228[1.360， 1.665]DN121.0150.210[0.882，1.149]
0.2071.4000.192[1.224， 1.577]111.2760.185[1.152，1.401]
0.50111.2760.185[1.152， 1.401]180.9650.164[0.883，1.046]
1.00230.9220.199[0.837， 1.008]120.7560.077[0.708，0.805]
1.50101.0220.216[0.868， 1.176]
2.00150.9010.108[0.842， 0.961]110.9320.117[0.854， 1.010]
0.50U4861.1870.043[1.183， 1.190]W4860.8990.067[0.893， 0.905]
1.001.2140.042[1.211， 1.218]1.2140.041[1.210， 1.218]

新窗口打开| 下载CSV


4 讨论

戈壁地表附近的沙尘运动由气流驱动。根据牛顿力学原理，沙尘颗粒的轨迹和运动速度应该能够被气流速度精确预测。目前，无论流沙地表，还是戈壁地表，瞬时输沙率都难以预测。主要原因可能在于地表附近湍流结构复杂多样、风沙流输运过程的间歇性特征、气流和沙粒在不同时空尺度上的相互作用。传统风沙输运模型中的摩阻速度代表了气流对地表的时均剪切力，不能描述瞬时气流场。本研究通过连续非零数据定义沙尘事件，克服了间歇性带来的困难，进而由Hilbert-Huang变换得到了特征物理量的能谱。符合幂次标度律的Hilbert边际谱意味着戈壁沙尘运动很可能像湍流那样建立统计理论，进而实现风沙运动的统计预测，克服瞬时输沙率预测的困难。各物理量的幂次律标度指数随高度变化则说明沙尘运动规律比壁湍流更复杂。本文的不足之处在于，风速和沙尘物理量的空间测量位置存在较大差异，难以从实验上测定戈壁气流和沙尘在采样间隔及更小时间尺度上的相互作用，进而揭示物理机制。利用更小尺寸与更高采样频率的传感器开展野外观测，并通过更先进的数学工具分析数据，可望克服这一缺陷。另外，弄清戈壁沙尘事件中风速、沙尘颗粒速度、浓度等物理量的概率密度是统计描述和预测的前提，亦有赖于将来的工作。

5 结论

本研究利用前人的风沙流数据集和自行获得的粉尘数据，通过Hilbert-Huang变换对戈壁沙尘事件的特征物理量进行非平稳时间序列分析，得到以下初步结论。

1 ）由连续非零观测数据定义的沙尘事件是探讨戈壁沙尘运动统计特征的可靠素材。

2 ）Hilbert-Huang变换能够成功提取沙尘运动特征物理量的变化趋势和具有特定频率的脉动细节。

3 ）沙尘事件中的风速、跃移沙粒动能与个数、粉尘浓度等物理量的Hilbert边际谱均符合幂次标度律，沙尘特征量和风速的标度指数范围分别为0.78~1.51和0.59~1.47。

参考文献

黄雨晖, 韩小元, 赵健, 等, 2019.

新疆戈壁地区风沙流结构及其粒径特征研究

[J]. 气象与减灾研究, 42(3):199-205.

[本文引用: 1]

蒋富强, 李荧, 李凯崇, 等, 2010.

兰新铁路百里风区风沙流结构特性研究

[J]. 铁道学报, 32(3):105-110.

[本文引用: 1]

李凯崇, 蒋富强, 薛春晓, 等, 2010.

兰新铁路十三间房段的戈壁风沙流特征分析

[J]. 铁道工程学报, 27(3):15-18.

[本文引用: 1]

研究目的：通过对兰新铁路“百里风区”风沙成因进行分析,研究戈壁砾漠地区风沙流密度沿高度变化的特点,不同高度风沙流密度与风速的关系及大风携沙的粒径分布特征。研究结论：戈壁地区大风携沙主要集中在3 m以下,约占总携沙量的87%左右;在同一高度条件下,风沙流密度随着风速的增加显现e的指数增加;而当风速条件一致时,风沙流密度随着高度的增加显现下降趋势。集沙颗粒分析则指出:在不大于3 m的高度范围内百里风区的风沙流运动形式是以跃移为主;超过3 m以后,则以悬移为主。

刘江, 王元, 2009.

利用静态小波变换的非定常阵风分析与模拟

[J]. 西安交通大学学报, 43(7):74-78.

[本文引用: 1]

刘旭阳, 宁文晓, 王振亭, 2020.

两种戈壁地表风沙流特征的野外观测

[J]. 干旱区研究, 37(4):1087-1 094.

[本文引用: 3]

柳本立, 张伟民, 刘小宁, 等, 2010.

莫高窟顶戈壁偏东风作用下输沙率变化的观测研究

[J]. 中国沙漠, 30(3):516-521.

[本文引用: 1]

通过一次沙尘暴天气过程对莫高窟顶不同观测点戈壁风沙流结构和输沙率的时空变化特征进行了系统同步观测。观测表明,偏东风条件下,当风速为10 m·s-1时,水平观测方向上存在一个各点输沙率平衡的阶段,即崖顶戈壁至鸣沙山边缘戈壁各观测点的输沙率相当;当风速大于10 m·s-1时,戈壁地表以风蚀搬运作用为主,出现了戈壁向鸣沙山的长距离风沙搬运;当风速小于10 m·s-1时,戈壁沙物质的长距离搬运作用不明显。不同风况及沙源条件下,风沙流结构和输沙特征会发生显著变化。借鉴Lattau输沙率方程,运用风速和输沙率推算各观测点输沙率的时空变化,计算结果表明,在风速为10 m·s-1左右时,窟顶戈壁带的输送效率约为60%。

屈建军, 黄宁, 拓万全, 等, 2005.

戈壁风沙流结构特性及其意义

[J]. 地球科学进展, 20(1):19-23.

DOI      [本文引用: 1]

通过砾质戈壁风沙流野外实测数据的分析以及风洞模拟实验研究发现：戈壁风沙流结构具有与沙漠风沙流完全不同的风沙流特征，戈壁风沙地表的粗糙度随风速的增大而增加，其表面风沙流输沙量高度分布表现出独特的"象鼻"效应，在一定高度处呈现最大值，并随风速的增加而增高。该"象鼻"效应导致戈壁风沙流结构特征值λ远大于 1，不论风速多大，风沙流都处于未饱和状态的非堆积搬运状态。这种特殊性质比较清楚地解释了敦煌莫高窟千年来不被沙山埋没的谜底。并且在风沙防治工程实践中，采用砾石压沙措施，构造类似与砾质戈壁的下垫面，人工促使风沙流结构呈现"象鼻"形状，可使防沙工程达到理想的输导作用。

屈建军, 张克存, 张伟民, 等, 2012.

几种典型戈壁床面风沙流特性比较

[J]. 中国沙漠, 32(2):285-290.

[本文引用: 1]

通过不同材料覆盖的戈壁床面风沙流特性风洞模拟实验,发现对于棱角状砾石戈壁床面,地表动力学粗糙度随风速的增加而增加;而卵石床面,地表动力学粗糙度随风速的增加呈减小趋势。戈壁床面风速随高度的分布同样满足对数规律,棱角状砾石床面对风速的减弱程度相对于卵石床面更趋于显著。沙粒与戈壁床面棱角状砾石发生碰撞时其起跳高度增大,引起含沙量随高度分布不再满足流沙地表的指数衰减规律,而呈“象鼻效应”,出现拐点。戈壁地表输沙率与风速服从幂函数关系,但其幂指数远大于流沙地表。输沙率与风速之间幂函数关系中幂指数的取值主要受控于地表粒度组成。

王涛, 陈广庭, 2008. 中国的沙漠·戈壁[M]. 上海: 上海科学技术文献出版社:216-224.

[本文引用: 1]

王旭, 2019.

金阿铁路戈壁风沙流输沙量研究

[J]. 科技资讯,2019, 17(13):180-182.

[本文引用: 1]

王志强, 黄晟敏, 于涛, 等, 2010.

戈壁地表风沙运动特征的野外观测研究

[J]. 新疆环境保护, 32(3):10-13.

[本文引用: 1]

尹永顺, 1990.

砾漠大风地区风沙流及铁路防沙技术研究

[J]. 铁道建筑, 30(5):5-8.

[本文引用: 1]

张克存, 屈建军, 俎瑞平, 等, 2005.

戈壁风沙流结构和风速廓线特征研究

[J]. 水土保持研究, 12(1):54-55.

[本文引用: 1]

张正偲, 董治宝, 赵爱国, 2011.

人工模拟戈壁风沙流与风程效应观测

[J]. 中国科学:地球科学, 41(10):1505-1 510.

[本文引用: 1]

邹学勇, 董光荣, 王周龙, 1995.

戈壁风沙流若干特征研究

[J]. 中国沙漠, 15(4):368-373.

[本文引用: 1]

本文从理论推导入手,阐述了戈壁风沙流中沙粒起跃角度、跃移高度和风沙流结构、强度、能量分布状况等特征,指出与流沙地表风沙流的显著差异,并得到风洞实验验证。

BAAS A C W, 2006.

Wavelet power spectra of aeolian sand transport by boundary layer turbulence

[J]. Geophysical Research Letters, 33(5):L05403. DOI: 10.1029/2005GL025547.

[本文引用: 2]

BAAS A C W, 2008.

Challenges in aeolian geomorphology: Investigating aeolian streamers

[J]. Geomorphology, 93(1/2):3-16. DOI: 10.1016/j.geomorph.2006.12.015.

[本文引用: 1]

BUTTERFIELD G R, 1998.

Transitional behaviour of saltation: Wind tunnel observations of unsteady winds

[J]. Journal of Arid Environments, 39(3):377-394. DOI: 10.1006/jare. 1997.0367.

[本文引用: 1]

CHAMECKI M, KOK J F, 2023.

Fundamental causes of model inaccuracies in predicting wind-blown sand fluxes

[J]. Geophysical Research Letters, 50(11):e2023GL103490. DOI: 10.1029/2023GL103490.

[本文引用: 1]

CHENG J J, JIANG F Q, XUE C X, et al, 2015.

Characteristics of the disastrous wind-sand environment along railways in the Gobi area of Xinjiang, China

[J]. Atmospheric Environment, 102:344-354. DOI: 10.1016/j.atmosenv.2014.12.018.

[本文引用: 2]

COLOMINAS M A, SCHLOTTHAUER G, TORRES M E, 2014.

Improved complete ensemble EMD: A suitable tool for biomedical signal processing

[J]. Biomedical Signal Processing and Control, 14:19-29. DOI: 10.1016/j.bspc.2014. 06.009.

[本文引用: 2]

DAVIDSON-ARNOTT R G D, BAUER B O, WALKER I J, et al, 2012.

High-frequency sediment transport responses on a vegetated foredune

[J]. Earth Surface Processes and Landforms, 37(11):1227-1 241. DOI: 10.1002/esp.3275.

[本文引用: 2]

DE SOUZA U B, ESCOLA J P L, DA CUNHA BRITO L, 2022.

A survey on Hilbert-Huang transform: Evolution, challenges and solutions

[J]. Digital Signal Processing,120:103292. DOI: 10.1016/j.dsp.2021.103292.

[本文引用: 4]

ELLIS J T, SHERMAN D J, 2023.

Cross-wavelet analysis of coherent wind and saltation events

[J]. Earth Surface Processes and Landforms, 48(2):406-414. DOI: 10.1002/esp.5493.

[本文引用: 1]

HUANG N E, SHEN Z, LONG S R, et al, 1998.

The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis

[J]. Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, 454(1971):903-995. DOI: 10.1098/rspa.1998.0193.

[本文引用: 1]

HUANG N E, WU Z H, 2008.

A review on Hilbert-Huang transform: Method and its applications to geophysical studies

[J]. Reviews of Geophysics, 46(2):RG2006. DOI: 10.1029/2007RG000228.

[本文引用: 1]

LIU B L, WANG Z Y, NIU B C, et al, 2021.

Large scale sand saltation over hard surface: A controlled experiment in still air

[J]. Journal of Arid Land, 13(6):599-611. DOI: 10.1007/s40333-021-0104-3.

[本文引用: 1]

Saltation is the major particle movement type in wind erosion process. Saltating sand grains can rebound up to tens of times larger in length and height over hard surface (such as gravel surface) than over loose sand surface. Gravels usually have different faces, causing distinct response of the impacting grains, but the effects of the grain and gravel-surface contact angle on grain rebound are not yet well quantified. We performed full-range controlled experiments of grain saltation using different contact angles, grain sizes and impact speeds in still air, to show that contact angle increases the height of representative saltation path but decreases particle travel length. The results were compared with outputs from the COMprehensive numerical model of SALTation (COMSALT). Large saltation height of 4.8 m and length of 9.0 m were recorded. The maximum and representative saltation height over the gravel surface were found to be about 4.9 times and 12.8 times those over the loose sandy surface, respectively. The maximum saltation length may be reduced by 58% and the representative saltation height may be increased by 77% as contact angle increases from 20° to 40°. We further showed that the collision inertia contributes 60% of the saltation length, and wind contributes to the other 40%. These quantitative findings have important implications for modeling saltation trajectory over gravel surface.

LIU J Q, KIMURA R, 2018.

Wind speed characteristics and blown sand flux over a gravel surface in a compact wind tunnel

[J]. Aeolian Research, 35:39-46. DOI: 10.1016/j.aeolia.2018.09.005.

[本文引用: 1]

MARUSIC I, MONTY J P, 2019.

Attached eddy model of wall turbulence

[J]. Annual Review of Fluid Mechanics, 51(1):49-74. DOI: 10.1146/annurev-fluid-010518-040427.

[本文引用: 1]

SHERMAN D J, LI B L, ELLIS J T, et al, 2018.

Intermittent aeolian saltation: A protocol for quantification

[J]. Geographical Review, 108(2):296-314. DOI: 10.1111/gere.12249.

[本文引用: 2]

SHU K, WANG W J, DING H H, et al, 2022.

Influence of sand transport rate on rolling wear and damage behaviors of wheel/rail in Gobi and desert windblown sand environments

[J]. Tribology International,172:107584. DOI: 10.1016/j.triboint.2022.107584.

[本文引用: 1]

SPIES P J, MCEWAN I K, BUTTERFIELD G R, 2000.

One-dimensional transitional behaviour in saltation

[J]. Earth Surface Processes and Landforms, 25(5):505-518.

[本文引用: 1]

STOUT J E, ZOBECK T M, 1997.

Intermittent saltation

[J]. Sedimentology, 44(5):959-970. DOI: 10.1046/j.1365-3091.1997.d01-55.x.

[本文引用: 2]

TAN L H, AN Z S, ZHANG K, et al, 2020.

Intermittent aeolian saltation over a gobi surface: Threshold, saltation layer height, and high-frequency variability

[J]. Journal of Geophysical Research: Earth Surface, 125(1):e2019JF005329. DOI: 10.1029/2019JF005329.

[本文引用: 9]

TAN L H, QU J J, WANG T, et al, 2021.

Field observation evidence for kink points in the vertical kinetic energy flux profiles of wind-blown sand over gobi and its significance

[J]. Geophysical Research Letters, 48(3):e2020GL091224. DOI: 10.1029/2020GL091224.

[本文引用: 1]

TAN L H, QU J J, WANG T, et al, 2022.

Vertical flux density and frequency profiles of wind-blown sand as a function of the grain size over gobi and implications for aeolian transport processes

[J]. Aeolian Research,55:100787. DOI: 10.1016/j.aeolia.2022.100787.

[本文引用: 1]

TAN L H, ZHANG W M, AN Z S, et al, 2019.

The effect of roughness density of gobi beds on the entrainment of sediment by wind: A wind tunnel study

[J]. Journal of Wind Engineering and Industrial Aerodynamics, 190:183-189. DOI: 10.1016/j.jweia.2019.05.002.

[本文引用: 1]

TAN L H, ZHANG W M, QU J J, et al, 2013.

Aeolian sand transport over gobi with different gravel coverages under limited sand supply: A mobile wind tunnel investigation

[J]. Aeolian Research, 11:67-74. DOI: 10.1016/j.aeolia.2013. 10.003.

[本文引用: 1]

TAN L H, ZHANG W M, QU J J, et al, 2016.

Aeolian sediment transport over gobi: Field studies atop the Mogao Grottoes, China

[J]. Aeolian Research, 21:53-60. DOI: 10.1016/j.aeolia.2016.03.002.

[本文引用: 1]

WANG P, ZHENG X J, 2015.

Unsteady saltation on Mars

[J]. Icarus, 260:161-166. DOI: 10.1016/j.icarus.2015.07.012.

[本文引用: 1]

WANG T, QU J J, TAN L H, 2023.

Aeolian sediment transport over sandy gobi: Field studies in the Nanhu gobi along the Hami-Lop Nor Railway

[J]. International Soil and Water Conservation Research, 11(1):125-134. DOI: 10.1016/j.iswcr.2022.03.009.

[本文引用: 5]

WANG T, QU J J, TAN L H, et al, 2022.

Aeolian sediment transport over the Gobi with high gravel coverage under extremely strong winds in the Hundred Miles windy area along the Lanzhou-Xinjiang High-Speed Railway

[J]. Journal of Wind Engineering and Industrial Aerodynamics,220:104857. DOI: 10.1016/j.jweia.2021.104857.

[本文引用: 4]

WANG Z T, ZHANG C L, CEN S B, 2018.

Unsteady aeolian saltation

[J]. The European Physical Journal E, 41(10):121. DOI: 10.1140/epje/i2018-11730-8.

[本文引用: 4]

WANG Z T, 2024.

Some statistical properties of aeolian saltation

[J]. AIP Advances, 14(2):025342. DOI: 10.1063/5.019 2219.

[本文引用: 6]

WANG Z T, ZHANG C L, WANG H T, 2014.

Intermittency of aeolian saltation

[J]. The European Physical Journal E, 37(12):126. DOI: 10.1140/epje/i2014-14126-x.

[本文引用: 2]

WEI W, ZHANG H S, SCHMITT F G, et al, 2017.

Investigation of turbulence behaviour in the stable boundary layer using arbitrary-order Hilbert spectra

[J]. Boundary-Layer Meteorology, 163(2):311-326. DOI: 10.1007/s10546-016-0227-9.

[本文引用: 1]

ZHANG C L, WEI G R, ZOU X Y, et al, 2022.

The varying fetch effect of aeolian sand transport above a gobi surface and its implication for gobi development process

[J]. International Soil and Water Conservation Research, 10(4):623-634. DOI: 10.1016/j.iswcr.2022.03.002.

[本文引用: 3]

/

〈 〉



陇ICP备09004376

Copyright © 2019 《干旱气象》 编辑部

地址： 甘肃省兰州市东岗东路2070号，中国气象局兰州干旱气象研究所 730020

电话： 0931-2402270、0931-2402775   Email：ghqx@iamcma.cn、ghs_ghqx@sina.com

技术支持: 北京玛格泰克科技发展有限公司

访问总数: 当日访问总数: 当前在线人数: