CN111914396B - 基于高分辨率dem数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法 - Google Patents
基于高分辨率dem数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法 Download PDFInfo
- Publication number
- CN111914396B CN111914396B CN202010612344.3A CN202010612344A CN111914396B CN 111914396 B CN111914396 B CN 111914396B CN 202010612344 A CN202010612344 A CN 202010612344A CN 111914396 B CN111914396 B CN 111914396B
- Authority
- CN
- China
- Prior art keywords
- solar
- grid point
- terrain
- angle
- mode
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000005855 radiation Effects 0.000 title claims abstract description 119
- 238000000034 method Methods 0.000 title claims abstract description 68
- 230000000694 effects Effects 0.000 title claims abstract description 25
- 230000004907 flux Effects 0.000 claims abstract description 81
- 238000012876 topography Methods 0.000 claims description 34
- 238000012937 correction Methods 0.000 claims description 31
- 238000012935 Averaging Methods 0.000 claims description 8
- 239000002131 composite material Substances 0.000 claims description 7
- 238000002834 transmittance Methods 0.000 claims description 7
- 230000005540 biological transmission Effects 0.000 abstract description 12
- 230000008901 benefit Effects 0.000 abstract description 9
- 230000007246 mechanism Effects 0.000 abstract description 9
- 230000000763 evoking effect Effects 0.000 abstract description 7
- 238000004364 calculation method Methods 0.000 description 28
- 230000008569 process Effects 0.000 description 11
- 238000004422 calculation algorithm Methods 0.000 description 8
- 230000001965 increasing effect Effects 0.000 description 8
- 238000003860 storage Methods 0.000 description 6
- 238000011160 research Methods 0.000 description 5
- 230000000295 complement effect Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000001556 precipitation Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 241000933095 Neotragus moschatus Species 0.000 description 1
- 229910018503 SF6 Inorganic materials 0.000 description 1
- 241000219094 Vitaceae Species 0.000 description 1
- 238000002679 ablation Methods 0.000 description 1
- 230000001154 acute effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000003851 biochemical process Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011438 discrete method Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000003344 environmental pollutant Substances 0.000 description 1
- 238000001704 evaporation Methods 0.000 description 1
- 230000008020 evaporation Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 235000021021 grapes Nutrition 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000000873 masking effect Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 230000029553 photosynthesis Effects 0.000 description 1
- 238000010672 photosynthesis Methods 0.000 description 1
- 231100000719 pollutant Toxicity 0.000 description 1
- 230000000191 radiation effect Effects 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- SFZCNBIFKDRMGX-UHFFFAOYSA-N sulfur hexafluoride Chemical compound FS(F)(F)(F)(F)F SFZCNBIFKDRMGX-UHFFFAOYSA-N 0.000 description 1
- 229960000909 sulfur hexafluoride Drugs 0.000 description 1
- 101150053276 sunI gene Proteins 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J1/00—Photometry, e.g. photographic exposure meter
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
- G06T17/205—Re-meshing
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Photometry And Measurement Of Optical Pulse Characteristics (AREA)
Abstract
本发明公开了基于高分辨率DEM数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法，包括如下步骤：在模式运行前提取太阳辐射强迫效应因子。利用高分辨率DEM数据计算模式次网格尺度的地形辐射强迫效应参量，再将这些次网格尺度的地形辐射强迫参量聚合到模式网格尺度，得到模式网格尺度的地形辐射强迫效应参量；在模式运行时，修正太阳入射辐射各分量。利用在模式网格尺度的地形辐射强迫效应参量和模式本身计算的模式网格尺度变量，根据山地三维地表辐射传输理论模型计算得到考虑地形辐射强迫效应后的模式网格尺度太阳直接辐射通量、太阳散射辐射通量和周围地形反射太阳辐射通量。该方法具有输入简单、物理机制明确清晰、运算量经济和可移植性强等优点。
Description
技术领域
本发明涉及一种基于DEM数据的次网格地形太阳山地辐射强迫效应快速参数化方法，尤其涉及一种基于高分辨率DEM数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法。
背景技术
入射太阳辐射通量对近地面温度、土壤墒情、蒸发、冰雪消融、光合作用、污染物生消等物理生物化学过程具有十分重要的作用，因此对入射太阳辐射通量的研究一直是气候、农业、生态、水文和环境等学科的重要研究内容。不同于气温、降水等气象资料的相对易于获取，入射太阳辐射的直接观测在地形复杂地区一直处于缺乏状态，辐射测量站点的空间密度和观测频次都不能很好地满足科学研究和社会生产的需求。除辐射测量仪器布设成本较高和日常维护难度较大等原因外，更为重要的原因是很难利用站点资料通过插值方法得到令人满意的山地地表入射太阳辐射格点资料。对于理想平坦地表，入射太阳辐射通量在大范围内保持均匀，单点测量值可以代表整个区域。地形是重要性仅次于云的辐射强迫因子，地形高度、坡度、坡向、天空可见程度和遮蔽等特征因子对入射太阳辐射通量具有显著作用，入射太阳辐射通量水平梯度在山地显著增加，山地辐射站点空间代表性很差。这就形成了矛盾：一方面是山区地表辐射通量复杂多变难以测量，其时空分布难以很好地展现；另一方面是山区地表辐射通量对山区地表其它过程错综复杂的影响具有重要研究意义，亟需辐射资料支撑。为了解决这个矛盾，数值模式模拟和卫星资料反演逐渐成为山地地表辐射资料的两大重要来源。
自20世纪50年代以来，国内外涌现出大量计算山地辐射通量各分量的计算方法。其中既有基于统计的经验算法，也有基于二维无限长地形假设的模型算法，还有基于三维地表辐射传输理论的模型算法。统计经验算法虽然普遍计算简单，但物理过程机制模糊，往往只适用于特定条件下的某一地区，算法的可移植性非常差。基于二维无限长地形假设的算法仅考虑了坡度、坡向对地表太阳入射角度的影响，即计算目标点自身地形对太阳入射辐射通量的影响。基于三维地表辐射传输理论的模型算法虽然计算较为复杂，但通常考虑了地形坡度、坡向、遮蔽、天空可见程度和地形配置等因素对辐射通量的综合影响，具有物理机制明晰和可移植性强的优点，目前在卫星资料反演山地辐射资料工作中已取得了广泛应用。
研究表明在大气数值模式中增加对山地辐射强迫效应对改善模式性能的重要性随着模式分辨率提高而提高。一些研究将山地辐射强迫效应引入到天气/气候数值模式(如MM5、WRF、GRAPES、RegCM4等)中，有效改善了这些数值模式在地形复杂地区及其周边区域的地表降水、气温和风速等气象要素的模拟性能。虽然地形辐射强迫效应对天气/气候过程及其数值模拟的重要性已经充分彰显，但绝大多数天气/气候数值模式普遍忽略了这一物理机制。究其原因，很大程度上是由于缺少物理机制明确、运算经济并且可移植性强的地形辐射强迫效应参数化方案。
现有山地辐射效应参数化方法存在一些缺点：(1)有的方法中物理机制描述不够明晰充分，如采用山地辐射经验统计公式、或采用无限长地形假设仅考虑地形对太阳入射角的影响。(2)有的方法虽然物理机制明晰，但运算过程十分复杂，客观上加重了天气/气候数值模式运算和存储负担，并且在模式间的移植过程复杂，如采用蒙特卡洛光子追踪技术。(3)有的方法采用中、低分辨率高程数据描述模式次网格地形，在天气/气候数值模式运行时计算考虑地形影响因子的次网格尺度地表辐射通量后再平均得到模式网格尺度地表辐射通量，在模式分辨率较低的情况下或许可行。但目前各类天气/气候数值模式的水平分辨率普遍提高到100～102km，需要采用高分辨率高程数据描述次网格地形。数字高程模型(Digital Elevation Model，DEM)数据的水平分辨率已经达到区域100m、全球范围101m，可以很好的表现目前天气/气候数值模式的网格尺度地形，然而在模式运行过程中计算水平分辨率100～101m的次网格地表辐射通量显然是不可行的。
发明内容
本发明所要解决的技术问题是在模式运行过程中计算高分辨率次网格地表辐射通量。
为了解决这一技术问题，本发明提供了一种基于高分辨率DEM数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法。包含如下步骤：在模式运行前，利用高分辨率DEM数据计算次网格点i坡度αi和坡向βi；利用次网格点坡度αi和坡向βi计算次网格点的cosαi、sinαicosβi和sinαisinβi；将模式网格点P在X轴和Y轴方向前后设定距离范围内所有次网格点i的cosαi、sinαicosβi和sinαisinβi分别求平均获得模式网格点P的地形修正量up、vp和wp；利用高分辨率DEM数据、次网格点坡度αi和坡向βi计算次网格点天空可见因子SVFi，并计算次网格点的SVFi+SVFi cosαi，将模式网格点P在X轴和Y轴方向前后设定距离范围内所有次网格点的SVFi+SVFi cosαi求平均获得模式网格点P的天空可见因子和坡度的复合修正量(SVF+SVFcosα)p；利用高分辨率DEM数据计算次网格点i的k个搜索方位中所有方位的最大地形遮蔽角∠εi，k的正弦值sinεi，k，组成所有次网格点i在不同搜索方位上的最大地形遮蔽角∠εi，k的正弦值数据集，模式网格点P在X轴和Y轴方向前后设定距离范围内所有次网格点i的所有搜索方位最大地形遮蔽角∠εi.k的正弦值数据分别与同搜索方位太阳高度角正弦值sin h(太阳高度角正弦值sin h等于太阳天顶角余弦值cos Z)进行比对判定，筛选出所有搜索方位中未被遮蔽的次网格点i，将任一搜索方位中未被遮蔽的次网格点i数和同搜索方位次网格点i总数的比值作为模式网格点P在该搜索方位的地形遮蔽系数SFp，将所有模式网格点P在k个搜索方位和设定的所有太阳高度角正弦值sin h下地形遮蔽系数SFp组成四维数组R(a，b，k，m)，a和b分别表示P点在X轴方向和Y轴方向的序号，k表示搜索方位的序号，m表示设定的太阳高度角正弦值的序号；在模式运行时，利用模式网格点P的太阳赤纬σp、纬度太阳时角ωp和太阳天顶角Zp计算模式网格点P的太阳方位角θP；利用模式网格点P的地形修正量up、vp和wp以及太阳方位角θP和太阳天顶角Zp计算模式网格点P的太阳入射角余弦值(cosγ)p；利用模式网格点P的太阳方位角θP和太阳天顶角余弦值cos Zp分别确定k和m序号数，从四维数组R(α，b，k，m)读取当前时刻模式网格点P的遮蔽系数SFP；利用太阳常数Eac、日地距离修正因子/>模式网格点P的大气透射率/>该时刻对应的地形遮蔽系数SFp和太阳入射角余弦值(cosγ)p计算获得该时刻模式网格点P考虑地形作用的太阳直接辐射通量/>利用太阳常数Eac、模式网格点P不考虑地形作用的太阳直接辐射通量/>及太阳散射辐射通量/>考虑地形作用的太阳直接辐射通量/>天空可见因子与坡度的复合修正量(SVF+SVFcosα)p和该时刻对应的太阳入射角余弦值(cosγ)p计算获得该时刻模式网格点P考虑地形作用的太阳散射辐射通量/>利用模式网格点P不考虑地形作用的太阳直接辐射通量/>太阳散射辐射通量/>地表反照率ρp和地形修正量up计算获得模式网格点P考虑地形作用的周围地形反射的辐射通量/>
上述技术方案中，所述
其中H为高分辨率DEM数据的高程值。
上述技术方案中，所述运算符< >i→p表示将模式网格P区域内次网格点的参量平均到模式网格点P。
上述技术方案中，所述
上述技术方案中，所述(cosγ)p＝up cos Zp+vp sin Zp cosθp+wp sin Zp sinθp。
上述技术方案中，所述
其中φ分别为次网格点i四周地形的方位角，/>是次网格点i与其φk方位角上最大遮蔽地形点间连线和铅直方向的夹角。
上述技术方案中，所述(SVF+SVFcosα)p＝＜SVFi+SVFi cosαi＞i→p。
上述技术方案中，所述
otherwise
上述技术方案中，所述
上述技术方案中，所述
本方法基于高分辨率DEM数据和山地三维地表辐射传输理论，在模式运行前计算模式网格尺度的地形辐射强迫效应修正量，在模式运行时通过这些修正量修正模式运算的平地地表入射太阳辐射通量各分量。本方法实现了在不影响运算效率和不增加存储负担的前提下，向各类天气/气候数值模式引入描述地形三维地表太阳辐射强迫效应的目的，具有输入简单、物理机制明确清晰、运算量经济和可移植性强等优点。
附图说明
图1.坡地入射太阳辐射示意图；
图2.三阶差分方法示意图；
图3.模式网格点P和次网格点i示意图；
图4.天空可见因子计算示意图；
图5.入射太阳辐射通量各各分量示意图；
图6.SFp(cosγ)p与＜SFicosγi＞i→p之间相对差异的频率分布图；
图7.本发明参数化方法实例计算流程图。
具体实施方式
1.模式运行前地形辐射强迫特征参量的计算
1.1坡度和坡向计算。
如图1所示，坡度α是坡面与水平面的夹角，也是坡面方向量与Z轴的夹角，表征了坡面的陡峭程度，变化范围为0°到90°。坡向β是坡面的朝向，定义为坡面的法向量/>在水平面上的投影/>与正北方向N的夹角，沿顺时针方向为正，变化范围为0°到360°，0°或360°是北坡，90°是东坡，180°是南坡，270°是西坡。坡度和坡向可以通过海拔高度H的在水平方向的梯度/>和/>求解得到：
以正东方向和正北方向分别作为高分辨DEM数据的x轴和y轴正方向，利用(1)式和(2)式计算坡度α、坡向β在于对高分辨率DEM数据在x轴方向高程梯度阳y轴方向高程梯度的离散化。通常采用二阶差分方法或者三阶差分方法实现，差分方法的差异对坡度、坡向的影响仅在误差大小上。这里简要介绍高程梯度的三阶差分方法。如图2所示，采用3×3虚拟移动窗格，利用z5点周围8个格点高程计算z5点处高程在x轴、y轴方向的梯度：
在高分辨率DEM数据矩阵中连续移动3×3虚拟窗格，遍历DEM数据网格点后完成坡度、坡向的计算。
1.2山地太阳入射角计算。
如图1所示，山地某点处的太阳入射角γ定义为太阳入射光线I与坡面法向量之间的夹角(锐角)。IZ和/>分别是入射光线I在z坐标轴和坡面法向量/>上的投影，Isploe是到达坡面的太阳光线辐射量：
Isploe＝I cosγ (4)
Isploe＝I[sin h cosα+sinαcos h cos(θ-β)] (6)
其中h为太阳高度角，θ为太阳方位角，α为坡度，β为坡向。这里定义太阳方位角θ和坡向β都以正北方向为0°，沿顺时针方向增大。
cosγ＝sin h cosα+sinαcos h cos(θ-β) (7)
由于太阳天顶角Z和太阳高度角h互为余角，展开(7)式得到：
cosγ＝cos Z cosα+sinαcosβsin Z cosθ+sinαsinβsin Z sinθ (8)
如图3所示，模式网格点记为P，模式网格水平分辨率记为Res，在P点x和y方向前后0.5分辨率距离范围内共有n个次网格点i。将n个次网格点i的太阳入射角余弦cosγi平均得到网格点P的太阳入射角余弦值(cosγ)p：
(cosγ)p＝＜cosγi＞i→p (9)
其中运算符< >i→p表示将次网格点i的参量平均到模式网格点P。
根据(8)式有cosγi＝cos Zi cosαi+sinαicosβisin Zi cosθi+sinαisinβisin Zisinθi，代入(9)式得到：
(cosγ)p＝＜cos Zi cosαi+sinαi cosβi sin Zi cosθi+sinαi sinβi sin Zi sinθi＞i→p (10)
在模式水平网格分辨率尺度范围内(Res量级100～102km)，模式网格P区域内的次网格点i的太阳天顶角Zi(太阳方位角θi)和模式网格点P的太阳天顶角Zp(太阳方位角θp)之间的差异可以忽略不计，即和/>因此(10)式可以变换为：
(cosγ)p＝up cos Zp+vp sin Zp cosθp+wp sin Zp sinθp (11)
其中
天气/气候数值模式运行过程中都会计算模式网格点P的太阳天顶角ZP，但通常并不计算太阳方位角θP。在模式运算过程中，可以根据(13)式计算模式网格点P的太阳方位角θP：
其中σp、ωp和Zp分别是模式网格点P的太阳赤纬、纬度、太阳时角和太阳天顶角，均为模式可以提供的参量，计算得到的太阳方位角以正北方向为0°，顺时针增大。
模式运行前，利用DEM数据计算模式次网格尺度的cosαi、sinαicosβi和sinαisinβi，并将它们平均到模式网格尺度得到up、vp和wp。模式运行时，根据(13)式计算模式网格点的太阳方位角θp后，再根据(11)式计算得到模式网格点的太阳入射角余弦值(cosγ)p。
1.3天空可见因子计算
天空可见因子(Sky View Factor，SVF)定义为研究目标点上方半球可见部分占整个半球的比例，是具有三维地表结构的陆地表面的重要表面几何参数，它常被用来计算地表太阳散射辐射通量和大气向下长波辐射通量。利用DEM数据计算SVF的离散化方法较多，比较多种SVF计算方法后，本发明采用如下SVF计算方法：
其中α、β和φ分别为目标点O的纬度、坡度、坡向和其四周地形的方位角。如图4所示，/>是目标点与φk地形方位角上最大遮蔽地形点间连线和水平面的夹角，即最大遮蔽地形仰角。/>是目标点与φk方位角上最大遮蔽地形点间连线和铅直方向的夹角，/>和/>互为余角。SVF的计算精度与搜索方向数目ktotal和最大搜索半径Rd有关，ktotal和Rd的数值越大，SVF精确越高，但耗费的运算资源也越多。最大搜索半径Rd建议设置为20～30km，搜索方向数目ktotal建议设置为200至360个。
在模式运行前利用高分辨率DEM数据计算各次网格点i的天空可见因子SVFi，并计算次网格点的SVFi+SVFicosαi，平均得到模式网格点P的天空可见因子和坡度的复合修正量：(SVF+SVFcosα)p＝＜SVFi+SVFicosαi＞i→p。
1.4遮蔽系数的计算
遮蔽系数(Shading Factor，SF)表征目标点是否被太阳方位角θ方向地形所遮蔽，如图4所示，当目标点与太阳方位角方向的最大遮蔽地形仰角大于太阳高度角h，即 时，目标点被地形遮蔽，遮蔽系数定义为0，否则目标点没有被地形遮蔽，遮蔽系数定义为1。遮蔽系数同太阳高度角和太阳方位角都有关系，是时间的变量。目前天气/气候数值模式中辐射传输模块更新频率通常是每小时数次，如果在模式次网格上计算对应时次的遮蔽系数，会造成计算和存储冗余。因此不能在模式次网格尺度上直接计算遮蔽系数，而应该计算蕴含次网格地理信息的模式网格尺度的遮蔽系数。
通过如下方法获取在模式格点上不随时间变化的地形遮蔽情况基础数据集。
在模式运行前，首先利用高分辨率DEM数据计算各次网格点i上以该点为中心第k个搜索方位上最大地形遮蔽角∠εi，k的正弦值sinεi，k。搜索方位如图4所示：搜索的方位角φi，k与太阳方位角定义方式一致：同以正北方向为0°，顺时针旋转增加。对所有方位进行计算，这样我们可以得到所有次网格点在不同搜索方位上的最大地形遮蔽角∠εi，k的正弦值数据集，ktotal为搜索方位总数。如ktotal＝360，则以每1°为一个方位进行搜索，搜索的方位角与太阳方位角定义方式一致：同以正北方向为0°，顺时针旋转增加。搜索的方位总数ktotal越多、最大搜索半径Rd越大，最后计算的模式网格点遮蔽系数SFP精度越高，但耗费的计算资源也随之增长。建议最大搜索半径Rd设置为20～30km。搜索方位数ktotal的最小值设置为天气/气候模式中1天(24小时)内太阳辐射通量计算频次的2倍，假设太阳辐射通量在模式中每Rt分钟计算一次，则其中nint()为四舍五入取整函数。
进而通过统计获取地形遮蔽数据集：设模式格点Pa，b中有n个次网格点，a为x方向序号为和b为y方向序号。由于最大地形遮蔽角正弦值介于0.0～1.0之间的，将0.0～1.0这个范围等分为101份。根据存储条件，可以是1001份，10001份....，分级越多则计算越精确。因而每份长度为0.01，这样有0，0.01，0.02....，1共计101个等级。为了统计和存储方便，将sinεi，k扩大100倍取整，这样变成了0，1，2....100。基于上面获取的所有次网格点不同搜索方位上最大地形遮蔽角正弦值数据集，我们在第k个搜索方位上分别统计最大地形遮蔽角正弦值在某个界限下的次网格格点数，设n1，n2，...，n101分别对应nint(100sinεi，k)＝0，nint(100sinεi，k)≤1，...，nint(100sinεi，k)≤100(i＝1，2，...，n)的次网格格点数，将不同界限下统计得到的次网格格点数除以模式格点Pa，b中总次网格点数n：表示模式格点Pa，b上在太阳方位角等于或最接近于第k个搜索方位角下当nint(100sin h)＝0，nint(100sin h)＝1，...，nint(100sin h)＝100(sinh为太阳高度角正弦值，sinh＝cosZ，Z为太阳天顶角)时没有被地形遮蔽的次网格点数占该模式格点中总次网格点数的比值。将/> 存放在数组R(a，b，k，m)中(m＝1，2，...，101)，这样就可以得到模式格点Pa，b第k个搜索方位上当太阳方位角等于或最接近于该搜索方位角时不同太阳高度角下地形遮蔽情况的基础数据。对所有搜索方位进行类似统计，就可以得到模式格点Pa，b处所有搜索方位上当太阳方位角等于或最接近于该搜索方位角时不同太阳高度角下的地形遮蔽情况。采用同样方法对所有模式格点进行统计，这样我们就可以得到所有模式格点不同太阳方位角和太阳高度角下地形遮蔽情况基础数据，存放在四维数组R(Inum，Jnum，Ktotal，101)中，Inum和Jnum分别为模式在x和y方向的格点总算，ktotal为搜索方位总数，101为最大地形高度角正弦值分级数。
2.模式运行时山地入射太阳辐射通量各分量修正的举例说明
水平面上某点入射太阳辐射通量E＝Edir+Edif包含两个分量：太阳直接辐射通量Edir和太阳散射辐射通量Edif。如图5所示，山地入射太阳辐射通量Et＝Edir，t+Edif，t+Eref，t包含三个分量：山地太阳直接辐射通量Edir，t、山地太阳散射辐射通量Edif，t和周围地形反射辐射通量Eref，t。
以下利用《晴空下山区太阳辐射模拟》文中介绍的山地太阳直接辐射通量Fdir，t、山地太阳散射辐射通量Edif，t和周围地形反射太阳辐射通量Eref，t计算公式，以实例方式说明考虑地形辐射强迫效应后模式网格尺度的入射太阳辐射通量修正方法。
2.1山地太阳直接辐射通量修正的举例说明
不考虑地形作用，某点的太阳直接辐射通量可以表示为：
其中Z和τb分别是太阳常数、日地距离修正系数、太阳天顶角和大气透射率。Eac为常数。/>与模式网格点或次网格点的空间位置无关，只是时间的函数，模式已经计算得到。在模式水平网格分辨率尺度范围内(量级100～102km)，模式次网格点i的太阳天顶角Zi(大气透射率/>)和模式网格点P的太阳天顶角ZP(大气透射率/>)之间的差异可以忽略不计，即/>和/>所以当不考虑地形作用时，可以认为模式网格尺度和次网格尺度的太阳直接辐射通量没有差异，即/>
考虑地形作用后，次网格点i的太阳直接辐射通量可以表示为：
其中SFi、γi分别是次网格点i的遮蔽系数、坡面太阳入射角。则模式网格点P处的太阳直接辐射通量：
其中和SFp分别是太阳常数、日地距离修正系数、模式网格点P的大气透射率、太阳入射角余弦值和地形遮蔽系数。Eac为常数，/>仅是时间的函数且模式已经计算。某一时刻模式格点P处的太阳高度角正弦值sinh(太阳天顶角余弦cosZP和太阳方位角θp模式运算时已经计算得到，SFp从模式运行前计算统计得到的所有模式格点不同太阳方位角和太阳高度角下地形遮蔽情况基础数据R(Inum，Jnum，Ktotal，101)中读取，具体如下：数组标号a和b分别是模式网格点P在x轴和y轴方向的序号，对于格点P这两个序号是确定的。需要确定某个时刻太阳方位角θp下的方位序号k和太阳高度角正弦值sinh下的分级序号m；首先需确定当前时刻P点处与此时太阳方位θp，由公式13求得最为接近的方位序号k，显然/>θp的单位需要由弧度换算成角度(°)。然后根据太阳高度角正弦＝cosZp确定最大地形遮蔽角正弦分级序号m，显然m＝nint(100sin h)+1，这样就可以直接得到该模式格点处的遮蔽系数SFp＝R(a，b，k，m)。
(cosγ)p利用模式运行前计算的模式网格尺度的修正量uP、vP和wp等修正量根据公式(11)：(cosγ)p＝up cos Zp+vp sin Zp cosθp+wp sin Zp sinθp求得，其中某时刻模式某一格点的Zp(模式格点P处太阳天顶角，与太阳高度角互余)和θp(模式格点P处太阳天方位角)已在模式运行时计算得到。
为了验证将＜SFicosγi＞i→p转化为SFp(cosγ)p计算的可行性，利用随机数开展了差异分析。总共开展90000组试验，即90000个样本(相当于有90000个模式网格点)，每组试验包含随机数对(SFi，cosγi)10000个(相当于每个模式网格点拥有10000个次网格点，模式网格分辨率量级取101km，DEM分辨率量级取102m)。随机数对(SFi，cosγi)中的SFi被赋予0或1的随机数值，cosγi被赋予0.0～1.0(包含0.0和1.0)间的随机数值。由10000个点的SFi平均得到SFp、cosγi平均得到(cosγ)p、SFi cosγi平均得到＜SFicosγi＞i→p，计算得到的网格尺度SFp和(cosγ)p数值范围均为0～1。以＜SFicosγi＞i→p作为真值，将90000个样本SFp(cosγ)p与＜SFicosγi＞i→p间的相对差异以0.1％分度值进行频率分布统计后得到图6。可以看到两种计算方式间的相对差异呈正态分布，相对差异在±1％之间的样本数占样本总量的95.5％，相对差异在±0.5％之间的样本数占样本总量的85.0％。因此＜SFicosγi＞i→p转化为SFp(cosγ)p计算是可行的。
2.2山地太阳散射辐射通量修正的举例说明
不考虑地形作用，某点的太阳散射辐射通量可以表示为：
其中Eac、Z、和τd分别是太阳常数、日地距离修正系数、太阳天顶角和大气散射率。/>均与模式网格点或次网格点的空间位置无关。在模式水平网格分辨率尺度范围内(量级100～102km)，模式次网格点i的太阳天顶角Zi(大气散射率/>)和模式网格点P的太阳天顶角ZP(大气散射率/>)之间的差异可以忽略不计，即/>和/>所以当不考虑地形作用时，可以认为模式网格尺度和次网格尺度的太阳散射辐射通量没有差异，即
考虑地形作用后，次网格点i的太阳散射辐射通量可以表示为：
其中SFi、SVFi、γi、αi、Zi、和/>分别是次网格点i的遮蔽系数、天空可见因子、太阳入射角、坡度、太阳天顶角、不考虑地形作用时的太阳直接辐射通量和太阳散射辐射通量。由于/>取次网格点的太阳散射辐射通量平均值得到模式网格点P处的太阳散射辐射通量：
结合公式(17)：
/>
otherwise
可以得到：
其中Eac、及/>(cosγ)p和(SVF+SVFcosα)p分别是太阳常数、模式网格点P不考虑地形作用的太阳直接辐射通量及太阳散射辐射通量、考虑地形作用的太阳直接辐射通量、太阳入射角余弦值和天空可见因子与坡度的复合修正量。Eac为常数。模式网格点P的遮蔽系数SFp和天空可见因子与坡度的复合修正量(SVF+SVFcosα)p在模式运行前已经求取。模式网格点P不考虑地形作用的太阳直接辐射通量/>及太阳散射辐射通量/>在模式运行过程中模式已经计算，太阳入射角余弦值(cosγ)p和考虑地形作用的太阳直接辐射通量/>已经在2.1中求取。
2.3周围地形反射辐射通量修正的举例说明
山地次网格点i接收到周围地形反射的辐射通量可以表示为：
其中ρi、αi、和/>分别是次网格点i的地表反照率、坡度、不考虑地形作用时的太阳直接辐射通量和太阳散射辐射通量。在模式水平网格分辨率尺度范围内(量级100～102km)，模式网格尺度和次网格尺度的的地表反照率的差异可以忽略不计，即/>由于 和/>取次网格点的周围地形反射辐射通量的平均值得到模式网格点P处的周围地形反射辐射通量：
其中和ρp分别是模式计算的模式网格点P上不考虑地形作用的太阳直接辐射通量、太阳散射辐射通量和地表反照率。up为模式运行前根据(12)式计算得到的模式网格尺度修正量。
如上述实例所示(图7)，本参数化方法的实施的主要有两个步骤：
第一步主要在模式运行前开展，提取太阳辐射强迫效应因子提取。利用高分辨率DEM数据计算模式次网格尺度的地形辐射强迫效应参量，如坡度、坡向、太阳入射角余弦值修正量、天空可见因子、最大遮蔽地形仰角正弦值等。再将这些次网格尺度的地形辐射强迫参量聚合到模式网格尺度，得到模式网格尺度的地形辐射强迫效应参量。
第二步主要在模式运行时开展，修正太阳入射辐射各分量。利用模式本身计算的模式网格尺度变量(如大气透射率、日地距离修正因子、地表反照率、太阳天顶角、不考虑地形作用的平面地表上太阳直接辐射通量和太阳散射辐射通量等)和第一步中得到模式网格尺度的地形辐射强迫效应参量，根据山地三维地表辐射传输理论模型计算得到考虑地形辐射强迫效应后的模式网格尺度太阳直接辐射通量、太阳散射辐射通量和周围地形反射太阳辐射通量。
本参数化方法设计思路是在误差允许范围内忽略大气透射率、地表反照率、太阳高度角和太阳方位角等变量在模式网格尺度和次网格尺度的差异，把与时间无关的地形辐射强迫效应参量从山地三维地表辐射传输理论模型中剥离出来，与时间和空间均有关的“将次网格尺度入射太阳辐射通量各分量(和/>)平均到模式网格尺度”问题转化为仅与空间有关的“将次网格尺度地形辐射强迫效应参量平均到模式网格尺度”问题，从而次网格尺度的相关计算均在模式运行前完成，模式运行效率不受影响。
本参数化方法具有四大优点：
(1)具有输入简单易得的优点。唯一的输入数据是高分辨率DEM数据，目前全球范围/准全球范围高精度高分辨率DEM数据(如SRTM90、Aster GDEM、TanDEM-X90、ALOSWorld3D等)非常易于获得。
(2)具有节省运算资源和存储空间的优点。从高分辨率DEM数据中提取模式网格尺度的地形辐射强迫效应信息在天气/气候数值模式运行之前完成，避免在模式运行过程中计算次网格尺度山地太阳入射辐射各分量，不增加模式运算量或内存占用，不影响模式运行效率。
(3)具有物理机制明晰的优点。目前三维地表山地太阳辐射传输理论日趋完善，可供选用的地表入射太阳辐射通量各分量的计算方法众多。相对于二维无限长地形理论或经验性山地辐射计算公式，三维地表山地太阳辐射传输模型充分考虑了地形对地表辐射通量的各种强迫效应。
(4)具有可移植性强的优点。主要体现在三方面：一是并不限定于某种三维地表山地辐射传输模型，大多数三维地表山地辐射传输模型适用于本方法；二是各种高分辨率DEM数据均可作为本方法的输入数据；三是在天气/气候数值模式间的可移植性强，在天气/气候数值模式源代码中增加和修改少量语句(用于读取修正量和修正地表辐射通量各分量)即可实现，不需要大范围改动模式框架、辐射传输模块或陆面过程模块。
Claims (1)
1.基于高分辨率DEM数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法，其特征在于：在模式运行前，利用高分辨率DEM数据计算次网格点i坡度αi和坡向βi；利用次网格点坡度αi和坡向βi计算次网格点的cosαi、sinαi cosβi和sinαisinβi；将模式网格点P在X轴和Y轴方向前后设定距离范围内所有次网格点i的cosαi、sinαicosβi和sinαisinβi分别求平均获得模式网格点P的地形修正量up、vp和wp；利用高分辨率DEM数据、次网格点坡度αi和坡向βi计算次网格点天空可见因子SVFi，并计算次网格点的SVFi+SVFicosαi，将模式网格点P在X轴和Y轴方向前后设定距离范围内所有次网格点的SVFi+SVFicosαi求平均获得模式网格点P的天空可见因子和坡度的复合修正量(SVF+SVFcosα)p；利用高分辨率DEM数据计算次网格点i的k个搜索方位中所有方位的最大地形遮蔽角∠εi,k的正弦值sinεi,k，组成所有次网格点i在不同搜索方位上的最大地形遮蔽角∠εi,k的正弦值数据集，模式网格点P在X轴和Y轴方向前后设定距离范围内所有次网格点i的所有搜索方位最大地形遮蔽角∠εi,k的正弦值数据分别与同搜索方位太阳高度角正弦值进行比对判定，筛选出所有搜索方位中未被遮蔽的次网格点i，将任一搜索方位中未被遮蔽的次网格点i数和同搜索方位次网格点总数的比值作为模式网格点P在该搜索方位的地形遮蔽系数SFp，将所有模式网格点P在k个搜索方位和设定的所有太阳天顶角余弦值cosZp下地形遮蔽系数SFp组成四维数组R(a,b,k,m)，a和b分别表示P点在X轴方向和Y轴方向的序号，k表示搜索方位的序号，m表示设定的太阳天顶角余弦值的序号；在模式运行时，利用模式网格点P的太阳赤纬σp、纬度太阳时角ωp和太阳天顶角Zp计算模式网格点P的太阳方位角θp；利用模式网格点P的地形修正量up、vp和wp以及太阳方位角θp和太阳天顶角Zp计算模式网格点P的太阳入射角余弦值(cosγ)p；利用模式网格点P的太阳方位角θp和太阳天顶角余弦值cosZp分别确定k和m序号数，从四维数组R(a,b,k,m)读取当前时刻模式网格点P的遮蔽系数SFp；利用太阳常数Eac、日地距离修正因子模式网格点P的大气透射率τbp、该时刻对应的地形遮蔽系数SFp和太阳入射角余弦值(cosγ)p计算获得该时刻模式网格点P考虑地形作用的太阳直接辐射通量/>利用太阳常数Eac、模式网格点P不考虑地形作用的太阳直接辐射通量/>及太阳散射辐射通量考虑地形作用的太阳直接辐射通量/>天空可见因子与坡度的复合修正量(SVF+SVFcosα)p和该时刻对应的太阳入射角余弦值(cosγ)p计算获得该时刻模式网格点P考虑地形作用的太阳散射辐射通量/>利用模式网格点P不考虑地形作用的太阳直接辐射通量太阳散射辐射通量/>地表反照率ρp和地形修正量up计算获得模式网格点P考虑地形作用的周围地形反射的辐射通量/>所述
坡度：
坡向：其中H为高分辨率DEM数据的高程值；
所述地形修正量运算符＜＞i→p表示将P网格区域次网格点的参量平均到模式网格点P；
所述
所述(cosγ)p＝up cosZp+vpsinZp cosθp+wp sinZp sinθp；
所述
其中φ为次网格点i四周地形的方位角,是次网格点i与其φk方位角上最大遮蔽地形点间连线和铅直方向的夹角；所述
otherwise
所述
所述
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010612344.3A CN111914396B (zh) | 2020-06-30 | 2020-06-30 | 基于高分辨率dem数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010612344.3A CN111914396B (zh) | 2020-06-30 | 2020-06-30 | 基于高分辨率dem数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111914396A CN111914396A (zh) | 2020-11-10 |
| CN111914396B true CN111914396B (zh) | 2023-09-08 |
Family
ID=73226955
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010612344.3A Active CN111914396B (zh) | 2020-06-30 | 2020-06-30 | 基于高分辨率dem数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111914396B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113591297B (zh) * | 2021-07-27 | 2022-03-29 | 河海大学 | 一种短波辐射与气温降尺度方法、装置及存储介质 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104406686A (zh) * | 2014-12-10 | 2015-03-11 | 西北师范大学 | 复杂地形条件下太阳短波入射辐射估算方法 |
| CN107917881A (zh) * | 2017-11-13 | 2018-04-17 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于冠层内辐射传输机理的光学遥感影像地形校正方法 |
| EP3308442A1 (de) * | 2015-07-30 | 2018-04-18 | Siemens Aktiengesellschaft | Verfahren zur rechnergestützten parametrierung eines umrichters in einem stromnetz |
| CN109446599A (zh) * | 2018-10-12 | 2019-03-08 | 深圳市沃泰凯尔科技开发有限公司 | 应用于大区域水流的计算方法和计算系统 |
| CN109460532A (zh) * | 2018-10-24 | 2019-03-12 | 中国科学院地理科学与资源研究所 | 一种太阳直接辐射遥感计算方法和装置 |
| CN111027175A (zh) * | 2019-11-06 | 2020-04-17 | 中国地质大学（武汉） | 基于耦合模型集成模拟的洪水对社会经济影响的评估方法 |
| CN111260111A (zh) * | 2020-01-08 | 2020-06-09 | 中国科学院上海技术物理研究所苏州研究院 | 基于气象大数据的径流预报改进方法 |
-
2020
- 2020-06-30 CN CN202010612344.3A patent/CN111914396B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104406686A (zh) * | 2014-12-10 | 2015-03-11 | 西北师范大学 | 复杂地形条件下太阳短波入射辐射估算方法 |
| EP3308442A1 (de) * | 2015-07-30 | 2018-04-18 | Siemens Aktiengesellschaft | Verfahren zur rechnergestützten parametrierung eines umrichters in einem stromnetz |
| CN107917881A (zh) * | 2017-11-13 | 2018-04-17 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于冠层内辐射传输机理的光学遥感影像地形校正方法 |
| CN109446599A (zh) * | 2018-10-12 | 2019-03-08 | 深圳市沃泰凯尔科技开发有限公司 | 应用于大区域水流的计算方法和计算系统 |
| CN109460532A (zh) * | 2018-10-24 | 2019-03-12 | 中国科学院地理科学与资源研究所 | 一种太阳直接辐射遥感计算方法和装置 |
| CN111027175A (zh) * | 2019-11-06 | 2020-04-17 | 中国地质大学（武汉） | 基于耦合模型集成模拟的洪水对社会经济影响的评估方法 |
| CN111260111A (zh) * | 2020-01-08 | 2020-06-09 | 中国科学院上海技术物理研究所苏州研究院 | 基于气象大数据的径流预报改进方法 |
Non-Patent Citations (1)
| Title |
|---|
| 顾春雷 ; 方德贤 ; 周洋 ; 黄安宁.次网格坡地辐射参数化对RegCM4.1模式模拟东亚夏季气候的影响.气象科学.2018,585-595. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111914396A (zh) | 2020-11-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110516816B (zh) | 基于机器学习的全天候地表温度生成方法及装置 | |
| CN113297528B (zh) | 一种基于多源大数据的no2高分辨率时空分布计算方法 | |
| US7305983B1 (en) | Assessment of solar energy potential on existing buildings in a region | |
| Fang et al. | A dataset of daily near-surface air temperature in China from 1979 to 2018 | |
| CN109460532B (zh) | 一种太阳直接辐射遥感计算方法和装置 | |
| CN105069295B (zh) | 基于卡尔曼滤波的卫星以及地面降水测量值同化方法 | |
| CN105868529A (zh) | 一种基于遥感的近地表日均大气温度反演方法 | |
| KR101288016B1 (ko) | 지형효과를 적용한 태양복사모델의 구성 방법 및 장치 | |
| Xie et al. | An adjusted two‐leaf light use efficiency model for improving GPP simulations over mountainous areas | |
| Jiang et al. | Space-time mapping of ground-level PM 2.5 and NO 2 concentrations in heavily polluted northern China during winter using the Bayesian maximum entropy technique with satellite data | |
| CN101936777A (zh) | 一种基于热红外遥感反演近地层气温的方法 | |
| Wang et al. | Constructing a gridded direct normal irradiance dataset in China during 1981–2014 | |
| Huang et al. | Evaluation and comparison of MODIS collection 6.1 and collection 6 dark target aerosol optical depth over mainland China under various conditions including spatiotemporal distribution, haze effects, and underlying surface | |
| CN111914396B (zh) | 基于高分辨率dem数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法 | |
| Wang et al. | Evaluation of five planetary boundary layer schemes in WRF over China's largest semi-fixed desert | |
| Liu | Study of topographic effects on hydrological patterns and the implication on hydrological modeling and data interpolation | |
| CN111079835A (zh) | 一种基于深度全连接网络的Himawari-8大气气溶胶反演方法 | |
| Li et al. | Temporal and spatial variations of global solar radiation over the Qinghai–Tibetan Plateau during the past 40 years | |
| Zhang et al. | Integrating a novel irrigation approximation method with a process-based remote sensing model to estimate multi-years' winter wheat yield over the North China Plain | |
| Huang et al. | A novel global grid model for atmospheric weighted mean temperature in real-time GNSS precipitable water vapor sounding | |
| Skirvin et al. | Climate spatial variability and data resolution in a semi-arid watershed, south-eastern Arizona | |
| CN111199092A (zh) | 太阳辐射遥感估算方法、系统及数据处理装置 | |
| Huang et al. | A global grid model for the estimation of zenith tropospheric delay considering the variations at different altitudes | |
| Nemes | A clear sky irradiation assessment using the European Solar Radiation Atlas model and Shuttle Radar Topography Mission database: A case study for Romanian territory | |
| Wang et al. | Spatial Downscaling of Remote Sensing Precipitation Data in the Beijing-Tianjin-Hebei Region |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |