SHG代表什么 SHGYIZHAO什么意思
本文内容来源于《测绘学报》2021年第6期(审图号GS(2021)3565号)
病态乘性误差模型的加权最小二乘正则化迭代解法及精度评定
王乐洋1,陈涛,1 邹传义1,2
1. 东华理工大学测绘工程学院, 江西 南昌 330013;2. 武汉大学测绘学院, 湖北 武汉 430079
基金项目:国家自然科学基金(41874001;41664001);东华理工大学研究生创新专项(DHYC-202020)
摘要:针对乘性误差模型的病态问题,引入Tikhonov正则化方法,导出了病态乘性误差模型的加权最小二乘正则化解。顾及加权最小二乘正则化法在求解病态乘性误差模型时,参数估值与观测值之间存在复杂的非线性关系,本文利用一种无需求导、通过加权的方式便能够计算非线性函数的均值和均方误差阵的比例对称采样的无迹变换(scaled unscented transformation,SUT)法,对病态乘性误差模型进行精度评定。模拟算例和真实算例结果表明,本文提出的加权最小二乘正则化迭代解法可以有效减弱模型的病态性,基于SUT法的精度评定方法能够得到比已有方法更为合理的精度信息,具有较强的适用性。
关键词:病态乘性误差模型 Tikhonov正则化 L曲线法 精度评定 SUT法
引文格式:王乐洋, 陈涛, 邹传义. 病态乘性误差模型的加权最小二乘正则化迭代解法及精度评定[J]. 测绘学报,2021,50(5):589-599. DOI: 10.11947/j.AGCS.2021.20200126
WANG Leyang, CHEN Tao, ZOU Chuanyi. Weighted least squares regularization iteration solution and precision estimation for ill-posed multiplicative error model[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(5): 589-599. DOI: 10.11947/j.AGCS.2021.20200126
阅读全文:http://xb.sinomaps.com/article/2021/1001-1595/2021-5-589.htm
引 言
在大地测量数据处理领域中,若观测误差随观测量的大小或物理信号的强弱而变化,则该类误差为乘性误差[1]。例如,现代观测手段中的合成孔径雷达SAR观测值的随机误差表现为乘性误差[2-3],光电测距EDM(electronic distance measurement)和全球定位系统GPS的观测误差表现为乘性误差或者加乘性混合误差[4-5]。目前,在大地测量领域中关于乘性误差模型参数估计和精度评定的研究成果相对较少[6-9],且其中尚无文献针对乘性误差模型中存在的病态问题进行研究,因此如何对病态乘性误差模型进行参数估计和精度评定是值得研究的问题。
病态问题广泛存在于测绘数据处理中,模型的病态性会造成参数解的不稳定,甚至严重偏离真值[10]。已有的总体最小二乘病态问题[11]、大地测量反演的病态问题[12]和控制网平差的病态问题[13]等都是基于加性误差模型建立的,而对于乘性误差模型的病态问题尚未发现研究成果。乘性误差模型的最小二乘解法最早由文献[6]提出,文献[7]在文献[6]的基础上,总结了乘性误差模型的最小二乘、加权最小二乘和偏差改正加权最小二乘方法,并推导了参数估计的精度评定公式。若直接采用文献[7]中的3种方法处理病态乘性误差模型,而不考虑系数矩阵的病态性,会使参数估值产生偏差且不稳定。文献[4]为避免系数矩阵病态,将模拟点集中在数字地面模型(digital terrain model, DTM)中的峰值处,其他的点均匀分布,这仅可用来处理模拟数据,无法适用于真实的观测数据。系数矩阵病态会使得系数矩阵的列向量之间存在复共线性,进而引起法方程条件数过大,导致解不稳定[10]。目前,处理病态问题的方法主要有Tikhonov正则化法[14-16]和虚拟观测值解法[17]等。正则化参数的确定方法,主要有广义交叉核实法[18]、岭迹法[19]和L曲线法[20-22]等。广义交叉核实法难以获得最优解,适用性弱;岭迹法虽然计算简单,但有一定的主观性;L曲线法适用性强,并且能获得较为合理的正则化参数,被广泛用于正则化参数的选取,可以很好地用于本文算法中正则化参数的选取。
经分析发现,使用加权最小二乘正则化法求解病态乘性误差模型时,参数估值是观测值的非线性函数,而观测值的协因数阵是参数估值的非线性函数,并且在逐步迭代的过程中参数的每一步估值都具有随机性,使得最终的参数估值与观测值为复杂的非线性关系,只能使用数值方法迭代求解,无法得到解析解。因此本文引入无需求导、无需获取非线性函数具体表达式的SUT法[23-24]对病态乘性误差模型进行精度评定。
综上所述,针对已有文献未考虑乘性误差模型的病态性,本文首先利用Tikhonov正则化方法导出病态乘性误差模型的加权最小二乘正则化解,其次通过SUT法对病态乘性误差模型进行精度评定;最后通过两个模拟数据算例和一个真实数据算例验证本文方法的适用性及优势。
1 乘性误差模型
给定一组带有乘性误差的观测值,相应的乘性误差模型可以表示为[7]
(1)
随机模型为
(2)
式中,y∈Rn×1表示观测值向量;β∈Rt×1表示未知参数;fi(•)表示未知参数的函数;1∈Rn×1表示元素全为1的列向量;εm∈Rn×1表示乘性误差向量;Qm∈Rn×n表示乘性误差向量εm的协因数阵;σ02表示单位权方差。
进一步可将式(1)改写为如下形式
(3)
式中,⊙表示哈达玛乘积符号。
本文研究的f(•)是β的线性函数形式,即fi(β)=aiTβ,其中aiT为n×t维向量,则式(3)可以写为
(4)
式中,A=[a1,a2, …,an]T;β=[β1,β2, …, βt]T。
由式(4)可知,在乘性误差模型中,观测值精度与系数矩阵和参数估值的乘积有关:信号Aβ越强,观测值y的误差越大;反之,信号Aβ越弱,误差越小。这明显区别于传统的加性误差模型,即观测值y的误差与系数矩阵和参数估值无关。目前,已有关于乘性误差模型的研究均不考虑系数矩阵A病态[7-8],然而在大地测量实际数据处理中,系数矩阵A很容易出现病态,此时会引起法方程条件数过大,观测量微小的扰动将会引起参数较大的变化,导致解不稳定。例如,文献[25]针对SAR型噪声的特性采用乘性误差模型进行影像去噪,去噪效果达到70%~93%。此时,若系数矩阵A病态,会引起法方程条件数过大,导致解不稳定,进而影响影像去噪。因此需要针对病态乘性误差模型进行更深入的研究来丰富现代大地测量数据处理理论。
2 病态乘性误差模型参数估计的加权最小二乘正则化迭代解法及精度评定
2.1 病态乘性误差模型参数估计的加权最小二乘正则化解法
为了更方便地进行公式推导,将式(4)改写为如下形式
(5)
令
(6)
根据协因数传播率,得到式(6)中e和y的协因数阵为
(7)
根据最小二乘准则求解式(6),得
(8)
β的最小二乘估值为
(9)
根据前文中的分析可知,当法方程N=ATA病态时,由式(9)求得的最小二乘解不可靠。因此,引入正则化因子α,构建如下病态乘性误差模型的正则化准则
(10)
对式(10)中的β求偏导,并令其为0,可得
(11)
将式(6)代入式(11),可得病态乘性误差模型的加权最小二乘正则化解
(12)
式中,
进而可以得到
(13)
式中,为加权最小二乘正则化解;α为正则化参数;In为单位阵。
与文献[7]中的偏差改正加权最小二乘解相比,式(13)右端法方程求逆部分增加了αIn项,有效抑制了法方程的病态性,进而可以得到可靠、稳定的参数解。
2.2 L曲线法基本原理
由式(10)可知,α的函数,L曲线法确定病态乘性误差模型参数估计正则化参数的主要步骤是[20-21]:以为横坐标、为纵坐标作图,得到多组α值作为所求的正则化参数值α。
在确定本文病态乘性误差模型中的正则化参数时,由于观测值的协因数阵为参数估值的非线性函数,因此观测值的权阵也为参数估值的非线性函数,这使得正则化参数值产生动态变化。
2.3 病态乘性误差模型的精度评定
根据文献[8],病态乘性误差模型的单位权方差估值可表示为
(14)
式中,
;n为观测值个数;t为待求参数个数。
由式(13)可知,加权最小二乘正则化解与观测值是线性的关系,由协因数传播率可得的协因数阵为
(15)
由偏差的定义[26],对式(13)两边取期望,可得参数估值的偏差为
(16)
式中,
表示参数真值。
利用矩阵反演公式[27],将正则化解作为参数真值[28-29],进一步可将式(16)表示为
(17)
病态问题中,由于产生了偏差项,此时的均方误差不等于协方差不能作为一个良好的精度评判标准。虽然最小二乘估值是无偏,但已不是最优值[7]。根据均方误差的定义[30-31],可知均方误差的表达式为
(18)
式中,等式右边第1部分表示参数估值协方差的迹;第2部分表示参数估值偏差平方的迹。
2.4 病态乘性误差模型参数估计的加权最小二乘正则化迭代解法
在求解式(13)时,协因数阵Qe是参数估值的非线性函数,而求得的参数估值又是协因数阵Qe的非线性函数,只能使用数值逼近的方法进行迭代求得数值解,而无法得到参数估值的解析解。因此,在由式(13)求解参数估值时,需利用上一次计算得到的参数估值更新协因数阵Qe和正则化参数α,进而求解参数估值,迭代循环直至前后两次参数估值之差小于某一特定的阈值ε,则迭代终止;本文取ε为10-10。
式(13)的迭代公式如下所示
(19)
式中,i为第i次迭代的正则化参数。
综上所述,本文将病态乘性误差模型的加权最小二乘正则化法迭代解法命名为算法1,具体步骤如下:
(1) 采用最小二乘法求得参数估值的初始值
(20)
(2) 求得协因数阵Qei。
(3) 求得正则化参数αi。
(4) 根据式(19)迭代求得病态乘性误差模型加权最小二乘正则化迭代解。
(5) 重复步骤(2)至(4),当
ε时,迭代终止。
3 病态乘性误差模型精度评定的SUT法
由算法1可以看出,病态乘性误差模型的加权最小二乘正则化迭代解法中,每一步迭代参数估值都是观测值的函数,而观测值的协因数阵又是参数估值的函数,使得最终的参数估值与观测值为一个复杂的非线性函数。本文将这种非线性关系表示为[32-33]
(21)
式中,φ(•)表示参数估值与观测值y之间的非线性函数。
第1步参数迭代估值可以表示为[32-33]
(22)
式中,
表示由最小二乘法求得的迭代初值。
最终的参数估值表达式为[32-33]
(23)
式(23)表明,病态乘性误差模型的加权最小二乘正则化法的迭代过程使得参数估值和观测值之间的非线性关系非常复杂,且迭代的过程使参数的每一步的估值都具有随机性。为了避免复杂的公式运算,本文在病态乘性误差模型加权最小二乘正则化法中引入SUT采样法,旨在获得更精确的参数估值和精度信息。
本文将病态乘性误差模型精度评定的SUT法命名为算法2,步骤如下:
(1) 确定采样的先验统计均值和协方差阵。本文将观测值和对应的协方差阵作为采样的先验统计均值E(y)和协方差阵D(y)。
(2) 基于先验统计均值和协方差阵构造采样点集{χi}(i=0, 1, 2, …, 2n),生成2n 1个采样点
(24)
式中,
表示对协方差阵D(y)进行Cholesky分解生成的下三角矩阵;(
)i表示矩阵的第i列;λ=γ2(t κ)-t表示比例参数,取值范围为10-4≤γ≤1[24],本文中γ的取值为10-3;κ表示另一个比例参数,通常取值为3-t或0[24],在本文中取0。
(3) 将采样点集{χi}(i=0, 1, 2, …, 2n)代入式(19)中,得到2n 1个解
。
(4) 确定采样点集的权
(25)
式中,Wim和Wic分别用于均值和协方差矩阵的计算;对于正态分布,ξ=2为最优取值[24]。
(5) 计算参数估值的均值和均方误差阵
(26)
(27)
(6) 计算单位权方差
(28)
本文首先对病态乘性误差模型加权最小二乘正则化迭代解法进行改化,将病态乘性误差模型中参数估值与观测值表示为嵌套函数的形式,然后引入无需求导、无需了解非线性函数构造的SUT法求解病态乘性误差模型中参数估值的均值、均方误差矩阵和单位权方差。
4 算例及分析
为了比较不同方法下参数估值和精度信息的差异,分别使用文献[7]中的未考虑病态性的三种方法即最小二乘法、加权最小二乘法和偏差改正加权最小二乘法,以及本文提出算法1和算法2进行解算。5种方案及对应的方法见表 1。
表 1 5种方案及对应的方法
Tab. 1 Five schemes and corresponding methods
4.1 算例1
为了验证本文公式推导以及算法的有效性,算例1根据文献[34]改编得到,通过利用GPS测量某一区域道路中心线地面点的高程,假定该区域道路中心线地面点高程真值与某一坐标原点的距离除以100符合以下函数模型
(29)
式中,y为地面点高程的真值;x为高程点与某一坐标原点的距离除以100。
本算例中x的取值范围为0~300 m,在x的取值范围内等间距选取31个点,根据式(29)计算出高程点的真值
假设地面点高程真值被乘性误差干扰,其中乘性误差向量εm相互独立,且服从均值为0、标准差为0.1的正态分布,即εm~0, 0.12×I31,相应的乘性误差模型的观测方程为
(30)
式中,Y(x)为受乘性误差干扰的地面高程点观测值向量;1为元素全为1的31维列向量;εm为31维乘性误差向量。
由式(30)得到一组地面高程点的数据模拟值见表 2。为了说明高程点受乘性误差的影响程度,本文将地面高程点受乘性误差干扰前后的数据绘于图 1中。
表 2 数据模拟值
Tab. 2 The simulated data
图 1 地面高程点受乘性误差干扰前后分布
Fig. 1 GPS elevation points before and after disturbed by multiplicative error
本算例选取的单位权中误差σ0为0.3。使用表 1中的5种方案计算,将参数真值、5种方案计算的参数估值、参数估值与真值的2范数及单位权中误差估值列于表 3;5种方案计算的参数估值的均方根误差列于表 4。
表 3 参数估值、参数估值与真值之间的2范数和单位权中误差
Tab. 3 Parameter estimation, 2-norm between the parameter estimation and true value and the unit weight error
表 4 参数估值的均方根误差
Tab. 4 Root mean square error of parameter estimation
由上文分析,本文算法中的正则化参数值会产生动态变化,正则化参数α随迭代次数的变化如图 2所示。经计算,本算例法矩阵N=的条件数为9.896 2×104,具有病态性,使用本文方法求得的法矩阵条件数随迭代次数的变化如图 3所示。
图 2 正则化参数α随迭代次数的变化
Fig. 2 The regularization parameter with the number of iterations
图 3 法矩阵条件数随迭代次数的变化
Fig. 3 The condition number of the normal matrix with the number of iterations
由图 1可知,由于地面高程点受到乘性误差干扰,导致点位产生了严重偏离。由表 3中参数估值与真值的2范数可知,LS方法求得的2范数结果最大;WLS方法由于考虑了观测值的权得到的结果优于LS法,其2范数为6.058 4;bcWLS方法对WLS方法进行了偏差改正,2范数较WLS方法减小了0.119 4,为5.939 0,但由于未考虑模型的病态性,仍与真值偏离较大。本文算法1和算法2求得的2范数分别为0.978 4和0.833 7,参数估值更接近真值,说明本文算法对于降低病态性有一定的效果。其中,由于算法2考虑了正则化解法逐步迭代过程中每一步的随机性和非线性迭代的过程对参数估值的影响,求得的2范数较算法1减小了0.144 7,表明了本文使用SUT法的可行性。由表 3中的单位权中误差可知,由于文献[7]中的方法未考虑模型的病态性,导致求得的单位权中误差偏离真值,而本文算法考虑了这些影响,求得的单位权中误差更接近真值,进一步表明了本文算法的优势。
由表 4可得,文献[7]中3种方法求得的参数估值的均方根误差较大,这是因为法矩阵中存在接近于0的特征值,导致求逆不稳定,进而使得均方根误差结果较大,而本文算法1和算法2求得的均方根误差均比文献[7]中的3种方法小,说明了结果的可靠性。算法2考虑到迭代过程对参数估值的影响,使用SUT法进行精度评定,求得的参数估值的均方根误差最小,其中
4.2 算例2
算例2为一个DTM模型的算例,本文主要考虑模型为病态的情况,因此模拟一个病态DTM模型。文献[4]在模拟DTM模型时,为防止法矩阵病态情况出现,在模拟的DTM模型的每个峰值周围采样更多点,并使峰值点之间远离,虽然可以在一定程度上降低病态性,但仅可用来处理模拟数据,无法适用于真实的观测数据。本文依据文献[4]的思路,使用插值法来模拟DTM模型,插值法是模拟DTM模型的主要方法之一[35-36],本文模拟的DTM模型由下面的插值函数生成
(31)
式中,f(x,y)表示通过插值函数计算得到的高度值;x与y的区域均为0~80 m,间隔为2 m;4个参数估值的真值分别为-1.5,20,5和-4。
函数fi(x,y)(i=1, 2, 3, 4)分别为如下
(32)
(33)
(34)
(35)
模拟的DTM模型如图 4所示,相应的乘性误差模型的观测方程为
(36)
图 4 模拟的DTM模型
Fig. 4 Simulated DTM model
式中,h(x,y)为受乘性误差干扰的观测值向量;1为元素全为1的1681维列向量;εm为1681维乘性误差向量。
在本文中,假设误差向量εm相互独立、且服从均值为0,标准差为0.1的正态分布,即εm~0, 0.12×I1681。为了说明DTM模型受乘性误差的影响程度,本文将不含误差的观测值绘于图 4中,受乘性误差干扰的观测值绘于图 5中。
图 5 受乘性误差干扰的DTM模型
Fig. 5 DTM model disturbed by multiplicative error
本算例选取的单位权中误差σ0为0.3。使用表 1中的5种方案计算,将参数真值、5种方案计算的参数估值、参数估值与真值的2范数及单位权中误差估值列于表 5;5种方案计算的参数估值的均方根误差列于表 6。正则化参数α随迭代次数的变化如图 6所示。经计算,法矩阵N=的条件数为5.190 5×104,具有病态性,使用本文方法求得的法矩阵条件数随迭代次数的变化如图 7所示。
表 5 参数估值、参数估值与真值之间的2范数和单位权中误差
Tab. 5 Parameter estimation, 2-norm between parameter estimation and true value and the unit weight error
表 6 参数估值的均方根误差
Tab. 6 Root mean square error of parameter estimation
图 6 正则化参数α随迭代次数的变化
Fig. 6 The regularization parameter with the number of iterations
图 7 法矩阵条件数随迭代次数的变化
Fig. 7 The condition number of the normal matrix with the number of iterations
由图 4和图 5可以发现,尽管本文模拟的DTM模型加入的乘性误差标准差仅为0.1,但对高程产生了较大的影响,此时法矩阵的条件数为5.190 5×104,严重病态。因此,针对病态乘性误差模型数据处理理论需要进行更深入的研究。
由表 5可得,由于未考虑模型的病态性,文献[7]中的方法求得的参数估值严重偏离真值,而本文算法求得的结果更接近真值,参数估值的2范数分别为0.817 7和0.770 5,比bcWLS方法小了10.031 8和10.079 0,进一步说明本文算法对于降低病态性有一定的效果。其中,本文算法2考虑到非线性迭代的过程对参数估值的影响,求得的2范数最小,这与算例1得到的结果一致。由表 5中的单位权中误差可知,LS方法求得的单位权中误差0.700 6,结果严重偏离真值;WLS方法和bcWLS方法由于考虑了观测值的权,求得的单位权中误差比LS方法更接近真值,为0.293 0和0.294 4;本文算法考虑了病态性对单位权中误差的影响,求得的单位权中误差最接近真值。
由表 6可得,本文使用正则化方法得到的均方根误差远小于文献[7]中的方法,表明了结果的可靠性。在本文算法中,算法1求得的参数估值、
4.3 算例3
为进一步说明本文方法的可行性和适用性,算例3为一个病态数字高程模型真实数据算例。数字高程模型是一种承载地面高程信息的空间数据模型,是DTM模型的一个分支,本文主要利用SRTM(shuttle radar topography mission)提供的免费DEM数据集,通过在网站http://srtm.csi.cgiar.org/srtmdata/下载了2304个点的三维坐标数据来说明本文算法在实际应用中的优势。SRTM数据是由航天飞机雷达测量生成的LiDAR影像数据[37],且LiDAR数据的噪声已被证明具有乘性噪声性质[38],因此本文含有6个未知参数的数字高程DEM乘性误差模型如下
(37)
式中,(xi, yi)表示地面点经纬度坐标;H(xi, yi)表示该点对应的高程,单位为米;εm表示乘性误差向量。
F(xi, yi)的具体形式如下[7]
(38)
由于SRTM数据均是在等精度测量中获取,因此在定权时假设未知的单位权中误差和乘性误差相等,即乘性误差的权阵为单位阵。本文获取的受乘性误差干扰的数字高程模型如图 8所示,分别使用表 1中的5种方案计算,WLS方法和bcWLS方法由于受病态性的影响,系数矩阵奇异,参数估值不收敛;LS方法、本文算法1和2方法拟合的数字高程模型如图 9、图 10和图 11所示。将3种收敛方案计算的参数估值和单位权中误差估值列于表 7,参数估值的均方根误差列于表 8。
图 8 受乘性误差干扰的数字高程模型
Fig. 8 Digital elevation model disturbed by multiplicative error
图 9 LS方法求得的数字高程模型
Fig. 9 Digital elevation model obtained by the LS method
图 10 RWLS方法求得的数字高程模型
Fig. 10 Digital elevation model obtained by the RWLS method
图 11 SUTRWLS方法求得的数字高程模型
Fig. 11 Digital elevation model obtained by the SUTRWLS method
表 7 参数估值和单位权中误差
Tab. 7 Parameter estimation and the unit weight error
表 8 参数估值的均方根误差
Tab. 8 Root mean square error of parameter estimation
由图 8和图 9可以发现,LS方法求得拟合高程的结果严重偏离原始的数字高程模型,这是由于LS方法未考虑观测值的权且模型为病态性所导致的。本文算法考虑了病态性的影响,在处理实际数据时更具优势,从图 10和图 11来看,求得拟合高程的结果能够较好地拟合图 8的高程。说明现有文献未考虑病态性的方法无法处理实际数据,因此顾及模型病态性带来的影响是十分必要的。
由表 7可以看出,LS方法求得的参数估值和单位权中误差过大,不符合实际情况,进一步说明了LS方法无法处理实际数据,本文算法在处理实际数据的优势。由表 8也可以看出,LS方法求得参数估值的根均方误差最大,本文算法远小于LS方法。其中,在本文算法中,由于算法2考虑了非线性迭代过程对参数估计带来的影响,使用SUT法求得的参数估值的均方根误差小于算法1,这也验证了本文算法的适用性和优势。对比算法1和算法2在参数估值上并没有很大差别,这可能是因为参数估值较小;而参数估值的均方根误差却差别很大,这是因为算法1根据误差传播定律求得的均方根误差只能反映参数估值的一阶精度信息,而算法2采用的SUT法为二阶精度评定方法,求得参数估值的精度信息为二阶精度。
5 结论
为了完善乘性误差模型数据处理理论,本文在分析已有研究的基础上,将乘性误差模型扩展至病态乘性误差模型,推导了病态乘性误差模型参数估计的Tikhonov正则化迭代解公式。考虑到迭代过程中参数估值与观测值之间存在复杂的非线性关系,引入SUT法对病态乘性误差模型进行精度评定。最后,通过算例的验证与分析得到以下几点结论:
(1) 本文推导的病态乘性误差模型参数估计的Tikhonov正则化迭代解法可以得到比未考虑模型病态性的方法更准确的参数估值,且参数估值的精度更高。
(2) 通过将加权最小二乘正则化迭代过程表示成非线性函数,引入的SUT精度评定方法能够得到更为合理的精度信息。
(3) 本文提出的方法是基于函数f(•)为β的线性函数的情况下推导的,若将函数f(•)扩展至非线性函数,需要进一步研究。
作者简介
第一作者简介:王乐洋(1983-), 男, 博士, 教授, 研究方向为大地测量反演及大地测量数据处理。E-mail: wleyang@163.com
初审:张艳玲
复审:宋启凡
终审:金 君
《测绘学报》专刊征稿函 | 滑坡灾害监测预警
《测绘学报(英文版)》(JGGS)专刊征稿 | Call for Papers:空间人文与社会地理计算 (SHGSS)
重磅 | 第八届“测绘科学前沿技术论坛”会议通知(一号)
自然资源部办公厅关于印发《自然资源部2021年政务公开工作要点》的通知