引言

矿产资源权益金测算是矿产资源权益金制度得以顺利实施的基础性工作，也是推进中国权益金制度改革中一个亟待解决的关键问题。由于中国目前尚未建议统一规范的权益金测算标准和方法，为适应矿产资源“三资”化管理的要求，维护国家矿产资源国家权益，有必要开展权益金测算研究。权益金是矿产资源权益价值货币化体现形式^[1]，资源租金是资源所有者凭借其所有权而获得的一种收益，其理论基础是矿租理论。因此，资源禀赋不同是矿产资源权益金产生差异性的根本原因^[2-4]。权益金测算方法研究主要基于税收理论和产权理论两个视角开展。基于税收理论的权益金测算研究主要从宏观经济角度分析，侧重于社会收入分配，如Grobler J.^[5]以南非采矿业为例，通过分析采矿业对国内生产总值贡献率(8.8%)和权益金改革对采矿行业盈利能力影响度分析，从而提出最优权益金率(2%~5%)；Thomas^[6]运用拉弗曲线原理，在对金矿收益测算基础上提出最优权益金率为生产价值的3.5%。随着CGE模型在实践中的运用，Jing Chen等^[7-9]运用CGE模型研究矿产资源行业税费政策对宏观经济影响，在此基础上提出合理资源税费率。李兴国^[10]根据公开的上市矿业公司的经济数据，通过对比研究，用有效税率、经济风险概率等经济指标来定量考察矿业权益金率。基于产权理论的权益金测算研究主要从微观经济角度分析，侧重于资源价值差异化。如朱学义^[11]则以煤炭为例，采用先定级、后定价的方法，从矿产资源价值视角构建权益金评价模型；张彦平^[12]运用模糊评价法进行综合评价计算矿床禀赋优势系数；周怀峰等^[13]以煤炭资源为例，对矿产资源禀赋要素进行了分析研究，运用系统动力学的方法构建了矿产资源禀赋要素体系，用于矿产资源价值评估；李娜^[14]从资源禀赋差异视角提出了矿业权出让收益测算方法。

上述文献从不同角度、不同方法分析了矿产资源权益金评价方法。基于税收理论视角研究权益金，强调权益金分配效果，不可避免的出于“政策需要”等目的，过分重视资源权益价值在国民经济中的分配功能，而忽视了资源禀赋的差异性，不能客观、真实反映个体矿山权益金标准。基于产权理论视角测算权益金，现有研究中多采用模糊评价法、神经网络法和层次分析法，这些方法不能解决评价指标间信息冗余问题，对隶属函数的解释性也存在不足，而层次分析法在对定性变量定量化处理中，存在一定的主观性。因此，权益金评价结果实用性有待商榷。基于此，本文借鉴相关评价指标，采用定性与定量指标相结合的方法，运用多元回归分析，建立了基于矿产资源资源禀赋条件视角的权益金评价模型。

1 权益金测算模型构建 1.1 回归模型及假定检验

回归模型是定量描述相关关系的数学模型。设因变量为y，自变量中定量变量为d₁, d₂, …, d_k，定性变量为Q₁, Q₂, …, Q_k，将定性变量转化为k个虚拟变量组x=(x₁, x₂, …, x_j)′后，因变量和自变量的多元回归关系为：

$ y=A x+\sum\limits_{i=1}^{k} B_{i} d_{i}+\varepsilon $

(1)

其中：A=(α₁, α₂, …, α_j)，为定量变量的回归系数，d_i=(D_i1, D_i2, D_imj-1)′，为第i个定性变量Q_i对应的m_i-1个虚拟变量，B_i=(β_i1, β_i2，…, β_imj-1), 为Q_i对应的虚拟变量的回归系数，ε为误差项。

鉴于调研区域成矿条件类似，导致一些影响因素的水平重复率较高，本文在拟合多元回归时采用了逐步回归法剔除变量，依据AIC信息准则判断模型拟合的优劣。

设因变量g(y)的观测样本为y₁, y₂, …, y_n，总体概率密度为f(y; q)，θ^*为θ的极大似然估计，那么θ^*的对数似然函数为：

$ 1\left( {{\theta ^*}} \right) = \sum\limits_{i = 1}^n {\mathit{ln}} f\left( {{y_i};{\theta ^*}} \right) $

(2)

假设模型的误差服从独立正态分布，则有：

$ A I C=2 k-21\left(\theta^{*}\right) $

(3)

k为θ的维数，k越小则AIC越小，1(θ^*)越大则AIC越小。k小说明模型简洁，1(θ^*)大说明模型精确。AIC值越小，说明g(y)与f(y; q)之间越接近。

假定1：随机扰动项与解释变量不相关。检验：绘制自变量与残差的线性关系图，并观察如果自变量与残差之间存在明显非线性关系，那么需要对自变量做非线性转化，之后重新拟合模型。

假定2：正态性假定。检验：对残差$\hat u$排序，${\hat u_1} \le {\hat u_2} \le \cdots \le {\hat u_n}$，分位数$Q_{i}=\frac{\hat{u}_{i}-\overline{u}}{n}, t_{i}=\frac{i-0.5}{n}$，对照正态分布表查询t_i的分位数Q_i'，之后绘制Q-Q图检验残差的正态性分布，如果Q-Q对应点较为均匀的分布在y=x曲线附近，且全部落在90%区间，如果对应的点都落在可接受区间内，则残差符合正态性分布。

假定3：同方差假定。检验：用布罗特-帕甘检验(BP-test)验证模型的同方差性。BP检验的基本思想：模型的残差为$\hat u$，构造线性函数:

$ \hat{u}^{2}=\delta_{1} x_{1}+\delta_{2} x_{2}+\cdots+\delta_{k} x_{k}+v $

(4)

原假设：H₀再解析δ₁=δ₂=…=δ_k=0，备择假设：H₁：残差平方和与一个或多个自变量相关，即H₀不成立。对式(2.4)构造统计量：

$ F = \frac{{RSS/k}}{{ESS/(n - k - 1)}} = \frac{{R_{\hat u}^2/k}}{{\left( {1 - R_{\hat u}^2} \right)/(n - k - 1)}} $

(5)

其中RSS，ESS，$R_{\hat u}^2$分别为式(2.4)的回归平方和，残差平方和，拟合优度。若H₀成立，那么F~F(k, n-k-1)；如果F>F_α(k, n-k-1)，那么拒绝原假设，即存在异方差。

假定4：无自相关假定。检验：对残差做DW检验(Durbin-Watson检验)，构造统计量：

$ DW = \frac{{\sum\limits_j^n = 2{{\left( {{{\hat u}_i} - {{\hat u}_{i - 1}}} \right)}^2}}}{{\sum\limits_j^n = 1\hat u_j^2}} $

(6)

其中$\hat{u}_{i}=y_{i}-\overline{y}_{i}$，是ε的估计值。由自变量个数k查DW表，得到下限dL和上限dU，dU≤DW≤4-dU时，则残差之间互相独立。

假定5：无多重共线性假定。检验：为了检验自变量之间是否存在多重共线性，需要计算自变量方差膨胀因子(VIF)。

$ V I F_{j}=1 /\left(1-R_{j}^{2}\right) $

(7)

式中，R_j²是以x_j为因变量时对其它自变量回归的拟合优度。如果最大的VIF超过10，则表示存在明显多重相关性，需要剔除存在多重相关变量，否则将严重影响模型估计效果。

1.2 样本选择与指标选取

(1) 样本选择。2016年，安徽省登记在产的铁矿矿山107座，其中大型18座，中型13座，小型63座，零星小矿14座。根据研究需要确定样本选取原则：①矿产地质数据完整；②2011—2016年连续生产；③财务报表齐完整真实。在这107座矿山中，选取了38座铁矿矿山作为调查样本，样本占总矿山数比例为35.51%，样本矿山产量占全省的84.63%，地质储量占全省铁矿矿山地质储量的65.54%(资料来源《安徽省矿产资源储量表 2016》《矿产资源开发利用年报(2011—2016)》。

(2) 指标选取。矿产资源禀赋条件差异是产生级差矿租决定性因素，权益金大小取决于矿产资源禀赋条件的优劣。本文在参考袁怀雨^[15]、郑重^[16]、赵志勇^[17]等建立的评价指标体系基础上，以蒋健明^[18]提出的权益金评价指标体系为依据，进行指标选取(见表 1)。

1.3 虚拟变量设置与筛选

权益金评价指标体系中有定量指标，也有定性指标。对定量变量直接引入原始值，对定性变量，采用虚拟变量法。虚拟变量设置原则：如果变量φ有t个类别，虚拟变量设置t-1个。W₁, W₂, …, W_t-1的关系为：最多只能有一个为1，其余必须为0。W₁为1代表φ为类别1，W₂为1代表φ为类别2，依次类推。W₁, W₂, …, W_t-1全部为0时，则表示φ为类别t。

根据研究样本，对于定量指标，如地质品位等5项指标，直接引入回归模型计算。对于定性指标，将9定性变量分别转化为好、中、差3个水平等级，剔除定性变量中水平重复率达到90%以上的指标，然后把保留下来的定性变量转化成虚拟变量引入模型计算。如图 1，本研究以安徽省铁矿为例，其地理位置、气候条件、产业配套条件类似，重复率超过90%，又由于铁矿的形成有其自然地理特征，并不是在安徽省范围内平均分布，其地质特征类似，故地形地貌、矿石类型等项指标样本重复率也超过90%，因此只能将其余3个定性变量(矿体形态、水文地质条件和有害元素)引入模型(见图 1)。

由于矿体形态有好、中、差3个形态，对其设置2个虚拟变量W₁₁和W₁₂。W₁₁代表矿体形态为差，W₁₂=1代表矿体形态为中。W₁₁和W₁₂同时为0代表矿体形态为好。同理，对水文地质设置虚拟变量，W₂₁=1为水文地质条件差，W₂₂为水文地质条件中，W₂₁=W₂₂=0为水文地质好；对有害元素及伴生元素设置虚拟变量W₃₁为有害元素及伴生元素差，W₃₂为有害元素及伴生元素中，W₃₁=W₃₂=0为有害元素及伴生元素好。据此，对变量进行定义(见表 2)。

1.4 权益金测算方法与参数确定

(1) 权益金测算方法。矿产资源权益金是资源租金的货币化体现形式^[1]，其实质是资源未来租金的资本化，即：权益金=资源租金/资本化率。由于矿产资源在未开发之前总是埋藏于地下，无法直接估算，一般用剩余法进行测算^[19-21]。即权益金是矿产品销售收入扣除生产成本、矿业行业平均利润、税收和矿业风险投入回报之后的剩余，即因开发矿产资源而产生的超额利润，其测算公式为：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{V}{\rm{ = }}\\ \sum\limits_{i = 1}^n {\frac{{\left[ {s - {c_a} \cdot \left( {1 \pm \sum\limits_{i = 1}^r {{k_i}} } \right) - \lambda - f - m - {p_a} - {T_a} - {I_a}} \right]}}{{{{(1 + r)}^n}}}} \end{array} $

(8)

$ \pi=\frac{V}{s} $

(9)

其中，V：权益金；s:矿产品销售收入；c_a：矿产品直接生产标准成本；K₁：资源禀赋条件修正系数(包括地质储量/生产规模、地质品位、平均开采深度等)；λ：销售费用；f：财务费用；m：管理费用；T_a：综合税负；p_a：平均利润；I_a：矿业风险投入回报；r：资本化率；n：开采年限；π：权益金率。

(2) 相关参数确定。根据公式(8)，权益金的测算需要确定相关参数，其依据和方法分别是：①直接生产标准成本，首先通过对调查的38个样本财务报表进行分析，由于各矿山的投资结构、生产技术与管理水平以及铁矿资源禀赋差异等原因，铁精矿直接生产成本差别较大，需要进行成本修正，修正假设条件参照条件为：用平均资产负债率指标和采选业成本费用平均利润率表示平均管理水平(资料来源：中国冶金矿山财务资料汇编(2013—2017)中57个铁矿企业平均资产负债率57%)；用人均产量代替生产力水平(资料来源：安徽省矿产资源开发利用情况年报(2013—2017)中产量与从业人员比值计算，人均生产矿石量1 430 t)。②销售费用、管理费用、财务费用和综合税负分别为1.5%、9%、4%和3.5%(资料来源：根据2013—2017年中国冶金矿山财务资料汇编中57个铁矿企业平均相关费用计算)。③矿业企业平均利润率与矿业风险投入回报。矿业是高风险行业，矿业企业平均利润率取值参照中国(1995—2016年)采矿业企业平均销售利润率，并结合《中国矿业权评估准则》中采矿权折现率取值标准，按销售收入的8%计算(资料来源：《中国统计年鉴》，1995—2016年中国规模以上工业企业和采矿业企业平均销售利润率；《安徽统计年鉴》安徽省1998—2016规模以上工业企业及四大类矿产采选业平均销售利润率)。

1.5 回归模型拟合

对该模型进行逐步回归拟合，根据AIC准则逐步剔除变量(表 3)。

一般而言，矿山生产规模也是矿山企业决定投资建设的一项重要指标，因而是相对重要的，但我们在选取样本时发现，由于2010—2016年铁矿石价格波动较大，国际铁矿石价格从138美元跌至40美元，铁矿石价格的连续下跌导致生产成本高居的中小型矿山企业被迫停产，同时也由于环保的要求，从2012年起，安徽省政府要求年产3万t以下的铁矿矿山一律关闭，20万t以下的小型铁矿矿山进行整合(事实上从2012年起，小型及以下规模矿山将不能取得矿山开采许可证)，导致小型以下矿山样本选择空间有限。从而导致矿山规模在回归分析中并不显著。基于尊重统计数据原则，不得已，我们舍去了规模因素(回归结果见表 4)。根据表 4，可得矿山直接生产成本和指标之间的回归方程：

$ \begin{array}{l}{c=541.90-290.71 \times L+0.068 \times D-} \\ {0.73 \times T+51.13 \times W_{21}+22.09 \times W_{22}}\end{array} $

(10)

2 检验分析与测算结果 2.1 模型检验

(1) 线性假设检验。先对地质水平，开采深度和平均厚度3个连续型自变量检测其线性假设，对上述拟合模型，调用R语言自带的成分残差图函数绘制各自变量与成本残差的关系，绘制方法为：对level_i, i=1, 2, 3, L 38的不同值，成分残差图散点的横、纵轴坐标为(level_i, $\cos {\hat t_i}$-cost_i，其中cost_i为第i个观测值，$\cos {\hat t_i}$为根据上述模型得到的第i个拟合值。可以得到图 2，自变量与残差之间没有明显的非线性关系。

(2) 正态性检验。对因变量的残差绘制Q-Q图，图 3显示，Q-Q点分布和y=x线大致重合，说明模型正态性良好。

(3) 同方差性检验。用BP检验验证残差的同方差性，卡方检验值：1.114 5；p-value：0.291 1>0.05。不存在异方差性。

(4) 误差独立性检验。残差的DW检验结果：滞后阶数：1；自相关系数：0.122 9；D-W统计量：1.657 9；p-value：0.224>0.05。误差之间相互独立。

(5) 多重共线性。各个自变量的方差膨大因子结果分别是：地质品位(L)：1.345 5；开采深度(D)：1.475 9；矿体厚度(T)：1.082 6；水文地质条件W₂₁：2.029 6，W₂₂：1.885 4。由于各自变量方差膨大因子较小，可以认为不存在多重共线性。

由于λ, f, m, p_a, t_a, I_a都可以用s, c_a的线性组合表示，结合公式(8)和(9)，可以得出：

$ \pi = 1 - {\rm{ }}\sum\limits_{i\mathitÓ\left\{ {\lambda ,f,m,{p_a},{T_a},{I_a}} \right\}} {{\Lambda _i} \times {{\left( {S,{C_i}} \right)}^\prime }/S} {\rm{ }} $

(11)

其中Λ_i为将C_i表示为s, c_a的线性组合时的系数。经过合并同类项后，得出：

$ \pi=60.1 \%-0.83 \times c_{a} / s $

(12)

将公式(10)拟合结果代入公式(12)可以得到给定矿山的权益金公式：

$ \begin{array}{l} \pi = 60.1\% - 0.83 \times (541.90 - 290.71 \times L + 0.068 \times D - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;0.73 \times T + 51.13 \times {W_{21}} + 22.09 \times {W_{22}})/S \end{array} $

(13)

根据公式(13)，在铁矿石资源禀赋条件确定和每吨销售价格已知的情况下，可计算出相应的权益金率。测算当铁矿石价格为600~900元/t时的样本矿山权益金率，用测算出来的权益金率实际值进行拟合，最终提出铁矿矿山基准权益金率修正系数。

2.2 权益金测算结果

以全部因素的中间水平作为基准，得出权益金率的基准值π_b，在此基础上计算不同因素水平下各水平对应的调整系数k_i(i=1, 2, 3, 4)，其中，k₁：地质品位调整系数，k₂：平均开采深度调整系数，k₃：矿体平均厚度调整系数，k₄：水文条件调整系数。用测算表测算矿山权益金率时，算法是查得矿山因素水平对应的分组，找到对应的调整系数，并且把各调整系数与基准相加，见式(14):

$ \pi=\pi_{b}+\sum\limits_{i=1}^{4} k_{i} $

(14)

按照相关规范和标准，为便于区分，将影响因素都被分为3个水平。对于平均开采深度(D)，相对于安徽省地下开采的铁矿赋存深度，条件较好的分隔点为D_u=300，条件较差的分割点为D_d=500，中等水平为D_u < =depth < =D_d，计算基准时用中等水平两个端点的均值代表中等水平，D^*=(D_u+D_d)/2=400。当D < 300时，用min(D)和D_u的均值代表该水平，[min(D)+D_u]/2=165，min(D)为平均开采深度的样本最小观测值；当D≥500时，用min(D)和D_d的均值代表该水平，[max(D)+D_d]/2=525，max(D)为平均开采深度的样本最大观测值。地质品位、矿体平均厚度的计算同理。根据《中华人民共和国地质矿产行业标准(DZ/T 0200—2002)》，铁矿一般工业指标是边界品位20%，工业品位为25%，但现实情况中，由于安徽铁矿地下开采深度平均在200~400 m左右，生产成本较高，平均地质品位低于25%矿山在现有的技术条件和矿产品价格下，不具备开采价值(露采矿山除外)，在实际中无法选取到平均地质品位低于25%的样本矿山，所以平均地质品位低于25%的调整值不做处理。对于水文条件，将虚拟变量对应赋值为0或1。在此基础上，控制变量，计算单一变量的值变动时，对应权益金率变动值(测算表见表 5)。

根据公式(13)和表 5，求得铁精矿(62%)价格为700~900元/t本矿山权益金率实际值与拟合值比较表(图 4)。

图 4显示了各样本矿山在铁精矿(62%)价格分别为600、700、800、900元/t条件下的实际值与拟合值对比情况，柱状图的每根柱子(黑色代表实际值，蓝色代表拟合值)分别按顺序对应上表中不同矿山在不同价格水平下权益金率。可以看出，在不同价格水平下，实际值与拟合值的相对走势接近，正负性基本一致。

图 5显示了四个样本矿山在铁精矿(62%)价格分别为600、700、800、900元/t条件下权益金率实际值与拟合值对比情况(实线代表实际值，虚线代表拟合值)，拟合度较好，限于篇幅的原因，本图只选取了四个矿山。

2.3 配对样本T检验

为了评估模型拟合效果，需要从统计学角度证明实际值与拟合值之间有无显著性差异。现对每个矿山拟合前后的权益金率之差构造t统计量，采用配对样本t检验。

先检验不同价格水平下权益金率实际值和拟合值之间之间差值是否服从正态分布。由于样本为小样本，采用S-W法(Shapiro-Wilk)，原假设H₀：权益金率的实际值和拟合值之间的差值服从正态分布。经检验，p值为0.064，说明不能拒绝原假设，因此铁精矿(62%)价格水平在600~900元/t下，权益金率实际值和拟合值之间差值是服从正态分布的(见表 6)。

再对不同价格水平下权益金率实际值和拟合值做配对样本T检验，原假设H₀为：权益金率实际值和拟合值的均值无显著差异。经检验，p>0.05，说明不能拒绝原假设，因此权益金率的实际值和拟合值的均值之间无显著差异，拟合权益金率与实际权益金率线性关系显著(见表 7)。因此，模型可以很好的拟合权益金率。

2.4 计算实例

现以姑山铁矿为例，介绍测算过程。姑山铁矿位于中国安徽省东部马鞍山地区，该矿累计查明储量为1.28亿t，矿床类型是与火山作用有关的矿床，铁矿类型是磁铁矿(TFe)，平均地质TFe品位43.94%，S平均品位0.07%，P平均品位0.51%。矿体平均厚度为60 m，矿体倾斜度40°~60°，水文地质条件为中等。开采方式为地下开采，采矿方法是分层充填采矿法，平均开采深度为-150 m，年产矿石量130万t。

首先根据表 2确定姑山铁矿该矿的赋存水平等级，然后再对照铁矿矿山权益金系数修正表(表 5)，结合矿产品市场价格，确定的相应等级修正系数。假如当前铁精矿(62%)市场价格为600元(含税价)，姑山铁矿的地质品位修正系数为4.64%，平均开采深度修正系数为2.22%，矿体厚度修正系数为2.29%，水文地质条件修正系数为0，基准系数为-6.10%，根据公式(5.16)，姑山铁矿的权益金率为：-6.10%+4.64%+2.29%+2.22%+0=3.05%，每吨矿产品权益金为：600×3.05%=18.3元。

3 结论

为了实现对矿产资源权益金精确评估，本文从矿产资源价值视角，运用多元回归模型和虚拟变量法，构建了基于矿产资源禀赋视角的矿产资源权益金评价模型，并以中国安徽省铁矿资源为例对权益金评价模型的适用性进行了验证。本文得到如下研究结论：

(1) 基于矿产资源禀赋视角的矿产资源权益金评价模型，弥补了因单纯从收入分配视角研究权益金而造成“价值失真”方法的不足；从评价指标的处理上，更注重定量化处理，避免了在权重处理上的主观性。为政府和企业客观测算矿产资源权益金提供了简单、可靠的评估模型。

(2) 基于矿产资源禀赋视角的矿产资源权益金评价模型，可以直观准确的反应出矿产资源禀赋条件对权益金大小的影响。在诸多影响因素中，铁矿资源的地质平均品位、平均开采深度、矿体平均厚度和水文地质条件等4个因素对其权益金标准影响尤为显著。

(3) 矿产资源权益金与矿产品价格在一定区间内呈正相关关系，当矿产品价格低到一定程度时，权益金率会是负数，当然这种情况将导致该矿产资源无人开采，当矿产品价格无限上升，权益金率的变化并不是无限增大，会无限逼近一个常数。

本文在中国权益金制度改革背景下，提出了基于矿产资源禀赋视角的矿产资源权益金评价模型，为权益金制度落实提供了便捷的测算方法，同时，该模型结果客观可靠，能够动态反应权益金变化情况。矿产资源种类繁多，各矿种的特点和用途大相径庭，本文虽仅以铁矿为例对权益金进行了研究，但这一新的研究思路为其他矿种权益金评估、矿业税费的测算提供了参考。本研究的主要局限在于数据获得的地域性限制，仅以安徽省域内的铁矿资源禀赋条件为基础对权益金进行测算，如果从全国范围内开展权益金研究，则需要更加完整的数据加以诠释，这也是今后需要补充和完善的内容。